Something went wrong. Try again.
Distances on Directed Graphs in R
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952
#include "run_sp.h"#include "flows.h"
#include "dgraph.h"#include "heaps/heap_lib.h"
#include <unordered_map>
template <typename T>void inst_graph (std::shared_ptr<DGraph> g, size_t nedges, const std::map <std::string, size_t>& vert_map, const std::vector <std::string>& from, const std::vector <std::string>& to, const std::vector <T>& dist, const std::vector <T>& wt){ for (size_t i = 0; i < nedges; ++i) { size_t fromi = vert_map.at(from [i]); size_t toi = vert_map.at(to [i]); g->addNewEdge (fromi, toi, dist [i], wt [i], i); }}
struct OneAggregate : public RcppParallel::Worker{ RcppParallel::RVector <int> dp_fromi; const std::vector <size_t> toi; const RcppParallel::RMatrix <double> flows; const std::vector <std::string> vert_name; const std::unordered_map <std::string, size_t> verts_to_edge_map; size_t nverts; // can't be const because of reinterpret cast size_t nedges; const bool norm_sums; const double tol; const std::string heap_type; std::shared_ptr <DGraph> g;
std::vector <double> output;
// Constructor 1: The main constructor OneAggregate ( const RcppParallel::RVector <int> fromi, const std::vector <size_t> toi_in, const RcppParallel::RMatrix <double> flows_in, const std::vector <std::string> vert_name_in, const std::unordered_map <std::string, size_t> verts_to_edge_map_in, const size_t nverts_in, const size_t nedges_in, const bool norm_sums_in, const double tol_in, const std::string &heap_type_in, const std::shared_ptr <DGraph> g_in) : dp_fromi (fromi), toi (toi_in), flows (flows_in), vert_name (vert_name_in), verts_to_edge_map (verts_to_edge_map_in), nverts (nverts_in), nedges (nedges_in), norm_sums (norm_sums_in), tol (tol_in), heap_type (heap_type_in), g (g_in), output () { output.resize (nedges, 0.0); }
// Constructor 2: The Split constructor OneAggregate ( const OneAggregate& oneAggregate, RcppParallel::Split) : dp_fromi (oneAggregate.dp_fromi), toi (oneAggregate.toi), flows (oneAggregate.flows), vert_name (oneAggregate.vert_name), verts_to_edge_map (oneAggregate.verts_to_edge_map), nverts (oneAggregate.nverts), nedges (oneAggregate.nedges), norm_sums (oneAggregate.norm_sums), tol (oneAggregate.tol), heap_type (oneAggregate.heap_type), g (oneAggregate.g), output () { output.resize (nedges, 0.0); }
// Parallel function operator void operator() (size_t begin, size_t end) { std::shared_ptr<PF::PathFinder> pathfinder = std::make_shared <PF::PathFinder> (nverts, *run_sp::getHeapImpl (heap_type), g); std::vector <double> w (nverts); std::vector <double> d (nverts); std::vector <long int> prev (nverts);
for (size_t i = begin; i < end; i++) { //if (RcppThread::isInterrupted (i % static_cast<int>(100) == 0)) //if (RcppThread::isInterrupted ()) // return;
// These have to be reserved within the parallel operator function! std::fill (w.begin (), w.end (), INFINITE_DOUBLE); std::fill (d.begin (), d.end (), INFINITE_DOUBLE); std::fill (prev.begin (), prev.end (), INFINITE_INT);
size_t from_i = static_cast <size_t> (dp_fromi [i]); d [from_i] = w [from_i] = 0.0;
// reduce toi to only those within tolerance limt double flim = 0.0; size_t nto = toi.size (); if (tol > 0.0) { double fmax = 0.0; for (size_t j = 0; j < static_cast <size_t> (flows.ncol ()); j++) if (flows (i, j) > fmax) fmax = flows (i, j); flim = fmax * tol; size_t nto = 0; for (size_t j = 0; j < static_cast <size_t> (flows.ncol ()); j++) if (flows (i, j) > flim) nto++;
if (nto == 0) continue; }
// toi_index is into cols of flow matrix // toi_reduced is into the vertex vectors (d, w, prev) std::vector <size_t> toi_reduced, toi_index; toi_reduced.reserve (nto); toi_index.reserve (nto); for (size_t j = 0; j < static_cast <size_t> (flows.ncol ()); j++) if (flows (i, j) > flim) { toi_index.push_back (static_cast <size_t> (j)); toi_reduced.push_back (toi [j]); }
pathfinder->Dijkstra (d, w, prev, from_i, toi_reduced); for (size_t j = 0; j < toi_reduced.size (); j++) { if (from_i != toi_reduced [j]) // Exclude self-flows { double flow_ij = flows (i, toi_index [j]); if (w [toi_reduced [j]] < INFINITE_DOUBLE && flow_ij > 0.0) { // count how long the path is, so flows on // each edge can be divided by this length int path_len = 1; if (norm_sums) { path_len = 0; long int target_t = static_cast <long int> (toi_reduced [j]); size_t from_t = static_cast <size_t> (dp_fromi [i]); while (target_t < INFINITE_INT) { path_len++; size_t target_size_t = static_cast <size_t> (target_t); target_t = prev [target_size_t]; if (target_t < 0 || target_size_t == from_t) break; } }
long int target = static_cast <int> (toi_reduced [j]); // can equal -1 while (target < INFINITE_INT) { size_t stt = static_cast <size_t> (target); if (prev [stt] >= 0 && prev [stt] < INFINITE_INT) { std::string v2 = "f" + vert_name [static_cast <size_t> (prev [stt])] + "t" + vert_name [stt]; output [verts_to_edge_map.at (v2)] += flow_ij / static_cast <double> (path_len); }
target = prev [stt]; // Only allocate that flow from origin vertex v to all // previous vertices up until the target vi if (target < 0L || target == dp_fromi [i]) { break; } } } } } } // end for i } // end parallel function operator
void join (const OneAggregate &rhs) { for (size_t i = 0; i < output.size (); i++) output [i] += rhs.output [i]; }};
struct OneAggregatePaired : public RcppParallel::Worker{ RcppParallel::RVector <int> dp_fromtoi; const RcppParallel::RVector <double> flows; const std::vector <std::string> vert_name; const std::unordered_map <std::string, size_t> verts_to_edge_map; size_t nfrom; size_t nverts; // can't be const because of reinterpret cast size_t nedges; const bool norm_sums; const double tol; const std::string heap_type; std::shared_ptr <DGraph> g;
std::vector <double> output;
// Constructor 1: The main constructor OneAggregatePaired ( const RcppParallel::RVector <int> fromtoi, const RcppParallel::RVector <double> flows_in, const std::vector <std::string> vert_name_in, const std::unordered_map <std::string, size_t> verts_to_edge_map_in, const size_t nfrom_in, const size_t nverts_in, const size_t nedges_in, const bool norm_sums_in, const double tol_in, const std::string &heap_type_in, const std::shared_ptr <DGraph> g_in ) : dp_fromtoi (fromtoi), flows (flows_in), vert_name (vert_name_in), verts_to_edge_map (verts_to_edge_map_in), nfrom (nfrom_in), nverts (nverts_in), nedges (nedges_in), norm_sums (norm_sums_in), tol (tol_in), heap_type (heap_type_in), g (g_in), output () { output.resize (nedges, 0.0); }
// Constructor 2: The Split constructor OneAggregatePaired ( const OneAggregatePaired& oneAggregatePaired, RcppParallel::Split) : dp_fromtoi (oneAggregatePaired.dp_fromtoi), flows (oneAggregatePaired.flows), vert_name (oneAggregatePaired.vert_name), verts_to_edge_map (oneAggregatePaired.verts_to_edge_map), nfrom (oneAggregatePaired.nfrom), nverts (oneAggregatePaired.nverts), nedges (oneAggregatePaired.nedges), norm_sums (oneAggregatePaired.norm_sums), tol (oneAggregatePaired.tol), heap_type (oneAggregatePaired.heap_type), g (oneAggregatePaired.g), output () { output.resize (nedges, 0.0); }
// Parallel function operator void operator() (size_t begin, size_t end) { std::shared_ptr<PF::PathFinder> pathfinder = std::make_shared <PF::PathFinder> (nverts, *run_sp::getHeapImpl (heap_type), g); std::vector <double> w (nverts); std::vector <double> d (nverts); std::vector <long int> prev (nverts);
for (size_t i = begin; i < end; i++) { //if (RcppThread::isInterrupted (i % static_cast<int>(100) == 0)) //if (RcppThread::isInterrupted ()) // return;
const size_t from_i = static_cast <size_t> (dp_fromtoi [i]); // to_i has to be a vector for Dijkstra algo, but has only 1 // element. const std::vector <size_t> to_i = {static_cast <size_t> (dp_fromtoi [nfrom + i])};
// These have to be reserved within the parallel operator function! std::fill (w.begin (), w.end (), INFINITE_DOUBLE); std::fill (d.begin (), d.end (), INFINITE_DOUBLE); std::fill (prev.begin (), prev.end (), INFINITE_INT);
d [from_i] = w [from_i] = 0.0;
pathfinder->Dijkstra (d, w, prev, from_i, to_i); // loop has always only one element here: for (size_t j = 0; j < to_i.size (); j++) { if (w [to_i [j]] < INFINITE_DOUBLE && flows [i] > 0.0) { // count how long the path is, so flows on // each edge can be divided by this length int path_len = 1; if (norm_sums) { path_len = 0; long int target_t = static_cast <long int> (to_i [j]); while (target_t < INFINITE_INT) { path_len++; size_t target_size_t = static_cast <size_t> (target_t); target_t = prev [target_size_t]; if (target_t < 0 || target_size_t == from_i) break; } }
long int target = static_cast <int> (to_i [j]); // can equal -1 while (target < INFINITE_INT) { size_t stt = static_cast <size_t> (target); if (prev [stt] >= 0 && prev [stt] < INFINITE_INT) { std::string v2 = "f" + vert_name [static_cast <size_t> (prev [stt])] + "t" + vert_name [stt]; output [verts_to_edge_map.at (v2)] += flows [i] / static_cast <double> (path_len); }
target = prev [stt]; // Only allocate that flow from origin vertex v to all // previous vertices up until the target vi if (target < 0L || target == static_cast <long int> (from_i)) { break; } } } } } // end for i } // end parallel function operator
void join (const OneAggregatePaired &rhs) { for (size_t i = 0; i < output.size (); i++) output [i] += rhs.output [i]; }};
struct OneDisperse : public RcppParallel::Worker{ RcppParallel::RVector <int> dp_fromi; const RcppParallel::RVector <double> dens; const std::vector <std::string> vert_name; const std::unordered_map <std::string, size_t> verts_to_edge_map; size_t nverts; // can't be const because of reinterpret cast size_t nedges; const RcppParallel::RVector <double> kfrom; const double tol; const std::string heap_type; std::shared_ptr <DGraph> g;
std::vector <double> output;
// Constructor 1: The main constructor OneDisperse ( const RcppParallel::RVector <int> fromi, const RcppParallel::RVector <double> dens_in, const std::vector <std::string> vert_name_in, const std::unordered_map <std::string, size_t> verts_to_edge_map_in, const size_t nverts_in, const size_t nedges_in, const RcppParallel::RVector <double> kfrom_in, const double tol_in, const std::string &heap_type_in, const std::shared_ptr <DGraph> g_in) : dp_fromi (fromi), dens (dens_in), vert_name (vert_name_in), verts_to_edge_map (verts_to_edge_map_in), nverts (nverts_in), nedges (nedges_in), kfrom (kfrom_in), tol (tol_in), heap_type (heap_type_in), g (g_in), output () { const R_xlen_t nfrom = static_cast <R_xlen_t> (dens.size ()); const R_xlen_t nk = static_cast <R_xlen_t> (kfrom.size ()) / nfrom; const size_t out_size = nedges * static_cast <size_t> (nk); output.resize (out_size, 0.0); }
// Constructor 2: The Split constructor OneDisperse ( const OneDisperse& oneDisperse, RcppParallel::Split) : dp_fromi (oneDisperse.dp_fromi), dens (oneDisperse.dens), vert_name (oneDisperse.vert_name), verts_to_edge_map (oneDisperse.verts_to_edge_map), nverts (oneDisperse.nverts), nedges (oneDisperse.nedges), kfrom (oneDisperse.kfrom), tol (oneDisperse.tol), heap_type (oneDisperse.heap_type), g (oneDisperse.g), output () { const R_xlen_t nfrom = static_cast <R_xlen_t> (dens.size ()); const R_xlen_t nk = static_cast <R_xlen_t> (kfrom.size ()) / nfrom; size_t out_size = nedges * static_cast <size_t> (nk); output.resize (out_size, 0.0); }
// Parallel function operator void operator() (size_t begin, size_t end) { std::shared_ptr<PF::PathFinder> pathfinder = std::make_shared <PF::PathFinder> (nverts, *run_sp::getHeapImpl (heap_type), g); std::vector <double> w (nverts); std::vector <double> d (nverts); std::vector <long int> prev (nverts);
const size_t nfrom = dens.size (); const size_t nk = kfrom.size () / nfrom;
for (size_t i = begin; i < end; i++) // over the from vertices { //if (RcppThread::isInterrupted (i % static_cast<int>(10) == 0)) if (RcppThread::isInterrupted ()) return;
// translate k-value to distance limit based on tol // exp(-d / k) = tol -> d = -k * log (tol) // k_from holds nk vectors of different k-values, each of length // nedges. double dlim = 0.0; for (size_t k = 0; k < nk; k++) if (kfrom [i + k * nfrom] > dlim) dlim = kfrom [i + k * nfrom]; // dlim is max k-value dlim = -dlim * log (tol); // converted to actual dist limit.
std::fill (w.begin (), w.end (), INFINITE_DOUBLE); std::fill (d.begin (), d.end (), INFINITE_DOUBLE); std::fill (prev.begin (), prev.end (), INFINITE_INT);
const size_t from_i = static_cast <size_t> (dp_fromi [i]); d [from_i] = w [from_i] = 0.0;
pathfinder->DijkstraLimit (d, w, prev, from_i, dlim);
std::vector <double> flows_i (nedges * nk, 0.0); std::vector <double> expsum (nk, 0.0);
for (size_t j = 0; j < nverts; j++) { if (prev [j] >= 0 && prev [j] < INFINITE_INT) { const std::string vert_to = vert_name [j], vert_from = vert_name [static_cast <size_t> (prev [j])]; const std::string two_verts = "f" + vert_from + "t" + vert_to;
size_t index = verts_to_edge_map.at (two_verts); if (d [j] < INFINITE_DOUBLE) { for (size_t k = 0; k < nk; k++) { double exp_jk; if (kfrom [i + k * nfrom] > 0.0) exp_jk = exp (-d [j] / kfrom [i + k * nfrom]); else { // standard logistic polygonal for UK cycling // models double lp = -3.894 + (-0.5872 * d [j]) + (1.832 * sqrt (d [j])) + (0.007956 * d [j] * d [j]); exp_jk = exp (lp) / (1.0 + exp (lp)); } const size_t k_st = static_cast <size_t> (k); expsum [k_st] += exp_jk; flows_i [index + k_st * nedges] += dens [i] * exp_jk; } } } } // end for j for (size_t k = 0; k < nk; k++) if (expsum [k] > tol) for (size_t j = 0; j < nedges; j++) output [j + k * nedges] += flows_i [j + k * nedges] / expsum [k]; } // end for i } // end parallel function operator
void join (const OneDisperse &rhs) { for (size_t i = 0; i < output.size (); i++) output [i] += rhs.output [i]; }};
struct OneSI : public RcppParallel::Worker{ RcppParallel::RVector <int> dp_fromi; const std::vector <size_t> toi; const RcppParallel::RVector <double> k_from; const RcppParallel::RVector <double> dens_from; const RcppParallel::RVector <double> dens_to; const std::vector <std::string> vert_name; const std::unordered_map <std::string, size_t> verts_to_edge_map; size_t nverts; // can't be const because of reinterpret cast size_t nedges; const bool norm_sums; const double tol; const std::string heap_type; std::shared_ptr <DGraph> g;
std::vector <double> output;
// Constructor 1: The main constructor OneSI ( const RcppParallel::RVector <int> fromi, const std::vector <size_t> toi_in, const RcppParallel::RVector <double> k_from_in, const RcppParallel::RVector <double> dens_from_in, const RcppParallel::RVector <double> dens_to_in, const std::vector <std::string> vert_name_in, const std::unordered_map <std::string, size_t> verts_to_edge_map_in, const size_t nverts_in, const size_t nedges_in, const bool norm_sums_in, const double tol_in, const std::string &heap_type_in, const std::shared_ptr <DGraph> g_in) : dp_fromi (fromi), toi (toi_in), k_from (k_from_in), dens_from (dens_from_in), dens_to (dens_to_in), vert_name (vert_name_in), verts_to_edge_map (verts_to_edge_map_in), nverts (nverts_in), nedges (nedges_in), norm_sums (norm_sums_in), tol (tol_in), heap_type (heap_type_in), g (g_in), output () { const size_t nfrom = dens_from.size (); const size_t nk = k_from.size () / nfrom; const size_t out_size = nedges * nk; output.resize (out_size, 0.0); }
// Constructor 2: The Split constructor OneSI ( const OneSI& oneSI, RcppParallel::Split) : dp_fromi (oneSI.dp_fromi), toi (oneSI.toi), k_from (oneSI.k_from), dens_from (oneSI.dens_from), dens_to (oneSI.dens_to), vert_name (oneSI.vert_name), verts_to_edge_map (oneSI.verts_to_edge_map), nverts (oneSI.nverts), nedges (oneSI.nedges), norm_sums (oneSI.norm_sums), tol (oneSI.tol), heap_type (oneSI.heap_type), g (oneSI.g), output () { const size_t nfrom = dens_from.size (); const size_t nk = k_from.size () / nfrom; const size_t out_size = nedges * nk; output.resize (out_size, 0.0); }
// Parallel function operator void operator() (size_t begin, size_t end) { std::shared_ptr<PF::PathFinder> pathfinder = std::make_shared <PF::PathFinder> (nverts, *run_sp::getHeapImpl (heap_type), g); std::vector <double> w (nverts); std::vector <double> d (nverts); std::vector <long int> prev (nverts);
// k_from can have multiple vectors of k-values, each equal in length to // the number of 'from' points. The output is then a single vector of // 'nedges' wrapped 'nk' times. const size_t nfrom = dens_from.size (); const size_t nk = k_from.size () / nfrom;
for (size_t i = begin; i < end; i++) { //if (RcppThread::isInterrupted (i % static_cast<int>(100) == 0)) if (RcppThread::isInterrupted ()) return;
// These have to be reserved within the parallel operator function! std::fill (w.begin (), w.end (), INFINITE_DOUBLE); std::fill (d.begin (), d.end (), INFINITE_DOUBLE); std::fill (prev.begin (), prev.end (), INFINITE_INT);
size_t from_i = static_cast <size_t> (dp_fromi [i]); d [from_i] = w [from_i] = 0.0;
// translate k-value to distance limit based on tol // exp(-d / k) = tol -> d = -k * log (tol) // const double dlim = -k_from [i] * log (tol); // k_from holds nk vectors of different k-values, each of length // nedges. double dlim = 0.0; for (size_t k = 0; k < nk; k++) if (k_from [i + k * nfrom] > dlim) dlim = k_from [i + k * nfrom]; dlim = -dlim * log (tol);
pathfinder->DijkstraLimit (d, w, prev, from_i, dlim);
std::vector <double> flows_i (nedges * nk, 0.0); std::vector <double> expsum (nk, 0.0); for (size_t j = 0; j < toi.size (); j++) { if (from_i != toi [j]) // Exclude self-flows { if (d [toi [j]] < INFINITE_DOUBLE) { /* Flow between i and j is based on d (i, j), but is * subsequently allocated to all intermediate edges. * The `norm_sums` option re-scales this intermediate * edge allocation such that SI value to one terminal * vertex from the origin vertex of: * N_i N_j exp (-d_{ij}) / sum_j (N_{ij} exp (-d_{ij})) * is divided by the number of intermediate edges * between i and j. This is achieved by initially only * allocating the values of the numerator divided by * numbers of intermediate edges, while accumulating the * denominator independent of that value in order to * divide all final values. */ std::vector <double> exp_jk (nk, 0.0), flow_ijk (nk, 0.0); for (size_t k = 0; k < nk; k++) { exp_jk [k] = dens_to [j] * exp (-d [toi [j]] / k_from [i + k * nfrom]); flow_ijk [k] = dens_from [i] * exp_jk [k]; expsum [k] += exp_jk [k]; }
// need to count how long the path is, so flows on // each edge can be divided by this length int path_len = 1; if (norm_sums) { path_len = 0; long int target_t = static_cast <long int> (toi [j]); size_t from_t = static_cast <size_t> (dp_fromi [i]); while (target_t < INFINITE_INT) { size_t target_size_t = static_cast <size_t> (target_t); path_len++; target_t = prev [target_size_t]; if (target_t < 0 || target_size_t == from_t) break; } }
long int target = static_cast <long int> (toi [j]); // can equal -1 while (target < INFINITE_INT) { size_t stt = static_cast <size_t> (target); if (prev [stt] >= 0 && prev [stt] < INFINITE_INT) { std::string v2 = "f" + vert_name [static_cast <size_t> (prev [stt])] + "t" + vert_name [stt]; // multiple flows can aggregate to same edge, so // this has to be +=, not just =! size_t index = verts_to_edge_map.at (v2); for (size_t k = 0; k < nk; k++) { flows_i [index + k * nedges] += flow_ijk [k] / static_cast <double> (path_len); } }
target = static_cast <long int> (prev [stt]); // Only allocate that flow from origin vertex v to all // previous vertices up until the target vi if (target < 0L || target == dp_fromi [i]) { break; } } } } } // end for j for (size_t k = 0; k < nk; k++) if (expsum [k] > tol) for (size_t j = 0; j < nedges; j++) output [j + k * nedges] += flows_i [j + k * nedges] / expsum [k]; } // end for i } // end parallel function operator
void join (const OneSI& rhs) { for (size_t i = 0; i < output.size (); i++) output [i] += rhs.output [i]; }};
//' rcpp_flows_aggregate_par//'//' @param graph The data.frame holding the graph edges//' @param vert_map_in map from <std::string> vertex ID to (0-indexed) integer//' index of vertices//' @param fromi Index into vert_map_in of vertex numbers//' @param toi Index into vert_map_in of vertex numbers//' @param tol Relative tolerance in terms of flows below which targets//' (to-vertices) are not considered.//'//' @note The parallelisation is achieved by dumping the results of each thread//' to a file, with aggregation performed at the end by simply reading back and//' aggregating all files. There is no way to aggregate into a single vector//' because threads have to be independent. The only danger with this approach//' is that multiple threads may generate the same file names, but with names 10//' characters long, that chance should be 1 / 62 ^ 10.//'//' @noRd// [[Rcpp::export]]Rcpp::NumericVector rcpp_flows_aggregate_par (const Rcpp::DataFrame graph, const Rcpp::DataFrame vert_map_in, Rcpp::IntegerVector fromi, Rcpp::IntegerVector toi_in, Rcpp::NumericMatrix flows, const bool norm_sums, const double tol, const std::string heap_type){ std::vector <size_t> toi = Rcpp::as <std::vector <size_t> > (toi_in); const size_t nfrom = static_cast <size_t> (fromi.size ());
const std::vector <std::string> from = graph ["from"]; const std::vector <std::string> to = graph ["to"]; const std::vector <double> dist = graph ["d"]; const std::vector <double> wt = graph ["d_weighted"];
const size_t nedges = static_cast <size_t> (graph.nrow ()); const std::vector <std::string> vert_name = vert_map_in ["vert"]; const std::vector <size_t> vert_indx = vert_map_in ["id"]; // Make map from vertex name to integer index std::map <std::string, size_t> vert_map_i; const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_name, vert_indx, vert_map_i);
std::unordered_map <std::string, size_t> verts_to_edge_map; std::unordered_map <std::string, double> verts_to_dist_map; run_sp::make_vert_to_edge_maps (from, to, wt, verts_to_edge_map, verts_to_dist_map);
std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts); inst_graph (g, nedges, vert_map_i, from, to, dist, wt);
// Create parallel worker OneAggregate oneAggregate (RcppParallel::RVector <int> (fromi), toi, RcppParallel::RMatrix <double> (flows), vert_name, verts_to_edge_map, nverts, nedges, norm_sums, tol, heap_type, g);
size_t chunk_size = run_sp::get_chunk_size (nfrom); RcppParallel::parallelReduce (0, nfrom, oneAggregate, chunk_size);
return Rcpp::wrap (oneAggregate.output);}
//' rcpp_flows_aggregate_pairwise//'//' Pairwise version of flows_aggregate_par//'//' @param graph The data.frame holding the graph edges//' @param vert_map_in map from <std::string> vertex ID to (0-indexed) integer//' index of vertices//' @param fromi Index into vert_map_in of vertex numbers//' @param toi Index into vert_map_in of vertex numbers//' @param tol Relative tolerance in terms of flows below which targets//' (to-vertices) are not considered.//'//' @note The parallelisation is achieved by dumping the results of each thread//' to a file, with aggregation performed at the end by simply reading back and//' aggregating all files. There is no way to aggregate into a single vector//' because threads have to be independent. The only danger with this approach//' is that multiple threads may generate the same file names, but with names 10//' characters long, that chance should be 1 / 62 ^ 10.//'//' @noRd// [[Rcpp::export]]Rcpp::NumericVector rcpp_flows_aggregate_pairwise (const Rcpp::DataFrame graph, const Rcpp::DataFrame vert_map_in, Rcpp::IntegerVector fromi, Rcpp::IntegerVector toi, Rcpp::NumericVector flows, const bool norm_sums, const double tol, const std::string heap_type){ if (fromi.size () != toi.size ()) Rcpp::stop ("pairwise dists must have from.size == to.size"); long int n = fromi.size (); size_t n_st = static_cast <size_t> (n);
const std::vector <std::string> from = graph ["from"]; const std::vector <std::string> to = graph ["to"]; const std::vector <double> dist = graph ["d"]; const std::vector <double> wt = graph ["d_weighted"];
const size_t nedges = static_cast <size_t> (graph.nrow ()); const std::vector <std::string> vert_name = vert_map_in ["vert"]; const std::vector <size_t> vert_indx = vert_map_in ["id"]; // Make map from vertex name to integer index std::map <std::string, size_t> vert_map_i; const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_name, vert_indx, vert_map_i);
std::unordered_map <std::string, size_t> verts_to_edge_map; std::unordered_map <std::string, double> verts_to_dist_map; run_sp::make_vert_to_edge_maps (from, to, wt, verts_to_edge_map, verts_to_dist_map);
std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts); inst_graph (g, nedges, vert_map_i, from, to, dist, wt);
// Paired (fromi, toi) in a single vector Rcpp::IntegerVector fromto (2 * n_st); for (int i = 0; i < n; i++) { size_t i_t = static_cast <size_t> (i); fromto [i] = fromi (i_t); fromto [i + n] = toi (i_t); }
// Create parallel worker OneAggregatePaired oneAggregatePaired (RcppParallel::RVector <int> (fromto), RcppParallel::RVector <double> (flows), vert_name, verts_to_edge_map, n_st, nverts, nedges, norm_sums, tol, heap_type, g);
size_t chunk_size = run_sp::get_chunk_size (n_st); RcppParallel::parallelReduce (0, n_st, oneAggregatePaired, chunk_size);
return Rcpp::wrap (oneAggregatePaired.output);}
//' rcpp_flows_disperse_par//'//' Modified version of \code{rcpp_flows_aggregate} that aggregates flows to all//' destinations from given set of origins, with flows attenuated by distance//' from those origins.//'//' @param graph The data.frame holding the graph edges//' @param vert_map_in map from <std::string> vertex ID to (0-indexed) integer//' index of vertices//' @param fromi Index into vert_map_in of vertex numbers//' @param k Coefficient of (current proof-of-principle-only) exponential//' distance decay function. If value of \code{k<0} is given, a standard//' logistic polynomial will be used.//'//' @note The flow data to be used for aggregation is a matrix mapping flows//' between each pair of from and to points.//'//' @noRd// [[Rcpp::export]]Rcpp::NumericVector rcpp_flows_disperse_par (const Rcpp::DataFrame graph, const Rcpp::DataFrame vert_map_in, Rcpp::IntegerVector fromi, Rcpp::NumericVector k, Rcpp::NumericVector dens, const double &tol, std::string heap_type){ const size_t nfrom = static_cast <size_t> (fromi.size ());
std::vector <std::string> from = graph ["from"]; std::vector <std::string> to = graph ["to"]; std::vector <double> dist = graph ["d"]; std::vector <double> wt = graph ["d_weighted"];
size_t nedges = static_cast <size_t> (graph.nrow ()); std::vector <std::string> vert_name = vert_map_in ["vert"]; std::vector <size_t> vert_indx = vert_map_in ["id"]; // Make map from vertex name to integer index std::map <std::string, size_t> vert_map_i; size_t nverts = run_sp::make_vert_map (vert_map_in, vert_name, vert_indx, vert_map_i);
std::unordered_map <std::string, size_t> verts_to_edge_map; std::unordered_map <std::string, double> verts_to_dist_map; run_sp::make_vert_to_edge_maps (from, to, wt, verts_to_edge_map, verts_to_dist_map);
std::shared_ptr<DGraph> g = std::make_shared<DGraph> (nverts); inst_graph (g, nedges, vert_map_i, from, to, dist, wt);
// Create parallel worker OneDisperse oneDisperse (RcppParallel::RVector <int> (fromi), RcppParallel::RVector <double> (dens), vert_name, verts_to_edge_map, nverts, nedges, RcppParallel::RVector <double> (k), tol, heap_type, g);
size_t chunk_size = run_sp::get_chunk_size (nfrom); RcppParallel::parallelReduce (0, nfrom, oneDisperse, chunk_size);
return Rcpp::wrap (oneDisperse.output);}
//' rcpp_flows_si//'//' @param graph The data.frame holding the graph edges//' @param vert_map_in map from <std::string> vertex ID to (0-indexed) integer//' index of vertices//' @param fromi Index into vert_map_in of vertex numbers//' @param toi Index into vert_map_in of vertex numbers//' @param kvec Vector of k-values for each fromi//' @param nvec Vector of density-values for each fromi//' @param tol Relative tolerance in terms of flows below which targets//' (to-vertices) are not considered.//'//' @noRd// [[Rcpp::export]]Rcpp::NumericVector rcpp_flows_si (const Rcpp::DataFrame graph, const Rcpp::DataFrame vert_map_in, Rcpp::IntegerVector fromi, Rcpp::IntegerVector toi_in, Rcpp::NumericVector kvec, Rcpp::NumericVector dens_from, Rcpp::NumericVector dens_to, const bool norm_sums, const double tol, const std::string heap_type){ std::vector <size_t> toi = Rcpp::as <std::vector <size_t> > ( toi_in); const size_t nfrom = static_cast <size_t> (fromi.size ());
const std::vector <std::string> from = graph ["from"]; const std::vector <std::string> to = graph ["to"]; const std::vector <double> dist = graph ["d"]; const std::vector <double> wt = graph ["d_weighted"];
const size_t nedges = static_cast <size_t> (graph.nrow ()); const std::vector <std::string> vert_name = vert_map_in ["vert"]; const std::vector <size_t> vert_indx = vert_map_in ["id"]; // Make map from vertex name to integer index std::map <std::string, size_t> vert_map_i; const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_name, vert_indx, vert_map_i);
std::unordered_map <std::string, size_t> verts_to_edge_map; std::unordered_map <std::string, double> verts_to_dist_map; run_sp::make_vert_to_edge_maps (from, to, wt, verts_to_edge_map, verts_to_dist_map);
std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts); inst_graph (g, nedges, vert_map_i, from, to, dist, wt);
// Create parallel worker OneSI oneSI (RcppParallel::RVector <int> (fromi), toi, RcppParallel::RVector <double> (kvec), RcppParallel::RVector <double> (dens_from), RcppParallel::RVector <double> (dens_to), vert_name, verts_to_edge_map, nverts, nedges, norm_sums, tol, heap_type, g);
size_t chunk_size = run_sp::get_chunk_size (nfrom); RcppParallel::parallelReduce (0, nfrom, oneSI, chunk_size);
return Rcpp::wrap (oneSI.output);}