diff --git a/apps/cadecon/src/App.tsx b/apps/cadecon/src/App.tsx
index 9da7e76..542af43 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 672a446..b82a4e1 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 0000000..a06d204
--- /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 6188710..82a7634 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,19 @@ 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 settingShort = (nc: boolean) => (nc ? 'noise-constr.' : 'standard');
+ const deconvLabel = () => `Deconv (${settingShort(noiseConstrained())})`;
+ const compareLabel = () => `Deconv (${settingShort(!noiseConstrained())})`;
+
// Zoom window state
const totalDuration = createMemo(() => {
const raw = fullRawTrace();
@@ -311,7 +328,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 +395,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 +437,7 @@ export function TraceInspector(): JSX.Element {
dsResid,
dsGTCalcium,
dsGTSpikes,
+ dsCompare,
];
});
@@ -424,10 +453,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 +494,7 @@ export function TraceInspector(): JSX.Element {
{
key: 'deconv',
color: TRACE_COLORS.deconv,
- label: 'Deconv',
+ label: deconvLabel(),
visible: showDeconv,
setVisible: setShowDeconv,
},
@@ -470,6 +506,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 +548,10 @@ 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;
+ };
// 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 d14c0f9..bb96c60 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 4f86ba4..4f76261 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 af68d46..1a510eb 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 34e8e63..79e131d 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 88cfc58..09ce4ad 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 4103cb7..06014a0 100644
--- a/apps/cadecon/src/lib/viz-store.ts
+++ b/apps/cadecon/src/lib/viz-store.ts
@@ -19,6 +19,12 @@ const [showResidual, setShowResidual] = createSignal(false);
const [showGTCalcium, setShowGTCalcium] = createSignal(true);
const [showGTSpikes, setShowGTSpikes] = createSignal(true);
+// 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 {
viewedIteration,
setViewedIteration,
@@ -38,4 +44,6 @@ export {
setShowGTCalcium,
showGTSpikes,
setShowGTSpikes,
+ showSparsityCompare,
+ setShowSparsityCompare,
};
diff --git a/apps/cadecon/src/styles/layout.css b/apps/cadecon/src/styles/layout.css
index 7ca4962..c0f9c80 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 4041000..20a3ddd 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 cd55fcb..860b49c 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 +316,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 +329,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 noise-constrained
+/// threshold selection.
+#[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 +412,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.
@@ -265,6 +440,7 @@ pub fn solve_trace(
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.
@@ -325,7 +501,7 @@ pub fn solve_trace(
threshold,
pve,
..
- } = threshold_search(
+ } = threshold_search_opts(
s_norm_slice,
&working_trace,
&banded,
@@ -333,12 +509,33 @@ pub fn solve_trace(
fs_up,
upsample_factor,
f64::INFINITY,
+ 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,
@@ -351,7 +548,7 @@ pub fn solve_trace(
}
// 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;
}
@@ -775,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/js_indeca.rs b/crates/solver/src/js_indeca.rs
index 001bd4f..036bbdf 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,7 @@ pub fn indeca_solve_trace(
hp_enabled,
lp_enabled,
lambda,
+ 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 6e74438..30e00a5 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))]
+#[allow(clippy::too_many_arguments)]
fn py_indeca_solve_trace<'py>(
py: Python<'py>,
trace: PyReadonlyArray1,
@@ -437,6 +438,7 @@ fn py_indeca_solve_trace<'py>(
lp_enabled: bool,
warm_counts: Option>,
lambda_: f64,
+ noise_constrained: bool,
) -> PyResult<(
Bound<'py, PyArray1>, // s_counts
f64, // alpha
@@ -449,7 +451,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 +463,7 @@ fn py_indeca_solve_trace<'py>(
hp_enabled,
lp_enabled,
lambda_,
+ indeca::SolveOptions { noise_constrained },
);
Ok((
diff --git a/crates/solver/src/threshold.rs b/crates/solver/src/threshold.rs
index eda6052..b782eec 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,111 @@ 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,
+ 1,
+ 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,
+ 1,
+ 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);
@@ -213,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],
@@ -229,17 +295,89 @@ 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
}
+/// 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();
+ // 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 `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 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 {
+ candidates.push(thr);
+ }
+ }
+
+ // 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 {
+ return thr;
+ }
+ if sse < best_effort_sse {
+ best_effort_sse = sse;
+ best_effort_thr = thr;
+ }
+ }
+ best_effort_thr
+}
+
/// 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
@@ -430,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/docs/CHANGELOG.md b/docs/CHANGELOG.md
index a498406..761c7ff 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 eeb094d..713d278 100644
--- a/python/docs/guides/cadecon.md
+++ b/python/docs/guides/cadecon.md
@@ -205,22 +205,24 @@ calab.solve_trace(
lp_enabled: bool = False,
warm_counts: np.ndarray | None = None,
lambda_: float = 0.0,
+ noise_constrained: bool = False,
) -> 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`.
diff --git a/python/src/calab/_compute.py b/python/src/calab/_compute.py
index 0a3f714..713b557 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 f82b3e7..369ddeb 100644
--- a/python/tests/test_solve_trace.py
+++ b/python/tests/test_solve_trace.py
@@ -86,6 +86,31 @@ 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
+
+ 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