mpadge.tngl.sh
Distances on Directed Graphs in R
1// Modified functions from run_sp.cpp to calculate categorical distances along
2// defined kinds of ways (issue #144)
3
4#include "run_sp.h"
5
6#include "dgraph.h"
7#include "heaps/heap_lib.h"
8#include <unordered_set>
9
10// # nocov start
11template <typename T>
12void inst_graph (std::shared_ptr<DGraph> g, size_t nedges,
13 const std::map <std::string, size_t>& vert_map,
14 const std::vector <std::string>& from,
15 const std::vector <std::string>& to,
16 const std::vector <size_t>& edge_type,
17 const std::vector <T>& dist,
18 const std::vector <T>& wt)
19{
20 for (size_t i = 0; i < nedges; ++i)
21 {
22 size_t fromi = vert_map.at(from [i]);
23 size_t toi = vert_map.at(to [i]);
24 g->addNewEdge (fromi, toi, dist [i], wt [i], edge_type [i]);
25 }
26}
27// # nocov end
28
29struct OneCategoricalDist : public RcppParallel::Worker
30{
31 RcppParallel::RVector <int> dp_fromi;
32 const std::vector <size_t> toi;
33 const std::vector <size_t> edge_type;
34 const size_t nverts;
35 const std::vector <double> vx;
36 const std::vector <double> vy;
37 const std::shared_ptr <DGraph> g;
38 const std::string heap_type;
39 const size_t num_edge_types;
40
41 RcppParallel::RMatrix <double> dout;
42
43 // constructor
44 OneCategoricalDist (
45 const RcppParallel::RVector <int> fromi,
46 const std::vector <size_t> toi_in,
47 const std::vector <size_t> edge_type_in,
48 const size_t nverts_in,
49 const std::vector <double> vx_in,
50 const std::vector <double> vy_in,
51 const std::shared_ptr <DGraph> g_in,
52 const std::string & heap_type_in,
53 const size_t & num_edge_types_in,
54 RcppParallel::RMatrix <double> dout_in) :
55 dp_fromi (fromi), toi (toi_in), edge_type (edge_type_in),
56 nverts (nverts_in), vx (vx_in), vy (vy_in),
57 g (g_in), heap_type (heap_type_in),
58 num_edge_types (num_edge_types_in),
59 dout (dout_in)
60 {
61 }
62
63 // Parallel function operator
64 void operator() (std::size_t begin, std::size_t end)
65 {
66 for (std::size_t i = begin; i < end; i++)
67 {
68 std::shared_ptr<PF::PathFinder> pathfinder =
69 std::make_shared <PF::PathFinder> (nverts,
70 *run_sp::getHeapImpl (heap_type), g);
71 std::vector <double> w (nverts);
72 std::vector <double> d (nverts * (num_edge_types + 1));
73 std::vector <long int> prev (nverts);
74
75 std::vector <double> heuristic (nverts, 0.0);
76
77 size_t from_i = static_cast <size_t> (dp_fromi [i]);
78
79 // only implemented for spatial graphs
80 for (size_t j = 0; j < nverts; j++)
81 {
82 const double dx = vx [j] - vx [from_i],
83 dy = vy [j] - vy [from_i];
84 heuristic [j] = sqrt (dx * dx + dy * dy);
85 }
86 pathfinder->AStarEdgeType (d, w, prev, heuristic, from_i, toi);
87
88 for (size_t j = 0; j < toi.size (); j++)
89 {
90 for (size_t k = 0; k <= num_edge_types; k++) {
91 const double dto = d [toi [j] + k * nverts];
92 if (dto < INFINITE_DOUBLE)
93 {
94 dout (i, j + k * toi.size ()) = dto;
95 }
96 }
97 }
98 }
99 }
100
101};
102
103struct OnePairedCategoricalDist : public RcppParallel::Worker
104{
105 RcppParallel::RVector <int> dp_fromi;
106 const std::vector <size_t> toi;
107 const std::vector <size_t> edge_type;
108 const size_t nverts;
109 const std::vector <double> vx;
110 const std::vector <double> vy;
111 const std::shared_ptr <DGraph> g;
112 const std::string heap_type;
113 const size_t num_edge_types;
114
115 RcppParallel::RMatrix <double> dout;
116
117 // constructor
118 OnePairedCategoricalDist (
119 const RcppParallel::RVector <int> fromi,
120 const std::vector <size_t> toi_in,
121 const std::vector <size_t> edge_type_in,
122 const size_t nverts_in,
123 const std::vector <double> vx_in,
124 const std::vector <double> vy_in,
125 const std::shared_ptr <DGraph> g_in,
126 const std::string & heap_type_in,
127 const size_t & num_edge_types_in,
128 RcppParallel::RMatrix <double> dout_in) :
129 dp_fromi (fromi), toi (toi_in), edge_type (edge_type_in),
130 nverts (nverts_in), vx (vx_in), vy (vy_in),
131 g (g_in), heap_type (heap_type_in),
132 num_edge_types (num_edge_types_in),
133 dout (dout_in)
134 {
135 }
136
137 // Parallel function operator
138 void operator() (std::size_t begin, std::size_t end)
139 {
140 for (std::size_t i = begin; i < end; i++)
141 {
142 std::shared_ptr<PF::PathFinder> pathfinder =
143 std::make_shared <PF::PathFinder> (nverts,
144 *run_sp::getHeapImpl (heap_type), g);
145 std::vector <double> w (nverts);
146 std::vector <double> d (nverts * (num_edge_types + 1));
147 std::vector <long int> prev (nverts);
148
149 std::vector <double> heuristic (nverts, 0.0);
150
151 const size_t from_i = static_cast <size_t> (dp_fromi [i]);
152 const std::vector <size_t> to_i = { toi [i] };
153
154 // only implemented for spatial graphs
155 for (size_t j = 0; j < nverts; j++)
156 {
157 const double dx = vx [j] - vx [from_i],
158 dy = vy [j] - vy [from_i];
159 heuristic [j] = sqrt (dx * dx + dy * dy);
160 }
161 pathfinder->AStarEdgeType (d, w, prev, heuristic, from_i, to_i);
162
163 for (size_t k = 0; k <= num_edge_types; k++)
164 {
165 const double dto = d [to_i [0] + k * nverts];
166 if (dto < INFINITE_DOUBLE)
167 dout (i, k) = dto;
168 }
169 }
170 }
171};
172
173// Modified version of OneCategoricalDist to aggregate only a vector of
174// distances, one value for each `from`.
175struct OneCategory : public RcppParallel::Worker
176{
177 RcppParallel::RVector <int> dp_fromi;
178 const std::vector <size_t> toi;
179 const std::vector <size_t> edge_type;
180 const size_t nverts;
181 const std::vector <double> vx;
182 const std::vector <double> vy;
183 const std::shared_ptr <DGraph> g;
184 const std::string heap_type;
185 const size_t num_edge_types;
186
187 RcppParallel::RMatrix <double> dout;
188
189 // constructor
190 OneCategory (
191 const RcppParallel::RVector <int> fromi,
192 const std::vector <size_t> toi_in,
193 const std::vector <size_t> edge_type_in,
194 const size_t nverts_in,
195 const std::vector <double> vx_in,
196 const std::vector <double> vy_in,
197 const std::shared_ptr <DGraph> g_in,
198 const std::string & heap_type_in,
199 const size_t & num_edge_types_in,
200 RcppParallel::RMatrix <double> dout_in) :
201 dp_fromi (fromi), toi (toi_in), edge_type (edge_type_in),
202 nverts (nverts_in), vx (vx_in), vy (vy_in),
203 g (g_in), heap_type (heap_type_in),
204 num_edge_types (num_edge_types_in),
205 dout (dout_in)
206 {
207 }
208
209 // Parallel function operator
210 void operator() (std::size_t begin, std::size_t end)
211 {
212 for (std::size_t i = begin; i < end; i++)
213 {
214 std::shared_ptr<PF::PathFinder> pathfinder =
215 std::make_shared <PF::PathFinder> (nverts,
216 *run_sp::getHeapImpl (heap_type), g);
217 std::vector <double> w (nverts);
218 std::vector <double> d (nverts * (num_edge_types + 1));
219 std::vector <long int> prev (nverts);
220
221 std::vector <double> heuristic (nverts, 0.0);
222
223 size_t from_i = static_cast <size_t> (dp_fromi [i]);
224
225 // only implemented for spatial graphs
226 for (size_t j = 0; j < nverts; j++)
227 {
228 const double dx = vx [j] - vx [from_i],
229 dy = vy [j] - vy [from_i];
230 heuristic [j] = sqrt (dx * dx + dy * dy);
231 }
232 pathfinder->AStarEdgeType (d, w, prev, heuristic, from_i, toi);
233
234 for (size_t j = 0; j < toi.size (); j++)
235 {
236 for (size_t k = 0; k <= num_edge_types; k++)
237 {
238 const double dto = d [toi [j] + k * nverts];
239 if (dto < INFINITE_DOUBLE) {
240 if (Rcpp::NumericMatrix::is_na (dout (i, k)))
241 dout (i, k) = dto;
242 else
243 dout (i, k) += dto;
244 }
245 }
246 }
247 }
248 }
249
250};
251
252// Modified version of OneCategory to aggregate vector of categorical
253// distances out to a specified threshold, one value for each `from`.
254struct OneCatThreshold : public RcppParallel::Worker
255{
256 RcppParallel::RVector <int> dp_fromi;
257 const std::vector <size_t> edge_type;
258 const size_t nverts;
259 const std::shared_ptr <DGraph> g;
260 const std::string heap_type;
261 const size_t num_edge_types;
262 const double dlimit;
263
264 RcppParallel::RMatrix <double> dout;
265
266 // constructor
267 OneCatThreshold (
268 const RcppParallel::RVector <int> fromi,
269 const std::vector <size_t> edge_type_in,
270 const size_t nverts_in,
271 const std::shared_ptr <DGraph> g_in,
272 const std::string & heap_type_in,
273 const size_t & num_edge_types_in,
274 const double & dlimit_in,
275 RcppParallel::RMatrix <double> dout_in) :
276 dp_fromi (fromi), edge_type (edge_type_in), nverts (nverts_in),
277 g (g_in), heap_type (heap_type_in),
278 num_edge_types (num_edge_types_in), dlimit (dlimit_in),
279 dout (dout_in)
280 {
281 }
282
283 // Parallel function operator
284 void operator() (std::size_t begin, std::size_t end)
285 {
286 for (std::size_t i = begin; i < end; i++)
287 {
288 std::shared_ptr<PF::PathFinder> pathfinder =
289 std::make_shared <PF::PathFinder> (nverts,
290 *run_sp::getHeapImpl (heap_type), g);
291 std::vector <double> w (nverts);
292 std::vector <double> d (nverts * (num_edge_types + 1));
293 std::vector <long int> prev (nverts);
294
295 size_t from_i = static_cast <size_t> (dp_fromi [i]);
296
297 pathfinder->DijkstraLimitEdgeType (d, w, prev, from_i, dlimit);
298
299 // Get the set of terminal vertices: those with w < dlimit_max but
300 // with no previous (outward-going) nodes
301 std::unordered_set <size_t> terminal_verts;
302 for (size_t j = 0; j < nverts; j++)
303 {
304 if (prev [j] == INFINITE_INT && w [j] < dlimit)
305 {
306 terminal_verts.emplace (j);
307 }
308 }
309
310 for (size_t j = 0; j < nverts; j++)
311 {
312 if (w [j] < INFINITE_DOUBLE)
313 {
314 size_t st_prev = static_cast <size_t> (prev [j]);
315
316 bool is_terminal =
317 terminal_verts.find (j) != terminal_verts.end () ||
318 (d [j] > dlimit && d [st_prev] < dlimit);
319
320 if (is_terminal) {
321
322 for (size_t k = 0; k <= num_edge_types; k++)
323 {
324 const double dto = d [j + k * nverts];
325 if (dto < INFINITE_DOUBLE) {
326 if (Rcpp::NumericMatrix::is_na (dout (i, k)))
327 dout (i, k) = dto;
328 else
329 dout (i, k) += dto;
330 }
331 }
332 }
333 }
334 }
335 }
336 }
337
338};
339
340size_t categorical::get_num_edge_types (const std::vector <size_t> &edge_type)
341{
342 std::unordered_set <size_t> type_set;
343 for (auto e: edge_type)
344 if (e != 0L)
345 type_set.emplace (e);
346
347 return type_set.size ();
348}
349
350void PF::PathFinder::AStarEdgeType (std::vector<double>& d,
351 std::vector<double>& w,
352 std::vector<long int>& prev,
353 const std::vector<double>& heur,
354 const size_t v0,
355 const std::vector <size_t> &to_index)
356{
357 // d is already nverts * num_edge_types, and init_arrays only fills without
358 // resizing.
359 const DGraphEdge *edge;
360
361 const size_t nverts = m_graph->nVertices();
362 const std::vector<DGraphVertex>& vertices = m_graph->vertices();
363
364 PF::PathFinder::init_arrays (d, w, prev, m_open, m_closed, v0, nverts);
365 m_heap->insert (v0, heur [v0]);
366 // initialise v0 values for all edge_types
367 const size_t num_edge_types = d.size () / nverts - 1L;
368 for (size_t i = 1; i <= num_edge_types; i++)
369 d [v0 + i * nverts] = 0.0;
370
371 size_t n_reached = 0;
372 const size_t n_targets = to_index.size ();
373 bool *is_target = new bool [nverts];
374 std::fill (is_target, is_target + nverts, false);
375 for (auto t: to_index)
376 is_target [t] = true;
377
378 while (m_heap->nItems() > 0) {
379 size_t v = m_heap->deleteMin();
380
381 m_closed [v] = true;
382 m_open [v] = false;
383
384 edge = vertices [v].outHead;
385 scan_edge_types_heur (edge, d, w, prev, m_open, m_closed, v, heur);
386
387 if (is_target [v])
388 n_reached++;
389 if (n_reached == n_targets)
390 break;
391 } // end while m_heap->nItems
392 delete [] is_target;
393}
394
395
396void PF::PathFinder::scan_edge_types_heur (const DGraphEdge *edge,
397 std::vector<double>& d,
398 std::vector<double>& w,
399 std::vector<long int>& prev,
400 bool *m_open_vec,
401 const bool *m_closed_vec,
402 const size_t &v0,
403 const std::vector<double> &heur) // heuristic for A*
404{
405 const size_t nverts = w.size ();
406 const size_t num_edge_types = d.size () / nverts - 1L;
407
408 while (edge) {
409
410 const size_t et = edge->target;
411 const size_t edge_id = edge->edge_id;
412
413 if (!m_closed_vec [et])
414 {
415 double wt = w [v0] + edge->wt;
416 if (wt < w [et]) {
417 d [et] = d [v0] + edge->dist;
418 // iterate over rest of matrix for each edge_type
419 for (size_t i = 1; i <= num_edge_types; i++) {
420 if (i == edge_id)
421 d [et + i * nverts] = d [v0 + i * nverts] + edge->dist;
422 else // carry over without incrementing dist
423 d [et + i * nverts] = d [v0 + i * nverts];
424 }
425
426 w [et] = wt;
427 prev [et] = static_cast <int> (v0);
428
429 if (m_open_vec [et]) {
430 m_heap->decreaseKey(et, wt + heur [et] - heur [v0]);
431 }
432 else {
433 m_heap->insert (et, wt + heur [et] - heur [v0]);
434 m_open_vec [et] = true;
435 }
436 } else
437 m_closed [et] = true;
438 }
439 edge = edge->nextOut;
440 }
441}
442
443// Modified pathfinder only out to specified distance limit. Done as separate
444// routine to avoid costs of the `if` clause in general non-limited case.
445void PF::PathFinder::DijkstraLimitEdgeType (
446 std::vector<double>& d,
447 std::vector<double>& w,
448 std::vector<long int>& prev,
449 const size_t v0,
450 const double &dlim)
451{
452 const DGraphEdge *edge;
453
454 const size_t n = m_graph->nVertices();
455 const std::vector<DGraphVertex>& vertices = m_graph->vertices();
456
457 PF::PathFinder::init_arrays (d, w, prev, m_open, m_closed, v0, n);
458 m_heap->insert (v0, 0.0);
459 // initialise v0 values for all edge_types
460 const size_t num_edge_types = d.size () / n - 1L;
461 const size_t nverts = w.size ();
462 for (size_t i = 1; i <= num_edge_types; i++)
463 d [v0 + i * nverts] = 0.0;
464
465 while (m_heap->nItems() > 0) {
466 size_t v = m_heap->deleteMin();
467
468 m_closed [v] = true;
469 m_open [v] = false;
470
471 // explore the OUT set of v only if distances are < threshold
472 bool explore = false;
473 edge = vertices [v].outHead;
474 while (edge) {
475 if ((d [v] + edge->dist) <= dlim)
476 {
477 explore = true;
478 break;
479 }
480 edge = edge->nextOut;
481 }
482
483 if (explore)
484 {
485 edge = vertices [v].outHead;
486 scan_edge_types (edge, d, w, prev, m_open, m_closed, v);
487 }
488 } // end while nItems > 0
489}
490
491void PF::PathFinder::scan_edge_types (const DGraphEdge *edge,
492 std::vector<double>& d,
493 std::vector<double>& w,
494 std::vector<long int>& prev,
495 bool *m_open_vec,
496 const bool *m_closed_vec,
497 const size_t &v0)
498{
499 const size_t n = w.size ();
500 const size_t num_edge_types = d.size () / n - 1L;
501
502 while (edge) {
503
504 const size_t et = edge->target;
505 const size_t edge_id = edge->edge_id;
506
507 if (!m_closed_vec [et])
508 {
509 double wt = w [v0] + edge->wt;
510 if (wt < w [et]) {
511 d [et] = d [v0] + edge->dist;
512 // iterate over rest of matrix for each edge_type
513 for (size_t i = 1; i <= num_edge_types; i++) {
514 if (i == edge_id)
515 d [et + edge_id * n] = d [v0 + edge_id * n] + edge->dist;
516 else // carry over without incrementing dist
517 d [et + i * n] = d [v0 + i * n];
518 }
519
520 w [et] = wt;
521 prev [et] = static_cast <int> (v0);
522
523 if (m_open_vec [et]) {
524 m_heap->decreaseKey(et, wt);
525 }
526 else {
527 m_heap->insert (et, wt);
528 m_open_vec [et] = true;
529 }
530 } else
531 m_closed [et] = true;
532 }
533 edge = edge->nextOut;
534 }
535}
536
537//' rcpp_get_sp_dists_categorical
538//'
539//' The `graph` must have an `edge_type` column of non-negative integers,
540//' with 0 denoting edges which are not aggregated, and all other values
541//' defining aggregation categories.
542//'
543//' Implemented in parallal form only; no single-threaded version, and
544//' only for AStar (so graphs must be spatial).
545//' @noRd
546// [[Rcpp::export]]
547Rcpp::NumericMatrix rcpp_get_sp_dists_categorical (const Rcpp::DataFrame graph,
548 const Rcpp::DataFrame vert_map_in,
549 Rcpp::IntegerVector fromi,
550 Rcpp::IntegerVector toi_in,
551 const std::string& heap_type,
552 const bool proportions_only)
553{
554 std::vector <size_t> toi =
555 Rcpp::as <std::vector <size_t> > ( toi_in);
556
557 size_t nfrom = static_cast <size_t> (fromi.size ());
558 size_t nto = static_cast <size_t> (toi.size ());
559
560 const std::vector <std::string> from = graph ["from"];
561 const std::vector <std::string> to = graph ["to"];
562 const std::vector <double> dist = graph ["d"];
563 const std::vector <double> wt = graph ["d_weighted"];
564 const std::vector <size_t> edge_type = graph ["edge_type"];
565
566 const size_t num_edge_types = categorical::get_num_edge_types (edge_type);
567
568 const size_t nedges = static_cast <size_t> (graph.nrow ());
569 std::map <std::string, size_t> vert_map;
570 std::vector <std::string> vert_map_id = vert_map_in ["vert"];
571 std::vector <size_t> vert_map_n = vert_map_in ["id"];
572 const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
573 vert_map_n, vert_map);
574
575 std::vector <double> vx (nverts), vy (nverts);
576 vx = Rcpp::as <std::vector <double> > (vert_map_in ["x"]);
577 vy = Rcpp::as <std::vector <double> > (vert_map_in ["y"]);
578
579 std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
580 inst_graph (g, nedges, vert_map, from, to, edge_type, dist, wt);
581
582 size_t ncol;
583 if (proportions_only)
584 ncol = num_edge_types + 1L;
585 else
586 ncol = nto * (num_edge_types + 1L);
587
588 Rcpp::NumericVector na_vec = Rcpp::NumericVector (nfrom * ncol,
589 Rcpp::NumericVector::get_na ());
590 Rcpp::NumericMatrix dout (static_cast <int> (nfrom),
591 static_cast <int> (ncol));
592
593 // Create parallel worker
594 size_t chunk_size = run_sp::get_chunk_size (nfrom);
595 if (proportions_only)
596 {
597 OneCategory one_dist (RcppParallel::RVector <int> (fromi), toi,
598 edge_type, nverts, vx, vy,
599 g, heap_type, num_edge_types,
600 RcppParallel::RMatrix <double> (dout));
601 RcppParallel::parallelFor (0, nfrom, one_dist, chunk_size);
602 } else
603 {
604 OneCategoricalDist one_dist (RcppParallel::RVector <int> (fromi),
605 toi, edge_type, nverts, vx, vy,
606 g, heap_type, num_edge_types,
607 RcppParallel::RMatrix <double> (dout));
608 RcppParallel::parallelFor (0, nfrom, one_dist, chunk_size);
609 }
610
611 return (dout);
612}
613
614//' rcpp_get_sp_dists_categ_paired
615//'
616//' Pairwise version of 'get_sp_dists_categorical'. The `graph` must have an
617//'`edge_type` column of non-negative integers, with 0 denoting edges which are
618//' not aggregated, and all other values defining aggregation categories.
619//'
620//' Implemented in parallal form only; no single-threaded version, and
621//' only for AStar (so graphs must be spatial).
622//' @noRd
623// [[Rcpp::export]]
624Rcpp::NumericMatrix rcpp_get_sp_dists_categ_paired (const Rcpp::DataFrame graph,
625 const Rcpp::DataFrame vert_map_in,
626 Rcpp::IntegerVector fromi,
627 Rcpp::IntegerVector toi_in,
628 const std::string& heap_type)
629{
630 std::vector <size_t> toi =
631 Rcpp::as <std::vector <size_t> > ( toi_in);
632
633 size_t nfrom = static_cast <size_t> (fromi.size ());
634 size_t nto = static_cast <size_t> (toi.size ());
635 if (nfrom != nto)
636 {
637 Rcpp::stop ("Pairwise categorical dists require equal numbers of 'from' and 'to' points.");
638 }
639
640 const std::vector <std::string> from = graph ["from"];
641 const std::vector <std::string> to = graph ["to"];
642 const std::vector <double> dist = graph ["d"];
643 const std::vector <double> wt = graph ["d_weighted"];
644 const std::vector <size_t> edge_type = graph ["edge_type"];
645
646 const size_t num_types = categorical::get_num_edge_types (edge_type);
647
648 const size_t nedges = static_cast <size_t> (graph.nrow ());
649 std::map <std::string, size_t> vert_map;
650 std::vector <std::string> vert_map_id = vert_map_in ["vert"];
651 std::vector <size_t> vert_map_n = vert_map_in ["id"];
652 const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
653 vert_map_n, vert_map);
654
655 std::vector <double> vx (nverts), vy (nverts);
656 vx = Rcpp::as <std::vector <double> > (vert_map_in ["x"]);
657 vy = Rcpp::as <std::vector <double> > (vert_map_in ["y"]);
658
659 std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
660 inst_graph (g, nedges, vert_map, from, to, edge_type, dist, wt);
661
662 const size_t ncol = num_types + 1L;
663
664 Rcpp::NumericVector na_vec = Rcpp::NumericVector (nfrom * ncol,
665 Rcpp::NumericVector::get_na ());
666 Rcpp::NumericMatrix dout (static_cast <int> (nfrom),
667 static_cast <int> (ncol));
668
669 // Create parallel worker
670 size_t chunk_size = run_sp::get_chunk_size (nfrom);
671 OnePairedCategoricalDist one_dist (RcppParallel::RVector <int> (fromi),
672 toi, edge_type, nverts, vx, vy,
673 g, heap_type, num_types,
674 RcppParallel::RMatrix <double> (dout));
675 RcppParallel::parallelFor (0, nfrom, one_dist, chunk_size);
676
677 return (dout);
678}
679
680//' rcpp_get_sp_dists_cat_threshold
681//'
682//' The `graph` must have an `edge_type` column of non-negative integers,
683//' with 0 denoting edges which are not aggregated, and all other values
684//' defining aggregation categories.
685//'
686//' Implemented in parallal form only; no single-threaded version, and
687//' only for AStar (so graphs must be spatial).
688//' @noRd
689// [[Rcpp::export]]
690Rcpp::NumericMatrix rcpp_get_sp_dists_cat_threshold (const Rcpp::DataFrame graph,
691 const Rcpp::DataFrame vert_map_in,
692 Rcpp::IntegerVector fromi,
693 const double dlimit,
694 const std::string& heap_type)
695{
696 size_t nfrom = static_cast <size_t> (fromi.size ());
697
698 const std::vector <std::string> from = graph ["from"];
699 const std::vector <std::string> to = graph ["to"];
700 const std::vector <double> dist = graph ["d"];
701 const std::vector <double> wt = graph ["d_weighted"];
702 const std::vector <size_t> edge_type = graph ["edge_type"];
703
704 const size_t num_types = categorical::get_num_edge_types (edge_type);
705
706 const size_t nedges = static_cast <size_t> (graph.nrow ());
707 std::map <std::string, size_t> vert_map;
708 std::vector <std::string> vert_map_id = vert_map_in ["vert"];
709 std::vector <size_t> vert_map_n = vert_map_in ["id"];
710 const size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
711 vert_map_n, vert_map);
712
713 std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
714 inst_graph (g, nedges, vert_map, from, to, edge_type, dist, wt);
715
716 const size_t ncol = num_types + 1L;
717
718 Rcpp::NumericVector na_vec = Rcpp::NumericVector (nfrom * ncol,
719 Rcpp::NumericVector::get_na ());
720 Rcpp::NumericMatrix dout (static_cast <int> (nfrom),
721 static_cast <int> (ncol));
722
723 // Create parallel worker
724 size_t chunk_size = run_sp::get_chunk_size (nfrom);
725 OneCatThreshold one_dist (RcppParallel::RVector <int> (fromi),
726 edge_type, nverts, g, heap_type,
727 num_types, dlimit,
728 RcppParallel::RMatrix <double> (dout));
729 RcppParallel::parallelFor (0, nfrom, one_dist, chunk_size);
730
731 return (dout);
732}