From eab3ece76212819a2759873b49dc563ddb7181d6 Mon Sep 17 00:00:00 2001 From: JeremyGelb Date: Sat, 11 May 2024 21:45:34 -0400 Subject: [PATCH] another try --- .github/workflows/rhub.yaml | 8 +- src/bw_selection_tnkde.cpp.bak | 1485 -------------------------------- src/matrices_functions.cpp | 12 - src/nkde_continuous.cpp.bak | 1238 -------------------------- 4 files changed, 4 insertions(+), 2739 deletions(-) delete mode 100644 src/bw_selection_tnkde.cpp.bak delete mode 100644 src/nkde_continuous.cpp.bak diff --git a/.github/workflows/rhub.yaml b/.github/workflows/rhub.yaml index a2471e6e..930b41bb 100644 --- a/.github/workflows/rhub.yaml +++ b/.github/workflows/rhub.yaml @@ -76,6 +76,10 @@ jobs: config: ${{ fromJson(needs.setup.outputs.platforms) }} steps: + - name: Query dependencies + run: | + install.packages('terra', repos='https://rspatial.r-universe.dev') + shell: Rscript {0} - uses: r-hub/actions/checkout@v1 - uses: r-hub/actions/setup-r@v1 with: @@ -85,10 +89,6 @@ jobs: with: token: ${{ secrets.RHUB_TOKEN }} job-config: ${{ matrix.config.job-config }} - - name: Query dependencies - run: | - install.packages('terra', repos='https://rspatial.r-universe.dev') - shell: Rscript {0} - uses: r-hub/actions/setup-deps@v1 with: job-config: ${{ matrix.config.job-config }} diff --git a/src/bw_selection_tnkde.cpp.bak b/src/bw_selection_tnkde.cpp.bak deleted file mode 100644 index ce95fdbc..00000000 --- a/src/bw_selection_tnkde.cpp.bak +++ /dev/null @@ -1,1485 +0,0 @@ -#include "spNetwork.h" -#include "base_kernel_funtions.h" -#include "matrices_functions.h" - -#include - - -// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -// THE FUNCTIONS TO CALCULATE BW SELECTION CV LOO TEMPORAL FOR SIMPLE TNKDE -// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -//' @title The worker function to calculate simple TNKDE likelihood cv -//' @name ess_kernel_loo_tnkde -//' @param kernel_func a cpp pointer function (selected with the kernel name) -//' @param edge_mat matrix, to find the id of each edge given two neighbours. -//' @param events a NumericVector indicating the nodes in the graph being events -//' @param time_events a NumericVector indicating the timestamp of each event -//' @param neighbour_list a List, giving for each node an IntegerVector with -//' its neighbours -//' @param v the actual node to consider (int) -//' @param v_time the time of v (double) -//' @param bws_net an arma::vec with the network bandwidths to consider -//' @param bws_time an arma::vec with the time bandwidths to consider -//' @param line_weights a vector with the length of the edges -//' @param depth the actual recursion depth -//' @param max_depth the maximum recursion depth -//' @return a cube with the impact of the event v on each other event for -//' each pair of bandwidths (cube(bws_net, bws_time, events)) -//' @keywords internal -arma::cube ess_kernel_loo_tnkde(fptros kernel_func, arma::sp_imat &edge_mat, - IntegerVector &events, - IntegerVector &events_wid, - NumericVector &time_events, - List &neighbour_list, - int v, int wid, double v_time, - arma::vec &bws_net, arma::vec &bws_time, - NumericVector &line_weights, int max_depth){ - - //step0 : generate the queue - int depth = 0; - std::queue data_holder; - arma::cube kvalues(bws_net.n_elem, bws_time.n_elem, events.length()); - kvalues.fill(0.0); - //step1 : generate the first case - List cas1 = List::create(Named("d")=0.0, - Named("v") = v, - Named("prev_node") = -999, - Named("depth") = 0 - ); - data_holder.push(cas1); - - double max_bw_net = arma::max(bws_net); - double bw_net, bw_time; - - // avant de lancer les iteration, nous devons vérifier s'il est necessaire - // de rajouter la densite a des evenement qui se situeraient sur le meme point - // de depart ! - std::vector index = get_all_indeces_int(events,v); - if(index.size() > 1){ - for(int ii = 0; ii < bws_net.n_elem ; ii++){ - bw_net = bws_net(ii); - double kernel_net = kernel_func(0,bw_net); - for(int j = 0 ; j < bws_time.n_elem; j ++){ - bw_time = bws_time(j); - // NOTE : we are doing the bw scaling here - for (int zz : index){ - if(time_events[zz] != v_time){ - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(ii,j,zz) = kvalues(ii,j,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - } - } - - //lancement des iterations - while(data_holder.empty()==FALSE){ - //unpacking (imagine some loop unrolling here with a function to deal with.) - List cas = data_holder.front(); - data_holder.pop(); - int v = cas["v"]; - int depth = cas["depth"]; - double d = cas["d"]; - int prev_node = cas["prev_node"]; - - //step1 : find all the neighbours - IntegerVector neighbours = neighbour_list[v-1]; - - //step2 : iterate over the neighbours - int cnt_n = neighbours.length(); - int new_depth; - if(cnt_n>2){ - new_depth = depth+1; - }else{ - new_depth = depth; - } - - //if we have only one neighbour, we must stop - if(cnt_n>1 or prev_node<=0){ - for(int i=0; i < cnt_n; ++i){ - int v2 = neighbours[i]; - //on ne veut pas revenir en arriere ! - if(v2!=prev_node){ - // find the edge between the two nodes - int edge_id = edge_mat(v,v2); - double d2 = line_weights[edge_id-1] + d; - std::vector index = get_all_indeces_int(events,v2); - if(index.size() >0 ){ - // il semble que v2 soit un noeud pour lequel au moins un evenement est present - for(int ii = 0; ii < bws_net.n_elem ; ii++){ - bw_net = bws_net(ii); - double kernel_net = kernel_func(d2,bw_net); - for(int j = 0 ; j < bws_time.n_elem; j ++){ - bw_time = bws_time(j); - // NOTE : we are doing the bw scaling here - for (int zz : index){ - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(ii,j,zz) = kvalues(ii,j,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - - } - } - } - //evaluating for the next move - if (d2< max_bw_net and new_depth data_holder; - arma::cube kvalues(bws_net.n_rows, bws_net.n_cols, events.size()); - kvalues.fill(0.0); - //step1 : generate the first case - List cas1 = List::create(Named("d")=0.0, - Named("v") = v, - Named("prev_node") = -999, - Named("depth") = 0 - ); - data_holder.push(cas1); - - // NOTE : POSSIBLE ENHANCEMENT HERE, I COULD RESET THE max_bw_net - // DURING THE ITERATIONS TO REDUCE A BIT THE CALCULATION TIME. WE WOULD - // AVOID SOME ITERATION THIS WAY. - double max_bw_net = bws_net.max(); - double bw_net, bw_time; - - // avant de lancer les iterations, il est important de verifier que - // certains evenement ne se trouvent pas sur la meme vertex d'origine. - // si c'est le cas, nous devons leur assigner la valeur de densite correspondante - std::vector index = get_all_indeces_int(events,v); - if(index.size() > 1){ - // Rcout << "See me ! I must add some density at the origin ! \n"; - for (int zz : index){ - if(time_events[zz] != v_time){ - // arma::mat slice_bw_net = bws_net.slice(zz); - // arma::mat slice_bw_time = bws_time.slice(zz); - //Rcout << "step8\n"; - for(int j = 0 ; j < bws_net.n_rows; j ++){ - for(int jj = 0 ; jj < bws_time.n_cols; jj ++){ - // NOTE : we are doing the bw scaling here - double bw_net = bws_net(j,jj); - double bw_time = bws_time(j,jj); - double kernel_net = kernel_func(0.0, bw_net); - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(j,jj,zz) = kvalues(j,jj,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - } - } - - - - //lancement des iterations - while(data_holder.empty()==FALSE){ - //unpacking (imagine some loop unrolling here with a function to deal with.) - List cas = data_holder.front(); - data_holder.pop(); - int v = cas["v"]; - int depth = cas["depth"]; - double d = cas["d"]; - int prev_node = cas["prev_node"]; - //step1 : find all the neighbours - IntegerVector neighbours = neighbour_list[v-1]; - - //step2 : iterate over the neighbours - int cnt_n = neighbours.length(); - int new_depth; - if(cnt_n>2){ - new_depth = depth+1; - }else{ - new_depth = depth; - } - //if we have only one neighbour, we must stop - if(cnt_n>1 or prev_node<=0){ - for(int i=0; i < cnt_n; ++i){ - int v2 = neighbours[i]; - // comment retrouver le potentiel wid du noeud ? - //on ne veut pas revenir en arriere ! - if(v2!=prev_node){ - //Rcout << "step4\n"; - // find the edge between the two nodes - int edge_id = edge_mat(v,v2); - double d2 = line_weights[edge_id-1] + d; - //Rcout << "step5\n"; - std::vector index = get_all_indeces_int(events,v2); - //Rcout << "step6\n"; - if(index.size() >0 ){ - // Rcout << "step7\n"; - // il semble que v2 soit un noeud pour lequel au moins un evenement est present - // NB : nrow pour les BW network et ncol pour le bw temporelles - for (int zz : index){ - //arma::mat slice_bw_net = bws_net.slice(zz); - //arma::mat slice_bw_time = bws_time.slice(zz); - //Rcout << "step8\n"; - for(int j = 0 ; j < bws_net.n_rows; j ++){ - for(int jj = 0 ; jj < bws_time.n_cols; jj ++){ - // NOTE : we are doing the bw scaling here - double bw_net = bws_net(j,jj); - double bw_time = bws_time(j,jj); - double kernel_net = kernel_func(d2, bw_net); - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - //Rcout << "step9\n"; - kvalues(j,jj,zz) = kvalues(j,jj,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - - } - //evaluating for the next move - if (d2< max_bw_net and new_depth data_holder; - arma::cube kvalues(bws_net.n_elem, bws_time.n_elem, events.length()); - kvalues.fill(0.0); - //step1 : generate the first case - List cas1 = List::create(Named("d")=0.0, - Named("v") = v, - Named("prev_node") = -999, - Named("depth") = 0, - Named("alpha") = 1 - ); - data_holder.push(cas1); - - double max_bw_net = arma::max(bws_net); - double bw_net, bw_time; - - // avant de lancer les iteration, nous devons vérifier s'il est necessaire - // de rajouter la densite a des evenement qui se situeraient sur le meme point - // de depart ! - std::vector index = get_all_indeces_int(events,v); - if(index.size() > 1){ - for(int ii = 0; ii < bws_net.n_elem ; ii++){ - bw_net = bws_net(ii); - double kernel_net = kernel_func(0,bw_net); - for(int j = 0 ; j < bws_time.n_elem; j ++){ - bw_time = bws_time(j); - // NOTE : we are doing the bw scaling here - for (int zz : index){ - if(time_events[zz] != v_time){ - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(ii,j,zz) = kvalues(ii,j,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - } - } - - - //lancement des iterations - while(data_holder.empty()==FALSE){ - //unpacking (imagine some loop unrolling here with a function to deal with.) - List cas = data_holder.front(); - data_holder.pop(); - - int v = cas["v"]; - int depth = cas["depth"]; - double d = cas["d"]; - double alpha = cas["alpha"]; - int prev_node = cas["prev_node"]; - - - //step1 : find all the neighbours - IntegerVector neighbours = neighbour_list[v-1]; - - //step2 : iterate over the neighbours - int cnt_n = neighbours.length(); - int new_depth; - if(cnt_n>2){ - new_depth = depth+1; - }else{ - new_depth = depth; - } - - //step 3 prepare the new alpha value - double new_alpha; - if((prev_node < 0) && (cnt_n > 2)){ - new_alpha = 2.0/(cnt_n); - }else if((prev_node < 0) && (cnt_n == 1)){ - new_alpha = 1.0; - }else{ - new_alpha = alpha * (1.0/(cnt_n-1.0)); - } - - //if we have only one neighbour, we must stop - if(cnt_n>1 or prev_node<=0){ - for(int i=0; i < cnt_n; ++i){ - int v2 = neighbours[i]; - //on ne veut pas revenir en arriere ! - if(v2!=prev_node){ - //find the edge between the two nodes - int edge_id = edge_mat(v,v2); - double d2 = line_weights[edge_id-1] + d; - - //est ce que v2 est un evenement pour lequel on doit garder la valeur - - std::vector index = get_all_indeces_int(events,v2); - if(index.size() >0 ){ - - for(int ii = 0; ii < bws_net.n_elem ; ii++){ - bw_net = bws_net(ii); - double kernel_net = kernel_func(d2,bw_net) * new_alpha; - for(int j = 0 ; j < bws_time.n_elem; j ++){ - bw_time = bws_time(j); - for (int zz : index){ - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(ii,j,zz) = kvalues(ii,j,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - } - - //evaluating for the next move - if (d2< max_bw_net and new_depth data_holder; - arma::cube kvalues(bws_net.n_rows, bws_time.n_cols, events.length()); - kvalues.fill(0.0); - //step1 : generate the first case - List cas1 = List::create(Named("d")=0.0, - Named("v") = v, - Named("prev_node") = -999, - Named("depth") = 0, - Named("alpha") = 1 - ); - data_holder.push(cas1); - - double max_bw_net = bws_net.max(); - double bw_net, bw_time; - - // avant de lancer les iterations, il est important de verifier que - // certains evenement ne se trouvent pas sur la meme vertex d'origine. - // si c'est le cas, nous devons leur assigner la valeur de densite correspondante - std::vector index = get_all_indeces_int(events,v); - if(index.size() > 1){ - // Rcout << "See me ! I must add some density at the origin ! \n"; - for (int zz : index){ - if(time_events[zz] != v_time){ - //arma::mat slice_bw_net = bws_net.slice(zz); - //arma::mat slice_bw_time = bws_time.slice(zz); - //Rcout << "step8\n"; - for(int j = 0 ; j < bws_net.n_rows; j ++){ - for(int jj = 0 ; jj < bws_time.n_cols; jj ++){ - // NOTE : we are doing the bw scaling here - double bw_net = bws_net(j,jj); - double bw_time = bws_time(j,jj); - double kernel_net = kernel_func(0.0, bw_net); - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(j,jj,zz) = kvalues(j,jj,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - } - } - - - - //lancement des iterations - while(data_holder.empty()==FALSE){ - //unpacking (imagine some loop unrolling here with a function to deal with.) - List cas = data_holder.front(); - data_holder.pop(); - int v = cas["v"]; - int depth = cas["depth"]; - double d = cas["d"]; - double alpha = cas["alpha"]; - int prev_node = cas["prev_node"]; - - - //step1 : find all the neighbours - IntegerVector neighbours = neighbour_list[v-1]; - - //step2 : iterate over the neighbours - int cnt_n = neighbours.length(); - int new_depth; - if(cnt_n>2){ - new_depth = depth+1; - }else{ - new_depth = depth; - } - - //step 3 prepare the new alpha value - double new_alpha; - if((prev_node < 0) && (cnt_n > 2)){ - new_alpha = 2.0/(cnt_n); - }else if((prev_node < 0) && (cnt_n == 1)){ - new_alpha = 1.0; - }else{ - new_alpha = alpha * (1.0/(cnt_n-1.0)); - } - - //if we have only one neighbour, we must stop - if(cnt_n>1 or prev_node<=0){ - for(int i=0; i < cnt_n; ++i){ - int v2 = neighbours[i]; - //on ne veut pas revenir en arriere ! - if(v2!=prev_node){ - //find the edge between the two nodes - int edge_id = edge_mat(v,v2); - double d2 = line_weights[edge_id-1] + d; - - //est ce que v2 est un evenement pour lequel on doit garder la valeur - - std::vector index = get_all_indeces_int(events,v2); - if(index.size() >0 ){ - for (int zz : index){ - //arma::mat slice_bw_net = bws_net.slice(zz); - //arma::mat slice_bw_time = bws_time.slice(zz); - for(int j = 0 ; j < bws_net.n_rows; j ++){ - for(int jj = 0 ; jj < bws_time.n_cols; jj ++){ - // NOTE : we are doing the bw scaling here - double bw_net = bws_net(j,jj); - double bw_time = bws_time(j,jj); - double kernel_net = kernel_func(d2,bw_net) * new_alpha; - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(j,jj,zz) = kvalues(j,jj,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - // for(int ii = 0; ii < bws_net.n_elem ; ii++){ - // bw_net = bws_net(ii); - // double kernel_net = kernel_func(d2,bw_net) * new_alpha; - // for(int j = 0 ; j < bws_time.n_elem; j ++){ - // bw_time = bws_time(j); - // for (int zz : index){ - // double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - // kvalues(ii,j,zz) = kvalues(ii,j,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - // } - // } - // } - } - - //evaluating for the next move - if (d2< max_bw_net and new_depth data_holder; - std::vector data_holder; - - // let us prepare the first cases - IntegerVector v_neighbours = neighbour_list[v-1]; - double alpha = 2.0/v_neighbours.length(); - - acase el = {v,-999,0,0.0,alpha}; - data_holder.push_back(el); - - // avant de lancer les iteration, nous devons vérifier s'il est necessaire - // de rajouter la densite a des evenement qui se situeraient sur le meme point - // de depart ! - std::vector index = get_all_indeces_int(events,v); - if(index.size() > 1){ - for(int ii = 0; ii < bws_net.n_elem ; ii++){ - bw_net = bws_net(ii); - double kernel_net = kernel_func(0,bw_net); - for(int j = 0 ; j < bws_time.n_elem; j ++){ - bw_time = bws_time(j); - // NOTE : we are doing the bw scaling here - for (int zz : index){ - if(time_events[zz] != v_time){ - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(ii,j,zz) = kvalues(ii,j,zz) + ((kernel_net * kernel_time) * (1.0/(bw_net*bw_time))); - } - } - } - } - } - - - - //lancement des iterations - while(data_holder.empty()==FALSE){ - - //unpacking (imagine some loop unrolling here with a function to deal with.) - acase cas = data_holder.back(); - data_holder.pop_back(); - // int v = cas.v; - // int prev_node = cas.prev_node; - // int depth = cas.depth; - // double d = cas.d; - // double alpha = cas.alpha; - - //Rcout << "v_old is : "<< prev_node <<"\n"; - //Rcout << "v is : "<< v <<"\n"; - //Rcout << "d is : "<< d<<"\n"; - //Rcout << "alpha is : "<< alpha <<"\n"; - //Rcout << "\n\n"; - - // we will update the densities on v - // but only if v is a vertex on wich I can find an event - - std::vector index = get_all_indeces_int(events,cas.v); - - // NOTE : WE SHOULD FIRST CALCULATE THE COMPLETE NKDE DENSITY CAUSED BY AN EVENT - // AND ONLY AFTER MULTIPLYING IT WITH THE TIME DENSITY - if(index.size() >0 ){ - for(int ii = 0; ii < bws_net.n_elem ; ii++){ - bw_net = bws_net(ii); - double kernel_net = kernel_func(cas.d,bw_net) * cas.alpha; - // NOTE : we are not doing the bw scaling here but later - for (int zz : index){ - net_k_values(zz,ii) = net_k_values(zz,ii) + kernel_net; - } - } - } - - // we can now prepare the next steps - if(max_bw_net >= cas.d){ - IntegerVector v_neighbours = neighbour_list[cas.v-1]; - int n = v_neighbours.length(); - int new_depth; - - //updating depth - if(n>2){ - new_depth = cas.depth+1; - }else{ - new_depth = cas.depth; - } - - if(n>1){ - // we must continue only if we are not too deep - if(new_depth <= max_depth){ - // and iterate over the neighbours - for(int j = 0; j < n ; j++){ - // ---- first, the simple case, we are not going back ---- - int v2 = v_neighbours[j]; - int l2 = edge_mat(cas.v,v2); - double d2 = cas.d + line_weights[l2-1]; - // first case, we must back fire - if(v2 == cas.prev_node){ - if(n>2){ - double p2 = (2.0-n)/n; - double new_alpha = cas.alpha * p2; - //data_holder.push_back(new_case); - data_holder.push_back((struct acase){v2,cas.v,new_depth,d2,new_alpha}); - } - }else{ - double new_alpha = cas.alpha * (2.0/n); - //data_holder.push_back(new_case); - data_holder.push_back((struct acase){v2,cas.v,new_depth,d2,new_alpha}); - } - } - } - } - } - } - - //Rcout << "HERE ARE THE CALCULATED NETWORK DENSITIES\n"; - //Rcout << net_k_values <<"\n\n"; - - // LET US CALCULATE THE TEMPORAL DENSITIES - for(int e = 0; e < events.length() ; e++){ - double d = std::abs(v_time - time_events(e)); - for(int t = 0; t < bws_time.n_elem; t++){ - time_k_values(e,t) = kernel_func(d,bws_time(t)); - } - } - //Rcout << "HERE ARE THE CALCULATED TIME DENSITIES\n"; - //Rcout << time_k_values <<"\n\n"; - - - // BEFORE THE SCALING HERE, I MUST APPLY THE TEMPORAL DENSITIES TOO - - for(int e = 0; e < events.length() ; e++){ - for(int n = 0; n < bws_net.n_elem; n++){ - for(int t = 0; t < bws_time.n_elem; t++){ - kvalues(n,t,e) = net_k_values(e,n) * time_k_values(e,t); - } - } - } - - //Rcout << "HERE ARE THE UNSCALED DENSITIES\n"; - //Rcout << kvalues <<"\n\n"; - - // and now we can apply the scaling - arma::mat scale_mat(bws_net.n_elem, bws_time.n_elem); - - for(int i = 0; i < bws_time.n_elem; i++){ - scale_mat.col(i) = 1.0/(bws_net * bws_time(i)); - } - - for(int i = 0; i < events.length(); i++){ - kvalues.slice(i) = kvalues.slice(i) % scale_mat; - } - - - return kvalues; -} - -//' @title The worker function to calculate continuous TNKDE likelihood cv (adaptive case) -//' @name esc_kernel_loo_tnkde_adpt -//' @description The worker function to calculate continuous TNKDE likelihood cv (INTERNAL) -//' @param kernel_func a cpp pointer function (selected with the kernel name) -//' @param edge_mat matrix, to find the id of each edge given two neighbours. -//' @param events a NumericVector indicating the nodes in the graph being events -//' @param time_events a NumericVector indicating the timestamp of each event -//' @param neighbour_list a List, giving for each node an IntegerVector with -//' its neighbours -//' @param v the actual node to consider (int) -//' @param v_time the time of v (double) -//' @param bws_net an arma::mat with the network bandwidths to consider -//' @param bws_time an arma::mat with the time bandwidths to consider -//' @param line_weights a vector with the length of the edges -//' @param depth the actual recursion depth -//' @param max_depth the maximum recursion depth -//' @return a cube with the impact of the event v on each other event for -//' each pair of bandwidths (cube(bws_net, bws_time, events)) -//' @keywords internal - arma::cube esc_kernel_loo_tnkde_adpt(fptros kernel_func, arma::sp_imat &edge_mat, - IntegerVector &events, - IntegerVector &events_wid, - NumericVector &time_events, - List &neighbour_list, - int v, int wid, double v_time, - arma::mat &bws_net, - arma::mat &bws_time, - NumericVector &line_weights, int max_depth){ - //step0 : generate the queue - int depth = 0; - struct acase{ - int v; - int prev_node; - int depth; - double d; - double alpha; - }; - - // prepare the data cube and matrix - int n_bw_net = bws_net.n_rows; - int n_bw_time = bws_time.n_cols; - arma::cube kvalues(bws_net.n_rows, bws_time.n_cols, events.length()); - arma::cube net_k_values(bws_net.n_rows, bws_time.n_cols, events.length()); - arma::cube time_k_values(bws_net.n_rows, bws_time.n_cols, events.length()); - kvalues.fill(0.0); - net_k_values.fill(0.0); - time_k_values.fill(0.0); - - double max_bw_net = bws_net.max(); - double bw_net, bw_time; - - //queue data_holder; - std::vector data_holder; - - // let us prepare the first cases - IntegerVector v_neighbours = neighbour_list[v-1]; - double alpha = 2.0/v_neighbours.length(); - - acase el = {v,-999,0,0.0,alpha}; - data_holder.push_back(el); - - // avant de lancer les iterations, il est important de verifier que - // certains evenement ne se trouvent pas sur la meme vertex d'origine. - // si c'est le cas, nous devons leur assigner la valeur de densite correspondante - std::vector index = get_all_indeces_int(events,v); - if(index.size() > 1){ - // Rcout << "See me ! I must add some density at the origin ! \n"; - for (int zz : index){ - if(time_events[zz] != v_time){ - //arma::mat slice_bw_net = bws_net.slice(zz); - //arma::mat slice_bw_time = bws_time.slice(zz); - //Rcout << "step8\n"; - for(int j = 0 ; j < bws_net.n_rows; j ++){ - for(int jj = 0 ; jj < bws_time.n_cols; jj ++){ - // NOTE : we are NOT doing the bw scaling here - double bw_net = bws_net(j,jj); - double bw_time = bws_time(j,jj); - double kernel_net = kernel_func(0.0, bw_net); - double kernel_time = kernel_func(std::abs(v_time-time_events[zz]), bw_time); - kvalues(j,jj,zz) = kvalues(j,jj,zz) + ((kernel_net * kernel_time)); - } - } - } - } - } - - - //lancement des iterations - while(data_holder.empty()==FALSE){ - - //unpacking (imagine some loop unrolling here with a function to deal with.) - acase cas = data_holder.back(); - data_holder.pop_back(); - int v = cas.v; - int prev_node = cas.prev_node; - int depth = cas.depth; - double d = cas.d; - double alpha = cas.alpha; - - // Rcout << "v_old is : "<< prev_node <<"\n"; - // Rcout << "v is : "<< v <<"\n"; - // Rcout << "d is : "<< d<<"\n"; - // Rcout << "alpha is : "<< alpha <<"\n"; - // Rcout << "\n\n"; - - // we will update the densities on v - // but only if v is a vertex on wich I can find an event - - std::vector index = get_all_indeces_int(events,v); - - // NOTE : WE SHOULD FIRST CALCULATE THE COMPLETE NKDE DENSITY CAUSED BY AN EVENT - // AND ONLY AFTER MULTIPLYING IT WITH THE TIME DENSITY - if(index.size() >0 ){ - for(int j = 0; j < bws_net.n_rows ; j++){ - for(int jj = 0 ; jj < bws_time.n_cols ; jj ++){ - bw_net = bws_net(j,jj); - // NOTE : we are not doing the bw scaling here but later - for (int zz : index){ - double kernel_net = kernel_func(d,bw_net) * alpha; - net_k_values(j,jj,zz) = net_k_values(j,jj,zz) + kernel_net; - } - } - } - } - - // we can now prepare the next steps - if(max_bw_net >= d){ - IntegerVector v_neighbours = neighbour_list[v-1]; - int n = v_neighbours.length(); - int new_depth; - - //updating depth - if(n>2){ - new_depth = depth+1; - }else{ - new_depth = depth; - } - - if(n>1){ - // we must continue only if we are not too deep - if(new_depth <= max_depth){ - // and iterate over the neighbours - for(int j = 0; j < n ; j++){ - // ---- first, the simple case, we are not going back ---- - int v2 = v_neighbours[j]; - int l2 = edge_mat(v,v2); - double d2 = d + line_weights[l2-1]; - // first case, we must back fire - if(v2 == prev_node){ - if(n>2){ - double p2 = (2.0-n)/n; - double new_alpha = alpha * p2; - acase new_case = {v2,v,new_depth,d2,new_alpha}; - //Rcout << " adding this cas : d="<= 0){ - k.slice(index).fill(0); - } - - // and summing the values at each iteration (% is the element wise product) - base_k += k; - }; - - return base_k; - // // and calculate the final values - // arma::mat result; - // - // arma::uvec neg_elems = arma::find(base_k <= 0); - // base_k.elem(neg_elems).fill(min_tol); - // - // - // result = arma::sum(arma::log(base_k),2); - // - // - // return result; -} - - -//' @title The exposed function to calculate TNKDE likelihood cv -//' @name tnkde_get_loo_values2 -//' @description The exposed function to calculate TNKDE likelihood cv (INTERNAL) when an adaptive bandwidth is used -//' @param method a string, one of "simple", "continuous", "discontinuous" -//' @param neighbour_list a List, giving for each node an IntegerVector with -//' its neighbours -//' @param sel_events a Numeric vector indicating the selected events (id of nodes) -//' @param sel_events_wid a Numeric Vector indicating the unique if of the selected events -//' @param sel_events_time a Numeric Vector indicating the time of the selected events -//' @param events a NumericVector indicating the nodes in the graph being events -//' @param events_wid a NumericVector indicating the unique id of all the events -//' @param events_time a NumericVector indicating the timestamp of each event -//' @param weights a cube with the weights associated with each event for each -//' bws_net and bws_time. -//' @param bws_net an arma::cube of three dimensions with the network bandwidths calculated for each observation for each global time and network bandwidths -//' @param bws_time an arma::cube of three dimensions with the time bandwidths calculated for each observation for each global time and network bandwidths -//' @param kernel_name a string with the name of the kernel to use -//' @param line_list a DataFrame describing the lines -//' @param max_depth the maximum recursion depth -//' @param min_tol a double indicating by how much 0 in density values must be replaced -//' @return a matrix with the CV score for each pair of global bandiwdths -//' @export -//' @examples -//' # no example provided, this is an internal function -// [[Rcpp::export]] - arma::cube tnkde_get_loo_values2(std::string method, List &neighbour_list, - IntegerVector &sel_events, - IntegerVector &sel_events_wid, - NumericVector &sel_events_time, - IntegerVector &events, - IntegerVector &events_wid, - NumericVector &events_time, - arma::cube &weights, - arma::cube &bws_net, - arma::cube &bws_time, - std::string kernel_name, - DataFrame &line_list, int max_depth, double min_tol){ - - //selecting the kernel function - //Rcout << "print stage 1 ... \n"; - fptros kernel_func = select_kernelos(kernel_name); - - //step0 extract the columns of the dataframe - NumericVector line_weights = line_list["weight"]; - - //step 1 : mettre toutes les valeurs a 0 (bw_net 0 et bw_time 0) - // NOTE WE calculate the values only for the events in sel_events - arma::cube base_k(bws_net.n_rows, bws_time.n_cols, sel_events.length()); - - //calculer la matrice des lignes - //arma::sp_mat edge_mat = make_matrix_sparse(line_list,neighbour_list); - arma::sp_imat edge_mat = make_imatrix_sparse(line_list, neighbour_list); - arma::cube k; - // Rcout << "Here is the neighbours wids : \n"; - // Rcout << events_wid << "\n\n"; - // Rcout << "Here is the neighbours time : \n"; - // Rcout << events_time << "\n\n"; - - // Rcout << "iterating to calculate the densities from each event \n"; - //step2 : iterer sur chaque event de la zone d'etude - int cnt_e = events.length()-1; - for(int i=0; i <= cnt_e; ++i){ - // Rcout << " "<< i<<"\n"; - //preparer les differentes valeurs de departs pour l'event y - int y = events[i]; - int wid = events_wid[i]; - double v_time = events_time[i]; - // launching recursion - // here we got the the influences of the vertex y on each other selected event in quadra - - if(method == "simple"){ - k = ess_kernel_loo_tnkde_adpt(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net.slice(i), bws_time.slice(i), - line_weights, max_depth); - - }else if (method == "discontinuous"){ - // NOTE TO MYSELF HERE : If we are considering the densities caused by - // the event y, then we only need the bandwidths of the event y - k = esd_kernel_loo_tnkde_adpt(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net.slice(i), bws_time.slice(i), - line_weights, max_depth); - }else{ - // NOTE TO MYSELF HERE : I MUST IMPLEMENT A VERSION OF THIS FUNCTION - // TO WORK WITH ADAPTIVE BANDWIDTH (SEE THE ONE PREPARED ABOVE FOR - // THE SIMPLE METHOD) - k = esc_kernel_loo_tnkde_adpt(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net.slice(i), bws_time.slice(i), - line_weights, max_depth); - } - // NOTE : the scaling by bws is applied above - arma::mat w = weights.slice(i); - - for(int ii = 0; ii < sel_events.length(); ii++){ - k.slice(ii) = k.slice(ii) % w; - } - - // if y was a selected event, its own weight must be set to 0 - int index = get_first_index_int(sel_events_wid,wid); - if(index >= 0){ - k.slice(index).fill(0); - } - - // and summing the values at each iteration (% is the element wise product) - base_k += k; - }; - - return base_k; - - } - - - -//// DEVELOPMENT ///// - - - -//' @title The exposed function to calculate adaptive bandwidth with space-time -//' interaction for TNKDE (INTERNAL) -//' @name adaptive_bw_tnkde_cpp -//' @param method a string, one of "simple", "continuous", "discontinuous" -//' @param neighbour_list a List, giving for each node an IntegerVector with -//' its neighbours -//' @param sel_events a Numeric vector indicating the selected events (id of nodes) -//' @param sel_events_wid a Numeric Vector indicating the unique if of the selected events -//' @param sel_events_time a Numeric Vector indicating the time of the selected events -//' @param events a NumericVector indicating the nodes in the graph being events -//' @param events_wid a NumericVector indicating the unique id of all the events -//' @param events_time a NumericVector indicating the timestamp of each event -//' @param weights a cube with the weights associated with each event for each -//' bws_net and bws_time. -//' @param bws_net an arma::vec with the network bandwidths to consider -//' @param bws_time an arma::vec with the time bandwidths to consider -//' @param kernel_name a string with the name of the kernel to use -//' @param line_list a DataFrame describing the lines -//' @param max_depth the maximum recursion depth -//' @param min_tol a double indicating by how much 0 in density values must be replaced -//' @return a vector witht the estimated density at each event location -//' @export -//' @examples -//' # no example provided, this is an internal function -// [[Rcpp::export]] -arma::rowvec adaptive_bw_tnkde_cpp(std::string method, - List &neighbour_list, - IntegerVector &sel_events, - IntegerVector &sel_events_wid, - NumericVector &sel_events_time, - IntegerVector &events, - IntegerVector &events_wid, - NumericVector &events_time, - arma::vec &weights, - arma::vec &bws_net, - arma::vec &bws_time, - std::string kernel_name, - DataFrame &line_list, - int max_depth, - double min_tol){ - - //selecting the kernel function - fptros kernel_func = select_kernelos(kernel_name); - - //step0 extract the columns of the dataframe - NumericVector line_weights = line_list["weight"]; - - //step 1 : mettre toutes les valeurs a 0 (bw_net 0 et bw_time 0) - // NOTE WE calculate the values only for the events in sel_events - arma::cube base_k(bws_time.n_elem, bws_net.n_elem, sel_events.length()); - - //calculer la matrice des lignes - //arma::sp_mat edge_mat = make_matrix_sparse(line_list,neighbour_list); - arma::sp_imat edge_mat = make_imatrix_sparse(line_list, neighbour_list); - arma::cube k; - //step2 : iterer sur chaque event de la zone d'etude - int cnt_e = events.length()-1; - for(int i=0; i <= cnt_e; ++i){ - //preparer les differentes valeurs de departs pour l'event y - int y = events[i]; - int wid = events_wid[i]; - double v_time = events_time[i]; - // launching recursion - // here we got the the influences of the vertex y on each other selected event in quadra - if(method == "simple"){ - k = ess_kernel_loo_tnkde(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net, bws_time, - line_weights, max_depth); - - }else if (method == "discontinuous"){ - k = esd_kernel_loo_tnkde(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net, bws_time, - line_weights, max_depth); - }else{ - k = esc_kernel_loo_tnkde(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net, bws_time, - line_weights, max_depth); - } - // NOTE : the scaling by bws is applied above - float w = weights(i); - for(int ii = 0; ii < sel_events.length(); ii++){ - k.slice(ii) = k.slice(ii) * w; - } - - // and summing the values at each iteration (% is the element wise product) - base_k += k; - }; - // and calculate the final values - - arma::uvec neg_elems = arma::find(base_k <= 0); - base_k.elem(neg_elems).fill(min_tol); - - arma::rowvec result = base_k(arma::span(0), arma::span(0), arma::span::all); - - return result; -} - - -//' @title The exposed function to calculate adaptive bandwidth with space-time -//' interaction for TNKDE (INTERNAL) -//' @name adaptive_bw_tnkde_cpp2 -//' @param method a string, one of "simple", "continuous", "discontinuous" -//' @param neighbour_list a List, giving for each node an IntegerVector with -//' its neighbours -//' @param sel_events a Numeric vector indicating the selected events (id of nodes) -//' @param sel_events_wid a Numeric Vector indicating the unique if of the selected events -//' @param sel_events_time a Numeric Vector indicating the time of the selected events -//' @param events a NumericVector indicating the nodes in the graph being events -//' @param events_wid a NumericVector indicating the unique id of all the events -//' @param events_time a NumericVector indicating the timestamp of each event -//' @param weights a cube with the weights associated with each event for each -//' bws_net and bws_time. -//' @param bws_net an arma::vec with the network bandwidths to consider -//' @param bws_time an arma::vec with the time bandwidths to consider -//' @param kernel_name a string with the name of the kernel to use -//' @param line_list a DataFrame describing the lines -//' @param max_depth the maximum recursion depth -//' @param min_tol a double indicating by how much 0 in density values must be replaced -//' @return a vector with the estimated density at each event location -//' @export -//' @examples -//' # no example provided, this is an internal function -// [[Rcpp::export]] - arma::cube adaptive_bw_tnkde_cpp2(std::string method, - List &neighbour_list, - IntegerVector &sel_events, - IntegerVector &sel_events_wid, - NumericVector &sel_events_time, - IntegerVector &events, - IntegerVector &events_wid, - NumericVector &events_time, - arma::vec &weights, - arma::vec &bws_net, - arma::vec &bws_time, - std::string kernel_name, - DataFrame &line_list, - int max_depth, - double min_tol){ - // Rcout << "step0\n"; - //selecting the kernel function - fptros kernel_func = select_kernelos(kernel_name); - // Rcout << "step1\n"; - //step0 extract the columns of the dataframe - NumericVector line_weights = line_list["weight"]; - // Rcout << "step2\n"; - //step 1 : mettre toutes les valeurs a 0 (bw_net 0 et bw_time 0) - // NOTE WE calculate the values only for the events in sel_events - arma::cube base_k(bws_net.n_elem, bws_time.n_elem, sel_events.length()); - // Rcout << "step3\n"; - //calculer la matrice des lignes - //arma::sp_mat edge_mat = make_matrix_sparse(line_list,neighbour_list); - arma::sp_imat edge_mat = make_imatrix_sparse(line_list, neighbour_list); - arma::cube k; - // Rcout << "step3\n"; - //step2 : iterer sur chaque event de la zone d'etude - int cnt_e = events.length()-1; - for(int i=0; i <= cnt_e; ++i){ - // Rcout << "Iterating on " << i<<"\n"; - //preparer les differentes valeurs de departs pour l'event y - int y = events[i]; - int wid = events_wid[i]; - double v_time = events_time[i]; - // launching recursion - // here we got the the influences of the vertex y on each other selected event in quadra - if(method == "simple"){ - k = ess_kernel_loo_tnkde(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net, bws_time, - line_weights, max_depth); - - }else if (method == "discontinuous"){ - k = esd_kernel_loo_tnkde(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net, bws_time, - line_weights, max_depth); - }else{ - k = esc_kernel_loo_tnkde(kernel_func, edge_mat, sel_events,sel_events_wid, - sel_events_time, neighbour_list, - y, wid, v_time, - bws_net, bws_time, - line_weights, max_depth); - } - // Rcout << "Here is k : \n"; - // Rcout << k <<"\n\n"; - // NOTE : the scaling by bws is applied above - float w = weights(i); - for(int ii = 0; ii < sel_events.length(); ii++){ - k.slice(ii) = k.slice(ii) * w; - } - // Rcout << "weight applied !\n"; - // and summing the values at each iteration (% is the element wise product) - base_k += k; - }; - // and calculate the final values - // Rcout << "end of iterations !\n"; - arma::uvec neg_elems = arma::find(base_k <= 0); - base_k.elem(neg_elems).fill(min_tol); - - return base_k; - } - - - - diff --git a/src/matrices_functions.cpp b/src/matrices_functions.cpp index 8ae2228e..ec5f5486 100644 --- a/src/matrices_functions.cpp +++ b/src/matrices_functions.cpp @@ -1,17 +1,5 @@ #include "spNetwork.h" -// a simple function to create a vector of values between a start and an end with defined step -std::vector seq_num(double start, double end, double step){ - - std::vector values; - double cumul = 0 - step; - while(cumul+step < end){ - cumul+=step; - values.push_back(cumul); - } - return values; -} - // a simple function to create a vector of values between a start and an end with defined step // [[Rcpp::export]] std::vector seq_num2(double start, double end, double step){ diff --git a/src/nkde_continuous.cpp.bak b/src/nkde_continuous.cpp.bak deleted file mode 100644 index 1226454e..00000000 --- a/src/nkde_continuous.cpp.bak +++ /dev/null @@ -1,1238 +0,0 @@ -#include "spNetwork.h" -#include "base_kernel_funtions.h" -#include "matrices_functions.h" - - -//##################################################################################### -// ###################### THE WORKER FUNCTIONS ###################################### -//##################################################################################### - - -//' @title The worker function to calculate continuous NKDE (with ARMADILLO and sparse matrix) -//' @name continuousWorker_sparse -//' @param kernel_func a cpp pointer function (selected with the kernel name) -//' @param samples_k a numeric vector of the actual kernel values, updates at -//' each recursion -//' @param neighbour_list a List, giving for each node an IntegerVector with -//' its neighbours -//' @param edge_mat matrix, to find the id of each edge given two neighbours. -//' @param v the actual node to consider for the recursion (int) -//' @param bw the kernel bandwidth -//' @param line_weights a vector with the length of the edges -//' @param samples_edgeid a vector associating each sample to an edge -//' @param samples_coords a matrix with the X and Y coordinates of the samples -//' @param nodes_coords a matrix with the X and Y coordinates of the nodes -//' @param depth the actual recursion depth -//' @param max_depth the maximum recursion depth -//' @return a vector with the kernel values calculated for each samples from -//' the first node given -//' @keywords internal -arma::vec esc_kernel_rcpp_arma_sparse(fptr kernel_func, List &neighbour_list, arma::sp_imat &edge_mat, int v, double bw, NumericVector &line_weights, arma::ivec &samples_edgeid, arma::mat &samples_coords, arma::mat &nodes_coords, int max_depth){ - - // let me create the vector of kernel values - arma::vec samples_k(samples_edgeid.n_elem); - //samples_k.fill(0.0); - - // instead of a recursion, we will iterate over cases stored in a queue - // a case will be an edge for which we want to recalculate the density - // of events on its - struct acase{ - int v1; - int v2; - int l; - int depth; - double d; - double alpha; - }; - - //queue data_holder; - std::vector data_holder; - - // let us prepare the first cases - IntegerVector v_neighbours = neighbour_list[v-1]; - double alpha = 2.0/v_neighbours.length(); - for(int j = 0; j < v_neighbours.length(); j++){ - int v2 = v_neighbours[j]; - int l = edge_mat(v,v2); - acase el = {v,v2,l,0,0.0,alpha}; - //Rcout << "Adding this cas : d="<<0.0<<", v1="<= cas.d){ - - // noice, now we can check the reachable nodes and create new cases - IntegerVector v2_neighbours = neighbour_list[cas.v2-1]; - int n = v2_neighbours.length(); - - int new_depth; - //updating depth - if(n>2){ - new_depth = cas.depth+1; - }else{ - new_depth = cas.depth; - } - - if(n>1){ - // we must continue only if we are not too deep - if(new_depth <= max_depth){ - // and iterate over the neighbours - for(int j = 0; j < n ; j++){ - - // ---- first, the simple case, we are not going back ---- - int v3 = v2_neighbours[j]; - // int l2 = edge_mat(cas.v2,v3); - // double d2 = cas.d + line_weights[cas.l-1]; - // first case, we must back fire - if(v3 == cas.v1){ - if(n>2){ - //acase new_case = {cas.v2,v3,l2,new_depth,d2,(cas.alpha * ((2.0-n)/n))}; - //Rcout << " adding this cas : d="< data_holder; -// std::vector data_holder; -// -// // let us prepare the first cases -// IntegerVector v_neighbours = neighbour_list[v-1]; -// double alpha = 2.0/v_neighbours.length(); -// for(int j = 0; j < v_neighbours.length(); j++){ -// int v2 = v_neighbours[j]; -// int l = edge_mat(v,v2); -// acase el = {v,v2,l,0,0.0,alpha}; -// //Rcout << "Adding this cas : d="<<0.0<<", v1="<= cas.d){ -// -// // noice, now we can check the reachable nodes and create new cases -// IntegerVector v2_neighbours = neighbour_list[cas.v2-1]; -// int n = v2_neighbours.length(); -// -// int new_depth; -// //updating depth -// if(n>2){ -// new_depth = cas.depth+1; -// }else{ -// new_depth = cas.depth; -// } -// -// if(n>1){ -// // we must continue only if we are not too deep -// if(new_depth <= max_depth){ -// // and iterate over the neighbours -// for(int j = 0; j < n ; j++){ -// -// // ---- first, the simple case, we are not going back ---- -// int v3 = v2_neighbours[j]; -// int l2 = edge_mat(cas.v2,v3); -// double d2 = cas.d + line_weights[cas.l-1]; -// // first case, we must back fire -// if(v3 == cas.v1){ -// if(n>2){ -// //acase new_case = {cas.v2,v3,l2,new_depth,d2,(cas.alpha * ((2.0-n)/n))}; -// //Rcout << " adding this cas : d="< data_holder; -// std::vector data_holder; -// -// // let us prepare the first cases -// IntegerVector v_neighbours = neighbour_list[v-1]; -// double alpha = 2.0/v_neighbours.length(); -// for(int j = 0; j < v_neighbours.length(); j++){ -// int v2 = v_neighbours[j]; -// int l = edge_mat(v,v2); -// acase el = {v,v2,l,0,0.0,alpha}; -// //Rcout << "Adding this cas : d="<<0.0<<", v1="<= d){ -// -// // noice, now we can check the reachable nodes and create new cases -// IntegerVector v2_neighbours = neighbour_list[v2-1]; -// int n = v2_neighbours.length(); -// -// int new_depth; -// //updating depth -// if(n>2){ -// new_depth = depth+1; -// }else{ -// new_depth = depth; -// } -// -// if(n>1){ -// // we must continue only if we are not too deep -// if(new_depth <= max_depth){ -// // and iterate over the neighbours -// for(int j = 0; j < n ; j++){ -// -// // ---- first, the simple case, we are not going back ---- -// int v3 = v2_neighbours[j]; -// int l2 = edge_mat(v2,v3); -// double d2 = d + line_weights[l-1]; -// // first case, we must back fire -// if(v3 == v1){ -// if(n>2){ -// double p2 = (2.0-n)/n; -// double new_alpha = alpha * p2; -// acase new_case = {v2,v3,l2,new_depth,d2,new_alpha}; -// //Rcout << " adding this cas : d="<d2) && (depth < max_depth)){ -// //on veut trouver toutes les lignes emannant de v (Lv) -// IntegerVector v_neighbours = neighbour_list[v1-1]; -// int n = v_neighbours.length(); -// int new_depth; -// //updating depth -// if(n>2){ -// new_depth = depth+1; -// }else{ -// new_depth = depth+0; -// } -// if(n>1){ -// for(int j=0; j < n; ++j){ -// //int li = Lv[j]; -// int vi = v_neighbours[j]; -// //int li = edge_dict[v1][vi]; -// int li = edge_mat(v1,vi); -// if(li==l1){ -// //we must backfire only if we have a true intersection -// if(n>2){ -// double p2 = (2.0-n)/n; -// double n_alpha = alpha * p2; -// samples_k = esc_kernel_rcpp_arma_sparse(kernel_func, samples_k, neighbour_list, edge_mat, v1,vi, li, d2, n_alpha, bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, new_depth, max_depth); -// } -// }else{ -// double n_alpha = alpha * (2.0/n); -// samples_k = esc_kernel_rcpp_arma_sparse(kernel_func, samples_k, neighbour_list, edge_mat, v1,vi, li, d2, n_alpha, bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, new_depth, max_depth); -// }; -// }; -// }; -// }; -// return samples_k; -// } - - - -//' @title The worker function to calculate continuous NKDE (with ARMADILLO and integer matrix) -//' @name continuousWorker -//' @param kernel_func a cpp pointer function (selected with the kernel name) -//' @param samples_k a numeric vector of the actual kernel values, updates at -//' each recursion -//' @param neighbour_list a List, giving for each node an IntegerVector with -//' its neighbours -//' @param edge_mat matrix, to find the id of each edge given two neighbours. -//' @param v the actual node to consider for the recursion (int) -//' @param bw the kernel bandwidth -//' @param line_weights a vector with the length of the edges -//' @param samples_edgeid a vector associating each sample to an edge -//' @param samples_x a vector with x coordinates of each sample -//' @param samples_y a vector with y coordinates of each sample -//' @param nodes_x a vector with x coordinates of each node -//' @param nodes_y a vector with y coordinates of each node -//' @param depth the actual recursion depth -//' @param max_depth the maximum recursion depth -//' @return a vector with the kernel values calculated for each samples from -//' the first node given -//' @keywords internal -arma::vec esc_kernel_rcpp_arma(fptr kernel_func, List &neighbour_list, IntegerMatrix &edge_mat, int v, double bw, NumericVector &line_weights, arma::ivec &samples_edgeid, arma::vec &samples_x, arma::vec &samples_y, arma::vec &nodes_x, arma::vec &nodes_y, int max_depth){ - - // let me create the vector of kernel values - arma::vec samples_k(samples_edgeid.n_elem); - //samples_k.fill(0.0); - - // instead of a recursion, we will iterate over cases stored in a queue - // a case will be an edge for which we want to recalculate the density - // of events on its - struct acase{ - int v1; - int v2; - int l; - int depth; - double d; - double alpha; - }; - - //queue data_holder; - std::vector data_holder; - - // let us prepare the first cases - IntegerVector v_neighbours = neighbour_list[v-1]; - double alpha = 2.0/v_neighbours.length(); - for(int j = 0; j < v_neighbours.length(); j++){ - int v2 = v_neighbours[j]; - int l = edge_mat(v,v2); - acase el = {v,v2,l,0,0.0,alpha}; - //Rcout << "Adding this cas : d="<<0.0<<", v1="< d){ - - // noice, now we can check the reachable nodes and create new cases - IntegerVector v2_neighbours = neighbour_list[v2-1]; - int n = v2_neighbours.length(); - - int new_depth; - //updating depth - if(n>2){ - new_depth = depth+1; - }else{ - new_depth = depth; - } - - if(n>1){ - // we must continue only if we are not too deep - if(new_depth <= max_depth){ - // and iterate over the neighbours - for(int j = 0; j < n ; j++){ - - // ---- first, the simple case, we are not going back ---- - int v3 = v2_neighbours[j]; - int l2 = edge_mat(v2,v3); - double d2 = d + line_weights[l-1]; - // first case, we must back fire - if(v3 == v1){ - if(n>2){ - double new_alpha = alpha * ((2.0-n)/n); - acase new_case = {v2,v3,l2,new_depth,d2,new_alpha}; - //Rcout << " adding this cas : d="<d2) && (depth < max_depth)){ -// //on veut trouver toutes les lignes emannant de v (Lv) -// IntegerVector v_neighbours = neighbour_list[v1-1]; -// int n = v_neighbours.length(); -// int new_depth; -// //updating depth -// if(n>2){ -// new_depth = depth+1; -// }else{ -// new_depth = depth+0; -// } -// if(n>1){ -// for(int j=0; j < n; ++j){ -// //int li = Lv[j]; -// int vi = v_neighbours[j]; -// //int li = edge_dict[v1][vi]; -// int li = edge_mat(v1,vi); -// if(li==l1){ -// //we must backfire only if we have a true intersection -// if(n>2){ -// double p2 = (2.0-n)/n; -// double n_alpha = alpha * p2; -// samples_k = esc_kernel_rcpp_arma(kernel_func, samples_k, neighbour_list, edge_mat, v1,vi, li, d2, n_alpha, bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, new_depth, max_depth); -// } -// }else{ -// double n_alpha = alpha * (2.0/n); -// samples_k = esc_kernel_rcpp_arma(kernel_func, samples_k, neighbour_list, edge_mat, v1,vi, li, d2, n_alpha, bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, new_depth, max_depth); -// }; -// }; -// }; -// }; -// return samples_k; -// } - - - - -//##################################################################################### -// ###################### THE CALLED FUNCTIONS ###################################### -//##################################################################################### - - -//' @title The main function to calculate continuous NKDE (with ARMADILO and sparse matrix) -//' @name continuousfunction2 -//' @param neighbour_list a list of the neighbours of each node -//' @param events a numeric vector of the node id of each event -//' @param weights a numeric vector of the weight of each event -//' @param samples a DataFrame of the samples (with spatial coordinates and belonging edge) -//' @param bws the kernel bandwidths for each event -//' @param kernel_name the name of the kernel to use -//' @param nodes a DataFrame representing the nodes of the graph (with spatial coordinates) -//' @param line_list a DataFrame representing the lines of the graph -//' @param mat_dist_samples a matrix with the following columns : sampleid, lineid, nodeid, dist. It will be used to not recalculate euclidean distance too much -//' @param max_depth the maximum recursion depth (after which recursion is stopped) -//' @param verbose a boolean indicating if the function must print its progress -//' @param div The divisor to use for the kernel. Must be "n" (the number of events within the radius around each sampling point), "bw" (the bandwidth) "none" (the simple sum). -//' @return a DataFrame with two columns : the kernel values (sum_k) and the number of events for each sample (n) -//' @export -//' @keywords internal -//' -// [[Rcpp::export]] - DataFrame continuous_nkde_cpp_arma_sparse(List neighbour_list, NumericVector events, - NumericVector weights, DataFrame samples, - NumericVector bws, std::string kernel_name, - DataFrame nodes, DataFrame line_list, - int max_depth, bool verbose, std::string div = "bw"){ - - //continuous_nkde_cpp_arma_sparse(neighbour_list,events$vertex_id, events$weight, samples@data, bws, kernel_name, nodes@data, graph_result$linelist, max_depth, tol, verbose) - //selecting the kernel function - fptr kernel_func = select_kernel(kernel_name); - - //step0 extract the columns of the dataframe - NumericVector line_weights = line_list["weight"]; - arma::ivec samples_edgeid = as(samples["edge_id"]); - - arma::mat samples_coords(samples.nrows(),2); - samples_coords.col(0) = as(samples["X_coords"]); - samples_coords.col(1) = as(samples["Y_coords"]); - - arma::mat nodes_coords(nodes.nrows(),2); - nodes_coords.col(0) = as(nodes["X_coords"]); - nodes_coords.col(1) = as(nodes["Y_coords"]); - - - //step 1 : mettre toutes les valeurs a 0 - arma::vec base_k(samples.nrow()); - arma::vec base_count(samples.nrow()); - arma::vec samples_k(samples.nrow()); - arma::vec count(samples.nrow()); - //base_k.fill(0.0); - //base_count.fill(0.0); - //count.fill(0.0); - //samples_k.fill(0.0); - - //calculer le dictionnaire des lignes - //IntegerMatrix edge_mat = make_matrix(line_list,neighbour_list); - arma::sp_imat edge_mat = make_imatrix_sparse(line_list, neighbour_list); - //step2 : iterer sur chaque event - int cnt_e = events.length()-1; - Progress p(cnt_e, verbose); - for(int i=0; i <= cnt_e; ++i){ - p.increment(); // update progress - //preparer les differentes valeurs de departs pour l'event y - int y = events[i]; - double w = weights[i]; - double bw = bws[i]; - - // on peut maintenant calculer la densite emanant de l'event y - samples_k = esc_kernel_rcpp_arma_sparse(kernel_func, neighbour_list, edge_mat, y, bw, line_weights, samples_edgeid, samples_coords,nodes_coords, max_depth); - if(div == "bw"){ - base_k += (samples_k*w) / bw; - }else{ - base_k += samples_k*w; - } - // calculating the new value - count.fill(0.0); - arma::uvec ids = arma::find(samples_k>0); - count.elem(ids).fill(1.0); - base_count += count*w; - - }; - NumericVector v1 = Rcpp::NumericVector(base_k.begin(), base_k.end()); - NumericVector v2 = Rcpp::NumericVector(base_count.begin(), base_count.end()); - DataFrame df = DataFrame::create( Named("sum_k") = v1 ,Named("n") = v2); - return df; - } - - - - -// DataFrame continuous_nkde_cpp_arma_sparse_old(List neighbour_list, NumericVector events, -// NumericVector weights, DataFrame samples, -// NumericVector bws, std::string kernel_name, -// DataFrame nodes, DataFrame line_list, -// int max_depth, bool verbose, std::string div = "bw"){ -// -// //continuous_nkde_cpp_arma_sparse(neighbour_list,events$vertex_id, events$weight, samples@data, bws, kernel_name, nodes@data, graph_result$linelist, max_depth, tol, verbose) -// //selecting the kernel function -// fptr kernel_func = select_kernel(kernel_name); -// -// //step0 extract the columns of the dataframe -// NumericVector line_weights = line_list["weight"]; -// arma::ivec samples_edgeid = as(samples["edge_id"]); -// arma::vec samples_x = as(samples["X_coords"]); -// arma::vec samples_y = as(samples["Y_coords"]); -// arma::vec nodes_x = as(nodes["X_coords"]); -// arma::vec nodes_y = as(nodes["Y_coords"]); -// -// //step 1 : mettre toutes les valeurs a 0 -// arma::vec base_k(samples.nrow()); -// arma::vec base_count(samples.nrow()); -// arma::vec samples_k(samples.nrow()); -// arma::vec count(samples.nrow()); -// //base_k.fill(0.0); -// //base_count.fill(0.0); -// //count.fill(0.0); -// //samples_k.fill(0.0); -// -// //calculer le dictionnaire des lignes -// //IntegerMatrix edge_mat = make_matrix(line_list,neighbour_list); -// arma::sp_mat edge_mat = make_matrix_sparse(line_list, neighbour_list); -// //step2 : iterer sur chaque event -// int cnt_e = events.length()-1; -// Progress p(cnt_e, verbose); -// for(int i=0; i <= cnt_e; ++i){ -// p.increment(); // update progress -// //preparer les differentes valeurs de departs pour l'event y -// int y = events[i]; -// double w = weights[i]; -// double bw = bws[i]; -// -// // on peut maintenant calculer la densite emanant de l'event y -// samples_k = esc_kernel_rcpp_arma_sparse(kernel_func, neighbour_list, edge_mat, y, bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, max_depth); -// if(div == "bw"){ -// base_k += (samples_k*w) / bw; -// }else{ -// base_k += samples_k*w; -// } -// // calculating the new value -// count.fill(0.0); -// arma::uvec ids = arma::find(samples_k>0); -// count.elem(ids).fill(1.0); -// base_count += count*w; -// -// }; -// NumericVector v1 = Rcpp::NumericVector(base_k.begin(), base_k.end()); -// NumericVector v2 = Rcpp::NumericVector(base_count.begin(), base_count.end()); -// DataFrame df = DataFrame::create( Named("sum_k") = v1 ,Named("n") = v2); -// return df; -// } - - -// THIS IS THE PREVIOUS VERSION USING RECURSION -// DataFrame continuous_nkde_cpp_arma_sparse(List neighbour_list, NumericVector events, NumericVector weights, DataFrame samples, NumericVector bws, std::string kernel_name, DataFrame nodes, DataFrame line_list, int max_depth, bool verbose){ -// -// //continuous_nkde_cpp_arma_sparse(neighbour_list,events$vertex_id, events$weight, samples@data, bws, kernel_name, nodes@data, graph_result$linelist, max_depth, tol, verbose) -// //selecting the kernel function -// fptr kernel_func = select_kernel(kernel_name); -// -// //step0 extract the columns of the dataframe -// NumericVector line_weights = line_list["weight"]; -// arma::vec samples_edgeid = as(samples["edge_id"]); -// arma::vec samples_x = as(samples["X_coords"]); -// arma::vec samples_y = as(samples["Y_coords"]); -// arma::vec nodes_x = as(nodes["X_coords"]); -// arma::vec nodes_y = as(nodes["Y_coords"]); -// -// //step 1 : mettre toutes les valeurs a 0 -// arma::vec base_k(samples.nrow()); -// arma::vec base_count(samples.nrow()); -// arma::vec samples_k(samples.nrow()); -// arma::vec count(samples.nrow()); -// base_k.fill(0.0); -// base_count.fill(0.0); -// count.fill(0.0); -// samples_k.fill(0.0); -// -// //calculer le dictionnaire des lignes -// //IntegerMatrix edge_mat = make_matrix(line_list,neighbour_list); -// arma::sp_mat edge_mat = make_matrix_sparse(line_list, neighbour_list); -// //step2 : iterer sur chaque event -// int cnt_e = events.length()-1; -// Progress p(cnt_e, verbose); -// for(int i=0; i <= cnt_e; ++i){ -// p.increment(); // update progress -// //preparer les differentes valeurs de departs pour l'event y -// int y = events[i]; -// double w = weights[i]; -// double bw = bws[i]; -// //on veut trouver toutes les voisins emannant de y -// IntegerVector y_neighbours = neighbour_list[y-1]; -// int cnt_y = y_neighbours.length(); -// double alpha = 2.0/cnt_y; -// for(int j=0; j < cnt_y; ++j){ -// //preparing all values for this loop -// int vi = y_neighbours[j]; -// //int li = edge_dict[y][vi]; -// int li = edge_mat(y,vi); -// double d = 0.0 ; -// int depth = 0 ; -// // launching recursion -// samples_k.fill(0.0); -// samples_k = esc_kernel_rcpp_arma_sparse(kernel_func, samples_k, neighbour_list, edge_mat ,y,vi,li,d,alpha,bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, depth,max_depth); -// base_k += samples_k*w; -// // calculating the new value -// count.fill(0.0); -// arma::uvec ids = arma::find(samples_k>0); -// count.elem(ids).fill(1.0); -// base_count += count*w; -// }; -// -// }; -// NumericVector v1 = Rcpp::NumericVector(base_k.begin(), base_k.end()); -// NumericVector v2 = Rcpp::NumericVector(base_count.begin(), base_count.end()); -// DataFrame df = DataFrame::create( Named("sum_k") = v1 ,Named("n") = v2); -// return df; -// } - - -//' @title The main function to calculate continuous NKDE (with ARMADILO and integer matrix) -//' @name continuousfunction -//' @param neighbour_list a list of the neighbours of each node -//' @param events a numeric vector of the node id of each event -//' @param weights a numeric vector of the weight of each event -//' @param samples a DataFrame of the samples (with spatial coordinates and belonging edge) -//' @param bws the kernel bandwidths for each event -//' @param kernel_name the name of the kernel to use -//' @param nodes a DataFrame representing the nodes of the graph (with spatial coordinates) -//' @param line_list a DataFrame representing the lines of the graph -//' @param max_depth the maximum recursion depth (after which recursion is stopped) -//' @param verbose a boolean indicating if the function must print its progress -//' @param div The divisor to use for the kernel. Must be "n" (the number of events within the radius around each sampling point), "bw" (the bandwidth) "none" (the simple sum). -//' @return a DataFrame with two columns : the kernel values (sum_k) and the number of events for each sample (n) -//' @export -//' @keywords internal -//' -// [[Rcpp::export]] -DataFrame continuous_nkde_cpp_arma(List neighbour_list, NumericVector events, - NumericVector weights, DataFrame samples, - NumericVector bws, std::string kernel_name, - DataFrame nodes, DataFrame line_list, - int max_depth, bool verbose, std::string div = "bw"){ - - //continuous_nkde_cpp_arma_sparse(neighbour_list,events$vertex_id, events$weight, samples@data, bws, kernel_name, nodes@data, graph_result$linelist, max_depth, tol, verbose) - //selecting the kernel function - fptr kernel_func = select_kernel(kernel_name); - - //step0 extract the columns of the dataframe - NumericVector line_weights = line_list["weight"]; - arma::ivec samples_edgeid = as(samples["edge_id"]); - arma::vec samples_x = as(samples["X_coords"]); - arma::vec samples_y = as(samples["Y_coords"]); - arma::vec nodes_x = as(nodes["X_coords"]); - arma::vec nodes_y = as(nodes["Y_coords"]); - - //step 1 : mettre toutes les valeurs a 0 - arma::vec base_k(samples.nrow()); - arma::vec base_count(samples.nrow()); - arma::vec samples_k(samples.nrow()); - arma::vec count(samples.nrow()); - // base_k.fill(0.0); - // base_count.fill(0.0); - // count.fill(0.0); - // samples_k.fill(0.0); - - //calculer le dictionnaire des lignes - //IntegerMatrix edge_mat = make_matrix(line_list,neighbour_list); - IntegerMatrix edge_mat = make_matrix(line_list, neighbour_list); - //step2 : iterer sur chaque event - int cnt_e = events.length()-1; - Progress p(cnt_e, verbose); - for(int i=0; i <= cnt_e; ++i){ - p.increment(); // update progress - //preparer les differentes valeurs de departs pour l'event y - int y = events[i]; - double w = weights[i]; - double bw = bws[i]; - - // on peut maintenant calculer la densite emanant de l'event y - samples_k = esc_kernel_rcpp_arma(kernel_func, neighbour_list, edge_mat ,y,bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, max_depth); - if(div == "bw"){ - base_k += (samples_k*w) / bw; - }else{ - base_k += samples_k*w; - } - // calculating the new value - count.fill(0.0); - arma::uvec ids = arma::find(samples_k>0); - count.elem(ids).fill(1.0); - base_count += count*w; - - }; - NumericVector v1 = Rcpp::NumericVector(base_k.begin(), base_k.end()); - NumericVector v2 = Rcpp::NumericVector(base_count.begin(), base_count.end()); - DataFrame df = DataFrame::create( Named("sum_k") = v1 ,Named("n") = v2); - return df; -} - - -// THIS IS THE PREVIOUS VERSION USING RECURSION -// DataFrame continuous_nkde_cpp_arma(List neighbour_list, NumericVector events, NumericVector weights, DataFrame samples, NumericVector bws, std::string kernel_name, DataFrame nodes, DataFrame line_list, int max_depth, bool verbose){ -// -// //selecting the kernel function -// fptr kernel_func = select_kernel(kernel_name); -// -// //step0 extract the columns of the dataframe -// NumericVector line_weights = line_list["weight"]; -// arma::vec samples_edgeid = as(samples["edge_id"]); -// arma::vec samples_x = as(samples["X_coords"]); -// arma::vec samples_y = as(samples["Y_coords"]); -// arma::vec nodes_x = as(nodes["X_coords"]); -// arma::vec nodes_y = as(nodes["Y_coords"]); -// -// //step 1 : mettre toutes les valeurs a 0 -// arma::vec base_k(samples.nrow()); -// arma::vec base_count(samples.nrow()); -// arma::vec samples_k(samples.nrow()); -// arma::vec count(samples.nrow()); -// base_k.fill(0.0); -// base_count.fill(0.0); -// count.fill(0.0); -// samples_k.fill(0.0); -// -// //calculer le dictionnaire des lignes -// IntegerMatrix edge_mat = make_matrix(line_list,neighbour_list); -// //step2 : iterer sur chaque event -// int cnt_e = events.length()-1; -// Progress p(cnt_e, verbose); -// for(int i=0; i <= cnt_e; ++i){ -// p.increment(); // update progress -// //preparer les differentes valeurs de departs pour l'event y -// int y = events[i]; -// double w = weights[i]; -// double bw = bws[i]; -// //on veut trouver toutes les voisins emannant de y -// IntegerVector y_neighbours = neighbour_list[y-1]; -// int cnt_y = y_neighbours.length(); -// double alpha = 2.0/cnt_y; -// for(int j=0; j < cnt_y; ++j){ -// //preparing all values for this loop -// int vi = y_neighbours[j]; -// //int li = edge_dict[y][vi]; -// int li = edge_mat(y,vi); -// double d = 0.0 ; -// int depth = 0 ; -// // launching recursion -// samples_k.fill(0.0); -// samples_k = esc_kernel_rcpp_arma(kernel_func, samples_k, neighbour_list, edge_mat ,y,vi,li,d,alpha,bw, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, depth,max_depth); -// base_k += samples_k*w; -// // calculating the new value -// count.fill(0.0); -// arma::uvec ids = arma::find(samples_k>0); -// count.elem(ids).fill(1.0); -// base_count += count*w; -// }; -// -// }; -// NumericVector v1 = Rcpp::NumericVector(base_k.begin(), base_k.end()); -// NumericVector v2 = Rcpp::NumericVector(base_count.begin(), base_count.end()); -// DataFrame df = DataFrame::create( Named("sum_k") = v1 ,Named("n") = v2); -// return df; -// } - - -//##################################################################################### -// ################# THE EXECUTION FUNCTIONS FOR TNKDE ############################## -//##################################################################################### - -//' @title The main function to calculate continuous TNKDE (with ARMADILO and sparse matrix) -//' @name tnkdecontinuousfunction -//' @param neighbour_list a list of the neighbours of each node -//' @param events a numeric vector of the node id of each event -//' @param events_time a numeric vector with the time for the events -//' @param weights a numeric vector of the weight of each event -//' @param samples a DataFrame of the samples (with spatial coordinates and belonging edge) -//' @param samples_time a NumericVector indicating when to do the samples -//' @param bws_net the network kernel bandwidths for each event -//' @param bws_time the time kernel bandwidths for each event -//' @param kernel_name the name of the kernel to use -//' @param nodes a DataFrame representing the nodes of the graph (with spatial coordinates) -//' @param line_list a DataFrame representing the lines of the graph -//' @param max_depth the maximum recursion depth (after which recursion is stopped) -//' @param verbose a boolean indicating if the function must print its progress -//' @param div a string indicating how to standardize the kernel values -//' @return a List with two matrices: the kernel values (sum_k) and the number of events for each sample (n) -//' @export -//' @keywords internal -//' -// [[Rcpp::export]] -List continuous_tnkde_cpp_arma_sparse(List neighbour_list, - IntegerVector events, NumericVector events_time,NumericVector weights, - DataFrame samples, arma::vec samples_time, - NumericVector bws_net, - NumericVector bws_time, - std::string kernel_name, DataFrame nodes, DataFrame line_list, - int max_depth, bool verbose, std::string div){ - - //selecting the kernel function - fptr kernel_func = select_kernel(kernel_name); - - //step0 extract the columns of the dataframe - NumericVector line_weights = line_list["weight"]; - arma::ivec samples_edgeid = as(samples["edge_id"]); - - arma::mat samples_coords(samples.nrows(),2); - samples_coords.col(0) = as(samples["X_coords"]); - samples_coords.col(1) = as(samples["Y_coords"]); - - arma::mat nodes_coords(nodes.nrows(),2); - nodes_coords.col(0) = as(nodes["X_coords"]); - nodes_coords.col(1) = as(nodes["Y_coords"]); - - //step 1 : set all values to 0 - // NOTE : we will produce matrices here (tnkde) - arma::mat base_k(samples.nrows(), samples_time.n_elem); - arma::mat base_count(samples.nrows(), samples_time.n_elem); - arma::mat count(samples.nrows(), samples_time.n_elem); - arma::vec samples_k(samples.nrows()); - - // NOTE It is possible that we must recalculate the density for the same vertex - // we will store them instead - std::map event_counter = count_values_intvec(events); - std::map saved_values; - - //calculer le dictionnaire des lignes - //IntegerMatrix edge_mat = make_matrix(line_list,neighbour_list); - arma::sp_imat edge_mat = make_imatrix_sparse(line_list, neighbour_list); - - //step2 : iterer sur chaque event - int cnt_e = events.length()-1; - Progress p(cnt_e, verbose); - for(int i=0; i <= cnt_e; ++i){ - p.increment(); // update progress - //preparer les differentes valeurs de departs pour l'event y - int y = events[i]; - double w = weights[i]; - double bw_net = bws_net[i]; - double bw_time = bws_time[i]; - double et = events_time[i]; - - // calculating the density value on the network - if(event_counter[y] > 1){ - // here we have a dupplicated event, we must either use the savec nkde or calculate it and save it latter - if(map_contains_key(saved_values,y)){ - // we already have the value ! - samples_k = saved_values[y]; - }else{ - // we have not yet the value - samples_k = esc_kernel_rcpp_arma_sparse(kernel_func, neighbour_list, edge_mat, y, bw_net, line_weights, samples_edgeid, samples_coords, nodes_coords, max_depth); - saved_values[y] = samples_k; - } - }else{ - // this is not a duplicate event, no need to store the value in memory - samples_k = esc_kernel_rcpp_arma_sparse(kernel_func, neighbour_list, edge_mat, y, bw_net, line_weights, samples_edgeid, samples_coords, nodes_coords, max_depth); - } - - // ok, now we must calculate the temporal densities - arma::vec temporal_density = kernel_func(arma::abs(samples_time - et), bw_time); - - // and standardize by bw if required - if(div == "bw"){ - temporal_density = temporal_density / bw_time; - samples_k = samples_k / bw_net; - } - - // and create an awesome spatio-temporal matrix - arma::mat k_mat(samples.nrows(), samples_time.n_elem); - - for(int j = 0; j < temporal_density.n_elem; j++){ - k_mat.col(j) = samples_k * temporal_density[j]; - } - - base_k += k_mat*w; - // calculating the new value - count.fill(0.0); - arma::umat ids = arma::find(samples_k>0); - count.elem(ids).fill(1.0); - base_count += count*w; - - }; - - // standardise the result if div = n - if(div == "n"){ - base_k = base_k / base_count; - } - - List results = List::create( - Named("sum_k") = base_k, - Named("n") = base_count ); - - return results; -} - - - -//' @title The main function to calculate continuous TNKDE (with ARMADILO and integer matrix) -//' @name tnkdecontinuousfunction -//' @param neighbour_list a list of the neighbours of each node -//' @param events a numeric vector of the node id of each event -//' @param events_time a numeric vector with the time for the events -//' @param weights a numeric vector of the weight of each event -//' @param samples a DataFrame of the samples (with spatial coordinates and belonging edge) -//' @param samples_time a NumericVector indicating when to do the samples -//' @param bws_net the network kernel bandwidths for each event -//' @param kernel_name the name of the kernel to use -//' @param nodes a DataFrame representing the nodes of the graph (with spatial coordinates) -//' @param line_list a DataFrame representing the lines of the graph -//' @param max_depth the maximum recursion depth (after which recursion is stopped) -//' @param verbose a boolean indicating if the function must print its progress -//' @param div a string indicating how to standardize the kernel values -//' @return a List with two matrices: the kernel values (sum_k) and the number of events for each sample (n) -//' @export -//' @keywords internal -//' -// [[Rcpp::export]] -List continuous_tnkde_cpp_arma(List neighbour_list, - IntegerVector events, NumericVector events_time,NumericVector weights, - DataFrame samples, arma::vec samples_time, - NumericVector bws_net, - NumericVector bws_time, - std::string kernel_name, DataFrame nodes, DataFrame line_list, - int max_depth, bool verbose, std::string div){ - - //selecting the kernel function - fptr kernel_func = select_kernel(kernel_name); - - //step0 extract the columns of the dataframe - NumericVector line_weights = line_list["weight"]; - arma::ivec samples_edgeid = as(samples["edge_id"]); - arma::vec samples_x = as(samples["X_coords"]); - arma::vec samples_y = as(samples["Y_coords"]); - arma::vec nodes_x = as(nodes["X_coords"]); - arma::vec nodes_y = as(nodes["Y_coords"]); - - //step 1 : set all values to 0 - // NOTE : we will produce matrices here (tnkde) - arma::mat base_k(samples_x.n_elem, samples_time.n_elem); - arma::mat base_count(samples_x.n_elem, samples_time.n_elem); - arma::mat count(samples_x.n_elem, samples_time.n_elem); - arma::vec samples_k(samples_x.n_elem); - - // NOTE It is possible that we must recalculate the density for the same vertex - // we will store them instead - std::map event_counter = count_values_intvec(events); - std::map saved_values; - - //calculer le dictionnaire des lignes - IntegerMatrix edge_mat = make_matrix(line_list,neighbour_list); - - //step2 : iterer sur chaque event - int cnt_e = events.length()-1; - Progress p(cnt_e, verbose); - for(int i=0; i <= cnt_e; ++i){ - p.increment(); // update progress - //preparer les differentes valeurs de departs pour l'event y - int y = events[i]; - double w = weights[i]; - double bw_net = bws_net[i]; - double bw_time = bws_time[i]; - double et = events_time[i]; - - // calculating the density value on the network - if(event_counter[y] > 1){ - // here we have a dupplicated event, we must either use the savec nkde or calculate it and save it latter - if(map_contains_key(saved_values,y)){ - // we already have the value ! - samples_k = saved_values[y]; - }else{ - // we have not yet the value - samples_k = esc_kernel_rcpp_arma(kernel_func, neighbour_list, edge_mat, y, bw_net, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, max_depth); - saved_values[y] = samples_k; - } - }else{ - // this is not a duplicate event, no need to store the value in memory - samples_k = esc_kernel_rcpp_arma(kernel_func, neighbour_list, edge_mat, y, bw_net, line_weights, samples_edgeid, samples_x, samples_y, nodes_x, nodes_y, max_depth); - } - - // ok, now we must calculate the temporal densities - arma::vec temporal_density = kernel_func(arma::abs(samples_time - et), bw_time); - - // and standardize by bw if required - if(div == "bw"){ - temporal_density = temporal_density / bw_time; - samples_k = samples_k / bw_net; - } - - // and create an awesome spatio-temporal matrix - arma::mat k_mat(samples_x.n_elem, samples_time.n_elem); - - for(int j = 0; j < temporal_density.n_elem; j++){ - k_mat.col(j) = samples_k * temporal_density[j]; - } - - base_k += k_mat*w; - // calculating the new value - count.fill(0.0); - arma::umat ids = arma::find(samples_k>0); - count.elem(ids).fill(1.0); - base_count += count*w; - - }; - - // standardise the result if div = n - if(div == "n"){ - base_k = base_k / base_count; - } - - List results = List::create( - Named("sum_k") = base_k, - Named("n") = base_count ); - return results; -} - -