|
| 1 | +//! Miscellaneous algorithms. |
| 2 | + |
| 3 | +/// A comparator on partially ordered elements, that panics if they are incomparable |
| 4 | +pub fn asserting_cmp<T: PartialOrd>(a: &T, b: &T) -> std::cmp::Ordering { |
| 5 | + a.partial_cmp(b).expect("Comparing incomparable elements") |
| 6 | +} |
| 7 | + |
| 8 | +/// Assuming slice is totally ordered and sorted, returns the minimum i for which |
| 9 | +/// slice[i] >= key, or slice.len() if no such i exists |
| 10 | +pub fn slice_lower_bound<T: PartialOrd>(slice: &[T], key: &T) -> usize { |
| 11 | + slice |
| 12 | + .binary_search_by(|x| asserting_cmp(x, key).then(std::cmp::Ordering::Greater)) |
| 13 | + .unwrap_err() |
| 14 | +} |
| 15 | + |
| 16 | +/// Assuming slice is totally ordered and sorted, returns the minimum i for which |
| 17 | +/// slice[i] > key, or slice.len() if no such i exists |
| 18 | +pub fn slice_upper_bound<T: PartialOrd>(slice: &[T], key: &T) -> usize { |
| 19 | + slice |
| 20 | + .binary_search_by(|x| asserting_cmp(x, key).then(std::cmp::Ordering::Less)) |
| 21 | + .unwrap_err() |
| 22 | +} |
| 23 | + |
| 24 | +/// A simple data structure for coordinate compression |
| 25 | +pub struct SparseIndex { |
| 26 | + coords: Vec<i64>, |
| 27 | +} |
| 28 | + |
| 29 | +impl SparseIndex { |
| 30 | + /// Build an index, given the full set of coordinates to compress. |
| 31 | + pub fn new(mut coords: Vec<i64>) -> Self { |
| 32 | + coords.sort_unstable(); |
| 33 | + coords.dedup(); |
| 34 | + Self { coords } |
| 35 | + } |
| 36 | + |
| 37 | + /// Return Ok(i) if the coordinate q appears at index i |
| 38 | + /// Return Err(i) if q appears between indices i-1 and i |
| 39 | + pub fn compress(&self, q: i64) -> Result<usize, usize> { |
| 40 | + self.coords.binary_search(&q) |
| 41 | + } |
| 42 | +} |
| 43 | + |
| 44 | +/// Represents a minimum (lower envelope) of a collection of linear functions of a variable, |
| 45 | +/// evaluated using the convex hull trick with square root decomposition. |
| 46 | +pub struct PiecewiseLinearFn { |
| 47 | + sorted_lines: Vec<(f64, f64)>, |
| 48 | + intersections: Vec<f64>, |
| 49 | + recent_lines: Vec<(f64, f64)>, |
| 50 | + merge_threshold: usize, |
| 51 | +} |
| 52 | + |
| 53 | +impl PiecewiseLinearFn { |
| 54 | + /// For N inserts interleaved with Q queries, a threshold of N/sqrt(Q) yields |
| 55 | + /// O(N sqrt Q + Q log N) time complexity. If all queries come after all inserts, |
| 56 | + /// any threshold less than N (e.g., 0) yields O(N + Q log N) time complexity. |
| 57 | + pub fn with_merge_threshold(merge_threshold: usize) -> Self { |
| 58 | + Self { |
| 59 | + sorted_lines: vec![], |
| 60 | + intersections: vec![], |
| 61 | + recent_lines: vec![], |
| 62 | + merge_threshold, |
| 63 | + } |
| 64 | + } |
| 65 | + |
| 66 | + /// Replaces the represented function with the minimum of itself and a provided line |
| 67 | + pub fn min_with(&mut self, slope: f64, intercept: f64) { |
| 68 | + self.recent_lines.push((slope, intercept)); |
| 69 | + } |
| 70 | + |
| 71 | + fn update_envelope(&mut self) { |
| 72 | + self.recent_lines.extend(self.sorted_lines.drain(..)); |
| 73 | + self.recent_lines.sort_unstable_by(asserting_cmp); |
| 74 | + self.intersections.clear(); |
| 75 | + |
| 76 | + for (new_m, new_b) in self.recent_lines.drain(..).rev() { |
| 77 | + while let Some(&(last_m, last_b)) = self.sorted_lines.last() { |
| 78 | + // If slopes are equal, get rid of the old line as its intercept is higher |
| 79 | + if (new_m - last_m).abs() > 1e-9 { |
| 80 | + let intr = (new_b - last_b) / (last_m - new_m); |
| 81 | + if self.intersections.last().map(|&x| x < intr).unwrap_or(true) { |
| 82 | + self.intersections.push(intr); |
| 83 | + break; |
| 84 | + } |
| 85 | + } |
| 86 | + self.intersections.pop(); |
| 87 | + self.sorted_lines.pop(); |
| 88 | + } |
| 89 | + self.sorted_lines.push((new_m, new_b)); |
| 90 | + } |
| 91 | + } |
| 92 | + |
| 93 | + fn eval_helper(&self, x: f64) -> f64 { |
| 94 | + let idx = slice_lower_bound(&self.intersections, &x); |
| 95 | + std::iter::once(self.sorted_lines.get(idx)) |
| 96 | + .flatten() |
| 97 | + .chain(self.recent_lines.iter()) |
| 98 | + .map(|&(m, b)| m * x + b) |
| 99 | + .min_by(asserting_cmp) |
| 100 | + .unwrap_or(1e18) |
| 101 | + } |
| 102 | + |
| 103 | + /// Evaluates the function at x |
| 104 | + pub fn evaluate(&mut self, x: f64) -> f64 { |
| 105 | + if self.recent_lines.len() > self.merge_threshold { |
| 106 | + self.update_envelope(); |
| 107 | + } |
| 108 | + self.eval_helper(x) |
| 109 | + } |
| 110 | +} |
| 111 | + |
| 112 | +#[cfg(test)] |
| 113 | +mod test { |
| 114 | + use super::*; |
| 115 | + |
| 116 | + #[test] |
| 117 | + fn test_bounds() { |
| 118 | + let mut vals = vec![16, 45, 45, 45, 82]; |
| 119 | + |
| 120 | + assert_eq!(slice_upper_bound(&vals, &44), 1); |
| 121 | + assert_eq!(slice_lower_bound(&vals, &45), 1); |
| 122 | + assert_eq!(slice_upper_bound(&vals, &45), 4); |
| 123 | + assert_eq!(slice_lower_bound(&vals, &46), 4); |
| 124 | + |
| 125 | + vals.dedup(); |
| 126 | + for (i, q) in vals.iter().enumerate() { |
| 127 | + assert_eq!(slice_lower_bound(&vals, q), i); |
| 128 | + assert_eq!(slice_upper_bound(&vals, q), i + 1); |
| 129 | + } |
| 130 | + } |
| 131 | + |
| 132 | + #[test] |
| 133 | + fn test_coord_compress() { |
| 134 | + let mut coords = vec![16, 99, 45, 18]; |
| 135 | + let index = SparseIndex::new(coords.clone()); |
| 136 | + |
| 137 | + coords.sort_unstable(); |
| 138 | + for (i, q) in coords.into_iter().enumerate() { |
| 139 | + assert_eq!(index.compress(q - 1), Err(i)); |
| 140 | + assert_eq!(index.compress(q), Ok(i)); |
| 141 | + assert_eq!(index.compress(q + 1), Err(i + 1)); |
| 142 | + } |
| 143 | + } |
| 144 | + |
| 145 | + #[test] |
| 146 | + fn test_range_compress() { |
| 147 | + let queries = vec![(0, 10), (10, 19), (20, 29)]; |
| 148 | + let coords = queries.iter().flat_map(|&(i, j)| vec![i, j + 1]).collect(); |
| 149 | + let index = SparseIndex::new(coords); |
| 150 | + |
| 151 | + assert_eq!(index.coords, vec![0, 10, 11, 20, 30]); |
| 152 | + } |
| 153 | + |
| 154 | + #[test] |
| 155 | + fn test_convex_hull_trick() { |
| 156 | + let lines = [(0, 3), (1, 0), (-1, 8), (2, -1), (-1, 4)]; |
| 157 | + let xs = [0, 1, 2, 3, 4, 5]; |
| 158 | + // results[i] consists of the expected y-coordinates after processing |
| 159 | + // the first i+1 lines. |
| 160 | + let results = [ |
| 161 | + [3, 3, 3, 3, 3, 3], |
| 162 | + [0, 1, 2, 3, 3, 3], |
| 163 | + [0, 1, 2, 3, 3, 3], |
| 164 | + [-1, 1, 2, 3, 3, 3], |
| 165 | + [-1, 1, 2, 1, 0, -1], |
| 166 | + ]; |
| 167 | + for threshold in 0..=lines.len() { |
| 168 | + let mut func = PiecewiseLinearFn::with_merge_threshold(threshold); |
| 169 | + assert_eq!(func.evaluate(0.0), 1e18); |
| 170 | + for (&(slope, intercept), expected) in lines.iter().zip(results.iter()) { |
| 171 | + func.min_with(slope as f64, intercept as f64); |
| 172 | + let ys: Vec<i64> = xs.iter().map(|&x| func.evaluate(x as f64) as i64).collect(); |
| 173 | + assert_eq!(expected, &ys[..]); |
| 174 | + } |
| 175 | + } |
| 176 | + } |
| 177 | +} |
0 commit comments