Skip to main content

solverforge_maps/routing/
matrix.rs

1//! Matrix computation for road networks.
2
3use rayon::prelude::*;
4use std::collections::HashMap;
5use tokio::sync::mpsc::{self, Sender};
6
7use super::algo::dijkstra;
8use super::coord::Coord;
9use super::graph::NodeIdx;
10use super::network::{RoadNetwork, SnappedCoord};
11use super::progress::RoutingProgress;
12
13pub const UNREACHABLE: i64 = i64::MAX;
14
15#[derive(Debug, Clone, Default)]
16pub struct TravelTimeMatrix {
17    data: Vec<i64>,
18    size: usize,
19    locations: Vec<SnappedCoord>,
20}
21
22impl TravelTimeMatrix {
23    pub(crate) fn new(data: Vec<i64>, size: usize, locations: Vec<SnappedCoord>) -> Self {
24        debug_assert_eq!(data.len(), size * size);
25        debug_assert_eq!(locations.len(), size);
26        Self {
27            data,
28            size,
29            locations,
30        }
31    }
32
33    /// Get the travel time from one location to another.
34    ///
35    /// Returns `None` if indices are out of bounds.
36    /// Returns `Some(UNREACHABLE)` if the pair is not reachable.
37    #[inline]
38    pub fn get(&self, from: usize, to: usize) -> Option<i64> {
39        if from < self.size && to < self.size {
40            Some(self.data[from * self.size + to])
41        } else {
42            None
43        }
44    }
45
46    #[inline]
47    pub fn is_reachable(&self, from: usize, to: usize) -> bool {
48        self.get(from, to)
49            .map(|v| v != UNREACHABLE)
50            .unwrap_or(false)
51    }
52
53    #[inline]
54    pub fn size(&self) -> usize {
55        self.size
56    }
57
58    #[inline]
59    pub fn locations(&self) -> &[SnappedCoord] {
60        &self.locations
61    }
62
63    #[inline]
64    pub fn row(&self, i: usize) -> Option<&[i64]> {
65        if i < self.size {
66            let start = i * self.size;
67            Some(&self.data[start..start + self.size])
68        } else {
69            None
70        }
71    }
72
73    pub fn min(&self) -> Option<i64> {
74        let mut min_val = None;
75        for i in 0..self.size {
76            for j in 0..self.size {
77                if i == j {
78                    continue;
79                }
80                let val = self.data[i * self.size + j];
81                if val != UNREACHABLE {
82                    min_val = Some(min_val.map(|m: i64| m.min(val)).unwrap_or(val));
83                }
84            }
85        }
86        min_val
87    }
88
89    pub fn max(&self) -> Option<i64> {
90        let mut max_val = None;
91        for i in 0..self.size {
92            for j in 0..self.size {
93                if i == j {
94                    continue;
95                }
96                let val = self.data[i * self.size + j];
97                if val != UNREACHABLE {
98                    max_val = Some(max_val.map(|m: i64| m.max(val)).unwrap_or(val));
99                }
100            }
101        }
102        max_val
103    }
104
105    pub fn mean(&self) -> Option<f64> {
106        let mut sum = 0i64;
107        let mut count = 0usize;
108        for i in 0..self.size {
109            for j in 0..self.size {
110                if i == j {
111                    continue;
112                }
113                let val = self.data[i * self.size + j];
114                if val != UNREACHABLE {
115                    sum = sum.saturating_add(val);
116                    count += 1;
117                }
118            }
119        }
120        if count > 0 {
121            Some(sum as f64 / count as f64)
122        } else {
123            None
124        }
125    }
126
127    pub fn reachability_ratio(&self) -> f64 {
128        if self.size <= 1 {
129            return 1.0;
130        }
131        let total_pairs = self.size * (self.size - 1);
132        let reachable = self.count_reachable();
133        reachable as f64 / total_pairs as f64
134    }
135
136    fn count_reachable(&self) -> usize {
137        let mut count = 0;
138        for i in 0..self.size {
139            for j in 0..self.size {
140                if i != j && self.data[i * self.size + j] != UNREACHABLE {
141                    count += 1;
142                }
143            }
144        }
145        count
146    }
147
148    pub fn unreachable_pairs(&self) -> Vec<(usize, usize)> {
149        let mut pairs = Vec::new();
150        for i in 0..self.size {
151            for j in 0..self.size {
152                if i != j && self.data[i * self.size + j] == UNREACHABLE {
153                    pairs.push((i, j));
154                }
155            }
156        }
157        pairs
158    }
159
160    pub fn as_slice(&self) -> &[i64] {
161        &self.data
162    }
163}
164
165impl RoadNetwork {
166    pub async fn compute_matrix(
167        &self,
168        locations: &[Coord],
169        progress: Option<&Sender<RoutingProgress>>,
170    ) -> TravelTimeMatrix {
171        let n = locations.len();
172        if n == 0 {
173            return TravelTimeMatrix::new(vec![], 0, vec![]);
174        }
175
176        let node_snapped: Vec<Option<SnappedCoord>> = locations
177            .iter()
178            .map(|&coord| self.snap_to_road_detailed(coord).ok())
179            .collect();
180
181        let snapped_locations: Vec<SnappedCoord> = locations
182            .iter()
183            .zip(node_snapped.iter())
184            .map(|(&coord, edge_snap)| {
185                if let Some(snapped) = edge_snap {
186                    *snapped
187                } else {
188                    SnappedCoord {
189                        original: coord,
190                        snapped: coord,
191                        snap_distance_m: f64::INFINITY,
192                        node_index: NodeIdx::new(0),
193                    }
194                }
195            })
196            .collect();
197
198        let row_progress = progress.map(|tx| {
199            let (progress_tx, mut progress_rx) = mpsc::unbounded_channel::<usize>();
200            let tx = tx.clone();
201            let handle = tokio::spawn(async move {
202                let mut completed_rows = 0usize;
203                while progress_rx.recv().await.is_some() {
204                    completed_rows += 1;
205                    let percent = 55 + ((completed_rows * 44) / n.max(1)) as u8;
206                    let _ = tx
207                        .send(RoutingProgress::ComputingMatrix {
208                            percent,
209                            row: completed_rows,
210                            total: n,
211                        })
212                        .await;
213                }
214            });
215            (progress_tx, handle)
216        });
217
218        // Compute rows in parallel via rayon - each row runs Dijkstra from source endpoints
219        let graph = &self.graph;
220        let progress_tx = row_progress.as_ref().map(|(tx, _)| tx.clone());
221        let rows: Vec<Vec<i64>> = (0..n)
222            .into_par_iter()
223            .map(|i| {
224                let mut row = vec![0i64; n];
225
226                let Some(from) = &node_snapped[i] else {
227                    for (j, cell) in row.iter_mut().enumerate() {
228                        if i != j {
229                            *cell = UNREACHABLE;
230                        }
231                    }
232                    if let Some(tx) = &progress_tx {
233                        let _ = tx.send(i);
234                    }
235                    return row;
236                };
237
238                let costs = dijkstra(graph, from.node_index, None, |e| e.travel_time_s);
239
240                for j in 0..n {
241                    if i == j {
242                        continue;
243                    }
244
245                    let Some(to) = &node_snapped[j] else {
246                        row[j] = UNREACHABLE;
247                        continue;
248                    };
249
250                    row[j] = if let Some(&best) = costs.get(&to.node_index) {
251                        best.round() as i64
252                    } else {
253                        UNREACHABLE
254                    };
255                }
256
257                if let Some(tx) = &progress_tx {
258                    let _ = tx.send(i);
259                }
260                row
261            })
262            .collect();
263
264        if let Some((progress_tx, handle)) = row_progress {
265            drop(progress_tx);
266            let _ = handle.await;
267        }
268
269        if let Some(tx) = progress {
270            let _ = tx.send(RoutingProgress::Complete).await;
271        }
272
273        let data: Vec<i64> = rows.into_iter().flatten().collect();
274        TravelTimeMatrix::new(data, n, snapped_locations)
275    }
276
277    pub async fn compute_geometries(
278        &self,
279        locations: &[Coord],
280        progress: Option<&Sender<RoutingProgress>>,
281    ) -> HashMap<(usize, usize), Vec<Coord>> {
282        let n = locations.len();
283        let total_pairs = n * (n - 1);
284        let mut geometries = HashMap::new();
285        let mut pair_count = 0usize;
286
287        for i in 0..n {
288            for j in 0..n {
289                if i == j {
290                    continue;
291                }
292                if let Ok(result) = self.route(locations[i], locations[j]) {
293                    geometries.insert((i, j), result.geometry);
294                }
295                pair_count += 1;
296
297                if let Some(tx) = progress {
298                    if pair_count.is_multiple_of(10) || pair_count == total_pairs {
299                        let percent = 80 + (pair_count * 15 / total_pairs.max(1)) as u8;
300                        let _ = tx
301                            .send(RoutingProgress::ComputingGeometries {
302                                percent,
303                                pair: pair_count,
304                                total: total_pairs,
305                            })
306                            .await;
307                    }
308                }
309            }
310        }
311
312        if let Some(tx) = progress {
313            let _ = tx.send(RoutingProgress::Complete).await;
314        }
315
316        geometries
317    }
318}