From 77546cda31743c0d642297b8bb0404676062cbd8 Mon Sep 17 00:00:00 2001 From: daharoni Date: Wed, 8 Jul 2026 21:23:31 -0700 Subject: [PATCH 1/4] feat(cadecon): noise-constrained sparsity for low-SNR deconvolution Add an optional, knob-free sparsity to the InDeCa spike inference that suppresses noise being fit as spurious activity. Instead of picking the binarization threshold that maximizes fit (which overfits into the noise), pick the sparsest spike support whose reconstruction residual still reaches the data-derived noise floor ("don't spend spikes explaining what's within the noise"). The effect concentrates at low SNR and is neutral at high SNR. Core (crates/solver): - threshold.rs: Selection::{MaxPve,NoiseFloor}; threshold_search_opts + select_noise_floor_threshold pick the highest (sparsest) threshold whose residual, evaluated at original-rate grid positions, stays within sigma^2*T. - indeca.rs: SolveOptions{noise_constrained, collapse_runs}; solve_trace_opts. Noise sigma is LP-robust: raw high-frequency-PSD estimate (high_band_sigma) rescaled by the empirically measured noise gain of the upsample+HP/LP chain (estimate_grid_noise_sigma), so it tracks the noise that survives the filters. - upsample.rs: collapse_runs (event-preserving reduction; not wired to UI). - js_indeca.rs / py_api.rs: expose the flag through the WASM + PyO3 bindings. App (apps/cadecon): - Thread noise_constrained through worker -> pool -> iteration-manager, read from a new algorithm-store signal. Default ON. - Teaching/impact overlay (opt-in, default OFF): when enabled, each trace is solved with BOTH sparsity settings per iteration and the opposite-setting counts are stored; the Trace Inspector overlays them and shows both spike counts. Precomputed in the run (no main-thread solving) so browsing stays smooth. Lives under a new "Display Options" panel (it is a visual-inspection aid, not an algorithm knob) and doubles run time, documented in-UI and in the "Reading Convergence & Results" tutorial. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01BDPugMWynodvqCBPH4p8ab --- apps/cadecon/src/App.tsx | 10 + .../components/controls/AlgorithmSettings.tsx | 10 + .../components/controls/DisplaySettings.tsx | 26 ++ .../src/components/traces/TraceInspector.tsx | 68 +++- apps/cadecon/src/lib/algorithm-store.ts | 12 + apps/cadecon/src/lib/cadecon-pool.ts | 4 + apps/cadecon/src/lib/iteration-manager.ts | 26 ++ apps/cadecon/src/lib/iteration-store.ts | 2 + .../lib/tutorial/content/04-interpreting.ts | 12 +- apps/cadecon/src/lib/viz-store.ts | 7 + apps/cadecon/src/styles/layout.css | 9 + apps/cadecon/src/workers/cadecon-types.ts | 9 + apps/cadecon/src/workers/cadecon-worker.ts | 25 ++ crates/solver/src/indeca.rs | 191 +++++++++++- crates/solver/src/js_indeca.rs | 7 +- crates/solver/src/py_api.rs | 11 +- crates/solver/src/threshold.rs | 292 +++++++++++++----- crates/solver/src/upsample.rs | 28 ++ 18 files changed, 667 insertions(+), 82 deletions(-) create mode 100644 apps/cadecon/src/components/controls/DisplaySettings.tsx diff --git a/apps/cadecon/src/App.tsx b/apps/cadecon/src/App.tsx index 9da7e764..542af43e 100644 --- a/apps/cadecon/src/App.tsx +++ b/apps/cadecon/src/App.tsx @@ -20,6 +20,7 @@ import { ImportOverlay } from './components/layout/ImportOverlay.tsx'; import { RasterOverview } from './components/raster/RasterOverview.tsx'; import { SubsetConfig } from './components/controls/SubsetConfig.tsx'; import { AlgorithmSettings } from './components/controls/AlgorithmSettings.tsx'; +import { DisplaySettings } from './components/controls/DisplaySettings.tsx'; import { RunControls } from './components/controls/RunControls.tsx'; import { ProgressBar } from './components/controls/ProgressBar.tsx'; import { ConvergencePanel } from './components/charts/ConvergencePanel.tsx'; @@ -163,6 +164,15 @@ const App: Component = () => { + +

Display Options

+ + +

Run Controls

diff --git a/apps/cadecon/src/components/controls/AlgorithmSettings.tsx b/apps/cadecon/src/components/controls/AlgorithmSettings.tsx index 672a4460..7c11bac8 100644 --- a/apps/cadecon/src/components/controls/AlgorithmSettings.tsx +++ b/apps/cadecon/src/components/controls/AlgorithmSettings.tsx @@ -9,6 +9,8 @@ import { setHpFilterEnabled, lpFilterEnabled, setLpFilterEnabled, + noiseConstrained, + setNoiseConstrained, maxIterations, setMaxIterations, convergenceTol, @@ -87,6 +89,14 @@ export function AlgorithmSettings(): JSX.Element { onChange={setLpFilterEnabled} disabled={isRunLocked()} /> + + ); diff --git a/apps/cadecon/src/components/controls/DisplaySettings.tsx b/apps/cadecon/src/components/controls/DisplaySettings.tsx new file mode 100644 index 00000000..a06d204e --- /dev/null +++ b/apps/cadecon/src/components/controls/DisplaySettings.tsx @@ -0,0 +1,26 @@ +import type { JSX } from 'solid-js'; +import { ToggleSwitch } from './ToggleSwitch.tsx'; +import { sparsityCompareEnabled, setSparsityCompareEnabled } from '../../lib/algorithm-store.ts'; +import { isRunLocked } from '../../lib/iteration-store.ts'; + +/** + * Display / diagnostic options. These do not change the deconvolution result — + * they add extra material for visual inspection. The sparsity-comparison overlay + * must be enabled BEFORE a run because it is computed during the run, so it lives + * with the pre-run controls rather than being a purely post-run view toggle. + */ +export function DisplaySettings(): JSX.Element { + return ( +
+
+ +
+
+ ); +} diff --git a/apps/cadecon/src/components/traces/TraceInspector.tsx b/apps/cadecon/src/components/traces/TraceInspector.tsx index 6188710a..178227d3 100644 --- a/apps/cadecon/src/components/traces/TraceInspector.tsx +++ b/apps/cadecon/src/components/traces/TraceInspector.tsx @@ -53,9 +53,11 @@ import { setShowGTCalcium, showGTSpikes, setShowGTSpikes, + showSparsityCompare, + setShowSparsityCompare, viewedIteration, } from '../../lib/viz-store.ts'; -import { upsampleFactor } from '../../lib/algorithm-store.ts'; +import { upsampleFactor, noiseConstrained } from '../../lib/algorithm-store.ts'; import { subsetRectangles, selectedSubsetIdx } from '../../lib/subset-store.ts'; import { createGroundTruthCalciumSeries, @@ -72,6 +74,8 @@ const RESID_GAP_FRAC = 0.05; const RESID_SCALE = 0.25; const TRANSIENT_TAU_MULTIPLIER = 2; const TRACE_INSPECTOR_ZOOM_WINDOW_S = 60; +// Distinct dashed color for the sparsity-comparison overlay (the OPPOSITE setting). +const COMPARE_COLOR = '#e040fb'; interface BandLayout { deconvTop: number; @@ -211,6 +215,18 @@ export function TraceInspector(): JSX.Element { (): Float32Array | null => effectiveResult()?.filteredTrace ?? null, ); + // Sparsity-comparison overlay: the opposite-setting spike counts, precomputed + // per iteration during the run (when the comparison option is enabled) and + // stored on the result. Read-only here — no solving on the browse path. + const comparisonDeconv = createMemo( + (): Float32Array | null => effectiveResult()?.comparisonSCounts ?? null, + ); + const hasComparison = createMemo(() => comparisonDeconv() != null); + + // Labels for the two deconv traces, by their actual setting. + const deconvLabel = () => (noiseConstrained() ? 'Deconv (noise-constr.)' : 'Deconv (standard)'); + const compareLabel = () => (noiseConstrained() ? 'Deconv (standard)' : 'Deconv (noise-constr.)'); + // Zoom window state const totalDuration = createMemo(() => { const raw = fullRawTrace(); @@ -311,7 +327,7 @@ export function TraceInspector(): JSX.Element { }); }; - const EMPTY_DATA: [number[], ...number[][]] = [[], [], [], [], [], [], [], []]; + const EMPTY_DATA: [number[], ...number[][]] = [[], [], [], [], [], [], [], [], []]; const DOWNSAMPLE_BUCKETS = 600; const zoomData = createMemo<[number[], ...number[][]]>(() => { @@ -378,6 +394,17 @@ export function TraceInspector(): JSX.Element { dsDeconv = new Array(dsX.length).fill(null) as number[]; } + // Sparsity comparison — opposite-setting deconv, same deconv band + const comp = comparisonDeconv(); + let dsCompare: number[]; + if (comp && comp.length >= endSample) { + const compSlice = comp.subarray(startSample, endSample); + const [, dsCompRaw] = downsampleMinMax(x, compSlice, DOWNSAMPLE_BUCKETS); + dsCompare = scaleToDeconvBand(dsCompRaw, rawMin, rawMax); + } else { + dsCompare = new Array(dsX.length).fill(null) as number[]; + } + // Residual — compute against the working trace (what the solver actually fit) const residSource = isFiltered ? (dsFiltered as number[]) : dsRaw; const dsResid = computeResiduals(residSource, dsFit, rawMin, rawMax, dsX.length); @@ -409,6 +436,7 @@ export function TraceInspector(): JSX.Element { dsResid, dsGTCalcium, dsGTSpikes, + dsCompare, ]; }); @@ -424,10 +452,17 @@ export function TraceInspector(): JSX.Element { { label: 'Raw', stroke: TRACE_COLORS.raw, width: 1, show: showRaw() }, { label: 'Filtered', stroke: TRACE_COLORS.filtered, width: 1.5, show: showFiltered() }, { label: 'Fit', stroke: TRACE_COLORS.fit, width: 1.5, show: showFit() }, - { label: 'Deconv', stroke: TRACE_COLORS.deconv, width: 1, show: showDeconv() }, + { label: deconvLabel(), stroke: TRACE_COLORS.deconv, width: 1, show: showDeconv() }, { label: 'Residual', stroke: TRACE_COLORS.resid, width: 1, show: showResidual() }, gtCaSeries, gtSpkSeries, + { + label: compareLabel(), + stroke: COMPARE_COLOR, + width: 1, + dash: [4, 3], + show: showSparsityCompare() && hasComparison(), + }, ]; }); @@ -458,7 +493,7 @@ export function TraceInspector(): JSX.Element { { key: 'deconv', color: TRACE_COLORS.deconv, - label: 'Deconv', + label: deconvLabel(), visible: showDeconv, setVisible: setShowDeconv, }, @@ -470,6 +505,16 @@ export function TraceInspector(): JSX.Element { setVisible: setShowResidual, }, ]; + if (hasComparison()) { + items.push({ + key: 'compare', + color: COMPARE_COLOR, + label: compareLabel(), + visible: showSparsityCompare, + setVisible: setShowSparsityCompare, + dashed: true, + }); + } if (gtVisible()) { items.push( { @@ -502,6 +547,11 @@ export function TraceInspector(): JSX.Element { const r = effectiveResult(); return r ? r.sCounts.reduce((s, v) => s + v, 0).toFixed(0) : '--'; }; + const comparisonSpikeCount = () => { + const c = comparisonDeconv(); + return c ? c.reduce((s, v) => s + v, 0).toFixed(0) : null; + }; + const settingShort = (nc: boolean) => (nc ? 'noise-constr.' : 'standard'); // Subset highlight zones for the minimap — show which time regions // the algorithm operates on for the currently selected cell. @@ -545,7 +595,15 @@ export function TraceInspector(): JSX.Element {
alpha: {alpha()} PVE: {pve()} - spikes: {spikeCount()} + spikes: {spikeCount()}} + > + + spikes: {spikeCount()} ({settingShort(noiseConstrained())}) · {comparisonSpikeCount()}{' '} + ({settingShort(!noiseConstrained())}) + +
diff --git a/apps/cadecon/src/lib/algorithm-store.ts b/apps/cadecon/src/lib/algorithm-store.ts index d14c0f93..bb96c609 100644 --- a/apps/cadecon/src/lib/algorithm-store.ts +++ b/apps/cadecon/src/lib/algorithm-store.ts @@ -7,6 +7,14 @@ import { samplingRate } from './data-store.ts'; const [upsampleTarget, setUpsampleTarget] = createSignal(300); const [hpFilterEnabled, setHpFilterEnabled] = createSignal(true); const [lpFilterEnabled, setLpFilterEnabled] = createSignal(true); +// Noise-constrained sparsity: pick the sparsest spike support whose residual +// still reaches the data-derived noise floor (no tuning knob). Suppresses noise +// fit as spurious spikes; benefit concentrates at low SNR. Off = current max-PVE. +const [noiseConstrained, setNoiseConstrained] = createSignal(true); +// Opt-in: during each iteration, also solve every trace with the OPPOSITE +// noise-constrained setting and store it, so the Trace Inspector can overlay the +// two instantly (no live solve). Doubles inference cost — off by default. +const [sparsityCompareEnabled, setSparsityCompareEnabled] = createSignal(false); const [maxIterations, setMaxIterations] = createSignal(20); // Convergence is tested in kernel SHAPE space (peak time + FWHM). convergenceTol @@ -53,6 +61,10 @@ export { setHpFilterEnabled, lpFilterEnabled, setLpFilterEnabled, + noiseConstrained, + setNoiseConstrained, + sparsityCompareEnabled, + setSparsityCompareEnabled, maxIterations, setMaxIterations, convergenceTol, diff --git a/apps/cadecon/src/lib/cadecon-pool.ts b/apps/cadecon/src/lib/cadecon-pool.ts index 4f86ba43..4f762610 100644 --- a/apps/cadecon/src/lib/cadecon-pool.ts +++ b/apps/cadecon/src/lib/cadecon-pool.ts @@ -26,6 +26,8 @@ interface TraceJobFields { hpEnabled: boolean; lpEnabled: boolean; lambda: number; + noiseConstrained: boolean; + computeComparison: boolean; warmCounts?: Float32Array; onComplete(result: TraceResult): void; } @@ -120,6 +122,8 @@ const caDeconRouter: MessageRouter = { hpEnabled: job.hpEnabled, lpEnabled: job.lpEnabled, lambda: job.lambda, + noiseConstrained: job.noiseConstrained, + computeComparison: job.computeComparison, warmCounts: warmCopy, }, transfers, diff --git a/apps/cadecon/src/lib/iteration-manager.ts b/apps/cadecon/src/lib/iteration-manager.ts index af68d46f..1a510ebc 100644 --- a/apps/cadecon/src/lib/iteration-manager.ts +++ b/apps/cadecon/src/lib/iteration-manager.ts @@ -43,6 +43,8 @@ import { finalSelectionWindow, hpFilterEnabled, lpFilterEnabled, + noiseConstrained, + sparsityCompareEnabled, traceFistaMaxIters, traceFistaTol, kernelFistaMaxIters, @@ -164,6 +166,8 @@ function dispatchTraceJobs( hpEnabled: boolean, lpEnabled: boolean, lambda: number, + noiseConstrained: boolean, + computeComparison: boolean, prevResults?: Map, ): Promise>> { return new Promise((resolve) => { @@ -211,6 +215,8 @@ function dispatchTraceJobs( hpEnabled, lpEnabled, lambda, + noiseConstrained, + computeComparison, warmCounts, onComplete(result: TraceResult) { results[subsetIdx].set(cell, result); @@ -447,6 +453,8 @@ export async function startRun(): Promise { const hpOn = hpFilterEnabled(); const lpOn = lpFilterEnabled(); const sparsityLambda = 0.0; + const noiseConstrainedOn = noiseConstrained(); + const computeComparison = sparsityCompareEnabled(); // Create pool pool = createCaDeconWorkerPool(); @@ -583,6 +591,8 @@ export async function startRun(): Promise { hpOn, lpOn, sparsityLambda, + noiseConstrainedOn, + computeComparison, prevTraceCounts, ); @@ -601,6 +611,8 @@ export async function startRun(): Promise { >(); // Map cell → full-length filtered trace (stitched from subset windows) const cellFiltered = new Map(); + // Map cell → full-length opposite-setting counts (comparison overlay) + const cellComparison = new Map(); const batchEntries: Record = {}; for (let si = 0; si < rects.length; si++) { const rect = rects[si]; @@ -621,6 +633,15 @@ export async function startRun(): Promise { } fullFilt.set(result.filteredTrace, rect.tStart); } + // Stitch comparison counts subset windows into full-length arrays + if (result.comparisonSCounts) { + let fullCmp = cellComparison.get(cell); + if (!fullCmp) { + fullCmp = new Float32Array(nTp); + cellComparison.set(cell, fullCmp); + } + fullCmp.set(result.comparisonSCounts, rect.tStart); + } cellScalars.set(cell, { alpha: result.alpha, baseline: result.baseline, @@ -638,6 +659,7 @@ export async function startRun(): Promise { baseline: result.baseline, threshold: result.threshold, pve: result.pve, + comparisonSCounts: result.comparisonSCounts, }; } } @@ -656,6 +678,7 @@ export async function startRun(): Promise { baseline: scalars.baseline, threshold: scalars.threshold, pve: scalars.pve, + comparisonSCounts: cellComparison.get(cell), }; } @@ -907,6 +930,8 @@ export async function startRun(): Promise { hpEnabled: hpOn, lpEnabled: lpOn, lambda: sparsityLambda, + noiseConstrained: noiseConstrainedOn, + computeComparison, warmCounts, onComplete(result: TraceResult) { batch(() => { @@ -919,6 +944,7 @@ export async function startRun(): Promise { baseline: result.baseline, threshold: result.threshold, pve: result.pve, + comparisonSCounts: result.comparisonSCounts, }); finCompleted++; setCompletedSubsetTraceJobs(finCompleted); diff --git a/apps/cadecon/src/lib/iteration-store.ts b/apps/cadecon/src/lib/iteration-store.ts index 34e8e633..79e131d3 100644 --- a/apps/cadecon/src/lib/iteration-store.ts +++ b/apps/cadecon/src/lib/iteration-store.ts @@ -59,6 +59,8 @@ export interface TraceResultEntry { baseline: number; threshold: number; pve: number; + /** Spike counts from the opposite noise-constrained setting (comparison overlay). */ + comparisonSCounts?: Float32Array; } function cellSubsetKey(cellIndex: number, subsetIdx: number): string { diff --git a/apps/cadecon/src/lib/tutorial/content/04-interpreting.ts b/apps/cadecon/src/lib/tutorial/content/04-interpreting.ts index 88cfc587..09ce4ad9 100644 --- a/apps/cadecon/src/lib/tutorial/content/04-interpreting.ts +++ b/apps/cadecon/src/lib/tutorial/content/04-interpreting.ts @@ -51,7 +51,17 @@ export const interpretingTutorial: Tutorial = { side: 'top', popoverClass: 'driver-popover--wide', }, - // Step 6: When to change algorithm settings + // Step 6: Sparsity comparison overlay (display-only diagnostic) + { + element: '[data-tutorial="display-options"]', + title: 'Seeing the Sparsity’s Impact', + description: + 'CaDecon runs noise-constrained sparsity by default: it keeps only the spikes needed to explain the trace down to its noise floor, which cleans up spurious activity — the effect is largest on low-SNR cells. To see how much it is doing, enable Sparsity Comparison Overlay here before a run. The Trace Inspector then overlays the deconvolution computed both with and without the sparsity, and shows both spike counts, so you can judge the impact per cell.

' + + 'This is a visual-inspection aid only — it does not change the result. Because it solves every trace twice, it roughly doubles run time, so it is off by default. Leave it off for normal runs.', + side: 'right', + popoverClass: 'driver-popover--wide', + }, + // Step 7: When to change algorithm settings { element: '[data-tutorial="algorithm-settings"]', title: 'When to Adjust Settings', diff --git a/apps/cadecon/src/lib/viz-store.ts b/apps/cadecon/src/lib/viz-store.ts index 4103cb76..df2d0495 100644 --- a/apps/cadecon/src/lib/viz-store.ts +++ b/apps/cadecon/src/lib/viz-store.ts @@ -19,6 +19,11 @@ const [showResidual, setShowResidual] = createSignal(false); const [showGTCalcium, setShowGTCalcium] = createSignal(true); const [showGTSpikes, setShowGTSpikes] = createSignal(true); +// Sparsity comparison overlay: when on, TraceInspector re-runs inference for the +// inspected cell with the OPPOSITE noise-constrained setting (same converged +// kernel) and overlays it, so the impact of noise-constrained sparsity is visible. +const [showSparsityCompare, setShowSparsityCompare] = createSignal(true); + export { viewedIteration, setViewedIteration, @@ -38,4 +43,6 @@ export { setShowGTCalcium, showGTSpikes, setShowGTSpikes, + showSparsityCompare, + setShowSparsityCompare, }; diff --git a/apps/cadecon/src/styles/layout.css b/apps/cadecon/src/styles/layout.css index 7ca4962c..c0f9c80c 100644 --- a/apps/cadecon/src/styles/layout.css +++ b/apps/cadecon/src/styles/layout.css @@ -1,3 +1,12 @@ +/* Controls sidebar: keep each panel at its natural height and let the sidebar + scroll, instead of letting flexbox compress the panels together when the + viewport is short or the browser text size is large. Without this, the flex + column shrinks every panel (default flex-shrink: 1), squishing and overlapping + their contents even though a scrollbar is present. */ +.viz-layout__sidebar > * { + flex-shrink: 0; +} + /* 2-column, 2-row resizable visualization grid */ .viz-grid { diff --git a/apps/cadecon/src/workers/cadecon-types.ts b/apps/cadecon/src/workers/cadecon-types.ts index 40410004..20a3ddd4 100644 --- a/apps/cadecon/src/workers/cadecon-types.ts +++ b/apps/cadecon/src/workers/cadecon-types.ts @@ -10,6 +10,9 @@ export interface TraceResult { pve: number; iterations: number; converged: boolean; + /** Spike counts from the OPPOSITE noise-constrained setting, for the comparison + * overlay. Present only when the trace job requested `computeComparison`. */ + comparisonSCounts?: Float32Array; } /** Results from peak-seeded spike detection on a single trace. */ @@ -54,6 +57,12 @@ export type CaDeconWorkerInbound = lpEnabled: boolean; /** L1 sparsity penalty on spike solution. */ lambda: number; + /** Noise-constrained threshold selection: choose the sparsest spike support + * whose residual meets the data-derived noise floor (no tuning knob). */ + noiseConstrained: boolean; + /** Also solve with the OPPOSITE noise-constrained setting and return it as + * comparisonSCounts (for the teaching/impact overlay). */ + computeComparison: boolean; /** Previous iteration's s_counts at original rate for warm-start. */ warmCounts?: Float32Array; } diff --git a/apps/cadecon/src/workers/cadecon-worker.ts b/apps/cadecon/src/workers/cadecon-worker.ts index cd55fcb3..860b49ca 100644 --- a/apps/cadecon/src/workers/cadecon-worker.ts +++ b/apps/cadecon/src/workers/cadecon-worker.ts @@ -37,6 +37,7 @@ function handleTraceJob(req: Extract f64 { + let n = raw_trace.len(); + if n < 8 { + return 0.0; + } + let mean = raw_trace.iter().map(|&v| v as f64).sum::() / n as f64; + let mut input: Vec = raw_trace.iter().map(|&v| v - mean as f32).collect(); + + let mut planner = RealFftPlanner::::new(); + let r2c = planner.plan_fft_forward(n); + let mut spectrum = r2c.make_output_vec(); + if r2c.process(&mut input, &mut spectrum).is_err() { + return 0.0; + } + + let nyq = spectrum.len(); // n/2 + 1 + let lo = nyq / 2; + let mut acc = 0.0_f64; + let mut cnt = 0usize; + for c in &spectrum[lo..nyq] { + acc += (c.re as f64 * c.re as f64 + c.im as f64 * c.im as f64) / n as f64; + cnt += 1; + } + if cnt == 0 { + return 0.0; + } + (acc / cnt as f64).max(0.0).sqrt() +} + +fn variance(x: &[f32]) -> f64 { + let n = x.len(); + if n == 0 { + return 0.0; + } + let mean = x.iter().map(|&v| v as f64).sum::() / n as f64; + x.iter().map(|&v| (v as f64 - mean).powi(2)).sum::() / n as f64 +} + +/// Per-sample noise std of the FILTERED trace at original-grid positions — the +/// quantity the residual budget needs. Estimates the raw noise std from the +/// unfiltered trace (LP-immune) and scales by the empirically-measured +/// noise-variance gain of the actual upsample→HP/LP chain: a deterministic +/// white probe is pushed through the same path and its grid-sample variance +/// ratio gives the gain. This makes the budget track whatever noise survives +/// the filters, regardless of where the (kernel-derived) LP cutoff falls. +/// (The rolling-baseline subtraction is a mild high-pass and is not replicated +/// in the probe; its effect on noise variance is negligible.) +fn estimate_grid_noise_sigma( + raw_trace: &[f32], + upsample_factor: usize, + tau_r: f64, + tau_d: f64, + fs_up: f64, + hp: bool, + lp: bool, +) -> f64 { + let sigma_raw = high_band_sigma(raw_trace); + if sigma_raw <= 0.0 { + return 0.0; + } + // No filtering: grid samples equal the raw samples, so gain is 1. + if !hp && !lp { + return sigma_raw; + } + + // Deterministic unit-ish white probe (LCG); absolute scale cancels in the ratio. + let n = raw_trace.len(); + let mut probe = vec![0.0_f32; n]; + let mut state: u64 = 0x9E3779B97F4A7C15; + for v in probe.iter_mut() { + state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + *v = (((state >> 33) as f64) / ((1u64 << 31) as f64) - 1.0) as f32; + } + + let up = upsample_trace(&probe, upsample_factor); + let mut solver = Solver::new(); + solver.set_conv_mode(ConvMode::BandedAR2); + solver.set_params(tau_r, tau_d, 0.0, fs_up); + solver.set_trace(&up); + solver.set_hp_filter_enabled(hp); + solver.set_lp_filter_enabled(lp); + solver.apply_filter(); + let filt = solver.get_trace(); + + let grid: Vec = filt + .iter() + .step_by(upsample_factor.max(1)) + .copied() + .collect(); + let var_probe = variance(&probe); + if var_probe <= 0.0 { + return sigma_raw; + } + let gain = (variance(&grid) / var_probe).max(0.0); + sigma_raw * gain.sqrt() +} #[cfg_attr(feature = "jsbindings", derive(serde::Serialize))] pub struct InDecaResult { @@ -195,6 +317,7 @@ fn interior_peak(s: &[f32], pad: usize) -> f32 { /// `warm_counts`: optional spike counts from a previous iteration at the **original** /// sampling rate. These are upsampled to a binary trace at the upsampled rate and /// used as FISTA warm-start, which typically reduces iterations by 30-60%. +#[allow(clippy::too_many_arguments)] pub fn solve_trace( trace: &[f32], tau_r: f64, @@ -207,6 +330,39 @@ pub fn solve_trace( hp_enabled: bool, lp_enabled: bool, lambda: f64, +) -> InDecaResult { + solve_trace_opts( + trace, + tau_r, + tau_d, + fs, + upsample_factor, + max_iters, + tol, + warm_counts, + hp_enabled, + lp_enabled, + lambda, + SolveOptions::default(), + ) +} + +/// See [`solve_trace`]; adds optional [`SolveOptions`] for the count-debiasing +/// experiments (noise-constrained threshold selection and run collapse). +#[allow(clippy::too_many_arguments)] +pub fn solve_trace_opts( + trace: &[f32], + tau_r: f64, + tau_d: f64, + fs: f64, + upsample_factor: usize, + max_iters: u32, + tol: f64, + warm_counts: Option<&[f32]>, + hp_enabled: bool, + lp_enabled: bool, + lambda: f64, + opts: SolveOptions, ) -> InDecaResult { let fs_up = fs * upsample_factor as f64; let upsampled = upsample_trace(trace, upsample_factor); @@ -257,6 +413,26 @@ pub fn solve_trace( let banded = BandedAR2::new(tau_r, tau_d, fs_up); + // Noise-constrained threshold selection needs the per-sample noise std of the + // filtered trace at the original grid. Estimated LP-cutoff-agnostically from + // the raw trace's high band scaled by the filter chain's measured noise gain. + // Fully data-derived — no knob. + let selection = if opts.noise_constrained { + Selection::NoiseFloor { + sigma: estimate_grid_noise_sigma( + trace, + upsample_factor, + tau_r, + tau_d, + fs_up, + hp_enabled, + lp_enabled, + ), + } + } else { + Selection::MaxPve + }; + // ── Step 3: Scale iteration loop ──────────────────────────────────── // Each round: prescale by alpha_est → Box[0,1] FISTA → threshold search // against the *original* trace → lstsq recovers alpha directly. @@ -325,7 +501,7 @@ pub fn solve_trace( threshold, pve, .. - } = threshold_search( + } = threshold_search_opts( s_norm_slice, &working_trace, &banded, @@ -333,6 +509,7 @@ pub fn solve_trace( fs_up, upsample_factor, f64::INFINITY, + selection, ); // Track the best result by PVE. @@ -369,6 +546,14 @@ pub fn solve_trace( (vec![0.0; wt_len], 0.0, 0.0, 0.0, 0.0, 0, false) }); + // Optionally collapse each smeared bump to a single event before summing, + // so the bin-sum reports events rather than run length. + let s_binary = if opts.collapse_runs { + collapse_runs(&s_binary) + } else { + s_binary + }; + // Downsample binary spike train to original rate using centered bins let s_counts = downsample_binary(&s_binary, upsample_factor); diff --git a/crates/solver/src/js_indeca.rs b/crates/solver/src/js_indeca.rs index 001bd4fb..dc9a0a74 100644 --- a/crates/solver/src/js_indeca.rs +++ b/crates/solver/src/js_indeca.rs @@ -34,6 +34,7 @@ pub fn indeca_solve_trace( lp_enabled: bool, warm_counts: &[f32], lambda: f64, + noise_constrained: bool, ) -> Result { if let Some(i) = crate::first_nonfinite(trace) { return Err(JsError::new(&format!( @@ -45,7 +46,7 @@ pub fn indeca_solve_trace( } else { Some(warm_counts) }; - let result = indeca::solve_trace( + let result = indeca::solve_trace_opts( trace, tau_r, tau_d, @@ -57,6 +58,10 @@ pub fn indeca_solve_trace( hp_enabled, lp_enabled, lambda, + indeca::SolveOptions { + noise_constrained, + collapse_runs: false, + }, ); Ok(serde_wasm_bindgen::to_value(&result).unwrap_or(JsValue::NULL)) } diff --git a/crates/solver/src/py_api.rs b/crates/solver/src/py_api.rs index 6e744380..838165ed 100644 --- a/crates/solver/src/py_api.rs +++ b/crates/solver/src/py_api.rs @@ -423,7 +423,8 @@ fn seed_kernel_estimate<'py>( /// /// Returns (s_counts, alpha, baseline, threshold, pve, iterations, converged). #[pyfunction] -#[pyo3(signature = (trace, tau_rise, tau_decay, fs, upsample_factor=1, max_iters=500, tol=1e-4, hp_enabled=false, lp_enabled=false, warm_counts=None, lambda_=0.0))] +#[pyo3(signature = (trace, tau_rise, tau_decay, fs, upsample_factor=1, max_iters=500, tol=1e-4, hp_enabled=false, lp_enabled=false, warm_counts=None, lambda_=0.0, noise_constrained=false, collapse_runs=false))] +#[allow(clippy::too_many_arguments)] fn py_indeca_solve_trace<'py>( py: Python<'py>, trace: PyReadonlyArray1, @@ -437,6 +438,8 @@ fn py_indeca_solve_trace<'py>( lp_enabled: bool, warm_counts: Option>, lambda_: f64, + noise_constrained: bool, + collapse_runs: bool, ) -> PyResult<( Bound<'py, PyArray1>, // s_counts f64, // alpha @@ -449,7 +452,7 @@ fn py_indeca_solve_trace<'py>( let trace_f32 = to_f32_vec(&trace)?; let warm = optional_to_f32_vec(warm_counts)?; - let result = indeca::solve_trace( + let result = indeca::solve_trace_opts( &trace_f32, tau_rise, tau_decay, @@ -461,6 +464,10 @@ fn py_indeca_solve_trace<'py>( hp_enabled, lp_enabled, lambda_, + indeca::SolveOptions { + noise_constrained, + collapse_runs, + }, ); Ok(( diff --git a/crates/solver/src/threshold.rs b/crates/solver/src/threshold.rs index eda60526..2433937d 100644 --- a/crates/solver/src/threshold.rs +++ b/crates/solver/src/threshold.rs @@ -18,6 +18,20 @@ pub struct ThresholdResult { pub error: f64, } +/// How the binarization threshold is chosen. +#[derive(Clone, Copy)] +pub enum Selection { + /// Current behavior: threshold that minimizes reconstruction error (max PVE). + /// This overfits — it drives the residual below the noise floor. + MaxPve, + /// Noise-constrained: the sparsest support (highest threshold) whose residual + /// at original-rate grid positions stays within the noise floor + /// `sigma^2 * n_grid_interior`. `sigma` is the noise std of the filtered trace + /// at original rate (data-derived, no tuning knob). Realizes an OASIS-style + /// "don't fit below the noise" sparsity at the stage where counts are decided. + NoiseFloor { sigma: f64 }, +} + /// Compute boundary padding for threshold search: ceil(2 * tau_d * fs_up). /// Used to exclude edge effects from the error computation. pub fn boundary_padding(tau_decay: f64, fs_up: f64) -> usize { @@ -39,6 +53,30 @@ pub fn threshold_search( fs_up: f64, upsample_factor: usize, max_alpha: f64, +) -> ThresholdResult { + threshold_search_opts( + s_relaxed, + y, + banded, + tau_decay, + fs_up, + upsample_factor, + max_alpha, + Selection::MaxPve, + ) +} + +/// Threshold search with a selectable criterion. See [`Selection`]. +#[allow(clippy::too_many_arguments)] +pub fn threshold_search_opts( + s_relaxed: &[f32], + y: &[f32], + banded: &BandedAR2, + tau_decay: f64, + fs_up: f64, + upsample_factor: usize, + max_alpha: f64, + selection: Selection, ) -> ThresholdResult { let n = s_relaxed.len(); let pad = boundary_padding(tau_decay, fs_up).min(n / 4); @@ -79,91 +117,109 @@ pub fn threshold_search( error: f64::INFINITY, }; - // Phase 1: Coarse search — ~50 evenly spaced thresholds - let coarse_n = 50.min(vals.len()); - let coarse_step = if vals.len() > 1 { - (vals.len() - 1) as f64 / (coarse_n - 1).max(1) as f64 - } else { - 1.0 - }; - - let mut coarse_thresholds: Vec = Vec::with_capacity(coarse_n); - for i in 0..coarse_n { - let idx = (i as f64 * coarse_step).round() as usize; - let idx = idx.min(vals.len() - 1); - coarse_thresholds.push(vals[idx] as f64); - } - coarse_thresholds.dedup_by(|a, b| (*a - *b).abs() < 1e-10); - - // Enforce minimum threshold floor - coarse_thresholds.retain(|&t| t >= min_threshold); - if coarse_thresholds.is_empty() { - // All candidates below minimum — use min_threshold as the only candidate - coarse_thresholds.push(min_threshold); - } - - let mut consecutive_increases = 0; - for &thresh in &coarse_thresholds { - let err = evaluate_threshold( + // Noise-constrained selection scans for the sparsest support within the + // noise floor; the max-PVE path keeps the original coarse→fine search. + if let Selection::NoiseFloor { sigma } = selection { + best.threshold = select_noise_floor_threshold( s_relaxed, y, banded, - thresh, pad, + upsample_factor, max_alpha, + sigma, + &vals, + min_threshold, &mut s_bin, &mut conv_buf, ); - if err < best.error { - best.error = err; - best.threshold = thresh; - consecutive_increases = 0; + } else { + // Phase 1: Coarse search — ~50 evenly spaced thresholds + let coarse_n = 50.min(vals.len()); + let coarse_step = if vals.len() > 1 { + (vals.len() - 1) as f64 / (coarse_n - 1).max(1) as f64 } else { - consecutive_increases += 1; - if consecutive_increases >= 10 { - break; - } + 1.0 + }; + + let mut coarse_thresholds: Vec = Vec::with_capacity(coarse_n); + for i in 0..coarse_n { + let idx = (i as f64 * coarse_step).round() as usize; + let idx = idx.min(vals.len() - 1); + coarse_thresholds.push(vals[idx] as f64); } - } + coarse_thresholds.dedup_by(|a, b| (*a - *b).abs() < 1e-10); - // Phase 2: Fine search — ~50 thresholds around the best coarse result - let spread = if vals.len() > 1 { - (vals[vals.len() - 1] - vals[0]) as f64 / coarse_n as f64 * 2.0 - } else { - best.threshold * 0.2 - }; - let fine_lo = (best.threshold - spread).max(min_threshold); - let fine_hi = best.threshold + spread; - let fine_n = 50; - let fine_step = (fine_hi - fine_lo) / (fine_n - 1).max(1) as f64; - - consecutive_increases = 0; - for i in 0..fine_n { - let thresh = fine_lo + i as f64 * fine_step; - if thresh < 0.0 { - continue; + // Enforce minimum threshold floor + coarse_thresholds.retain(|&t| t >= min_threshold); + if coarse_thresholds.is_empty() { + // All candidates below minimum — use min_threshold as the only candidate + coarse_thresholds.push(min_threshold); } - let err = evaluate_threshold( - s_relaxed, - y, - banded, - thresh, - pad, - max_alpha, - &mut s_bin, - &mut conv_buf, - ); - if err < best.error { - best.error = err; - best.threshold = thresh; - consecutive_increases = 0; + + let mut consecutive_increases = 0; + for &thresh in &coarse_thresholds { + let err = evaluate_threshold( + s_relaxed, + y, + banded, + thresh, + pad, + max_alpha, + &mut s_bin, + &mut conv_buf, + ); + if err < best.error { + best.error = err; + best.threshold = thresh; + consecutive_increases = 0; + } else { + consecutive_increases += 1; + if consecutive_increases >= 10 { + break; + } + } + } + + // Phase 2: Fine search — ~50 thresholds around the best coarse result + let spread = if vals.len() > 1 { + (vals[vals.len() - 1] - vals[0]) as f64 / coarse_n as f64 * 2.0 } else { - consecutive_increases += 1; - if consecutive_increases >= 10 { - break; + best.threshold * 0.2 + }; + let fine_lo = (best.threshold - spread).max(min_threshold); + let fine_hi = best.threshold + spread; + let fine_n = 50; + let fine_step = (fine_hi - fine_lo) / (fine_n - 1).max(1) as f64; + + consecutive_increases = 0; + for i in 0..fine_n { + let thresh = fine_lo + i as f64 * fine_step; + if thresh < 0.0 { + continue; + } + let err = evaluate_threshold( + s_relaxed, + y, + banded, + thresh, + pad, + max_alpha, + &mut s_bin, + &mut conv_buf, + ); + if err < best.error { + best.error = err; + best.threshold = thresh; + consecutive_increases = 0; + } else { + consecutive_increases += 1; + if consecutive_increases >= 10 { + break; + } } } - } + } // end max-PVE search branch // Final pass: compute full result at best threshold binarize(s_relaxed, best.threshold, &mut s_bin); @@ -240,6 +296,102 @@ fn evaluate_threshold( err } +/// Scan candidate thresholds and return the highest (sparsest) one whose +/// grid-residual stays within the noise budget `sigma^2 * n_grid_interior`. +/// Residual increases as the threshold rises (fewer spikes → worse fit), so the +/// highest feasible threshold is the sparsest support that still explains the +/// signal down to — but not below — the noise floor. Falls back to the +/// min-grid-residual threshold if none is feasible (budget tighter than the +/// best achievable fit). +#[allow(clippy::too_many_arguments)] +fn select_noise_floor_threshold( + s_relaxed: &[f32], + y: &[f32], + banded: &BandedAR2, + pad: usize, + upsample_factor: usize, + max_alpha: f64, + sigma: f64, + vals: &[f32], + min_threshold: f64, + s_bin: &mut [f32], + conv_buf: &mut [f32], +) -> f64 { + let n = y.len(); + let stride = upsample_factor.max(1); + let n_grid = (pad..n.saturating_sub(pad)).step_by(stride).count(); + let budget = sigma * sigma * n_grid as f64; + + // Up to 256 candidate thresholds, evenly spaced through the sorted values. + let cap = 256usize; + let step = if vals.len() > cap { + vals.len() as f64 / cap as f64 + } else { + 1.0 + }; + + let mut best_feasible = f64::NEG_INFINITY; // highest feasible threshold + let mut best_effort_thr = min_threshold; // min grid-SSE fallback + let mut best_effort_sse = f64::INFINITY; + + let mut idx = 0.0_f64; + while (idx as usize) < vals.len() { + let thr = vals[idx as usize] as f64; + idx += step; + if thr < min_threshold { + continue; + } + let sse = grid_sse_at_threshold( + s_relaxed, y, banded, thr, pad, stride, max_alpha, s_bin, conv_buf, + ); + if sse <= budget && thr > best_feasible { + best_feasible = thr; + } + if sse < best_effort_sse { + best_effort_sse = sse; + best_effort_thr = thr; + } + } + + if best_feasible.is_finite() { + best_feasible + } else { + best_effort_thr + } +} + +/// Residual SSE at original-rate grid positions only (stride = upsample_factor), +/// with alpha/baseline fit over the full interior. Evaluating the noise +/// constraint on the original grid — where the upsampled trace equals the real +/// samples — avoids the correlated, reduced-variance noise at interpolated +/// positions, so the budget can use a noise std estimated at the original rate. +#[allow(clippy::too_many_arguments)] +fn grid_sse_at_threshold( + s_relaxed: &[f32], + y: &[f32], + banded: &BandedAR2, + threshold: f64, + pad: usize, + stride: usize, + max_alpha: f64, + s_bin: &mut [f32], + conv_buf: &mut [f32], +) -> f64 { + binarize(s_relaxed, threshold, s_bin); + banded.convolve_forward(s_bin, conv_buf); + let (alpha, baseline) = lstsq_alpha_baseline(conv_buf, y, pad, max_alpha); + let n = y.len(); + let mut sse = 0.0_f64; + let mut i = pad; + while i < n.saturating_sub(pad) { + let pred = alpha * conv_buf[i] as f64 + baseline; + let d = y[i] as f64 - pred; + sse += d * d; + i += stride; + } + sse +} + /// Least-squares fit for alpha and baseline: y ≈ alpha * conv + baseline. /// Solves the 2x2 normal equations over the inner region [pad..n-pad]. /// Alpha is constrained to [0, max_alpha]. When max_alpha is f64::INFINITY diff --git a/crates/solver/src/upsample.rs b/crates/solver/src/upsample.rs index 7150d9bb..e3d96d5e 100644 --- a/crates/solver/src/upsample.rs +++ b/crates/solver/src/upsample.rs @@ -89,6 +89,34 @@ pub fn downsample_average(signal: &[f32], factor: usize) -> Vec { .collect() } +/// Collapse each maximal run of consecutive nonzero samples in an upsampled +/// binary spike train to a single spike at the run's center, zeroing the rest. +/// +/// The Box[0,1] FISTA relaxation spreads one true event into a smooth bump that +/// binarizes to a multi-sample run; `downsample_binary` would then SUM that run +/// and report one event as several spikes (count inflation ∝ run length). This +/// collapse makes each contiguous bump count as exactly one event while still +/// preserving genuinely-separated events (they are broken by zeros, so remain +/// distinct runs and are not merged). +pub fn collapse_runs(s_bin: &[f32]) -> Vec { + let n = s_bin.len(); + let mut out = vec![0.0_f32; n]; + let mut i = 0; + while i < n { + if s_bin[i] > 0.0 { + let start = i; + while i < n && s_bin[i] > 0.0 { + i += 1; + } + let center = start + (i - start - 1) / 2; + out[center] = 1.0; + } else { + i += 1; + } + } + out +} + /// Downsample a binary spike signal by bin-summing with centered bins. /// /// Each output bin is centered on the original sample position (`i * factor`) From b36a25004e587ab8b496e3df1bcd6ad66f4887bf Mon Sep 17 00:00:00 2001 From: daharoni Date: Wed, 8 Jul 2026 22:00:24 -0700 Subject: [PATCH 2/4] feat(cadecon): finish collapse-runs removal, wire noise_constrained into Python, backfill changelog Remove the last collapse-runs reference (stale doc comment), expose the noise_constrained option through the calab.solve_trace Python binding, and document it in the CaDecon guide. Restructure docs/CHANGELOG.md: split the consolidated post-2.0.6 blob into the real minor releases (2.1.0-2.4.0), add 2.5.0 (PRs #153-#167), and scope 2.6.0 to PR #168. Co-Authored-By: Claude Opus 4.8 --- crates/solver/src/indeca.rs | 21 +-- crates/solver/src/js_indeca.rs | 5 +- crates/solver/src/py_api.rs | 8 +- crates/solver/src/upsample.rs | 28 ---- docs/CHANGELOG.md | 214 +++++++++++++++++++++---------- python/docs/guides/cadecon.md | 2 + python/src/calab/_compute.py | 8 ++ python/tests/test_solve_trace.py | 9 ++ 8 files changed, 172 insertions(+), 123 deletions(-) diff --git a/crates/solver/src/indeca.rs b/crates/solver/src/indeca.rs index d4ddd645..70adb5c9 100644 --- a/crates/solver/src/indeca.rs +++ b/crates/solver/src/indeca.rs @@ -15,22 +15,19 @@ use crate::banded::BandedAR2; use crate::threshold::{threshold_search_opts, Selection, ThresholdResult}; use crate::upsample::{ - collapse_runs, downsample_average, downsample_binary, upsample_counts_to_binary, upsample_trace, + downsample_average, downsample_binary, upsample_counts_to_binary, upsample_trace, }; use crate::{Constraint, ConvMode, Solver}; use realfft::RealFftPlanner; /// Optional spike-inference behaviors, off by default (current shipping path). /// -/// Both target the ~2× spike overcount without changing the default output: -/// `noise_constrained` chooses the binarization threshold at the noise floor -/// instead of maximizing fit (suppresses low-SNR spurious events), and -/// `collapse_runs` collapses each smeared bump to one event before the -/// bin-summing downsample (removes the run-length count inflation). +/// `noise_constrained` chooses the binarization threshold at the data-derived +/// noise floor instead of maximizing fit, suppressing low-SNR spurious spikes +/// without changing the default (max-PVE) output. #[derive(Clone, Copy, Default)] pub struct SolveOptions { pub noise_constrained: bool, - pub collapse_runs: bool, } /// Raw measurement-noise std from the high-frequency band of the periodogram @@ -348,7 +345,7 @@ pub fn solve_trace( } /// See [`solve_trace`]; adds optional [`SolveOptions`] for the count-debiasing -/// experiments (noise-constrained threshold selection and run collapse). +/// experiment (noise-constrained threshold selection). #[allow(clippy::too_many_arguments)] pub fn solve_trace_opts( trace: &[f32], @@ -546,14 +543,6 @@ pub fn solve_trace_opts( (vec![0.0; wt_len], 0.0, 0.0, 0.0, 0.0, 0, false) }); - // Optionally collapse each smeared bump to a single event before summing, - // so the bin-sum reports events rather than run length. - let s_binary = if opts.collapse_runs { - collapse_runs(&s_binary) - } else { - s_binary - }; - // Downsample binary spike train to original rate using centered bins let s_counts = downsample_binary(&s_binary, upsample_factor); diff --git a/crates/solver/src/js_indeca.rs b/crates/solver/src/js_indeca.rs index dc9a0a74..036bbdf1 100644 --- a/crates/solver/src/js_indeca.rs +++ b/crates/solver/src/js_indeca.rs @@ -58,10 +58,7 @@ pub fn indeca_solve_trace( hp_enabled, lp_enabled, lambda, - indeca::SolveOptions { - noise_constrained, - collapse_runs: false, - }, + indeca::SolveOptions { noise_constrained }, ); Ok(serde_wasm_bindgen::to_value(&result).unwrap_or(JsValue::NULL)) } diff --git a/crates/solver/src/py_api.rs b/crates/solver/src/py_api.rs index 838165ed..30e00a51 100644 --- a/crates/solver/src/py_api.rs +++ b/crates/solver/src/py_api.rs @@ -423,7 +423,7 @@ fn seed_kernel_estimate<'py>( /// /// Returns (s_counts, alpha, baseline, threshold, pve, iterations, converged). #[pyfunction] -#[pyo3(signature = (trace, tau_rise, tau_decay, fs, upsample_factor=1, max_iters=500, tol=1e-4, hp_enabled=false, lp_enabled=false, warm_counts=None, lambda_=0.0, noise_constrained=false, collapse_runs=false))] +#[pyo3(signature = (trace, tau_rise, tau_decay, fs, upsample_factor=1, max_iters=500, tol=1e-4, hp_enabled=false, lp_enabled=false, warm_counts=None, lambda_=0.0, noise_constrained=false))] #[allow(clippy::too_many_arguments)] fn py_indeca_solve_trace<'py>( py: Python<'py>, @@ -439,7 +439,6 @@ fn py_indeca_solve_trace<'py>( warm_counts: Option>, lambda_: f64, noise_constrained: bool, - collapse_runs: bool, ) -> PyResult<( Bound<'py, PyArray1>, // s_counts f64, // alpha @@ -464,10 +463,7 @@ fn py_indeca_solve_trace<'py>( hp_enabled, lp_enabled, lambda_, - indeca::SolveOptions { - noise_constrained, - collapse_runs, - }, + indeca::SolveOptions { noise_constrained }, ); Ok(( diff --git a/crates/solver/src/upsample.rs b/crates/solver/src/upsample.rs index e3d96d5e..7150d9bb 100644 --- a/crates/solver/src/upsample.rs +++ b/crates/solver/src/upsample.rs @@ -89,34 +89,6 @@ pub fn downsample_average(signal: &[f32], factor: usize) -> Vec { .collect() } -/// Collapse each maximal run of consecutive nonzero samples in an upsampled -/// binary spike train to a single spike at the run's center, zeroing the rest. -/// -/// The Box[0,1] FISTA relaxation spreads one true event into a smooth bump that -/// binarizes to a multi-sample run; `downsample_binary` would then SUM that run -/// and report one event as several spikes (count inflation ∝ run length). This -/// collapse makes each contiguous bump count as exactly one event while still -/// preserving genuinely-separated events (they are broken by zeros, so remain -/// distinct runs and are not merged). -pub fn collapse_runs(s_bin: &[f32]) -> Vec { - let n = s_bin.len(); - let mut out = vec![0.0_f32; n]; - let mut i = 0; - while i < n { - if s_bin[i] > 0.0 { - let start = i; - while i < n && s_bin[i] > 0.0 { - i += 1; - } - let center = start + (i - start - 1) / 2; - out[center] = 1.0; - } else { - i += 1; - } - } - out -} - /// Downsample a binary spike signal by bin-summing with centered bins. /// /// Each output bin is centered on the original sample position (`i * factor`) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a4984065..761c7ffe 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,83 +3,92 @@ Repo-level changelog for the CaLab monorepo. Uses [Keep a Changelog](https://keepachangelog.com/) format. Versions correspond to git tags (`v*`) and apply to the entire monorepo. -## [Unreleased] +## [2.6.0] -> This section accumulates every monorepo change since `v2.0.6` (PRs #58–#165); -> there has been no version tag in this window. Entries for PRs #60–#155 were -> reconstructed from git history (the changelog had lapsed) and consolidate -> closely-related PRs into single bullets for readability. +> Unreleased. Covers every change since `v2.5.0` (PR #168). + +### Added + +- **CaDecon** noise-constrained sparsity — an optional `noise_constrained` + spike-inference mode that picks the binarization threshold as the sparsest + spike support whose residual still reaches the data-derived noise floor, + instead of the fit-maximizing threshold. Knob-free and off by default; + suppresses spurious low-SNR spikes. Exposed through the WASM solver, the + CaDecon UI, and the `calab.solve_trace` Python binding (PR #168) + +## [2.5.0] - 2026-07-08 + +> Covers PRs #153–#167 (all merged 2026-07-08). ### Added -- **CaDecon** — a new app for automated calcium deconvolution (the InDeCa - algorithm) that estimates the kernel and deconvolution parameters directly - from the data, no manual tuning required: app scaffold + data loading + - subset UI, the InDeCa compute engine with warm-start, visualization / QC - distributions / drill-down, community-database integration, ground-truth - overlay, and reaching functional parity with the reference InDeCa - implementation (PRs #85, #86, #87, #88, #90, #91, #94) -- CaDecon convergence redesign — converge in kernel **shape space** (peak time + - FWHM asymptote) with median-tail kernel selection (PR #154), plus an - **asymptote dashboard** charting the four convergence signals (PR #155); - earlier migrated the kernel parameterization from (tau_rise, tau_decay) to - (t_peak, FWHM) (PR #104) - **CaDecon** bi-exponential fit outcome surfaced as `FitMode` (`TwoComponent` / `SlowOnly` / `Degenerate` / `Empty`) on the kernel result; Python `fit_biexponential` now returns an 8-tuple (trailing `fit_mode` string) and `BiexpFitResult` gained a `fit_mode` field; KernelDisplay warns when subset fits are degenerate (PR #162) -- **`calab` Python package** greatly expanded — PyO3 bindings, CaImAn/Minian - loaders, browser bridge, and a CLI (PR #66); CaDecon Python bridge with - config, autorun, progress, and auto-export (PRs #108, #109); headless-browser - batch mode + InDeCa PyO3 bindings (PR #110) -- Shared Rust **simulation module** producing synthetic ground-truth traces, - exposed to both Python and WASM (PR #113) -- **Usage-analytics pipeline + admin dashboard** with breakdowns and bulk - moderation (PRs #62, #69), extended to track CaDecon submissions (PR #93) -- Shared **auth menu** in the header across all CaLab apps (PR #61) -- Community: highlight your own submissions in the scatter plot (PR #63); - DataSource tracking + bridge export button & heartbeat detection (PR #67) -- Solver: banded AR(2) O(T) convolution + box constraint (PR #78); peak-seeded - initial-kernel auto-estimation (PR #103); an independent fast component in the - bi-exponential fit (PR #105); a `skip` parameter for bi-exponential fitting - (PR #99) -- Dynamic worker-pool scaling with a URL override (PR #77) +- **CaDecon** convergence redesign — converge in kernel **shape space** (peak + time + FWHM asymptote) with median-tail kernel selection and both filters on + by default (PR #154), plus an **asymptote dashboard** charting the four + convergence signals (PR #155) - Shared uPlot chart primitives in `@calab/ui/chart` — colorblind-safe Okabe-Ito palette (`TRACE_COLORS`, `GROUND_TRUTH_COLORS`, `KERNEL_FIT_COLORS`, `METRIC_COLORS`, `subsetColor`), viridis colormap (`VIRIDIS_LUT`, `viridisRGB`/`viridisCss`), tick math (`niceTicks`), and axis/cursor/range helpers (`chartAxis`, `labeledAxis`, `syncCursor`, `safeRange`) (PRs #158, #159, #160) -- Chart/UX: transient-zone visual indicator (PR #81); draggable minimap edges on - the trace overview (PR #145) -- Tutorials: CaDecon tutorial set (PR #151); Python Package tutorial (PR #68); - Python syntax highlighting in code blocks (PR #75) -- Sphinx + ReadTheDocs documentation site for the Python package (PR #115) -- Comprehensive README files across all packages and apps (PR #58) ### Changed -- Moved the Rust solver to `crates/solver/` with dual WASM (`jsbindings`) / - PyO3 (`pybindings`) Cargo features (PR #65) -- Made `@calab/community` app-agnostic (PR #60) - **CaDecon** raster overview uses the shared viridis colormap and drops the intensity colorbar (activity is assumed to span 0→full; absolute values are not meaningful) (PR #159) -- Replaced the export-to-Python page with a dismissible modal (PR #79) -- CaDecon left-sidebar layout/UX (PR #89); convergence-UI improvements and - kernel-estimation groundwork (PR #94) -- CaTune: log-scale DualRangeSlider, card-grid fix, tutorial baseline docs (PR #96) -- Performance: FISTA pipeline (SIMD, loop fusion, Fenwick baseline) (PR #92); - CaDecon iteration hot paths (PR #107); solver cleanup/dedup/optimize (PR #106); - snappier Peak/FWHM slider drag (PR #134) - `calab-solver` tuning-constant hygiene: introduced `SeedConfig`, shared `baseline::DEFAULT_BASELINE_QUANTILE`, and a named `BASELINE_EMA_WEIGHT`; deduplicated the bi-exponential fast-component grid bounds so the grid search and golden-section refinement cannot drift (no behavior change) (PR #163) +- Tooling: ignore local Python virtualenvs `.venv*/` (PR #156) +- Documentation: reconciled repo docs with the CaDecon review series (PR #164), + aligned the CaDecon tutorials with it (PR #165), and backfilled the changelog + from git history (PR #167) + +### Fixed + +- `calab-solver` FFI boundaries (WASM and PyO3) reject non-finite (NaN/Inf) + input traces with an explicit error instead of returning garbage results + (PR #161) +- Solver: banded AR(2) forward model aligned via a one-sample source delay so + the reconvolution matches the double-exponential kernel (PR #157) +- CaDecon: correct per-subset kernel attribution + init/variance robustness + (PR #153) + +## [2.4.0] - 2026-03-20 + +> Covers the entire 2.4.x line (PRs #99–#152). Reconstructed from git history; +> closely-related PRs are consolidated into single bullets for readability. + +### Added + +- **`calab` Python package** — CaDecon Python bridge with config, autorun, + progress, and auto-export (PRs #108, #109); headless-browser batch mode + + InDeCa PyO3 bindings (PR #110) +- Shared Rust **simulation module** producing synthetic ground-truth traces, + exposed to both Python and WASM (PR #113) +- Solver: peak-seeded initial-kernel auto-estimation (PR #103); an independent + fast component in the bi-exponential fit (PR #105); a `skip` parameter for + bi-exponential fitting (PR #99) +- Migrated the kernel parameterization from (tau_rise, tau_decay) to + (t_peak, FWHM) (PR #104) +- CaDecon tutorial set (PR #151) +- Draggable minimap edges on the trace overview (PR #145) +- Sphinx + ReadTheDocs documentation site for the Python package (PR #115) + +### Changed + +- Performance: CaDecon iteration hot paths (PR #107); solver + cleanup/dedup/optimize (PR #106); snappier Peak/FWHM slider drag (PR #134) - Tooling: ESLint/Prettier/lint-surface cleanup (PR #120); prune unused exports and internalize test-only surface (PR #125); bump GHA for Node 24 and clear - reactivity lint (PR #133); gitignore the whole `.claude/` directory (PR #135); - ignore local Python virtualenvs `.venv*/` (PR #156) + reactivity lint (PR #133); gitignore the whole `.claude/` directory (PR #135) - CI: Rust + Python lint/type jobs, a build matrix, and SHA-pinned actions (PR #124) - Tests: smoke / export-roundtrip / sub-frame-timing / warm-start quick-wins @@ -88,29 +97,16 @@ Versions correspond to git tags (`v*`) and apply to the entire monorepo. RLS policy matrix (PR #130); bridge timeout & mid-run crash detection (PR #131) - Documentation: separated CaTune and CaDecon into dedicated guides (PR #117); promoted CaDecon to stable + root README update (PR #118); reviewed/improved - all Python docs (PR #116); improved tutorial terminology and scientific - accuracy (PR #59); reconciled repo docs with the review series (PR #164) and - aligned the CaDecon tutorials with it (PR #165) + all Python docs (PR #116) ### Fixed -- Codebase-wide quality sweep — 26 fixes (PR #64) - Address pre-merge audit findings — WASM drift, RLS PII, FFI panics, config, tests (PR #150) -- `calab-solver` FFI boundaries (WASM and PyO3) reject non-finite (NaN/Inf) - input traces with an explicit error instead of returning garbage results - (PR #161) - Solver: corrected a binning-induced time offset in iterative kernel fitting - (PR #102); golden-section refinement bug fix (PR #147); banded AR(2) forward - model aligned via a one-sample source delay so the reconvolution matches the - double-exponential kernel (PR #157) -- CaDecon: correct per-subset kernel attribution + init/variance robustness - (PR #153) -- CaTune: minimap no longer pushes the zoom window off-screen (PRs #70, #82); - clamp rise/decay sliders to prevent a negative kernel (PR #83); GT marker - alignment + spectrum/zoom-window perf sweep (PR #142); repair tutorial - highlighting after the Peak/FWHM migration (PR #143) -- Analytics: reliable session-duration tracking via heartbeat (PR #84) + (PR #102); golden-section refinement bug fix (PR #147) +- CaTune: GT marker alignment + spectrum/zoom-window perf sweep (PR #142); + repair tutorial highlighting after the Peak/FWHM migration (PR #143) - Headless: prevent resource leaks on browser start/close failures (PR #121) - Logic + UX polish — tau constraints, bridge errors, reactivity (PR #126) - Community: show bridge/training submissions and hide demo presets under @@ -123,6 +119,81 @@ Versions correspond to git tags (`v*`) and apply to the entire monorepo. geo-session edge function (PR #122) - Locked down analytics row-level security (PR #123) +## [2.3.0] - 2026-02-26 + +> Covers the entire 2.3.x line (PRs #85–#96). Reconstructed from git history. + +### Added + +- **CaDecon** — a new app for automated calcium deconvolution (the InDeCa + algorithm) that estimates the kernel and deconvolution parameters directly + from the data, no manual tuning required: app scaffold + data loading + + subset UI, the InDeCa compute engine with warm-start, visualization / QC + distributions / drill-down, community-database integration, and ground-truth + overlay (PRs #85, #86, #87, #88, #90, #91) +- Usage analytics extended to track CaDecon submissions (PR #93) + +### Changed + +- CaDecon left-sidebar layout/UX (PR #89); convergence-UI improvements and + kernel-estimation groundwork, including rise-time-collapse mitigation (PR #94) +- CaTune: log-scale DualRangeSlider, card-grid fix, tutorial baseline docs (PR #96) +- Performance: FISTA pipeline (SIMD, loop fusion, Fenwick baseline) (PR #92) + +### Fixed + +- Solver: alpha/PVE double-counting and energy-pooling correctness (PR #91) + +## [2.2.0] - 2026-02-23 + +> Covers the entire 2.2.x line (PRs #65–#84). Reconstructed from git history. + +### Added + +- **`calab` Python package** greatly expanded — PyO3 bindings, CaImAn/Minian + loaders, browser bridge, and a CLI (PR #66) +- Community: DataSource tracking + bridge export button & heartbeat detection + (PR #67) +- Solver: banded AR(2) O(T) convolution + box constraint (PR #78) +- Dynamic worker-pool scaling with a URL override (PR #77) +- Admin dashboard: analytics breakdowns and bulk moderation (PR #69) +- Chart/UX: transient-zone visual indicator (PR #81) +- Tutorials: Python Package tutorial (PR #68); Python syntax highlighting in + code blocks (PR #75) + +### Changed + +- Moved the Rust solver to `crates/solver/` with dual WASM (`jsbindings`) / + PyO3 (`pybindings`) Cargo features (PR #65) +- Replaced the export-to-Python page with a dismissible modal (PR #79) + +### Fixed + +- CaTune: minimap no longer pushes the zoom window off-screen (PRs #70, #82); + clamp rise/decay sliders to prevent a negative kernel (PR #83) +- Analytics: reliable session-duration tracking via heartbeat (PR #84) + +## [2.1.0] - 2026-02-20 + +> Covers the 2.0.8, 2.0.9, and 2.1.x patch line (PRs #58–#64). Reconstructed +> from git history. + +### Added + +- **Usage-analytics pipeline + admin dashboard** (PR #62) +- Shared **auth menu** in the header across all CaLab apps (PR #61) +- Community: highlight your own submissions in the scatter plot (PR #63) +- Comprehensive README files across all packages and apps (PR #58) + +### Changed + +- Made `@calab/community` app-agnostic (PR #60) +- Documentation: improved tutorial terminology and scientific accuracy (PR #59) + +### Fixed + +- Codebase-wide quality sweep — 26 fixes (PR #64) + ## [2.0.6] - 2026-02-19 ### Changed @@ -214,7 +285,12 @@ Major restructuring into a monorepo with reusable packages. - Renamed repo references from CaTune to CaLab - Stabilized tooling and codified conventions (Prettier, ESLint, CI) (PR #41) -[Unreleased]: https://github.com/miniscope/CaLab/compare/v2.0.6...HEAD +[2.6.0]: https://github.com/miniscope/CaLab/compare/v2.5.0...HEAD +[2.5.0]: https://github.com/miniscope/CaLab/compare/v2.4.10...v2.5.0 +[2.4.0]: https://github.com/miniscope/CaLab/compare/v2.3.8...v2.4.10 +[2.3.0]: https://github.com/miniscope/CaLab/compare/v2.2.7...v2.3.8 +[2.2.0]: https://github.com/miniscope/CaLab/compare/v2.1.2...v2.2.7 +[2.1.0]: https://github.com/miniscope/CaLab/compare/v2.0.6...v2.1.2 [2.0.6]: https://github.com/miniscope/CaLab/compare/v2.0.5...v2.0.6 [2.0.5]: https://github.com/miniscope/CaLab/compare/v2.0.4...v2.0.5 [2.0.4]: https://github.com/miniscope/CaLab/compare/v2.0.3...v2.0.4 diff --git a/python/docs/guides/cadecon.md b/python/docs/guides/cadecon.md index eeb094d0..6a80cc46 100644 --- a/python/docs/guides/cadecon.md +++ b/python/docs/guides/cadecon.md @@ -205,6 +205,7 @@ calab.solve_trace( lp_enabled: bool = False, warm_counts: np.ndarray | None = None, lambda_: float = 0.0, + noise_constrained: bool = False, ) -> SolveTraceResult ``` @@ -221,6 +222,7 @@ calab.solve_trace( | `lp_enabled` | Enable low-pass filtering. | | `warm_counts` | Spike counts from a previous iteration for warm-start. | | `lambda_` | L1 sparsity penalty (0 = auto). | +| `noise_constrained` | Pick the binarization threshold as the sparsest spike support whose residual still reaches the data-derived noise floor, instead of the fit-maximizing threshold. Knob-free; suppresses spurious low-SNR spikes. Default `False`. | Returns a `SolveTraceResult` namedtuple with fields: `s_counts`, `alpha`, `baseline`, `threshold`, `pve`, `iterations`, `converged`. diff --git a/python/src/calab/_compute.py b/python/src/calab/_compute.py index 0a3f7145..713b5571 100644 --- a/python/src/calab/_compute.py +++ b/python/src/calab/_compute.py @@ -355,6 +355,7 @@ def solve_trace( lp_enabled: bool = False, warm_counts: np.ndarray | None = None, lambda_: float = 0.0, + noise_constrained: bool = False, ) -> SolveTraceResult: """Run the InDeCa pipeline on a single trace. Delegates to Rust. @@ -378,6 +379,12 @@ def solve_trace( Spike counts from a previous iteration (at original rate) for warm-start. lambda_ : float L1 sparsity penalty. + noise_constrained : bool + Choose the binarization threshold as the sparsest spike support whose + residual still reaches the data-derived noise floor, instead of the one + that maximizes fit. Suppresses noise fit as spurious spikes; the effect + concentrates at low SNR. Knob-free (the noise floor is measured from the + trace). Default False. Returns ------- @@ -392,6 +399,7 @@ def solve_trace( trace_1d, tau_rise, tau_decay, fs, upsample_factor, max_iters, tol, hp_enabled, lp_enabled, warm, lambda_, + noise_constrained, ) return SolveTraceResult( s_counts=np.asarray(s_counts), diff --git a/python/tests/test_solve_trace.py b/python/tests/test_solve_trace.py index f82b3e75..ee9d8cba 100644 --- a/python/tests/test_solve_trace.py +++ b/python/tests/test_solve_trace.py @@ -86,6 +86,15 @@ def test_tuple_unpacking(self): assert isinstance(alpha, float) assert isinstance(converged, bool) + def test_noise_constrained_accepted(self): + # The noise_constrained knob is exposed through the binding and produces + # a valid result without changing output shape. + trace = _make_trace(0.02, 0.4, 30.0, 300, [30, 100, 200], alpha=10.0, baseline=2.0) + result = solve_trace(trace, 0.02, 0.4, 30.0, noise_constrained=True) + assert isinstance(result, SolveTraceResult) + assert result.s_counts.shape == (300,) + assert result.s_counts.sum() >= 0 + # --------------------------------------------------------------------------- # estimate_kernel From 55ff539e6d4f43946a6beb45da2c65515dea12b4 Mon Sep 17 00:00:00 2001 From: daharoni Date: Wed, 8 Jul 2026 22:33:30 -0700 Subject: [PATCH 3/4] fix(cadecon): mode-aware scale-loop selection under noise-constrained; review cleanups MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Council review caught a correctness bug: the outer scale-iteration loop always selected the best iterate by max PVE, even under noise-constrained selection — re-selecting the densest-fitting iteration after the inner threshold search had deliberately stopped at the noise floor, partially re-laundering the sparsity. Under NoiseFloor the outer loop now selects the best-calibrated prescale (min scale-convergence error), the loop's own fixed point; MaxPve is unchanged (historical output byte-identical). Also folds in the prior review cleanups: - unify grid_sse_at_threshold into evaluate_threshold via a stride param - top-down early-exit noise-floor scan (identical result, fewer evaluations) - Rust tests: budget-ordering invariant, infeasible fallback, high_band_sigma recovery, end-to-end noisy-trace recovery; Python noisy-trace test - doc rewording (drop "experiment"/"count-debiasing"); corrected viz-store overlay comment; consolidated TraceInspector label helpers - UI: clarify the toggle describes a stopping rule + SNR≈1 trim caveat Co-Authored-By: Claude Opus 4.8 --- .../components/controls/AlgorithmSettings.tsx | 2 +- .../src/components/traces/TraceInspector.tsx | 6 +- apps/cadecon/src/lib/viz-store.ts | 7 +- crates/solver/src/indeca.rs | 128 +++++++++- crates/solver/src/threshold.rs | 228 ++++++++++++++---- python/tests/test_solve_trace.py | 16 ++ 6 files changed, 321 insertions(+), 66 deletions(-) diff --git a/apps/cadecon/src/components/controls/AlgorithmSettings.tsx b/apps/cadecon/src/components/controls/AlgorithmSettings.tsx index 7c11bac8..b82a4e11 100644 --- a/apps/cadecon/src/components/controls/AlgorithmSettings.tsx +++ b/apps/cadecon/src/components/controls/AlgorithmSettings.tsx @@ -92,7 +92,7 @@ export function AlgorithmSettings(): JSX.Element { comparisonDeconv() != null); // Labels for the two deconv traces, by their actual setting. - const deconvLabel = () => (noiseConstrained() ? 'Deconv (noise-constr.)' : 'Deconv (standard)'); - const compareLabel = () => (noiseConstrained() ? 'Deconv (standard)' : 'Deconv (noise-constr.)'); + const settingShort = (nc: boolean) => (nc ? 'noise-constr.' : 'standard'); + const deconvLabel = () => `Deconv (${settingShort(noiseConstrained())})`; + const compareLabel = () => `Deconv (${settingShort(!noiseConstrained())})`; // Zoom window state const totalDuration = createMemo(() => { @@ -551,7 +552,6 @@ export function TraceInspector(): JSX.Element { const c = comparisonDeconv(); return c ? c.reduce((s, v) => s + v, 0).toFixed(0) : null; }; - const settingShort = (nc: boolean) => (nc ? 'noise-constr.' : 'standard'); // Subset highlight zones for the minimap — show which time regions // the algorithm operates on for the currently selected cell. diff --git a/apps/cadecon/src/lib/viz-store.ts b/apps/cadecon/src/lib/viz-store.ts index df2d0495..06014a0a 100644 --- a/apps/cadecon/src/lib/viz-store.ts +++ b/apps/cadecon/src/lib/viz-store.ts @@ -19,9 +19,10 @@ const [showResidual, setShowResidual] = createSignal(false); const [showGTCalcium, setShowGTCalcium] = createSignal(true); const [showGTSpikes, setShowGTSpikes] = createSignal(true); -// Sparsity comparison overlay: when on, TraceInspector re-runs inference for the -// inspected cell with the OPPOSITE noise-constrained setting (same converged -// kernel) and overlays it, so the impact of noise-constrained sparsity is visible. +// Sparsity comparison overlay (display toggle): when on, TraceInspector overlays +// the opposite-noise-constrained deconvolution that was precomputed during the +// run (requires `sparsityCompareEnabled`), so the impact of noise-constrained +// sparsity is visible. Read-only on the browse path — no live solve here. const [showSparsityCompare, setShowSparsityCompare] = createSignal(true); export { diff --git a/crates/solver/src/indeca.rs b/crates/solver/src/indeca.rs index 70adb5c9..be9462eb 100644 --- a/crates/solver/src/indeca.rs +++ b/crates/solver/src/indeca.rs @@ -20,7 +20,9 @@ use crate::upsample::{ use crate::{Constraint, ConvMode, Solver}; use realfft::RealFftPlanner; -/// Optional spike-inference behaviors, off by default (current shipping path). +/// Optional spike-inference behaviors. The library default is off (`MaxPve`, +/// preserving the historical output); the CaDecon app enables `noise_constrained` +/// by default. /// /// `noise_constrained` chooses the binarization threshold at the data-derived /// noise floor instead of maximizing fit, suppressing low-SNR spurious spikes @@ -344,8 +346,8 @@ pub fn solve_trace( ) } -/// See [`solve_trace`]; adds optional [`SolveOptions`] for the count-debiasing -/// experiment (noise-constrained threshold selection). +/// See [`solve_trace`]; adds optional [`SolveOptions`] for noise-constrained +/// threshold selection. #[allow(clippy::too_many_arguments)] pub fn solve_trace_opts( trace: &[f32], @@ -438,6 +440,7 @@ pub fn solve_trace_opts( const SCALE_RTOL: f64 = 0.05; let mut best_pve = f64::NEG_INFINITY; + let mut best_scale_err = f64::INFINITY; let mut best_result: Option<(Vec, f64, f64, f64, f64, u32, bool)> = None; // Pre-allocate scratch buffers reused across scale iterations. @@ -509,10 +512,30 @@ pub fn solve_trace_opts( selection, ); - // Track the best result by PVE. - // alpha_lstsq is already the true alpha (fit against original trace). - if pve > best_pve { + // Scale-loop convergence error: how close the lstsq-recovered alpha is to + // the prescale used this round. This is the loop's own objective. + let scale_err = if alpha_est > 1e-10 { + (alpha_lstsq / alpha_est - 1.0).abs() + } else { + f64::INFINITY + }; + + // Select the best iterate across scale rounds. alpha_lstsq is already the + // true alpha (fit against the original trace). + // + // MaxPve keeps the highest-PVE iterate (historical behavior). Under + // NoiseFloor, ranking by PVE would defeat the criterion — the inner search + // deliberately stops at the noise floor (below max PVE), so a max-PVE outer + // pick would re-select the densest-fitting iteration and re-launder the + // sparsity. Instead select the best-calibrated prescale (smallest scale + // error) — the scale loop's own fixed point, which is criterion-neutral. + let is_better = match selection { + Selection::MaxPve => pve > best_pve, + Selection::NoiseFloor { .. } => scale_err < best_scale_err, + }; + if is_better { best_pve = pve; + best_scale_err = scale_err; best_result = Some(( s_binary, alpha_lstsq, @@ -525,7 +548,7 @@ pub fn solve_trace_opts( } // Converged: alpha_lstsq ≈ alpha_est means the prescale was correct. - if alpha_est > 1e-10 && (alpha_lstsq / alpha_est - 1.0).abs() < SCALE_RTOL { + if scale_err < SCALE_RTOL { break; } @@ -949,4 +972,95 @@ mod tests { result.pve ); } + + /// Deterministic white-ish noise in [-amp, amp) (variance ≈ amp²/3). + fn lcg_noise(n: usize, amp: f32, seed: u64) -> Vec { + let mut state = seed; + (0..n) + .map(|_| { + state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + let u = ((state >> 32) as f64) / ((1u64 << 31) as f64) - 1.0; + (u as f32) * amp + }) + .collect() + } + + #[test] + fn high_band_sigma_recovers_white_noise_std() { + // On pure white noise the high-band periodogram mean estimates the noise + // variance, so its sqrt should recover the injected std. + let n = 4096; + let amp = 0.3_f32; + let noise = lcg_noise(n, amp, 0x1234_5678); + let sigma_true = (amp as f64) / 3.0_f64.sqrt(); + let sigma_est = high_band_sigma(&noise); + let rel = (sigma_est - sigma_true).abs() / sigma_true; + assert!( + rel < 0.15, + "estimated sigma {:.4} should match injected {:.4} (rel err {:.3})", + sigma_est, + sigma_true, + rel + ); + } + + #[test] + fn noise_constrained_recovers_events_on_noisy_trace() { + // Exercise the noise-floor selection path end-to-end (solve_trace_opts → + // estimate_grid_noise_sigma → Selection::NoiseFloor) on a genuinely noisy + // trace. It should recover the real events with a reasonable fit and not + // wildly overcount. The per-trace count is NOT guaranteed to be below the + // max-PVE default — the two criteria use different search strategies — so + // the sparsity ordering is asserted deterministically in the threshold + // unit test `noise_floor_larger_budget_is_sparser` instead. + let spike_positions = [30usize, 100, 200, 260]; + let alpha_true = 6.0_f32; + let baseline_true = 2.0_f32; + let kernel = build_kernel(0.02, 0.4, 30.0); + let n = 300; + let mut trace = vec![baseline_true; n]; + for &pos in &spike_positions { + for (k, &kv) in kernel.iter().enumerate() { + if pos + k < n { + trace[pos + k] += alpha_true * kv; + } + } + } + let noise = lcg_noise(n, 0.4, 0x0BADC0DE); + for (t, &e) in trace.iter_mut().zip(&noise) { + *t += e; + } + + let constrained = solve_trace_opts( + &trace, + 0.02, + 0.4, + 30.0, + 1, + 500, + 1e-4, + None, + false, + false, + 0.0, + SolveOptions { + noise_constrained: true, + }, + ); + + assert_eq!(constrained.s_counts.len(), n); + let count: f32 = constrained.s_counts.iter().sum(); + assert!( + (1.0..=(spike_positions.len() as f32 * 2.0)).contains(&count), + "should recover the events without gross overcounting, got {}", + count + ); + assert!( + constrained.pve > 0.5, + "fit should be reasonable, pve {}", + constrained.pve + ); + } } diff --git a/crates/solver/src/threshold.rs b/crates/solver/src/threshold.rs index 2433937d..b782eec5 100644 --- a/crates/solver/src/threshold.rs +++ b/crates/solver/src/threshold.rs @@ -165,6 +165,7 @@ pub fn threshold_search_opts( banded, thresh, pad, + 1, max_alpha, &mut s_bin, &mut conv_buf, @@ -204,6 +205,7 @@ pub fn threshold_search_opts( banded, thresh, pad, + 1, max_alpha, &mut s_bin, &mut conv_buf, @@ -269,13 +271,21 @@ fn binarize(s: &[f32], threshold: f64, s_bin: &mut [f32]) { } } -/// Evaluate a single threshold: binarize → convolve → lstsq → error. +/// Evaluate a single threshold: binarize → convolve → lstsq → residual SSE. +/// +/// Alpha/baseline are always fit over the full interior; the SSE is accumulated +/// over the interior sampling every `stride`-th position (`stride >= 1`). +/// `stride == 1` sums the whole interior — the max-PVE search. A larger stride +/// restricts the residual to original-rate grid positions (`stride == upsample_ +/// factor`), where the upsampled trace equals the real samples and the noise is +/// uncorrelated — the quantity the noise-floor budget is calibrated against. fn evaluate_threshold( s_relaxed: &[f32], y: &[f32], banded: &BandedAR2, threshold: f64, pad: usize, + stride: usize, max_alpha: f64, s_bin: &mut [f32], conv_buf: &mut [f32], @@ -285,13 +295,14 @@ fn evaluate_threshold( let (alpha, baseline) = lstsq_alpha_baseline(conv_buf, y, pad, max_alpha); - // Error over the interior (excluding boundary padding) let n = y.len(); let mut err = 0.0_f64; - for i in pad..n.saturating_sub(pad) { + let mut i = pad; + while i < n.saturating_sub(pad) { let pred = alpha * conv_buf[i] as f64 + baseline; let d = y[i] as f64 - pred; err += d * d; + i += stride.max(1); } err } @@ -318,78 +329,53 @@ fn select_noise_floor_threshold( conv_buf: &mut [f32], ) -> f64 { let n = y.len(); + // Evaluate the noise constraint on original-rate grid positions only + // (stride = upsample_factor): there the upsampled trace equals the real + // samples, avoiding the correlated, reduced-variance noise at interpolated + // positions, so the budget can use a noise std estimated at the original rate. let stride = upsample_factor.max(1); let n_grid = (pad..n.saturating_sub(pad)).step_by(stride).count(); let budget = sigma * sigma * n_grid as f64; - // Up to 256 candidate thresholds, evenly spaced through the sorted values. + // Up to `cap` candidate thresholds, evenly spaced through the sorted values. + // `cap` is a fixed search resolution (mirroring the coarse/fine counts in the + // max-PVE path), not a result-affecting tuning knob. let cap = 256usize; let step = if vals.len() > cap { vals.len() as f64 / cap as f64 } else { 1.0 }; - - let mut best_feasible = f64::NEG_INFINITY; // highest feasible threshold - let mut best_effort_thr = min_threshold; // min grid-SSE fallback - let mut best_effort_sse = f64::INFINITY; - + let mut candidates: Vec = Vec::with_capacity(cap.min(vals.len())); let mut idx = 0.0_f64; while (idx as usize) < vals.len() { let thr = vals[idx as usize] as f64; idx += step; - if thr < min_threshold { - continue; + if thr >= min_threshold { + candidates.push(thr); } - let sse = grid_sse_at_threshold( + } + + // Grid residual rises as the threshold rises (fewer spikes → worse fit), so + // the sparsest feasible support is the *highest* threshold within budget. + // Scan top-down and return the first feasible one. If none is feasible + // (budget tighter than the best achievable fit), fall back to the + // min-residual threshold found over the same scan. + let mut best_effort_thr = min_threshold; + let mut best_effort_sse = f64::INFINITY; + for &thr in candidates.iter().rev() { + let sse = evaluate_threshold( s_relaxed, y, banded, thr, pad, stride, max_alpha, s_bin, conv_buf, ); - if sse <= budget && thr > best_feasible { - best_feasible = thr; + if sse <= budget { + return thr; } if sse < best_effort_sse { best_effort_sse = sse; best_effort_thr = thr; } } - - if best_feasible.is_finite() { - best_feasible - } else { - best_effort_thr - } -} - -/// Residual SSE at original-rate grid positions only (stride = upsample_factor), -/// with alpha/baseline fit over the full interior. Evaluating the noise -/// constraint on the original grid — where the upsampled trace equals the real -/// samples — avoids the correlated, reduced-variance noise at interpolated -/// positions, so the budget can use a noise std estimated at the original rate. -#[allow(clippy::too_many_arguments)] -fn grid_sse_at_threshold( - s_relaxed: &[f32], - y: &[f32], - banded: &BandedAR2, - threshold: f64, - pad: usize, - stride: usize, - max_alpha: f64, - s_bin: &mut [f32], - conv_buf: &mut [f32], -) -> f64 { - binarize(s_relaxed, threshold, s_bin); - banded.convolve_forward(s_bin, conv_buf); - let (alpha, baseline) = lstsq_alpha_baseline(conv_buf, y, pad, max_alpha); - let n = y.len(); - let mut sse = 0.0_f64; - let mut i = pad; - while i < n.saturating_sub(pad) { - let pred = alpha * conv_buf[i] as f64 + baseline; - let d = y[i] as f64 - pred; - sse += d * d; - i += stride; - } - sse + best_effort_thr } /// Least-squares fit for alpha and baseline: y ≈ alpha * conv + baseline. @@ -582,4 +568,142 @@ mod tests { assert_eq!(boundary_padding(0.2, 100.0), 40); assert_eq!(boundary_padding(1.0, 10.0), 20); } + + /// Deterministic white-ish noise in [-amp, amp) via an LCG (tests must not + /// use real RNG). Zero-mean in expectation; variance ≈ amp²/3. + fn lcg_noise(n: usize, amp: f32, seed: u64) -> Vec { + let mut state = seed; + (0..n) + .map(|_| { + state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + let u = ((state >> 32) as f64) / ((1u64 << 31) as f64) - 1.0; + (u as f32) * amp + }) + .collect() + } + + /// Build a graded relaxed solution + noisy observed trace with three true + /// spikes whose relaxed values differ (1.0 / 0.85 / 0.7) plus weaker spurious + /// bumps (0.55), so the budget genuinely controls how many survive. + fn graded_case(banded: &BandedAR2, n: usize) -> (Vec, Vec, f64) { + let true_pos = [80usize, 250, 430]; + let true_val = [1.0_f32, 0.85, 0.7]; + let mut s_bin_true = vec![0.0_f32; n]; + let mut s_relaxed = vec![0.0_f32; n]; + for (&p, &v) in true_pos.iter().zip(&true_val) { + s_bin_true[p] = 1.0; + s_relaxed[p] = v; + } + for &p in &[40usize, 150, 320, 500, 560] { + s_relaxed[p] = 0.55; + } + let mut conv = vec![0.0_f32; n]; + banded.convolve_forward(&s_bin_true, &mut conv); + let alpha = 8.0_f64; + let baseline = 1.0_f64; + let amp = 0.25_f32; + let noise = lcg_noise(n, amp, 0xABCDEF01); + let y: Vec = conv + .iter() + .zip(&noise) + .map(|(&c, &e)| (alpha * c as f64 + baseline) as f32 + e) + .collect(); + let noise_std = (amp as f64) / 3.0_f64.sqrt(); + (s_relaxed, y, noise_std) + } + + #[test] + fn noise_floor_larger_budget_is_sparser() { + // A larger noise budget admits sparser (higher) thresholds, so it should + // never yield more spikes than a tight budget. With a huge budget the + // sparsest support (highest candidate) is immediately feasible; with a + // budget just above the noise floor, all three true spikes are required. + let banded = BandedAR2::new(0.02, 0.4, 30.0); + let n = 600; + let (s_relaxed, y, noise_std) = graded_case(&banded, n); + + let tight = threshold_search_opts( + &s_relaxed, + &y, + &banded, + 0.4, + 30.0, + 1, + f64::INFINITY, + Selection::NoiseFloor { + sigma: 1.5 * noise_std, + }, + ); + let loose = threshold_search_opts( + &s_relaxed, + &y, + &banded, + 0.4, + 30.0, + 1, + f64::INFINITY, + Selection::NoiseFloor { sigma: 100.0 }, + ); + + let tight_count: f32 = tight.s_binary.iter().sum(); + let loose_count: f32 = loose.s_binary.iter().sum(); + assert!( + loose.threshold >= tight.threshold, + "looser budget should pick a higher threshold ({} vs {})", + loose.threshold, + tight.threshold + ); + assert!( + loose_count <= tight_count, + "looser budget should not add spikes ({} vs {})", + loose_count, + tight_count + ); + assert!( + (tight_count - 3.0).abs() < 0.5, + "tight budget should keep all three true spikes, got {}", + tight_count + ); + assert!( + loose_count < tight_count, + "huge budget should drop the weaker spikes ({} vs {})", + loose_count, + tight_count + ); + } + + #[test] + fn noise_floor_falls_back_when_infeasible() { + // sigma → 0 makes the budget unreachable, so no threshold is feasible and + // the scan must fall back to the minimum-residual threshold — which on + // this signal recovers the true spikes rather than returning garbage. + let banded = BandedAR2::new(0.02, 0.4, 30.0); + let n = 600; + let (s_relaxed, y, _) = graded_case(&banded, n); + + let result = threshold_search_opts( + &s_relaxed, + &y, + &banded, + 0.4, + 30.0, + 1, + f64::INFINITY, + Selection::NoiseFloor { sigma: 1e-9 }, + ); + assert!(result.threshold.is_finite()); + let count: f32 = result.s_binary.iter().sum(); + assert!( + (count - 3.0).abs() < 0.5, + "fallback should recover the three true spikes, got {}", + count + ); + assert!( + result.pve > 0.9, + "fallback fit should explain the signal, pve {}", + result.pve + ); + } } diff --git a/python/tests/test_solve_trace.py b/python/tests/test_solve_trace.py index ee9d8cba..369ddeb0 100644 --- a/python/tests/test_solve_trace.py +++ b/python/tests/test_solve_trace.py @@ -95,6 +95,22 @@ def test_noise_constrained_accepted(self): assert result.s_counts.shape == (300,) assert result.s_counts.sum() >= 0 + def test_noise_constrained_detects_events_with_noise(self): + # Exercise the noise-floor selection path end-to-end on a genuinely noisy + # trace (a noiseless trace only hits the fallback branch). The option + # should still recover the real events and produce a valid fit. Note the + # per-trace spike count is not guaranteed to be <= the default: the two + # criteria use different search strategies (see the Rust + # `noise_floor_larger_budget_is_sparser` test for the sparsity invariant). + rng = np.random.default_rng(0) + clean = _make_trace(0.02, 0.4, 30.0, 400, [40, 150, 300], alpha=4.0, baseline=2.0) + trace = clean + rng.normal(0.0, 0.4, size=clean.shape) + + result = solve_trace(trace, 0.02, 0.4, 30.0, noise_constrained=True) + assert result.s_counts.shape == (400,) + assert result.s_counts.sum() >= 1 + assert result.pve > 0.5 + # --------------------------------------------------------------------------- # estimate_kernel From 83d8bc7067569b37e825d7c3b007dbfd0cf723a6 Mon Sep 17 00:00:00 2001 From: daharoni Date: Wed, 8 Jul 2026 22:53:22 -0700 Subject: [PATCH 4/4] docs(cadecon): fix prettier formatting in solve_trace guide Co-Authored-By: Claude Opus 4.8 --- python/docs/guides/cadecon.md | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/python/docs/guides/cadecon.md b/python/docs/guides/cadecon.md index 6a80cc46..713d2782 100644 --- a/python/docs/guides/cadecon.md +++ b/python/docs/guides/cadecon.md @@ -209,19 +209,19 @@ calab.solve_trace( ) -> SolveTraceResult ``` -| Parameter | Description | -| ----------------- | ------------------------------------------------------ | -| `trace` | 1-D calcium trace. | -| `tau_rise` | Rise time constant in seconds. | -| `tau_decay` | Decay time constant in seconds. | -| `fs` | Sampling rate in Hz. | -| `upsample_factor` | Upsampling multiplier (1 = no upsampling). | -| `max_iters` | Maximum FISTA iterations. | -| `tol` | Convergence tolerance. | -| `hp_enabled` | Enable high-pass filtering. | -| `lp_enabled` | Enable low-pass filtering. | -| `warm_counts` | Spike counts from a previous iteration for warm-start. | -| `lambda_` | L1 sparsity penalty (0 = auto). | +| Parameter | Description | +| ------------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `trace` | 1-D calcium trace. | +| `tau_rise` | Rise time constant in seconds. | +| `tau_decay` | Decay time constant in seconds. | +| `fs` | Sampling rate in Hz. | +| `upsample_factor` | Upsampling multiplier (1 = no upsampling). | +| `max_iters` | Maximum FISTA iterations. | +| `tol` | Convergence tolerance. | +| `hp_enabled` | Enable high-pass filtering. | +| `lp_enabled` | Enable low-pass filtering. | +| `warm_counts` | Spike counts from a previous iteration for warm-start. | +| `lambda_` | L1 sparsity penalty (0 = auto). | | `noise_constrained` | Pick the binarization threshold as the sparsest spike support whose residual still reaches the data-derived noise floor, instead of the fit-maximizing threshold. Knob-free; suppresses spurious low-SNR spikes. Default `False`. | Returns a `SolveTraceResult` namedtuple with fields: `s_counts`, `alpha`, `baseline`, `threshold`, `pve`, `iterations`, `converged`.