mpadge.tngl.sh
many-to-many multi-modal routing aggregator
1
2#include "run_sp.h"
3
4#include "dgraph.h"
5#include "bheap.h"
6
7// # nocov start
8template <typename T>
9void inst_graph (std::shared_ptr<DGraph> g, size_t nedges,
10 const std::map <std::string, size_t>& vert_map,
11 const std::vector <std::string>& from,
12 const std::vector <std::string>& to,
13 const std::vector <T>& dist,
14 const std::vector <T>& wt)
15{
16 for (size_t i = 0; i < nedges; ++i)
17 {
18 size_t fromi = vert_map.at(from [i]);
19 size_t toi = vert_map.at(to [i]);
20 g->addNewEdge (fromi, toi, dist [i], wt [i], i);
21 }
22}
23// # nocov end
24
25// RcppParallel jobs can be chunked to a specified "grain size"; see
26// https://rcppcore.github.io/RcppParallel/#grain_size
27// This function determines chunk size such that there are at least 100 chunks
28// for a given `nfrom`.
29size_t run_sp::get_chunk_size (const size_t nfrom)
30{
31 size_t chunk_size;
32
33 if (nfrom > 1000)
34 chunk_size = 100;
35 else if (nfrom > 100)
36 chunk_size = 10;
37 else
38 chunk_size = 1;
39
40 return chunk_size;
41}
42
43
44std::shared_ptr <HeapDesc> run_sp::getHeapImpl(const std::string& heap_type)
45{
46 return std::make_shared <HeapD<BHeap> >();
47}
48
49
50struct OneDist : public RcppParallel::Worker
51{
52 RcppParallel::RVector <int> dp_fromi;
53 const std::vector <size_t> toi;
54 const size_t nverts;
55 const std::vector <double> vx;
56 const std::vector <double> vy;
57 const std::shared_ptr <DGraph> g;
58
59 RcppParallel::RMatrix <double> dout;
60
61 // constructor
62 OneDist (
63 const RcppParallel::RVector <int> fromi,
64 const std::vector <size_t> toi_in,
65 const size_t nverts_in,
66 const std::vector <double> vx_in,
67 const std::vector <double> vy_in,
68 const std::shared_ptr <DGraph> g_in,
69 RcppParallel::RMatrix <double> dout_in) :
70 dp_fromi (fromi), toi (toi_in), nverts (nverts_in),
71 vx (vx_in), vy (vy_in), g (g_in),
72 dout (dout_in)
73 {
74 }
75
76 // Parallel function operator
77 void operator() (std::size_t begin, std::size_t end)
78 {
79 const std::string heap_type = "BHeap";
80 std::shared_ptr<PF::PathFinder> pathfinder =
81 std::make_shared <PF::PathFinder> (nverts,
82 *run_sp::getHeapImpl (heap_type), g);
83 std::vector <double> w (nverts);
84 std::vector <double> d (nverts);
85 std::vector <long int> prev (nverts);
86
87 std::vector <double> heuristic (nverts, 0.0);
88
89 for (std::size_t i = begin; i < end; i++)
90 {
91 size_t from_i = static_cast <size_t> (dp_fromi [i]);
92
93 for (size_t j = 0; j < nverts; j++)
94 {
95 const double dx = vx [j] - vx [from_i],
96 dy = vy [j] - vy [from_i];
97 heuristic [j] = sqrt (dx * dx + dy * dy);
98 }
99 pathfinder->AStar (d, w, prev, heuristic, from_i, toi);
100
101 for (size_t j = 0; j < toi.size (); j++)
102 {
103 if (w [toi [j]] < INFINITE_DOUBLE)
104 {
105 dout (i, j) = d [toi [j]];
106 }
107 }
108 }
109 }
110
111};
112
113struct OneWeightedDist : public RcppParallel::Worker
114{
115 RcppParallel::RVector <int> dp_fromi;
116 const std::vector <size_t> toi;
117 const std::vector <double> wti;
118 const size_t nverts;
119 const double k;
120 const double dlim;
121 const std::vector <double> vx;
122 const std::vector <double> vy;
123 const std::shared_ptr <DGraph> g;
124
125 RcppParallel::RMatrix <double> dout;
126
127 // constructor
128 OneWeightedDist (
129 const RcppParallel::RVector <int> fromi,
130 const std::vector <size_t> toi_in,
131 const std::vector <double> wti_in,
132 const size_t nverts_in,
133 const double k_in,
134 const double dlim_in,
135 const std::vector <double> vx_in,
136 const std::vector <double> vy_in,
137 const std::shared_ptr <DGraph> g_in,
138 RcppParallel::RMatrix <double> dout_in) :
139 dp_fromi (fromi), toi (toi_in), wti (wti_in),
140 nverts (nverts_in), k (k_in), dlim (dlim_in),
141 vx (vx_in), vy (vy_in), g (g_in),
142 dout (dout_in)
143 {
144 }
145
146 // Parallel function operator
147 void operator() (std::size_t begin, std::size_t end)
148 {
149 const std::string heap_type = "BHeap";
150 std::shared_ptr<PF::PathFinder> pathfinder =
151 std::make_shared <PF::PathFinder> (nverts,
152 *run_sp::getHeapImpl (heap_type), g);
153 std::vector <double> w (nverts);
154 std::vector <double> d (nverts);
155 std::vector <long int> prev (nverts);
156
157 std::vector <double> heuristic (nverts, 0.0);
158
159 for (std::size_t i = begin; i < end; i++)
160 {
161 size_t from_i = static_cast <size_t> (dp_fromi [i]);
162
163 for (size_t j = 0; j < nverts; j++)
164 {
165 const double dx = vx [j] - vx [from_i],
166 dy = vy [j] - vy [from_i];
167 heuristic [j] = sqrt (dx * dx + dy * dy);
168 }
169 pathfinder->DijkstraLimit (d, w, prev, from_i, dlim);
170
171 for (size_t j = 0; j < toi.size (); j++)
172 {
173 if (d [toi [j]] < INFINITE_DOUBLE)
174 {
175 dout (i, 0) += exp (-d [toi [j]] / k);
176 dout (i, 1) += dout (i, 0) * wti [j];
177 }
178 }
179 }
180 }
181
182};
183
184struct OneDistNTargets : public RcppParallel::Worker
185{
186 RcppParallel::RVector <int> dp_fromi;
187 const std::vector <size_t> toi;
188 const size_t nverts;
189 const size_t n_targets;
190 const std::vector <double> vx;
191 const std::vector <double> vy;
192 const std::shared_ptr <DGraph> g;
193
194 RcppParallel::RMatrix <double> dout;
195
196 // constructor
197 OneDistNTargets (
198 const RcppParallel::RVector <int> fromi,
199 const std::vector <size_t> toi_in,
200 const size_t nverts_in,
201 const size_t n_targets_in,
202 const std::vector <double> vx_in,
203 const std::vector <double> vy_in,
204 const std::shared_ptr <DGraph> g_in,
205 RcppParallel::RMatrix <double> dout_in) :
206 dp_fromi (fromi), toi (toi_in),
207 nverts (nverts_in), n_targets (n_targets_in),
208 vx (vx_in), vy (vy_in), g (g_in),
209 dout (dout_in)
210 {
211 }
212
213 // Parallel function operator
214 void operator() (std::size_t begin, std::size_t end)
215 {
216 const std::string heap_type = "BHeap";
217 std::shared_ptr<PF::PathFinder> pathfinder =
218 std::make_shared <PF::PathFinder> (nverts,
219 *run_sp::getHeapImpl (heap_type), g);
220 std::vector <double> w (nverts);
221 std::vector <double> d (nverts);
222 std::vector <long int> prev (nverts);
223
224 std::vector <double> heuristic (nverts, 0.0);
225
226 for (std::size_t i = begin; i < end; i++)
227 {
228 std::fill (d.begin (), d.end (), INFINITE_DOUBLE);
229 std::fill (w.begin (), w.end (), INFINITE_DOUBLE);
230 std::fill (heuristic.begin (), heuristic.end (), 0.0);
231 std::fill (prev.begin (), prev.end (), INFINITE_INT);
232
233 size_t from_i = static_cast <size_t> (dp_fromi [i]);
234
235 for (size_t j = 0; j < nverts; j++)
236 {
237 const double dx = vx [j] - vx [from_i],
238 dy = vy [j] - vy [from_i];
239 heuristic [j] = sqrt (dx * dx + dy * dy);
240 }
241 pathfinder->DijkstraNTargets (d, w, prev, from_i, toi, n_targets);
242
243 // The "NTargets" scanner only accumulates distances up to the
244 // closest "n_targets" values, but still includes the full "toi"
245 // values, most of which are NA. The loop to find the closest ones
246 // therefore has to examine the whole "toi" vector and use an index.
247 size_t count = 0;
248 for (size_t j = 0; j < toi.size (); j++)
249 {
250 if (d [toi [j]] < INFINITE_DOUBLE)
251 {
252 dout (i, count) = d [toi [j]];
253 dout (i, n_targets + count) = static_cast <double> (toi [j]);
254 count++;
255 }
256 if (count >= n_targets)
257 {
258 break;
259 }
260 }
261 }
262 }
263
264};
265
266
267struct SaveOneDist : public RcppParallel::Worker
268{
269 RcppParallel::RVector <int> dp_fromi;
270 const std::vector <std::string> from_names;
271 const std::vector <size_t> toi;
272 const size_t nverts;
273 const std::vector <double> vx;
274 const std::vector <double> vy;
275 const std::string path;
276 const std::shared_ptr <DGraph> g;
277
278 // constructor
279 SaveOneDist (
280 const RcppParallel::RVector <int> fromi,
281 const std::vector <std::string> from_names_in,
282 const std::vector <size_t> toi_in,
283 const size_t nverts_in,
284 const std::vector <double> vx_in,
285 const std::vector <double> vy_in,
286 const std::string path_in,
287 const std::shared_ptr <DGraph> g_in) :
288 dp_fromi (fromi), from_names (from_names_in), toi (toi_in),
289 nverts (nverts_in), vx (vx_in), vy (vy_in),
290 path (path_in), g (g_in)
291 {
292 }
293
294 // Parallel function operator
295 void operator() (std::size_t begin, std::size_t end)
296 {
297 const std::string heap_type = "BHeap";
298 std::shared_ptr<PF::PathFinder> pathfinder =
299 std::make_shared <PF::PathFinder> (nverts,
300 *run_sp::getHeapImpl (heap_type), g);
301 std::vector <double> w (nverts);
302 std::vector <double> d (nverts);
303 std::vector <long int> prev (nverts);
304
305 std::vector <double> heuristic (nverts, 0.0);
306
307 for (std::size_t i = begin; i < end; i++)
308 {
309 size_t from_i = static_cast <size_t> (dp_fromi [i]);
310
311 for (size_t j = 0; j < nverts; j++)
312 {
313 const double dx = vx [j] - vx [from_i],
314 dy = vy [j] - vy [from_i];
315 heuristic [j] = sqrt (dx * dx + dy * dy);
316 }
317 pathfinder->AStar (d, w, prev, heuristic, from_i, toi);
318
319 const std::string fname = path + "m4ra_from_" + from_names [i];
320 std::ofstream out_file;
321 out_file.open (fname.c_str (), std::ofstream::out);
322
323 for (size_t j = 0; j < toi.size (); j++)
324 {
325 double dtemp = -1.0;
326 if (w [toi [j]] < INFINITE_DOUBLE)
327 {
328 dtemp = d [toi [j]];
329 }
330 out_file << std::setprecision (10) << dtemp << std::endl;
331 }
332 out_file.close ();
333 }
334 }
335
336};
337
338
339size_t run_sp::make_vert_map (const Rcpp::DataFrame &vert_map_in,
340 const std::vector <std::string> &vert_map_id,
341 const std::vector <size_t> &vert_map_n,
342 std::map <std::string, size_t> &vert_map)
343{
344 for (size_t i = 0;
345 i < static_cast <size_t> (vert_map_in.nrow ()); ++i)
346 {
347 vert_map.emplace (vert_map_id [i], vert_map_n [i]);
348 }
349 size_t nverts = static_cast <size_t> (vert_map.size ());
350 return (nverts);
351}
352
353//' rcpp_weighted_dists
354//'
355//' @noRd
356// [[Rcpp::export]]
357Rcpp::NumericMatrix rcpp_weighted_dists (const Rcpp::DataFrame graph,
358 const Rcpp::DataFrame vert_map_in,
359 Rcpp::IntegerVector fromi,
360 Rcpp::IntegerVector toi_in,
361 Rcpp::NumericVector weights,
362 const double dlim,
363 const double k)
364{
365 std::vector <size_t> toi =
366 Rcpp::as <std::vector <size_t> > ( toi_in);
367
368 size_t nfrom = static_cast <size_t> (fromi.size ());
369
370 const std::vector <std::string> from = graph [".vx0"];
371 const std::vector <std::string> to = graph [".vx1"];
372 const std::vector <double> dist = graph ["d"];
373 const std::vector <double> wt = graph ["d_weighted"];
374
375 const size_t nedges = static_cast <size_t> (graph.nrow ());
376 std::map <std::string, size_t> vert_map;
377 std::vector <std::string> vert_map_id = vert_map_in ["vert"];
378 std::vector <size_t> vert_map_n = vert_map_in ["id"];
379 const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
380 vert_map_n, vert_map);
381
382 std::vector <double> vx (nverts), vy (nverts), wts;
383 vx = Rcpp::as <std::vector <double> > (vert_map_in ["x"]);
384 vy = Rcpp::as <std::vector <double> > (vert_map_in ["y"]);
385 wts = Rcpp::as <std::vector <double> > (weights);
386
387 std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
388 inst_graph (g, nedges, vert_map, from, to, dist, wt);
389
390 Rcpp::NumericMatrix dout (static_cast <int> (nfrom), 2);
391
392 // Create parallel worker
393 OneWeightedDist one_dist (RcppParallel::RVector <int> (fromi), toi,
394 wts, nverts, k, dlim, vx, vy, g,
395 RcppParallel::RMatrix <double> (dout));
396
397 size_t chunk_size = run_sp::get_chunk_size (nfrom);
398 RcppParallel::parallelFor (0, nfrom, one_dist, chunk_size);
399
400 return (dout);
401}
402
403//' rcpp_get_sp_dists_par
404//'
405//' @noRd
406// [[Rcpp::export]]
407Rcpp::NumericMatrix rcpp_get_sp_dists_par (const Rcpp::DataFrame graph,
408 const Rcpp::DataFrame vert_map_in,
409 Rcpp::IntegerVector fromi,
410 Rcpp::IntegerVector toi_in)
411{
412 std::vector <size_t> toi =
413 Rcpp::as <std::vector <size_t> > ( toi_in);
414
415 size_t nfrom = static_cast <size_t> (fromi.size ());
416 size_t nto = static_cast <size_t> (toi.size ());
417
418 const std::vector <std::string> from = graph [".vx0"];
419 const std::vector <std::string> to = graph [".vx1"];
420 const std::vector <double> dist = graph ["d"];
421 const std::vector <double> wt = graph ["d_weighted"];
422
423 const size_t nedges = static_cast <size_t> (graph.nrow ());
424 std::map <std::string, size_t> vert_map;
425 std::vector <std::string> vert_map_id = vert_map_in ["vert"];
426 std::vector <size_t> vert_map_n = vert_map_in ["id"];
427 const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
428 vert_map_n, vert_map);
429
430 std::vector <double> vx (nverts), vy (nverts);
431 vx = Rcpp::as <std::vector <double> > (vert_map_in ["x"]);
432 vy = Rcpp::as <std::vector <double> > (vert_map_in ["y"]);
433
434 std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
435 inst_graph (g, nedges, vert_map, from, to, dist, wt);
436
437 Rcpp::NumericVector na_vec = Rcpp::NumericVector (nfrom * nto,
438 Rcpp::NumericVector::get_na ());
439 Rcpp::NumericMatrix dout (static_cast <int> (nfrom),
440 static_cast <int> (nto), na_vec.begin ());
441
442 // Create parallel worker
443 OneDist one_dist (RcppParallel::RVector <int> (fromi), toi,
444 nverts, vx, vy, g,
445 RcppParallel::RMatrix <double> (dout));
446
447 size_t chunk_size = run_sp::get_chunk_size (nfrom);
448 RcppParallel::parallelFor (0, nfrom, one_dist, chunk_size);
449
450 return (dout);
451}
452
453//' rcpp_dists_to_n_targets
454//'
455//' @noRd
456// [[Rcpp::export]]
457Rcpp::NumericMatrix rcpp_dists_to_n_targets (const Rcpp::DataFrame graph,
458 const Rcpp::DataFrame vert_map_in,
459 Rcpp::IntegerVector fromi,
460 Rcpp::IntegerVector toi_in,
461 const int n_targets)
462{
463 std::vector <size_t> toi =
464 Rcpp::as <std::vector <size_t> > ( toi_in);
465
466 size_t nfrom = static_cast <size_t> (fromi.size ());
467
468 const std::vector <std::string> from = graph [".vx0"];
469 const std::vector <std::string> to = graph [".vx1"];
470 const std::vector <double> dist = graph ["d"];
471 const std::vector <double> wt = graph ["d_weighted"];
472
473 const size_t nedges = static_cast <size_t> (graph.nrow ());
474 std::map <std::string, size_t> vert_map;
475 std::vector <std::string> vert_map_id = vert_map_in ["vert"];
476 std::vector <size_t> vert_map_n = vert_map_in ["id"];
477 const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
478 vert_map_n, vert_map);
479
480 std::vector <double> vx (nverts), vy (nverts);
481 vx = Rcpp::as <std::vector <double> > (vert_map_in ["x"]);
482 vy = Rcpp::as <std::vector <double> > (vert_map_in ["y"]);
483
484 std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
485 inst_graph (g, nedges, vert_map, from, to, dist, wt);
486
487 Rcpp::NumericVector na_vec = Rcpp::NumericVector (nfrom * 2 * static_cast <size_t> (n_targets),
488 Rcpp::NumericVector::get_na ());
489 Rcpp::NumericMatrix dout (static_cast <int> (nfrom),
490 static_cast <int> (2 * n_targets), na_vec.begin ());
491
492 // Create parallel worker
493 OneDistNTargets one_dist_n_targets (RcppParallel::RVector <int> (fromi), toi,
494 nverts, static_cast <size_t> (n_targets),
495 vx, vy, g,
496 RcppParallel::RMatrix <double> (dout));
497
498 const size_t chunk_size = run_sp::get_chunk_size (nfrom);
499
500 RcppParallel::parallelFor (0, nfrom, one_dist_n_targets, chunk_size);
501
502 return (dout);
503}
504
505//' rcpp_save_sp_dists_par
506//'
507//' This function doesn't return anything useful. The R function which calls it
508//' ultimately returns the paths to all files created here, but that is not
509//' really possible with Rcpp, because of difficulties getting CharacterVectors
510//' to play nicely with RcppParallel::RVector. The file name matching is easily
511//' done on the R side anyway.
512//'
513//' @noRd
514// [[Rcpp::export]]
515bool rcpp_save_sp_dists_par (const Rcpp::DataFrame graph,
516 const Rcpp::DataFrame vert_map_in,
517 Rcpp::IntegerVector from_index,
518 Rcpp::CharacterVector from_names_in,
519 Rcpp::IntegerVector toi_in,
520 const std::string path)
521{
522 std::vector <size_t> toi =
523 Rcpp::as <std::vector <size_t> > ( toi_in);
524 std::vector <std::string> from_names =
525 Rcpp::as <std::vector <std::string> > (from_names_in);
526
527 size_t nfrom = static_cast <size_t> (from_index.size ());
528
529 const std::vector <std::string> from = graph [".vx0"];
530 const std::vector <std::string> to = graph [".vx1"];
531 const std::vector <double> dist = graph ["d"];
532 const std::vector <double> wt = graph ["d_weighted"];
533
534 const size_t nedges = static_cast <size_t> (graph.nrow ());
535 std::map <std::string, size_t> vert_map;
536 std::vector <std::string> vert_map_id = vert_map_in ["vert"];
537 std::vector <size_t> vert_map_n = vert_map_in ["id"];
538 const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
539 vert_map_n, vert_map);
540
541 std::vector <double> vx (nverts), vy (nverts);
542 vx = Rcpp::as <std::vector <double> > (vert_map_in ["x"]);
543 vy = Rcpp::as <std::vector <double> > (vert_map_in ["y"]);
544
545 std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
546 inst_graph (g, nedges, vert_map, from, to, dist, wt);
547
548 // Create parallel worker
549 SaveOneDist one_dist (RcppParallel::RVector <int> (from_index), from_names,
550 toi, nverts, vx, vy, path, g);
551
552 size_t chunk_size = run_sp::get_chunk_size (nfrom);
553 RcppParallel::parallelFor (0, nfrom, one_dist, chunk_size);
554
555 return (true);
556}