mpadge.tngl.sh
Routing and analysis engine for GTFS (General Transit Feed Specification) data
1#include "transfers.h"
2
3//' Haversine for variable x and y
4//'
5//' @return single distance
6//'
7//' @note The sxd and syd values could be calculated in arrays, each value of
8//' which could be determined with only n operations, rather than the n2 used
9//' here. Doing so, however, requires very large C arrays which are often
10//' problematic, so this is safer.
11//'
12//' @noRd
13double transfers::one_haversine (const double &x1, const double &y1,
14 const double &x2, const double &y2,
15 const double &cosy1, const double &cosy2)
16{
17 double sxd = sin ((x2 - x1) * M_PI / 360.0);
18 double syd = sin ((y2 - y1) * M_PI / 360.0);
19 double d = syd * syd + cosy1 * cosy2 * sxd * sxd;
20 d = 2.0 * earth * asin (sqrt (d));
21 return (d);
22}
23
24
25// [[Rcpp::export]]
26Rcpp::List rcpp_transfer_nbs (Rcpp::DataFrame stops,
27 const double dlim)
28{
29 const size_t n = static_cast <size_t> (stops.nrow ());
30
31 const std::vector <double> stop_x = stops ["stop_lon"];
32 const std::vector <double> stop_y = stops ["stop_lat"];
33
34 Rcpp::List res (n * 2);
35 // Initialize res:
36 for (size_t i = 0; i < n; i++)
37 {
38 res (i) = std::vector <size_t> ();
39 res (i + n) = std::vector <double> ();
40 }
41
42 for (size_t i = 0; i < (n - 1); i++)
43 {
44 const double cosy1 = cos (stop_y [i] * pi / 180.0);
45
46 std::vector <size_t> index;
47 std::vector <double> dist;
48 for (size_t j = (i + 1); j < n; j++)
49 {
50 const double cosy2 = cos (stop_y [j] * pi / 180.0);
51 const double d_j = transfers::one_haversine (stop_x [i], stop_y [i],
52 stop_x [j], stop_y [j], cosy1, cosy2);
53 if (d_j <= dlim)
54 {
55 index.push_back (j);
56 dist.push_back (d_j);
57 }
58 } // end for j
59
60 // Each vector in 'res' needs to be extended by these new values:
61 std::vector <size_t> index_i = res (i);
62 std::vector <double> dist_i = res (n + i);
63 for (size_t j = 0; j < index.size (); j++)
64 {
65 index_i.push_back (index [j]);
66 dist_i.push_back (dist [j]);
67 }
68
69 res (i) = index_i;
70 res (n + i) = dist_i;
71
72 // But the 'j' values are only incrementally obtained, so have to be
73 // expanded each time. Each element of 'index' has to be mapped back to
74 // the 'res' entry, which needs an 'i' appended. Each element of 'dist'
75 // maps back to 'res(index[j])`, and has to have 'd(i,j)' appended.
76 for (size_t j = 0; j < index.size (); j++)
77 {
78 std::vector <size_t> vec_j = res (index [j]);
79 vec_j.push_back (i);
80 res (index [j]) = vec_j;
81
82 std::vector <double> dist_j = res (index [j] + n);
83 dist_j.push_back (dist [j]);
84 res (index [j] + n) = dist_j;
85 }
86
87 index.clear ();
88 dist.clear ();
89 }
90
91 // Then increment all 'index' values by 1 for 1-based R:
92 for (size_t i = 0; i < n; i++)
93 {
94 std::vector <size_t> vec_i = res (i);
95 for (size_t j = 0; j < vec_i.size (); j++)
96 {
97 vec_i [j] = vec_i [j] + 1;
98 }
99 res (i) = vec_i;
100 }
101
102 return res;
103}