subetha_cxc/
periodicity_sensor.rs1use std::collections::VecDeque;
18
19const BIN_MS: u64 = 500;
22const MAX_BINS: usize = 120;
25const LEO_PERIOD_MIN_S: f64 = 4.0;
27const LEO_PERIOD_MAX_S: f64 = 20.0;
28const CONF_FLOOR: f64 = 0.40;
30
31#[derive(Debug, Clone)]
33pub struct PeriodicitySensor {
34 bins: VecDeque<f64>,
35 cur_bin_start_us: u64,
36 cur_sum: f64,
37 cur_n: u32,
38 started: bool,
39}
40
41impl Default for PeriodicitySensor {
42 fn default() -> Self {
43 Self::new()
44 }
45}
46
47impl PeriodicitySensor {
48 pub fn new() -> Self {
49 Self {
50 bins: VecDeque::with_capacity(MAX_BINS),
51 cur_bin_start_us: 0,
52 cur_sum: 0.0,
53 cur_n: 0,
54 started: false,
55 }
56 }
57
58 pub fn observe(&mut self, owd_us: f64, t_us: u64) {
62 if !self.started {
63 self.cur_bin_start_us = t_us;
64 self.started = true;
65 }
66 while t_us >= self.cur_bin_start_us + BIN_MS * 1000 {
69 let v = if self.cur_n > 0 {
70 self.cur_sum / self.cur_n as f64
71 } else {
72 *self.bins.back().unwrap_or(&owd_us)
73 };
74 self.push_bin(v);
75 self.cur_bin_start_us += BIN_MS * 1000;
76 self.cur_sum = 0.0;
77 self.cur_n = 0;
78 }
79 self.cur_sum += owd_us;
80 self.cur_n += 1;
81 }
82
83 fn push_bin(&mut self, v: f64) {
84 if self.bins.len() == MAX_BINS {
85 self.bins.pop_front();
86 }
87 self.bins.push_back(v);
88 }
89
90 pub fn detected_period(&self) -> Option<(f64, f64)> {
94 let n = self.bins.len();
95 let min_bins = (2.0 * LEO_PERIOD_MIN_S * 1000.0 / BIN_MS as f64) as usize;
97 if n < min_bins {
98 return None;
99 }
100 let mean = self.bins.iter().sum::<f64>() / n as f64;
101 let var: f64 = self.bins.iter().map(|x| (x - mean).powi(2)).sum();
102 if var <= 0.0 {
103 return None;
104 }
105 let lag_min = (LEO_PERIOD_MIN_S * 1000.0 / BIN_MS as f64).round() as usize;
106 let lag_max = ((LEO_PERIOD_MAX_S * 1000.0 / BIN_MS as f64).round() as usize).min(n / 2);
107 let xs: Vec<f64> = self.bins.iter().copied().collect();
108 let (mut best_lag, mut best_r) = (0usize, 0.0f64);
109 for lag in lag_min..=lag_max {
110 let mut acc = 0.0;
111 for i in 0..(n - lag) {
112 acc += (xs[i] - mean) * (xs[i + lag] - mean);
113 }
114 let r = acc / var;
115 if r > best_r {
116 best_r = r;
117 best_lag = lag;
118 }
119 }
120 if best_lag == 0 || best_r < CONF_FLOOR {
121 return None;
122 }
123 Some((best_lag as f64 * BIN_MS as f64 / 1000.0, best_r))
124 }
125
126 pub fn secs_to_next_spike(&self) -> Option<f64> {
129 let (period_s, _) = self.detected_period()?;
130 let period_bins = (period_s * 1000.0 / BIN_MS as f64).round() as usize;
131 if period_bins == 0 {
132 return None;
133 }
134 let n = self.bins.len();
136 let start = n.saturating_sub(period_bins);
137 let xs: Vec<f64> = self.bins.iter().copied().collect();
138 let peak_off = (start..n)
139 .max_by(|&a, &b| xs[a].partial_cmp(&xs[b]).unwrap())
140 .unwrap_or(n - 1);
141 let since = (n - 1).saturating_sub(peak_off);
143 let to_next = period_bins.saturating_sub(since);
144 Some(to_next as f64 * BIN_MS as f64 / 1000.0)
145 }
146}
147
148#[cfg(test)]
149mod tests {
150 use super::*;
151
152 fn feed_periodic(s: &mut PeriodicitySensor, p_s: f64, cycles: usize) {
155 let dt_us = 100_000u64; let n = (cycles as f64 * p_s * 1e6 / dt_us as f64) as u64;
157 for i in 0..n {
158 let t = i * dt_us;
159 let phase = (t as f64 / 1e6) / p_s * std::f64::consts::TAU;
160 let owd = 40_000.0 + 60_000.0 * (phase.sin().max(0.0)).powi(4);
162 s.observe(owd, t);
163 }
164 }
165
166 #[test]
167 fn detects_a_clean_periodic_cadence() {
168 let mut s = PeriodicitySensor::new();
169 feed_periodic(&mut s, 15.0, 4);
170 let (period, conf) = s.detected_period().expect("a period is detected");
171 assert!((period - 15.0).abs() <= 1.0, "period ~ 15 s, got {period}");
172 assert!(conf > CONF_FLOOR, "confident, got {conf}");
173 }
174
175 #[test]
176 fn detects_a_shorter_cadence_too() {
177 let mut s = PeriodicitySensor::new();
178 feed_periodic(&mut s, 6.0, 6);
179 let (period, _) = s.detected_period().expect("a period is detected");
180 assert!((period - 6.0).abs() <= 1.0, "period ~ 6 s, got {period}");
181 }
182
183 #[test]
184 fn a_flat_trace_has_no_period() {
185 let mut s = PeriodicitySensor::new();
186 for i in 0..400 {
187 s.observe(40_000.0, i * 100_000); }
189 assert!(s.detected_period().is_none(), "a flat OWD has no cadence");
190 }
191
192 #[test]
193 fn predicts_the_next_spike_within_a_period() {
194 let mut s = PeriodicitySensor::new();
195 feed_periodic(&mut s, 10.0, 5);
196 let to_next = s.secs_to_next_spike().expect("a spike is predicted");
197 assert!(
198 (0.0..=10.0).contains(&to_next),
199 "next spike within one 10 s period, got {to_next}"
200 );
201 }
202}