1use std::collections::VecDeque;
30
31pub const ENVELOPE_SLOT_NS: u64 = 50_000_000;
33pub const ENVELOPE_BIN_NS: u64 = 2_000_000_000;
35pub const ENVELOPE_SIGMA_BINS: u32 = 3;
37pub const ENVELOPE_LOCK_SPAN_NS: u64 = 40_000_000_000;
39pub const ENVELOPE_SLOPE_SPAN_NS: u64 = 2_000_000_000;
43
44#[derive(Debug, Clone)]
47pub struct EnvelopeFit {
48 window_ns: u64,
49 pts: VecDeque<(u64, i128)>,
51 hull: Vec<(u64, i128)>,
53 model: Option<Line>,
54}
55
56#[derive(Debug, Clone, Copy)]
58struct Line {
59 x0: u64,
60 y0: i128,
61 s: f64,
62 sigma_ns: i64,
63}
64
65impl Line {
66 fn rel(&self, host_ns: u64) -> f64 {
67 self.s * (host_ns as i128 - self.x0 as i128) as f64
68 }
69}
70
71impl EnvelopeFit {
72 pub fn new(window_ns: u64) -> Self {
74 EnvelopeFit {
75 window_ns,
76 pts: VecDeque::new(),
77 hull: Vec::new(),
78 model: None,
79 }
80 }
81
82 pub fn sixty_seconds() -> Self {
84 Self::new(60_000_000_000)
85 }
86
87 pub fn push(&mut self, host_ns: u64, source_ns: i64) {
89 let y = source_ns as i128 - host_ns as i128;
90 let mut at = self.pts.len();
92 while at > 0 && self.pts[at - 1].0 > host_ns {
93 at -= 1;
94 }
95 let slot = host_ns / ENVELOPE_SLOT_NS;
97 let same = [at.checked_sub(1), Some(at)]
98 .into_iter()
99 .flatten()
100 .find(|&i| {
101 self.pts
102 .get(i)
103 .is_some_and(|p| p.0 / ENVELOPE_SLOT_NS == slot)
104 });
105 if let Some(i) = same {
106 if self.pts[i].1 >= y {
107 return;
108 }
109 self.pts.remove(i);
110 if i < at {
111 at -= 1;
112 }
113 }
114 self.pts.insert(at, (host_ns, y));
115 let newest = self.pts.back().map_or(host_ns, |p| p.0);
116 while let Some(&(h, _)) = self.pts.front() {
117 if newest.saturating_sub(h) > self.window_ns {
118 self.pts.pop_front();
119 } else {
120 break;
121 }
122 }
123 self.refit();
124 }
125
126 pub fn reset(&mut self) {
128 self.pts.clear();
129 self.model = None;
130 }
131
132 pub fn len(&self) -> usize {
134 self.pts.len()
135 }
136
137 pub fn is_empty(&self) -> bool {
139 self.pts.is_empty()
140 }
141
142 pub fn span_ns(&self) -> u64 {
144 match (self.pts.front(), self.pts.back()) {
145 (Some(a), Some(b)) => b.0 - a.0,
146 _ => 0,
147 }
148 }
149
150 fn refit(&mut self) {
151 let (Some(&(oldest, _)), Some(&(newest, _))) = (self.pts.front(), self.pts.back()) else {
152 self.model = None;
153 return;
154 };
155 if newest - oldest < ENVELOPE_SLOPE_SPAN_NS {
156 let (x0, y0) = self
157 .pts
158 .iter()
159 .copied()
160 .max_by_key(|p| p.1)
161 .unwrap_or((oldest, 0));
162 let mut line = Line {
163 x0,
164 y0,
165 s: 0.0,
166 sigma_ns: -1,
167 };
168 line.sigma_ns = self.sigma(&line);
169 self.model = Some(line);
170 return;
171 }
172 self.hull.clear();
174 for &(x, y) in &self.pts {
175 if let Some(&(lx, ly)) = self.hull.last() {
176 if lx == x {
177 if y <= ly {
178 continue;
179 }
180 self.hull.pop();
181 }
182 }
183 while let [.., (ax, ay), (bx, by)] = self.hull[..] {
184 let cross = (bx - ax) as i128 * (y - ay) - (by - ay) * (x - ax) as i128;
186 if cross >= 0 {
187 self.hull.pop();
188 } else {
189 break;
190 }
191 }
192 self.hull.push((x, y));
193 }
194 let Some(&first) = self.hull.first() else {
195 self.model = None;
196 return;
197 };
198 let n = self.pts.len() as u128;
199 let sum: u128 = self.pts.iter().map(|&(x, _)| (x - first.0) as u128).sum();
200 let mean_x = first.0 + (sum / n) as u64;
201 let (a, b) = self
202 .hull
203 .windows(2)
204 .map(|w| (w[0], w[1]))
205 .find(|&(_, b)| b.0 >= mean_x)
206 .unwrap_or((first, first));
207 let s = if b.0 > a.0 {
208 (b.1 - a.1) as f64 / (b.0 - a.0) as f64
209 } else {
210 0.0
211 };
212 let mut line = Line {
213 x0: a.0,
214 y0: a.1,
215 s,
216 sigma_ns: -1,
217 };
218 line.sigma_ns = self.sigma(&line);
219 self.model = Some(line);
220 }
221
222 fn sigma(&self, line: &Line) -> i64 {
226 let (Some(&(oldest, _)), Some(&(newest, _))) = (self.pts.front(), self.pts.back()) else {
227 return -1;
228 };
229 let (first_bin, open_bin) = (oldest / ENVELOPE_BIN_NS, newest / ENVELOPE_BIN_NS);
230 let mut bin = u64::MAX;
231 let mut min = f64::INFINITY;
232 let (mut sum_sq, mut bins) = (0.0f64, 0u32);
233 for &(x, y) in &self.pts {
234 let b = x / ENVELOPE_BIN_NS;
235 if b == first_bin {
236 continue;
237 }
238 if b == open_bin {
239 break;
240 }
241 if b != bin {
242 if min.is_finite() {
243 sum_sq += min * min;
244 bins += 1;
245 }
246 bin = b;
247 min = f64::INFINITY;
248 }
249 let late = ((line.y0 - y) as f64 + line.rel(x)).max(0.0);
250 min = min.min(late);
251 }
252 if min.is_finite() {
253 sum_sq += min * min;
254 bins += 1;
255 }
256 if bins < ENVELOPE_SIGMA_BINS {
257 return -1;
258 }
259 (sum_sq / f64::from(bins)).sqrt().round() as i64
260 }
261
262 pub fn offset_i64(&self, host_ns: u64) -> Option<i64> {
264 let m = self.model?;
265 Some((m.y0 + m.rel(host_ns).round() as i128) as i64)
266 }
267
268 pub fn to_source(&self, host_ns: u64) -> Option<i64> {
270 self.offset_i64(host_ns)
271 .map(|o| (host_ns as i128 + o as i128) as i64)
272 }
273
274 pub fn skew_ppm(&self) -> f64 {
276 self.model.map_or(0.0, |m| m.s * 1e6)
277 }
278
279 pub fn sigma_ns(&self) -> i64 {
282 self.model.map_or(-1, |m| m.sigma_ns)
283 }
284}
285
286#[cfg(test)]
287mod tests {
288 use super::*;
289
290 const MS: u64 = 1_000_000;
291
292 #[test]
293 fn a_noiseless_line_is_recovered_exactly() {
294 let mut f = EnvelopeFit::sixty_seconds();
295 for i in 0..100u64 {
296 let h = 1_000 * MS + i * 100 * MS;
297 f.push(h, (h as f64 * (1.0 + 40e-6)) as i64 - 37 * MS as i64);
299 }
300 let h = 1_000 * MS + 99 * 100 * MS;
301 let want = (h as f64 * (1.0 + 40e-6)) as i64 - 37 * MS as i64;
302 assert!(f.to_source(h).unwrap().abs_diff(want) <= 2);
303 assert!((f.skew_ppm() - 40.0).abs() < 0.01);
304 assert!(f.sigma_ns() <= 2, "sigma {}", f.sigma_ns());
305 }
306
307 #[test]
308 fn late_pairs_do_not_move_the_line() {
309 let mut f = EnvelopeFit::sixty_seconds();
310 for i in 0..50u64 {
311 let h = i * 100 * MS;
312 f.push(h, h as i64);
313 }
314 let before = f.offset_i64(49 * 100 * MS);
315 f.push(4_950 * MS, (4_750 * MS) as i64);
318 f.push(4_920 * MS, (4_720 * MS) as i64);
319 assert_eq!(f.offset_i64(49 * 100 * MS), before);
320 assert_eq!(f.len(), 51);
321 }
322
323 #[test]
324 fn each_slot_keeps_its_least_late_pair() {
325 let mut f = EnvelopeFit::sixty_seconds();
326 for k in 0..70 * 30u64 {
328 for cam in 0..4u64 {
329 let h = k * 33_333_333 + cam * 7 * MS;
330 f.push(h, h as i64 - (cam * 3 * MS) as i64);
332 }
333 }
334 assert!(f.len() <= (60_000 / 50) as usize + 1, "{} pairs", f.len());
335 assert!(f.span_ns() >= 59_900 * MS, "span {}", f.span_ns());
336 assert_eq!(f.offset_i64(69 * 1_000 * MS), Some(0));
338 let h = 70_000 * MS + 10 * MS;
340 f.push(h, h as i64 - 20 * MS as i64);
341 f.push(h + 5 * MS, (h + 5 * MS) as i64);
342 assert_eq!(f.pts.back(), Some(&(h + 5 * MS, 0)));
343 }
344
345 #[test]
346 fn a_young_fit_is_flat_through_the_least_late_pair() {
347 let mut f = EnvelopeFit::sixty_seconds();
348 f.push(1_000 * MS, (800 * MS) as i64);
350 f.push(1_050 * MS, (1_050 * MS) as i64);
351 assert_eq!(f.skew_ppm(), 0.0);
352 assert_eq!(f.offset_i64(3_000 * MS), Some(0));
353 }
354
355 #[test]
356 fn sigma_waits_for_three_whole_bins() {
357 let mut f = EnvelopeFit::sixty_seconds();
358 let mut h = 0;
359 while h < 4 * ENVELOPE_BIN_NS {
360 f.push(h, h as i64);
361 h += 100 * MS;
362 }
363 assert_eq!(
364 f.sigma_ns(),
365 -1,
366 "the oldest bin is left out, the fourth is open"
367 );
368 f.push(h, h as i64);
369 assert_eq!(f.sigma_ns(), 0);
370 }
371}