stats_claw/streaming/p2.rs
1//! The P² streaming-quantile estimator (Jain & Chlamtac, 1985).
2//
3// `indexing_slicing` is allowed in this module only: every index is a
4// compile-time constant or a loop bound in `0..5` / `1..4` indexing the
5// fixed-length `[_; 5]` marker arrays, so the access is provably in bounds.
6// Rewriting each as `.get(..).unwrap_or(..)` would obscure the dense P²
7// arithmetic and add unreachable fallbacks, so a scoped allow is the clearer
8// choice here (no `unwrap`/`as`/`unsafe` is introduced).
9#![allow(clippy::indexing_slicing)]
10
11use super::position_as_f64;
12
13/// Online estimator of a single p-quantile using the P² algorithm.
14///
15/// Tracks five markers — the running minimum, the p/2, p, and (1+p)/2 quantiles,
16/// and the running maximum — and adjusts their heights with a piecewise-parabolic
17/// (hence "P²") interpolation as each value arrives. It estimates the quantile in
18/// O(1) time and O(1) memory per update without storing or sorting the stream,
19/// trading a small approximation error for unbounded scalability.
20///
21/// # Invariants
22///
23/// All state lives in fixed-length arrays (five markers), so
24/// `size_of::<P2Quantile>()` is constant regardless of how many values have been
25/// consumed — the bounded-memory guarantee.
26///
27/// # Examples
28///
29/// ```
30/// use stats_claw::streaming::P2Quantile;
31///
32/// let mut q = P2Quantile::new(0.5);
33/// for x in 1..=99 {
34/// q.update(f64::from(x));
35/// }
36/// // Median of 1..=99 is 50; P² lands within a couple of integers on this ramp.
37/// assert!((q.value() - 50.0).abs() < 2.0);
38/// ```
39#[derive(Debug, Clone)]
40pub struct P2Quantile {
41 /// Number of values consumed so far.
42 count: u64,
43 /// Marker heights (the running quantile estimates), ascending.
44 heights: [f64; 5],
45 /// Integer marker positions (1-based, as in the original paper).
46 positions: [i64; 5],
47 /// Desired (real-valued) marker positions.
48 desired: [f64; 5],
49 /// Per-marker increments to the desired positions on each new value.
50 increments: [f64; 5],
51}
52
53impl P2Quantile {
54 /// Creates an estimator for the `p`-quantile.
55 ///
56 /// # Arguments
57 ///
58 /// * `p` — the quantile probability, expected in `[0, 1]` (e.g. `0.5` for the
59 /// median, `0.99` for the 99th percentile). Values outside the range are not
60 /// rejected here — the estimate is simply meaningless, mirroring the
61 /// panic-free convention of the rest of the crate.
62 ///
63 /// # Returns
64 ///
65 /// A fresh estimator with no values consumed; [`Self::value`] is `0.0` until
66 /// the first [`Self::update`].
67 #[must_use]
68 pub fn new(p: f64) -> Self {
69 Self {
70 count: 0,
71 heights: [0.0; 5],
72 positions: [1, 2, 3, 4, 5],
73 desired: [
74 1.0,
75 2.0_f64.mul_add(p, 1.0),
76 4.0_f64.mul_add(p, 1.0),
77 2.0_f64.mul_add(p, 3.0),
78 5.0,
79 ],
80 increments: [0.0, p / 2.0, p, f64::midpoint(1.0, p), 1.0],
81 }
82 }
83
84 /// Folds one observation into the estimator.
85 ///
86 /// # Arguments
87 ///
88 /// * `x` — the next value of the stream.
89 pub fn update(&mut self, x: f64) {
90 if self.count < 5 {
91 self.bootstrap(x);
92 return;
93 }
94 let k = self.locate_cell(x);
95 for marker in (k + 1)..5 {
96 self.positions[marker] += 1;
97 }
98 for marker in 0..5 {
99 self.desired[marker] += self.increments[marker];
100 }
101 self.adjust_interior();
102 self.count += 1;
103 }
104
105 /// Returns the current quantile estimate: the middle marker height once five
106 /// values are in, the largest seen value during warm-up, or `0.0` if empty.
107 #[must_use]
108 pub fn value(&self) -> f64 {
109 if self.count == 0 {
110 return 0.0;
111 }
112 if self.count < 5 {
113 let last = usize::try_from(self.count).unwrap_or(0).saturating_sub(1);
114 return self.heights.get(last).copied().unwrap_or(0.0);
115 }
116 self.heights[2]
117 }
118
119 /// Returns the number of values consumed so far.
120 #[must_use]
121 pub const fn count(&self) -> u64 {
122 self.count
123 }
124
125 /// Fills the initial five markers, keeping them sorted, before P² proper runs.
126 fn bootstrap(&mut self, x: f64) {
127 let slot = usize::try_from(self.count).unwrap_or(0);
128 if let Some(h) = self.heights.get_mut(slot) {
129 *h = x;
130 }
131 self.count += 1;
132 if self.count == 5 {
133 self.heights
134 .sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
135 }
136 }
137
138 /// Finds the marker cell `x` falls into and clamps the running extrema.
139 ///
140 /// # Returns
141 ///
142 /// The index `k` in `0..4` such that `heights[k] <= x < heights[k+1]`, after
143 /// extending the min/max markers to cover `x`.
144 fn locate_cell(&mut self, x: f64) -> usize {
145 if x < self.heights[0] {
146 self.heights[0] = x;
147 return 0;
148 }
149 for k in 0..4 {
150 if x < self.heights[k + 1] {
151 return k;
152 }
153 }
154 self.heights[4] = x;
155 3
156 }
157
158 /// Adjusts the three interior markers toward their desired positions using the
159 /// parabolic formula, falling back to linear when the parabola would break the
160 /// ascending order.
161 fn adjust_interior(&mut self) {
162 for i in 1..4 {
163 let d = self.desired[i] - position_as_f64(self.positions[i]);
164 let gap_hi = self.positions[i + 1] - self.positions[i];
165 let gap_lo = self.positions[i] - self.positions[i - 1];
166 if (d >= 1.0 && gap_hi > 1) || (d <= -1.0 && gap_lo > 1) {
167 let dir = d.signum();
168 let candidate = self.parabolic(i, dir);
169 self.heights[i] =
170 if self.heights[i - 1] < candidate && candidate < self.heights[i + 1] {
171 candidate
172 } else {
173 self.linear(i, dir)
174 };
175 self.positions[i] += i64::from(dir >= 0.0) * 2 - 1;
176 }
177 }
178 }
179
180 /// Piecewise-parabolic prediction (P²) of marker `i`'s new height.
181 fn parabolic(&self, i: usize, dir: f64) -> f64 {
182 let pos_lo = position_as_f64(self.positions[i - 1]);
183 let pos_mid = position_as_f64(self.positions[i]);
184 let pos_hi = position_as_f64(self.positions[i + 1]);
185 let upper =
186 (pos_mid - pos_lo + dir) * (self.heights[i + 1] - self.heights[i]) / (pos_hi - pos_mid);
187 let lower =
188 (pos_hi - pos_mid - dir) * (self.heights[i] - self.heights[i - 1]) / (pos_mid - pos_lo);
189 (dir / (pos_hi - pos_lo)).mul_add(upper + lower, self.heights[i])
190 }
191
192 /// Linear fallback prediction of marker `i`'s new height.
193 fn linear(&self, i: usize, dir: f64) -> f64 {
194 let neighbor = if dir >= 0.0 { i + 1 } else { i - 1 };
195 let pos_i = position_as_f64(self.positions[i]);
196 let pos_neighbor = position_as_f64(self.positions[neighbor]);
197 self.heights[i] + dir * (self.heights[neighbor] - self.heights[i]) / (pos_neighbor - pos_i)
198 }
199}