From 5214e1edf76a3bb23c6738d30cd3247c215f2f9b Mon Sep 17 00:00:00 2001 From: daharoni Date: Wed, 8 Jul 2026 14:54:10 -0700 Subject: [PATCH] chore(solver): eliminate tuning magic-number drift hazards Numerical-hygiene pass on the solver crate for the pre-publication review. No behavior change: every default reproduces the previously hardcoded value. Biexp fast-grid dedup (biexp_fit.rs): - Hoist the fast-component grid bounds (0.25/2.0/0.5/8.0 x dt, 0.15 x tau_d) into shared named consts (TRF_*/TDF_*), used by both cold_grid_search and golden_section_refine so the two stages can no longer drift. - Fix a real drift: refinement floored tau_r_fast at 0.1 x dt while the grid floored at 0.25 x dt, letting refinement search a region the grid never saw. Refinement now clamps to the same shared bounds. - Resolve the two long-standing TODOs (sequential ceiling; dt-relative vs fixed-ms bounds) as documented design decisions with rationale. Seed/baseline/EMA config- ification: - Introduce SeedConfig (mad_multiplier, onset_fraction, max_walkback_s, min_peak_distance_s) with Default; find_seed_spikes reads from it instead of bare literals. FFI binding signatures (seed_trace, seed_kernel_estimate) unchanged. - Home the rolling-baseline percentile as baseline::DEFAULT_BASELINE_QUANTILE, used by both Solver::subtract_baseline and indeca::solve_trace (was the literal 0.2 duplicated across two modules). - Name the display-baseline EMA weight as BASELINE_EMA_WEIGHT. Tests: 131 Rust + 207 Python pass; clippy clean under jsbindings and pybindings. Co-Authored-By: Claude Fable 5 --- crates/solver/src/baseline.rs | 7 +++ crates/solver/src/biexp_fit.rs | 87 +++++++++++++++++++++------------- crates/solver/src/indeca.rs | 6 ++- crates/solver/src/lib.rs | 15 +++++- crates/solver/src/peak_seed.rs | 55 ++++++++++++++++----- 5 files changed, 122 insertions(+), 48 deletions(-) diff --git a/crates/solver/src/baseline.rs b/crates/solver/src/baseline.rs index 035344ed..b9d339c4 100644 --- a/crates/solver/src/baseline.rs +++ b/crates/solver/src/baseline.rs @@ -9,6 +9,13 @@ /// Uses a coordinate-compressed Fenwick tree (Binary Indexed Tree) for /// O(N log M) sliding-window k-th element queries, where M = distinct values. +/// Default rolling-baseline percentile (quantile of the causal window taken as +/// the floor estimate). Single source of truth for the callers that subtract a +/// rolling baseline (`Solver::subtract_baseline` and `indeca::solve_trace`), so +/// the two paths cannot drift; `subtract_rolling_baseline` still accepts any +/// quantile for callers that want to override it. +pub const DEFAULT_BASELINE_QUANTILE: f64 = 0.2; + /// Compute the rolling-baseline window size in samples. /// /// `5 * kernel_length` where `kernel_length = ceil(5 * tau_d * fs)`, diff --git a/crates/solver/src/biexp_fit.rs b/crates/solver/src/biexp_fit.rs index ade02152..04ee9090 100644 --- a/crates/solver/src/biexp_fit.rs +++ b/crates/solver/src/biexp_fit.rs @@ -71,16 +71,18 @@ /// • No external state: works identically on iteration 1 and iteration N, /// no dependence on warm-start or previous iteration history. /// -/// **TODO**: investigate whether this sequential ceiling is still necessary -/// now that the fast grid ranges have been tightened (tau_r_fast ≤ 2×dt, -/// tau_d_fast ≤ min(8×dt, tau_d×0.15), sub-sample floors). The original -/// slow/fast mixing problem may have been caused entirely by the fast grid -/// floors being too high (tau_r_fast ≥ dt, tau_d_fast ≥ 2×dt), which at -/// low sampling rates left almost no search range and forced wide fast -/// templates that overlapped heavily with the slow component. If removing -/// the ceiling with the current grid ranges does not reintroduce mixing, -/// the ceiling logic in eval_two_component can be simplified back to a -/// plain `bs >= bf` gate. +/// **Design decision — the ceiling is retained.** The fast grid ranges were +/// since tightened (tau_r_fast ≤ 2×dt, tau_d_fast ≤ min(8×dt, tau_d×0.15), +/// sub-sample floors), which addresses the *original* trigger of slow/fast +/// mixing: fast grid floors that were too high (tau_r_fast ≥ dt, +/// tau_d_fast ≥ 2×dt) left almost no search range at low sampling rates and +/// forced wide fast templates overlapping the slow component. Even so, the +/// sequential ceiling is kept as defence in depth rather than reverted to a +/// plain `bs >= bf` gate: it is O(1) (reuses inner products already computed), +/// self-calibrating (large artifact → generous ceiling), and independent of +/// the grid — so it still bounds redistribution for any (tau_r, tau_d) the +/// joint NNLS reaches, including warm-started points outside the cold grid. +/// The two guards are complementary, and the cost of keeping both is nil. /// Explicit outcome of a bi-exponential fit, so a degenerate / fallback fit is /// reported rather than silently inferred by the caller from `beta_fast == 0`. @@ -301,6 +303,21 @@ fn refine_candidate( /// search and the golden-section refinement so the two stages cannot drift. const TAU_D_HI: f64 = 5.0; +/// Fast-component grid bounds. Expressed as multipliers of `dt` (the sample +/// interval) except `TDF_REL_CAP`, which is relative to the slow `tau_d`. +/// +/// These are the single source of truth for the fast (noise-artifact) template's +/// search range: [`cold_grid_search`] uses them to lay down the initial grid, and +/// [`golden_section_refine`] uses the same values to clamp its narrowing search. +/// Keeping them here prevents the two stages from drifting apart (an earlier +/// version let refinement explore `tau_r_fast` down to `0.1 × dt` while the grid +/// floored at `0.25 × dt`, so refinement searched a region the grid never saw). +const TRF_LO_FACTOR: f64 = 0.25; // tau_r_fast floor = 0.25 × dt +const TRF_HI_FACTOR: f64 = 2.0; // tau_r_fast ceiling = 2 × dt +const TDF_LO_FACTOR: f64 = 0.5; // tau_d_fast floor = 0.5 × dt +const TDF_HI_FACTOR: f64 = 8.0; // tau_d_fast absolute ceiling = 8 × dt +const TDF_REL_CAP: f64 = 0.15; // tau_d_fast ≤ tau_d × 0.15 (relative ceiling) + /// Cold-start grid search. Returns (best_slow_only, best_two_component). fn cold_grid_search(h_free: &[f32], fs: f64, dt: f64, skip: usize) -> (BiexpResult, BiexpResult) { // Slow component grid ranges (in seconds). @@ -337,21 +354,23 @@ fn cold_grid_search(h_free: &[f32], fs: f64, dt: f64, skip: usize) -> (BiexpResu // - relative cap (tau_d × 0.15) prevents the fast template's decay // from extending into the slow component's domain // - // TODO: revisit whether these bounds should depend on dt at all, or - // whether fixed absolute limits in milliseconds would be more physically - // meaningful. The noise artifact's shape comes from the noise - // autocorrelation structure, which has a characteristic timescale - // independent of the sampling rate. dt-relative bounds were chosen - // pragmatically (the artifact occupies 1-3 discrete bins), but the - // underlying physics may justify fixed ms bounds (e.g. tau_r_fast ≤ 20ms, - // tau_d_fast ≤ 60ms) that would behave more consistently across sampling - // rates. + // Design decision — bounds are dt-relative (not fixed ms). One could argue + // for fixed absolute limits (e.g. tau_r_fast ≤ 20ms) on the grounds that the + // noise autocorrelation has a characteristic timescale independent of fs. + // We keep dt-relative bounds because the artifact this template absorbs is a + // *discretization* artifact: false-positive spikes imprint a feature onto + // h_free that occupies a small, fixed number of discrete bins (1-3) around + // the kernel's rising edge, so its extent in seconds scales with dt by + // construction. Fixed-ms bounds would, at high fs, span many bins (letting + // the fast template poach real slow-kernel structure) and, at low fs, span + // less than one bin (making it unrepresentable). The shared factors above + // (TRF_*/TDF_*) are the single knob if this ever needs revisiting. let trf_grid_n = 5; let tdf_grid_n = 8; - let trf_lo = 0.25 * dt; - let trf_hi = 2.0 * dt; - let tdf_lo = 0.5 * dt; - let tdf_abs_hi = 8.0 * dt; + let trf_lo = TRF_LO_FACTOR * dt; + let trf_hi = TRF_HI_FACTOR * dt; + let tdf_lo = TDF_LO_FACTOR * dt; + let tdf_abs_hi = TDF_HI_FACTOR * dt; let mut best_slow = BiexpResult::sentinel(); let mut best_two = BiexpResult::sentinel(); @@ -388,7 +407,7 @@ fn cold_grid_search(h_free: &[f32], fs: f64, dt: f64, skip: usize) -> (BiexpResu // Inner grid: scan independent (tau_r_fast, tau_d_fast) // Upper bound for tau_d_fast is the tighter of the absolute cap // (8×dt) and a relative cap (tau_d × 0.15) to prevent degeneracy. - let tdf_hi = tdf_abs_hi.min(tau_d * 0.15); + let tdf_hi = tdf_abs_hi.min(tau_d * TDF_REL_CAP); if tdf_hi <= tdf_lo { continue; // tau_d too small for a distinct fast component } @@ -634,9 +653,10 @@ fn golden_section_refine( } } 2 => { - // Refine tau_r_fast - let lo = (tau_r_fast * 0.5).max(0.1 * dt); - let hi = (tau_r_fast * 2.0).min((2.0 * dt).min(tau_d_fast * 0.99)); + // Refine tau_r_fast — clamped to the same grid bounds so + // refinement cannot wander outside the range the grid searched. + let lo = (tau_r_fast * 0.5).max(TRF_LO_FACTOR * dt); + let hi = (tau_r_fast * 2.0).min((TRF_HI_FACTOR * dt).min(tau_d_fast * 0.99)); if lo < hi { tau_r_fast = golden_bracket(lo, hi, |x| { eval_two_component(h_free, tau_r, tau_d, x, tau_d_fast, dt, skip).2 @@ -644,9 +664,12 @@ fn golden_section_refine( } } _ => { - // Refine tau_d_fast - let lo = (tau_d_fast * 0.5).max(tau_r_fast * 1.01); - let hi = (tau_d_fast * 2.0).min((8.0 * dt).min(tau_d * 0.15)); + // Refine tau_d_fast — floor clamped to the grid floor and the + // ordering constraint; ceiling to the shared absolute/relative caps. + let lo = (tau_d_fast * 0.5) + .max(tau_r_fast * 1.01) + .max(TDF_LO_FACTOR * dt); + let hi = (tau_d_fast * 2.0).min((TDF_HI_FACTOR * dt).min(tau_d * TDF_REL_CAP)); if lo < hi { tau_d_fast = golden_bracket(lo, hi, |x| { eval_two_component(h_free, tau_r, tau_d, tau_r_fast, x, dt, skip).2 @@ -978,7 +1001,7 @@ mod tests { tau_d, fs ); - let tdf_cap = (8.0 * dt).min(tau_d * 0.15); + let tdf_cap = (TDF_HI_FACTOR * dt).min(tau_d * TDF_REL_CAP); assert!( result.tau_decay_fast <= tdf_cap * 1.05, "tau_d_fast ({:.6}) should be ≤ cap ({:.6}) for (tau_r={}, tau_d={}, fs={})", @@ -1085,7 +1108,7 @@ mod tests { // If there's a fast component, it should be confined if result.tau_decay_fast > 0.0 { - let tdf_cap = (8.0 * dt).min(tau_d_true * 0.15); + let tdf_cap = (TDF_HI_FACTOR * dt).min(tau_d_true * TDF_REL_CAP); assert!( result.tau_decay_fast <= tdf_cap * 1.05, "tau_d_fast ({:.6}s = {:.1} bins) should be ≤ cap ({:.6}s) for large tau_d", diff --git a/crates/solver/src/indeca.rs b/crates/solver/src/indeca.rs index 503c0620..2028c291 100644 --- a/crates/solver/src/indeca.rs +++ b/crates/solver/src/indeca.rs @@ -233,7 +233,11 @@ pub fn solve_trace( // Rolling-percentile baseline subtraction: brings the floor to ~0. let bl_window = crate::baseline::baseline_window(tau_d, fs_up); - crate::baseline::subtract_rolling_baseline(&mut working_trace, bl_window, 0.2); + crate::baseline::subtract_rolling_baseline( + &mut working_trace, + bl_window, + crate::baseline::DEFAULT_BASELINE_QUANTILE, + ); // ── Step 2: Boundary padding + initial alpha estimate ─────────────── // Compute boundary padding: edge effects from AR2 convolution make the first diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 8aac9296..2618f095 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -43,6 +43,12 @@ pub(crate) fn first_nonfinite(data: &[f32]) -> Option { data.iter().position(|v| !v.is_finite()) } +/// EMA weight for the displayed scalar baseline. The per-iteration raw baseline +/// is smoothed as `ema = W·raw + (1-W)·ema` purely to keep the displayed value +/// steady across iterations; it does not enter the solve. Named here so the +/// smoothing strength is a single, documented knob rather than a bare literal. +const BASELINE_EMA_WEIGHT: f64 = 0.3; + /// Convolution mode for forward/adjoint operations in FISTA. #[derive(Clone, Copy, PartialEq, Eq)] #[cfg_attr(feature = "jsbindings", wasm_bindgen)] @@ -428,7 +434,8 @@ impl Solver { self.baseline_ema = raw_baseline; self.baseline_ema_init = true; } else { - self.baseline_ema = 0.3 * raw_baseline + 0.7 * self.baseline_ema; + self.baseline_ema = BASELINE_EMA_WEIGHT * raw_baseline + + (1.0 - BASELINE_EMA_WEIGHT) * self.baseline_ema; } } @@ -477,7 +484,11 @@ impl Solver { return; } let window = baseline::baseline_window(self.tau_decay, self.fs); - baseline::subtract_rolling_baseline(&mut self.trace[..n], window, 0.2); + baseline::subtract_rolling_baseline( + &mut self.trace[..n], + window, + baseline::DEFAULT_BASELINE_QUANTILE, + ); self.filtered = true; } diff --git a/crates/solver/src/peak_seed.rs b/crates/solver/src/peak_seed.rs index f496a3bc..6feef9ff 100644 --- a/crates/solver/src/peak_seed.rs +++ b/crates/solver/src/peak_seed.rs @@ -28,7 +28,7 @@ pub struct SeedTraceResult { pub fn seed_trace(trace: &[f32], fs: f64) -> SeedTraceResult { let n = trace.len(); let (bl, _) = median_and_mad(trace); - let onsets = find_seed_spikes(trace, fs, 5.0); + let onsets = find_seed_spikes(trace, fs, &SeedConfig::default()); let mut s_counts = vec![0.0_f32; n]; for &idx in &onsets { @@ -108,16 +108,45 @@ pub(crate) fn median_and_mad(data: &[f32]) -> (f32, f32) { (med, mad_val) } +/// Tunable parameters for peak-seeded spike detection. +/// +/// `Default` reproduces the values that were previously hardcoded in +/// `find_seed_spikes`, so behaviour is unchanged unless a caller overrides a +/// field. Every field shifts which onsets get seeded, so per the no-magic-number +/// rule none of them appears as a bare literal in the detection logic. +#[derive(Clone, Copy, Debug)] +pub struct SeedConfig { + /// Peak threshold = median + `mad_multiplier` × MAD. + pub mad_multiplier: f32, + /// Onset threshold = baseline + `onset_fraction` × (peak − baseline). + pub onset_fraction: f32, + /// Maximum walk-back from a peak to its onset, in seconds. + pub max_walkback_s: f64, + /// Minimum spacing between accepted peaks, in seconds. + pub min_peak_distance_s: f64, +} + +impl Default for SeedConfig { + fn default() -> Self { + SeedConfig { + mad_multiplier: 4.0, + onset_fraction: 0.10, + max_walkback_s: 1.0, + min_peak_distance_s: 5.0, + } + } +} + /// Find seed spike onset locations from a single trace. /// /// 1. Compute baseline = median(trace) -/// 2. Compute threshold = baseline + 4 * MAD -/// 3. Find local maxima above threshold, at least `min_peak_distance_s` seconds apart -/// 4. Walk back from each peak to onset (where trace drops to baseline + 10% of peak height) -/// with a max walk-back of 1 second +/// 2. Compute threshold = baseline + `cfg.mad_multiplier` * MAD +/// 3. Find local maxima above threshold, at least `cfg.min_peak_distance_s` seconds apart +/// 4. Walk back from each peak to onset (where trace drops to +/// baseline + `cfg.onset_fraction` * peak height), bounded by `cfg.max_walkback_s` /// /// Returns indices into the trace where spikes should be placed. -pub fn find_seed_spikes(trace: &[f32], fs: f64, min_peak_distance_s: f64) -> Vec { +pub fn find_seed_spikes(trace: &[f32], fs: f64, cfg: &SeedConfig) -> Vec { let n = trace.len(); if n < 3 { return Vec::new(); @@ -128,10 +157,10 @@ pub fn find_seed_spikes(trace: &[f32], fs: f64, min_peak_distance_s: f64) -> Vec if mad_val < 1e-10 { return Vec::new(); } - let threshold = baseline + 4.0 * mad_val; + let threshold = baseline + cfg.mad_multiplier * mad_val; - let min_peak_dist = (min_peak_distance_s * fs).round() as usize; - let max_walkback = (1.0 * fs).round() as usize; + let min_peak_dist = (cfg.min_peak_distance_s * fs).round() as usize; + let max_walkback = (cfg.max_walkback_s * fs).round() as usize; // Find local maxima above threshold let mut peaks: Vec<(usize, f32)> = Vec::new(); @@ -160,7 +189,7 @@ pub fn find_seed_spikes(trace: &[f32], fs: f64, min_peak_distance_s: f64) -> Vec let mut onsets: Vec = Vec::with_capacity(selected_peaks.len()); for &peak_idx in &selected_peaks { let peak_val = trace[peak_idx]; - let onset_threshold = baseline + 0.10 * (peak_val - baseline); + let onset_threshold = baseline + cfg.onset_fraction * (peak_val - baseline); let earliest = peak_idx.saturating_sub(max_walkback); @@ -202,7 +231,7 @@ pub fn seed_kernel_estimate( ) -> SeedKernelResult { let total_len: usize = trace_lengths.iter().sum(); - let min_peak_distance_s = 5.0; + let seed_cfg = SeedConfig::default(); // Kernel length: ~1.5 seconds at the given sampling rate let kernel_length = (1.5 * fs).ceil() as usize; @@ -236,7 +265,7 @@ pub fn seed_kernel_estimate( let trace = &traces_flat[offset..offset + len]; let (bl, _) = median_and_mad(trace); - let onsets = find_seed_spikes(trace, fs, min_peak_distance_s); + let onsets = find_seed_spikes(trace, fs, &seed_cfg); for &onset in &onsets { spike_trains[offset + onset] = 1.0; total_seed_spikes += 1; @@ -344,7 +373,7 @@ mod tests { } } - let onsets = find_seed_spikes(&trace, fs, 5.0); + let onsets = find_seed_spikes(&trace, fs, &SeedConfig::default()); assert!( onsets.len() >= 2 && onsets.len() <= 5,