many-to-many multi-modal routing aggregator
1
2#include "multi-modal.h"
3
4//' @param t_net_to_gtfs Reduced matrix of times from each "from" point to
5//' subset of GTFS points within nominated time limit.
6//' @param t_gtfs_to_gtfs Reduced subset of full GTFS travel times matrix, with
7//' numbers of rows reduced to only those stops reachable from each "from"
8//' point, but retaining all destination stops.
9//' @noRd
10struct OneTimesToEndStops : public RcppParallel::Worker
11{
12 const RcppParallel::RMatrix <int> t_net_to_gtfs;
13 const Rcpp::IntegerMatrix t_gtfs_to_gtfs;
14 const size_t nfrom;
15 const size_t n_gtfs_start;
16 const size_t n_gtfs_total;
17 // t_gtfs_to_gtfs is actually 3 matrices of (times, transfers, intervals)
18 // dim (n_gtfs_start, 3 * n_gtfs_total)
19
20 RcppParallel::RMatrix <int> tout;
21
22 // constructor
23 OneTimesToEndStops (
24 const RcppParallel::RMatrix <int> t_net_to_gtfs_in,
25 const Rcpp::IntegerMatrix t_gtfs_to_gtfs_in,
26 const size_t nfrom_in,
27 const size_t n_gtfs_start_in,
28 const size_t n_gtfs_total_in,
29 RcppParallel::RMatrix <int> tout_in) :
30 t_net_to_gtfs (t_net_to_gtfs_in), t_gtfs_to_gtfs (t_gtfs_to_gtfs_in),
31 nfrom (nfrom_in), n_gtfs_start (n_gtfs_start_in),
32 n_gtfs_total (n_gtfs_total_in), tout (tout_in)
33 {
34 }
35
36 // Parallel function operator.
37 //
38 // This post from JJ Aliare suggests that iteration can be implemented
39 // directly over rows of RMatrix objects, and not just over elements:
40 // https://gallery.rcpp.org/articles/parallel-distance-matrix/
41 void operator() (std::size_t begin, std::size_t end)
42 {
43 // i over all vertices in the input distance matrix
44 for (std::size_t i = begin; i < end; i++)
45 {
46 RcppParallel::RMatrix <int>::Row times_row = t_net_to_gtfs.row (i);
47 // Length of times_row is equal to number of GTFS stops reachable
48 // from each "from" point in < duration_max. Dim of 't_gtfs_to_gtfs'
49 // is also reduced to (n_reachable, n_total), so first dimension can
50 // be indexed directly with "j":
51
52 for (size_t j = 0; j < times_row.size (); j++)
53 {
54 // Time from network point to GTFS stop:
55 const int time_i_to_j = times_row [j];
56 if (time_i_to_j < 0) // NA's have been replaced with -ve
57 {
58 continue;
59 }
60
61 for (size_t k = 0; k < n_gtfs_total; k++)
62 {
63 const int time_j_to_k = t_gtfs_to_gtfs (j, k);
64 if (time_j_to_k < 0)
65 {
66 continue;
67 }
68
69 int time_i_to_k = time_i_to_j + time_j_to_k;
70 if (time_i_to_k < tout (i, k))
71 {
72 tout (i, k) = time_i_to_k;
73 tout (i, k + n_gtfs_total) = t_gtfs_to_gtfs (j, k + n_gtfs_total); // ntransfers
74 tout (i, k + 2 * n_gtfs_total) = t_gtfs_to_gtfs (j, k + 2 * n_gtfs_total); // interval
75 }
76 }
77 }
78 }
79 }
80};
81
82/* Both 'times_to_end_stops' and 'tout' are tripe-matrices, with successive
83 * thirds stored in columns for:
84 * 1. times;
85 * 2. transfers;
86 * 3. intervals
87 */
88struct AddTwoMatricesWorker : public RcppParallel::Worker
89{
90 const RcppParallel::RMatrix <int> times_to_end_stops;
91 const std::vector <std::vector <size_t> > gtfs_to_net_index_vec;
92 const std::vector <std::vector <double> > gtfs_to_net_dist_vec;
93 const size_t nfrom;
94
95 RcppParallel::RMatrix <int> tout;
96
97 // constructor
98 AddTwoMatricesWorker (
99 const RcppParallel::RMatrix <int> times_to_end_stops_in,
100 std::vector <std::vector <size_t> > gtfs_to_net_index_vec_in,
101 std::vector <std::vector <double> > gtfs_to_net_dist_vec_in,
102 const size_t nfrom_in,
103 RcppParallel::RMatrix <int> tout_in) :
104 times_to_end_stops (times_to_end_stops_in),
105 gtfs_to_net_index_vec (gtfs_to_net_index_vec_in),
106 gtfs_to_net_dist_vec (gtfs_to_net_dist_vec_in),
107 nfrom (nfrom_in),
108 tout (tout_in)
109 {
110 }
111
112 // Parallel function operator
113 void operator() (std::size_t begin, std::size_t end)
114 {
115 const size_t n_gtfs_total = gtfs_to_net_dist_vec.size ();
116 const size_t n_verts = static_cast <int> (tout.ncol () / 3);
117 // i over all vertices in the input distance matrix
118 for (std::size_t i = begin; i < end; i++)
119 {
120 for (size_t j = 0; j < n_gtfs_total; j++)
121 {
122 if (times_to_end_stops (i, j) == INFINITE_INT)
123 {
124 continue;
125 }
126
127 const size_t n_j = gtfs_to_net_index_vec [j].size ();
128
129 for (size_t k = 0; k < n_j; k++)
130 {
131 if (gtfs_to_net_dist_vec [j] [k] < 0)
132 {
133 continue;
134 }
135
136 const int dist_j_k = static_cast <int> (round (gtfs_to_net_dist_vec [j] [k]));
137 const int time_i_to_k = times_to_end_stops (i, j) + dist_j_k;
138 const size_t index_k = gtfs_to_net_index_vec [j] [k];
139
140 if (time_i_to_k < tout (i, index_k))
141 {
142 tout (i, index_k) = time_i_to_k;
143 tout (i, n_verts + index_k) = times_to_end_stops (i, j + n_gtfs_total); // transfers
144 tout (i, 2 * n_verts + index_k) = times_to_end_stops (i, j + 2 * n_gtfs_total); // interval
145 }
146 }
147 }
148 }
149 }
150};
151
152//' rcpp_add_net_to_gtfs
153//'
154//' Add times from selected start points to all GTFS stations to total GTFS '
155//' travel time matrix to generate fastest travel times to all GTFS end points.
156//' @noRd
157// [[Rcpp::export]]
158Rcpp::IntegerMatrix rcpp_add_net_to_gtfs (Rcpp::IntegerMatrix t_net_to_gtfs,
159 Rcpp::IntegerMatrix gtfs_times, Rcpp::List gtfs_to_net_index,
160 Rcpp::List gtfs_to_net_dist, const int nverts)
161{
162 const size_t nfrom = static_cast <size_t> (t_net_to_gtfs.nrow ());
163 const size_t n_gtfs_start = static_cast <size_t> (t_net_to_gtfs.ncol ());
164 const size_t n_gtfs_total = static_cast <size_t> (gtfs_times.ncol () / 3);
165
166 if (t_net_to_gtfs.ncol () != gtfs_times.nrow ())
167 {
168 Rcpp::stop ("GTFS and travel time matrices are not compatible");
169 }
170
171 // times_to_end_stops holds 3 matrices of:
172 // 1. Actual times;
173 // 2. Numbers of transfers;
174 // 3. Effective intervals to next service.
175 Rcpp::IntegerMatrix times_to_end_stops (
176 static_cast <int> (nfrom),
177 static_cast <int> (3 * n_gtfs_total)
178 );
179 std::fill (times_to_end_stops.begin (), times_to_end_stops.end (), INFINITE_INT);
180
181 // Add initial times to all closest GTFS stops to the times to all terminal
182 // GTFS stops:
183 OneTimesToEndStops one_times (RcppParallel::RMatrix <int> (t_net_to_gtfs),
184 gtfs_times, nfrom, n_gtfs_start, n_gtfs_total,
185 RcppParallel::RMatrix <int> (times_to_end_stops));
186
187 RcppParallel::parallelFor (0, nfrom, one_times);
188
189 // Those are then the fastest times to each terminal GTFS stop. Then use the
190 // lists of closest stops to each network point to calculate final fastest
191 // times to each terminal network point.
192
193 Rcpp::IntegerMatrix res (static_cast <int> (nfrom), 3 * nverts);
194 std::fill (res.begin (), res.end (), INFINITE_INT);
195
196 // convert the Rcpp::Lists into std::vecs:
197 std::vector <std::vector <size_t> > gtfs_to_net_index_vec (gtfs_to_net_index.size ());
198 std::vector <std::vector <double> > gtfs_to_net_dist_vec (gtfs_to_net_dist.size ());
199
200 for (size_t i = 0; i < gtfs_to_net_index.size (); i++)
201 {
202 Rcpp::IntegerVector index_i = gtfs_to_net_index (i);
203 gtfs_to_net_index_vec [i] = Rcpp::as <std::vector <size_t> > (index_i);
204
205 Rcpp::NumericVector d_i = gtfs_to_net_dist (i);
206 gtfs_to_net_dist_vec [i] = Rcpp::as <std::vector <double> > (d_i);
207 }
208
209 AddTwoMatricesWorker combine_two_mats (
210 RcppParallel::RMatrix <int> (times_to_end_stops),
211 gtfs_to_net_index_vec, gtfs_to_net_dist_vec, nfrom,
212 RcppParallel::RMatrix <int> (res));
213
214 RcppParallel::parallelFor (0, nfrom, combine_two_mats);
215
216 return res;
217}
218
219//' rcpp_min_from_two_matrices
220//'
221//' Iterate over all rows and columns of two identical matrices, and return
222//' minimal values from same indices in both.
223//' @noRd
224// [[Rcpp::export]]
225Rcpp::IntegerMatrix rcpp_min_from_two_matrices (Rcpp::IntegerMatrix mat1,
226 Rcpp::IntegerMatrix mat2)
227{
228 if (mat1.ncol () != mat2.ncol () || mat1.nrow () != mat2.nrow ())
229 {
230 Rcpp::stop ("Matrices must have identical dimensions.");
231 }
232
233 Rcpp::IntegerMatrix res (mat1.nrow (), mat1.ncol ());
234
235 for (size_t i = 0; i < mat1.nrow (); i++)
236 {
237 for (size_t j = 0; j < mat1.ncol (); j++)
238 {
239 res (i, j) = std::min (mat1 (i, j), mat2 (i, j));
240 }
241 }
242
243 return res;
244}