chore(solver): eliminate tuning magic-number drift hazards

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)`,

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,19 +653,23 @@ 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
                     });
                 }
             }
             _ => {
-                // 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",

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

crates/solver/src/lib.rs

@@ -43,6 +43,12 @@ pub(crate) fn first_nonfinite(data: &[f32]) -> Option<usize> {
     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;
     }
 

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<usize> {
+pub fn find_seed_spikes(trace: &[f32], fs: f64, cfg: &SeedConfig) -> Vec<usize> {
     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<usize> = 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,