many-to-many multi-modal routing aggregator
at main 244 lines 9.2 kB view raw
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}