diff --git a/apps/cadecon/src/components/controls/AlgorithmSettings.tsx b/apps/cadecon/src/components/controls/AlgorithmSettings.tsx
index b82a4e1..8ec6ccd 100644
--- a/apps/cadecon/src/components/controls/AlgorithmSettings.tsx
+++ b/apps/cadecon/src/components/controls/AlgorithmSettings.tsx
@@ -11,6 +11,8 @@ import {
setLpFilterEnabled,
noiseConstrained,
setNoiseConstrained,
+ massCount,
+ setMassCount,
maxIterations,
setMaxIterations,
convergenceTol,
@@ -97,6 +99,14 @@ export function AlgorithmSettings(): JSX.Element {
onChange={setNoiseConstrained}
disabled={isRunLocked()}
/>
+
+
);
diff --git a/apps/cadecon/src/lib/algorithm-store.ts b/apps/cadecon/src/lib/algorithm-store.ts
index bb96c60..1f44072 100644
--- a/apps/cadecon/src/lib/algorithm-store.ts
+++ b/apps/cadecon/src/lib/algorithm-store.ts
@@ -11,6 +11,12 @@ const [lpFilterEnabled, setLpFilterEnabled] = createSignal(true);
// 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);
+// Mass-based count readout: count events by relaxed mass (calibrated to the
+// single-spike mass) and refit alpha, instead of binarize->bin-sum. Undoes the
+// coherent-grid overcount that inflates spike counts and halves alpha at high
+// upsampling. Off by default (preserves historical output); see
+// docs/cadecon-mass-count.md.
+const [massCount, setMassCount] = createSignal(false);
// 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.
@@ -63,6 +69,8 @@ export {
setLpFilterEnabled,
noiseConstrained,
setNoiseConstrained,
+ massCount,
+ setMassCount,
sparsityCompareEnabled,
setSparsityCompareEnabled,
maxIterations,
diff --git a/apps/cadecon/src/lib/bridge-effects.ts b/apps/cadecon/src/lib/bridge-effects.ts
index cc63351..1722101 100644
--- a/apps/cadecon/src/lib/bridge-effects.ts
+++ b/apps/cadecon/src/lib/bridge-effects.ts
@@ -13,6 +13,7 @@ import {
setUpsampleTarget,
setHpFilterEnabled,
setLpFilterEnabled,
+ setMassCount,
setMaxIterations,
setConvergenceTol,
} from './algorithm-store.ts';
@@ -47,6 +48,7 @@ function applyConfig(config: BridgeConfig): void {
if (config.upsample_target != null) setUpsampleTarget(config.upsample_target);
if (config.hp_filter_enabled != null) setHpFilterEnabled(config.hp_filter_enabled);
if (config.lp_filter_enabled != null) setLpFilterEnabled(config.lp_filter_enabled);
+ if (config.mass_count != null) setMassCount(config.mass_count);
if (config.max_iterations != null) setMaxIterations(config.max_iterations);
if (config.convergence_tol != null) setConvergenceTol(config.convergence_tol);
if (config.num_subsets != null) setNumSubsets(config.num_subsets);
diff --git a/apps/cadecon/src/lib/cadecon-pool.ts b/apps/cadecon/src/lib/cadecon-pool.ts
index 4f76261..d615584 100644
--- a/apps/cadecon/src/lib/cadecon-pool.ts
+++ b/apps/cadecon/src/lib/cadecon-pool.ts
@@ -27,6 +27,7 @@ interface TraceJobFields {
lpEnabled: boolean;
lambda: number;
noiseConstrained: boolean;
+ massCount: boolean;
computeComparison: boolean;
warmCounts?: Float32Array;
onComplete(result: TraceResult): void;
@@ -123,6 +124,7 @@ const caDeconRouter: MessageRouter = {
lpEnabled: job.lpEnabled,
lambda: job.lambda,
noiseConstrained: job.noiseConstrained,
+ massCount: job.massCount,
computeComparison: job.computeComparison,
warmCounts: warmCopy,
},
diff --git a/apps/cadecon/src/lib/iteration-manager.ts b/apps/cadecon/src/lib/iteration-manager.ts
index 1a510eb..86e6f70 100644
--- a/apps/cadecon/src/lib/iteration-manager.ts
+++ b/apps/cadecon/src/lib/iteration-manager.ts
@@ -44,6 +44,7 @@ import {
hpFilterEnabled,
lpFilterEnabled,
noiseConstrained,
+ massCount,
sparsityCompareEnabled,
traceFistaMaxIters,
traceFistaTol,
@@ -167,6 +168,7 @@ function dispatchTraceJobs(
lpEnabled: boolean,
lambda: number,
noiseConstrained: boolean,
+ massCount: boolean,
computeComparison: boolean,
prevResults?: Map,
): Promise>> {
@@ -216,6 +218,7 @@ function dispatchTraceJobs(
lpEnabled,
lambda,
noiseConstrained,
+ massCount,
computeComparison,
warmCounts,
onComplete(result: TraceResult) {
@@ -454,6 +457,7 @@ export async function startRun(): Promise {
const lpOn = lpFilterEnabled();
const sparsityLambda = 0.0;
const noiseConstrainedOn = noiseConstrained();
+ const massCountOn = massCount();
const computeComparison = sparsityCompareEnabled();
// Create pool
@@ -592,6 +596,7 @@ export async function startRun(): Promise {
lpOn,
sparsityLambda,
noiseConstrainedOn,
+ massCountOn,
computeComparison,
prevTraceCounts,
);
@@ -931,6 +936,7 @@ export async function startRun(): Promise {
lpEnabled: lpOn,
lambda: sparsityLambda,
noiseConstrained: noiseConstrainedOn,
+ massCount: massCountOn,
computeComparison,
warmCounts,
onComplete(result: TraceResult) {
diff --git a/apps/cadecon/src/workers/cadecon-types.ts b/apps/cadecon/src/workers/cadecon-types.ts
index 20a3ddd..d9103e8 100644
--- a/apps/cadecon/src/workers/cadecon-types.ts
+++ b/apps/cadecon/src/workers/cadecon-types.ts
@@ -60,6 +60,9 @@ export type CaDeconWorkerInbound =
/** Noise-constrained threshold selection: choose the sparsest spike support
* whose residual meets the data-derived noise floor (no tuning knob). */
noiseConstrained: boolean;
+ /** Mass-based count readout: count events by relaxed mass and refit alpha,
+ * undoing the coherent-grid overcount that halves alpha (no tuning knob). */
+ massCount: boolean;
/** Also solve with the OPPOSITE noise-constrained setting and return it as
* comparisonSCounts (for the teaching/impact overlay). */
computeComparison: boolean;
diff --git a/apps/cadecon/src/workers/cadecon-worker.ts b/apps/cadecon/src/workers/cadecon-worker.ts
index 860b49c..3dd5554 100644
--- a/apps/cadecon/src/workers/cadecon-worker.ts
+++ b/apps/cadecon/src/workers/cadecon-worker.ts
@@ -38,6 +38,7 @@ function handleTraceJob(req: Extract,
+ /// Calibrated continuous firing-rate estimate (graded), on the same absolute scale
+ /// as `s_counts` but not rounded. Populated only by the `mass_count` path; empty
+ /// otherwise. See docs/masscount_R_metrics.pdf.
+ pub s_rate: Vec,
pub filtered_trace: Option>,
pub alpha: f64,
pub baseline: f64,
@@ -346,6 +360,179 @@ pub fn solve_trace(
)
}
+/// Mass of a single isolated spike, used to calibrate the
+/// mass-based count. Deconvolves a synthetic peak-normalized single-spike
+/// transient through the identical Box[0,1] FISTA path, normalizes the relaxed
+/// solution to unit interior peak, and integrates it above `mass_floor` (the fixed
+/// low floor also used for the events, so both masses cover the whole bump
+/// and are comparable).
+/// Depends only on (tau, fs_up, mass_floor)
+#[allow(clippy::too_many_arguments)]
+fn single_spike_mass(
+ solver: &mut Solver,
+ banded: &BandedAR2,
+ tau_r: f64,
+ tau_d: f64,
+ fs_up: f64,
+ max_iters: u32,
+ tol: f64,
+ lambda: f64,
+ mass_floor: f64,
+) -> f64 {
+ let n = ((8.0 * tau_d * fs_up).ceil() as usize).max(64);
+ let mut s = vec![0.0_f32; n];
+ s[n / 2] = 1.0;
+ let mut transient = vec![0.0_f32; n];
+ banded.convolve_forward(&s, &mut transient); // peak-normalized (peak = 1.0)
+
+ let (relaxed, _, _, _) = solve_upsampled(
+ solver,
+ &transient,
+ tau_r,
+ tau_d,
+ fs_up,
+ max_iters,
+ tol,
+ None,
+ false,
+ false,
+ Constraint::Box01,
+ true,
+ lambda,
+ );
+
+ let pad = crate::threshold::boundary_padding(tau_d, fs_up).min(n / 4);
+ let peak = interior_peak(&relaxed, pad);
+ if peak <= 1e-10 {
+ return 1.0;
+ }
+ let thr = mass_floor as f32;
+ let mass: f64 = relaxed
+ .iter()
+ .map(|&v| v / peak)
+ .filter(|&v| v >= thr)
+ .map(|v| v as f64)
+ .sum();
+ mass.max(1e-6)
+}
+
+/// Mass-based count readout. Each maximal contiguous run of the relaxed solution
+/// above the low floor `mass_floor` is one event (kept only if its peak clears the
+/// realness gate `theta`); its spike count is `round(run_mass / calibration_mass)`
+/// (at least 1), where `run_mass` is the integral of the relaxed solution over the
+/// run. The `count` spikes are placed within the run (single events at the run peak;
+/// bursts spread across the run span) and alpha/baseline are refit by least squares
+/// against the resulting event train, so alpha becomes the per-spike amplitude rather
+/// than the mass-conserving halved value. `s_rate` is the graded relaxed/calibration_mass
+/// over the kept runs — a continuous firing-rate estimate on the same absolute scale.
+///
+/// Returns `(s_counts, s_rate, alpha, baseline, pve)` at the original rate.
+#[allow(clippy::too_many_arguments)]
+fn mass_count_readout(
+ relaxed: &[f32],
+ theta: f64,
+ mass_floor: f64,
+ working_trace: &[f32],
+ banded: &BandedAR2,
+ tau_d: f64,
+ fs_up: f64,
+ upsample_factor: usize,
+ calibration_mass: f64,
+) -> (Vec, Vec, f64, f64, f64) {
+ let n = relaxed.len();
+ // Runs (events) are extracted down to the low floor mass_floor so the mass covers the
+ // whole bump; theta only gates realness (a run is kept iff its peak clears theta).
+ let flo = mass_floor as f32;
+ let mut s_events = vec![0.0_f32; n];
+ // Graded s_rate: the s_relaxed is scaled by the single-spike
+ // calibration mass. Each run integrates to
+ // mass/calibration_mass ≈ true count, so summed to the frame rate it is a
+ // continuous firing-rate estimate on the correct scale.
+ let mut s_soft = vec![0.0_f32; n];
+ let inv_calibration_mass = (1.0 / calibration_mass) as f32;
+
+ let mut i = 0;
+ while i < n {
+ if relaxed[i] >= flo {
+ let start = i;
+ let mut mass = 0.0_f64;
+ let mut peak = 0.0_f32;
+ let mut arg = i;
+ while i < n && relaxed[i] >= flo {
+ mass += relaxed[i] as f64;
+ if relaxed[i] > peak {
+ peak = relaxed[i];
+ arg = i;
+ }
+ i += 1;
+ }
+ let end = i; // exclusive
+ let span = end - start;
+ // Drop runs whose peak never clears the noise/selection floor.
+ if (peak as f64) < theta {
+ continue;
+ }
+ // Graded rate over the whole run.
+ for t in start..end {
+ s_soft[t] = relaxed[t] * inv_calibration_mass;
+ }
+ let count = (mass / calibration_mass).round().max(1.0) as usize;
+
+ if count <= 1 || span <= 1 {
+ s_events[arg] += 1.0; // single event at the run peak
+ } else {
+ // Spread the events across the run so the reconvolution matches a
+ // burst rather than a single tall spike.
+ let placeable = count.min(span);
+ for j in 0..placeable {
+ let idx = start + (j * span) / placeable + span / (2 * placeable);
+ s_events[idx.min(end - 1)] += 1.0;
+ }
+ // If there are more events than distinct bins, stack the remainder at the peak.
+ if count > placeable {
+ s_events[arg] += (count - placeable) as f32;
+ }
+ }
+ } else {
+ i += 1;
+ }
+ }
+
+ // Refit alpha + baseline against the event train.
+ let pad = crate::threshold::boundary_padding(tau_d, fs_up).min(n / 4);
+ let mut conv = vec![0.0_f32; n];
+ banded.convolve_forward(&s_events, &mut conv);
+ let (alpha, baseline) = lstsq_alpha_baseline(&conv, working_trace, pad, f64::INFINITY);
+
+ // PVE over the interior.
+ let lo = pad;
+ let hi = n.saturating_sub(pad);
+ let mut pve = 0.0;
+ if hi > lo {
+ let len = (hi - lo) as f64;
+ let mut y_sum = 0.0_f64;
+ for i in lo..hi {
+ y_sum += working_trace[i] as f64;
+ }
+ let y_mean = y_sum / len;
+ let mut ss_tot = 0.0_f64;
+ let mut ss_res = 0.0_f64;
+ for i in lo..hi {
+ let yi = working_trace[i] as f64;
+ let d = yi - y_mean;
+ ss_tot += d * d;
+ let pred = alpha * conv[i] as f64 + baseline;
+ let r = yi - pred;
+ ss_res += r * r;
+ }
+ pve = if ss_tot > 1e-20 { 1.0 - ss_res / ss_tot } else { 0.0 };
+ }
+
+ let s_counts = downsample_binary(&s_events, upsample_factor);
+ let s_rate = downsample_binary(&s_soft, upsample_factor); // sum graded values per frame
+ (s_counts, s_rate, alpha, baseline, pve)
+}
+
/// See [`solve_trace`]; adds optional [`SolveOptions`] for noise-constrained
/// threshold selection.
#[allow(clippy::too_many_arguments)]
@@ -442,6 +629,9 @@ pub fn solve_trace_opts(
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;
+ // Relaxed (normalized) solution of the selected iterate — only needed for
+ // the mass-count readout, so captured only when that mode is on.
+ let mut best_relaxed: Vec = Vec::new();
// Pre-allocate scratch buffers reused across scale iterations.
let wt_len = working_trace.len();
@@ -534,6 +724,9 @@ pub fn solve_trace_opts(
Selection::NoiseFloor { .. } => scale_err < best_scale_err,
};
if is_better {
+ if opts.mass_count {
+ best_relaxed = s_norm_slice.to_vec();
+ }
best_pve = pve;
best_scale_err = scale_err;
best_result = Some((
@@ -560,14 +753,46 @@ pub fn solve_trace_opts(
}
// ── Step 4: Extract best result ─────────────────────────────────────
- let (s_binary, alpha, baseline, threshold, pve, iterations, converged) = best_result
+ let (s_binary, mut alpha, mut baseline, threshold, mut pve, iterations, converged) = best_result
.unwrap_or_else(|| {
// Fallback: no valid result found (shouldn't happen)
(vec![0.0; wt_len], 0.0, 0.0, 0.0, 0.0, 0, false)
});
- // Downsample binary spike train to original rate using centered bins
- let s_counts = downsample_binary(&s_binary, upsample_factor);
+ // ── Step 4b (optional): mass-based count readout ────────────────────
+ // Replace the bin-summed binary count with a mass-based event count that
+ // undoes the coherent-grid overcount (see [`SolveOptions::mass_count`]).
+ // `s_rate` (mass_count only): calibrated continuous firing-rate estimate — graded and
+ // on the correct absolute scale, complementary to the integer `s_counts` (empty otherwise).
+ let (s_counts, s_rate) = if opts.mass_count && !best_relaxed.is_empty() && threshold > 1e-9 {
+ // Mass is integrated over the WHOLE bump down to a fixed low floor mass_floor,
+ // not the (possibly near-peak) selection threshold — integrating at the tip
+ // is ill-conditioned (see docs/cadecon-mass-count.md §"calibration floor").
+ // The selection threshold is used only as a realness gate (peak >= threshold).
+ let mass_floor = (0.5 / upsample_factor.max(1) as f64).min(0.15);
+ let calibration_mass = single_spike_mass(
+ &mut solver, &banded, tau_r, tau_d, fs_up, max_iters, tol, lambda, mass_floor,
+ );
+ let (s_counts_mc, s_rate_mc, alpha_mc, baseline_mc, pve_mc) = mass_count_readout(
+ &best_relaxed,
+ threshold, // realness gate (peak >= threshold)
+ mass_floor,
+ &working_trace,
+ &banded,
+ tau_d,
+ fs_up,
+ upsample_factor,
+ calibration_mass,
+ );
+ alpha = alpha_mc;
+ baseline = baseline_mc;
+ pve = pve_mc;
+ (s_counts_mc, s_rate_mc)
+ } else {
+ // Downsample binary spike train to original rate using centered bins.
+ // s_rate is only produced by the mass-count path.
+ (downsample_binary(&s_binary, upsample_factor), Vec::new())
+ };
// Downsample filtered trace to original rate directly from working_trace
// (working_trace is not modified after baseline subtraction).
@@ -575,6 +800,7 @@ pub fn solve_trace_opts(
InDecaResult {
s_counts,
+ s_rate,
filtered_trace,
alpha,
baseline,
@@ -857,6 +1083,143 @@ mod tests {
);
}
+ /// Fast guard for the readout logic (no FISTA): a run of mass ≈ calibration_mass
+ /// counts as one spike, a run of mass ≈ 2·calibration_mass as two, a sub-threshold
+ /// run is dropped, and s_rate is the graded relaxed/calibration_mass over kept runs.
+ #[test]
+ fn mass_count_readout_counts_and_gates() {
+ use crate::banded::BandedAR2;
+ let (tau_d, fs_up, factor) = (0.6, 300.0, 10usize);
+ let banded = BandedAR2::new(0.1, tau_d, fs_up);
+ let n = 300;
+ let calibration_mass = 3.0; // mass of the unit bump below
+ let (theta, mass_floor) = (0.5, 0.05);
+
+ let mut relaxed = vec![0.0_f32; n];
+ for (o, &v) in [0.3, 0.7, 1.0, 0.7, 0.3].iter().enumerate() {
+ relaxed[90 + o] = v; // event A: mass 3.0 = calibration_mass → count 1
+ }
+ for (o, &v) in [0.6, 1.0, 1.0, 1.0, 1.0, 1.0, 0.4].iter().enumerate() {
+ relaxed[130 + o] = v; // burst B: one run, mass 6.0 = 2·calibration_mass → count 2
+ }
+ relaxed[190] = 0.3; // sub-threshold C: peak 0.3 < theta → gated out
+ relaxed[191] = 0.3;
+
+ // Working trace correlated with the events so the alpha refit is well-posed.
+ let mut probe = vec![0.0_f32; n];
+ for &p in &[92usize, 131, 134] {
+ probe[p] = 1.0;
+ }
+ let mut wt = vec![0.0_f32; n];
+ banded.convolve_forward(&probe, &mut wt);
+ for v in wt.iter_mut() {
+ *v = 5.0 * *v + 2.0;
+ }
+
+ let (s_counts, s_rate, alpha, _b, _pve) = mass_count_readout(
+ &relaxed, theta, mass_floor, &wt, &banded, tau_d, fs_up, factor, calibration_mass,
+ );
+
+ // counts: A → 1, B → 2, C gated → total 3
+ assert!((s_counts.iter().sum::() - 3.0).abs() < 1e-6);
+ // s_rate integrates to (3+6)/calibration_mass = 3 over kept runs, graded (fractional)
+ assert!((s_rate.iter().sum::() - 3.0).abs() < 1e-2);
+ assert!(s_rate.iter().any(|&v| (v - v.round()).abs() > 0.05));
+ assert!(alpha.is_finite() && alpha > 0.0);
+ }
+
+ /// mass_count readout restores an unbiased alpha and spike count at the
+ /// default 10x upsample, where the bin-sum readout inflates the count ~k x
+ /// and halves alpha. Uses moderate (~1 Hz) firing where events are mostly
+ /// temporally resolvable, so both count and alpha should recover.
+ #[test]
+ #[ignore = "slow end-to-end regression (2 full + 2 calibration solves); run in CI with `cargo test -- --ignored`"]
+ fn mass_count_restores_alpha_and_count() {
+ let tau_r = 0.1;
+ let tau_d = 0.6;
+ let fs = 30.0;
+ let n = 1500;
+ let factor = 10;
+ let alpha_true = 5.0_f32;
+ let baseline = 2.0_f32;
+ let kernel = build_kernel(tau_r, tau_d, fs);
+
+ // Deterministic ~1 Hz spike train (xorshift; no rand/Date dependency).
+ let mut seed: u64 = 0x51A5_2025;
+ let mut next = || {
+ seed ^= seed << 13;
+ seed ^= seed >> 7;
+ seed ^= seed << 17;
+ (seed >> 11) as f64 / (1u64 << 53) as f64
+ };
+ let p = 1.0 / fs; // ~1 Hz
+ let mut s_true = vec![0.0_f32; n];
+ for v in s_true.iter_mut() {
+ if next() < p {
+ *v = 1.0;
+ }
+ }
+ let n_true: f32 = s_true.iter().sum();
+ let mut trace = vec![baseline; n];
+ for i in 0..n {
+ if s_true[i] > 0.5 {
+ for (k, &kv) in kernel.iter().enumerate() {
+ if i + k < n {
+ trace[i + k] += alpha_true * kv;
+ }
+ }
+ }
+ }
+ // SNR ~20 additive Gaussian noise.
+ let sigma = alpha_true / 20.0;
+ for v in trace.iter_mut() {
+ let u1 = next().max(1e-12);
+ let u2 = next();
+ let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
+ *v += (z as f32) * sigma;
+ }
+
+ let solve = |mass_count: bool| {
+ solve_trace_opts(
+ &trace, tau_r, tau_d, fs, factor, 1000, 1e-4, None, false, false, 0.0,
+ SolveOptions { noise_constrained: true, mass_count },
+ )
+ };
+
+ // Baseline (current behavior): alpha collapses, count inflates.
+ let off = solve(false);
+ let count_off: f32 = off.s_counts.iter().sum();
+ let alpha_ratio_off = off.alpha / alpha_true as f64;
+ assert!(
+ alpha_ratio_off < 0.4,
+ "sanity: without mass_count alpha should be badly under-estimated, got ratio {:.3}",
+ alpha_ratio_off
+ );
+ assert!(
+ count_off as f32 / n_true > 2.0,
+ "sanity: without mass_count count should be inflated, got ratio {:.2}",
+ count_off / n_true
+ );
+
+ // Fixed: alpha and count both recover.
+ let on = solve(true);
+ let count_on: f32 = on.s_counts.iter().sum();
+ let alpha_ratio_on = on.alpha / alpha_true as f64;
+ let count_ratio_on = count_on / n_true;
+ assert!(
+ (0.8..=1.2).contains(&alpha_ratio_on),
+ "mass_count alpha ratio should be ~1, got {:.3} (count_ratio {:.2})",
+ alpha_ratio_on,
+ count_ratio_on
+ );
+ assert!(
+ (0.75..=1.25).contains(&count_ratio_on),
+ "mass_count count ratio should be ~1 at 1 Hz, got {:.2} (alpha_ratio {:.3})",
+ count_ratio_on,
+ alpha_ratio_on
+ );
+ }
+
/// HP+LP filter path should produce valid results and return a filtered trace.
#[test]
fn filter_path_hp_lp() {
@@ -1047,6 +1410,7 @@ mod tests {
0.0,
SolveOptions {
noise_constrained: true,
+ mass_count: false,
},
);
diff --git a/crates/solver/src/js_indeca.rs b/crates/solver/src/js_indeca.rs
index 036bbdf..c1cfa84 100644
--- a/crates/solver/src/js_indeca.rs
+++ b/crates/solver/src/js_indeca.rs
@@ -35,6 +35,7 @@ pub fn indeca_solve_trace(
warm_counts: &[f32],
lambda: f64,
noise_constrained: bool,
+ mass_count: bool,
) -> Result {
if let Some(i) = crate::first_nonfinite(trace) {
return Err(JsError::new(&format!(
@@ -58,7 +59,10 @@ pub fn indeca_solve_trace(
hp_enabled,
lp_enabled,
lambda,
- indeca::SolveOptions { noise_constrained },
+ indeca::SolveOptions {
+ noise_constrained,
+ mass_count,
+ },
);
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 30e00a5..516b79e 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))]
+#[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, mass_count=false))]
#[allow(clippy::too_many_arguments)]
fn py_indeca_solve_trace<'py>(
py: Python<'py>,
@@ -439,8 +439,10 @@ fn py_indeca_solve_trace<'py>(
warm_counts: Option>,
lambda_: f64,
noise_constrained: bool,
+ mass_count: bool,
) -> PyResult<(
Bound<'py, PyArray1>, // s_counts
+ Bound<'py, PyArray1>, // s_rate (graded calibrated rate; empty unless mass_count)
f64, // alpha
f64, // baseline
f64, // threshold
@@ -463,11 +465,15 @@ fn py_indeca_solve_trace<'py>(
hp_enabled,
lp_enabled,
lambda_,
- indeca::SolveOptions { noise_constrained },
+ indeca::SolveOptions {
+ noise_constrained,
+ mass_count,
+ },
);
Ok((
PyArray1::from_vec(py, result.s_counts),
+ PyArray1::from_vec(py, result.s_rate),
result.alpha,
result.baseline,
result.threshold,
diff --git a/crates/solver/src/threshold.rs b/crates/solver/src/threshold.rs
index b782eec..09452f5 100644
--- a/crates/solver/src/threshold.rs
+++ b/crates/solver/src/threshold.rs
@@ -383,7 +383,7 @@ fn select_noise_floor_threshold(
/// Alpha is constrained to [0, max_alpha]. When max_alpha is f64::INFINITY
/// (the default from solve_trace), alpha is effectively uncapped — the
/// free-solve phase calibrates the prescale so alpha_lstsq lands near 1.0.
-fn lstsq_alpha_baseline(conv: &[f32], y: &[f32], pad: usize, max_alpha: f64) -> (f64, f64) {
+pub(crate) fn lstsq_alpha_baseline(conv: &[f32], y: &[f32], pad: usize, max_alpha: f64) -> (f64, f64) {
let n = y.len();
let lo = pad;
let hi = n.saturating_sub(pad);
diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 761c7ff..16e8536 100644
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -15,6 +15,21 @@ Versions correspond to git tags (`v*`) and apply to the entire monorepo.
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)
+- **CaDecon** calibrated continuous rate — when `mass_count` is on, `solve_trace`
+ additionally returns `s_rate`, a graded firing-rate estimate on the correct absolute
+ scale (the relaxed solution scaled by the single-spike mass and gated by realness).
+ Unlike the integer `s_counts`, it preserves the smoothed rate envelope, so it matches
+ the uncorrected bin-sum's correlation-based recovery while also being correctly scaled.
+ Exposed via the WASM solver result and the `calab.solve_trace` Python binding;
+ empty unless `mass_count=True`
+- **CaDecon** mass-based count readout — an optional `mass_count` mode that
+ counts each contiguous event by its deconvolved mass (calibrated to the
+ single-spike mass) and refits alpha, replacing the per-bin binarize→bin-sum.
+ Corrects the coherent-grid spike overcount that inflates counts and halves
+ alpha at high upsampling, while preserving multiplicity for temporally
+ resolvable bursts. Knob-free and off by default; exposed through the WASM
+ solver, the CaDecon UI, the bridge config, and the `calab.solve_trace` /
+ `calab.decon` Python bindings.
## [2.5.0] - 2026-07-08
diff --git a/packages/io/src/__fixtures__/decon-config-full.json b/packages/io/src/__fixtures__/decon-config-full.json
index 5123c4b..4a475a0 100644
--- a/packages/io/src/__fixtures__/decon-config-full.json
+++ b/packages/io/src/__fixtures__/decon-config-full.json
@@ -3,6 +3,7 @@
"upsample_target": 300,
"hp_filter_enabled": true,
"lp_filter_enabled": false,
+ "mass_count": true,
"max_iterations": 20,
"convergence_tol": 0.01,
"num_subsets": 4,
diff --git a/packages/io/src/bridge.ts b/packages/io/src/bridge.ts
index bf85030..0add485 100644
--- a/packages/io/src/bridge.ts
+++ b/packages/io/src/bridge.ts
@@ -22,6 +22,7 @@ export interface BridgeConfig {
upsample_target?: number;
hp_filter_enabled?: boolean;
lp_filter_enabled?: boolean;
+ mass_count?: boolean;
max_iterations?: number;
convergence_tol?: number;
num_subsets?: number;
@@ -36,6 +37,7 @@ export const BRIDGE_CONFIG_KEYS: readonly string[] = [
'upsample_target',
'hp_filter_enabled',
'lp_filter_enabled',
+ 'mass_count',
'max_iterations',
'convergence_tol',
'num_subsets',
diff --git a/python/docs/guides/cadecon.md b/python/docs/guides/cadecon.md
index 713d278..3fa0978 100644
--- a/python/docs/guides/cadecon.md
+++ b/python/docs/guides/cadecon.md
@@ -206,6 +206,7 @@ calab.solve_trace(
warm_counts: np.ndarray | None = None,
lambda_: float = 0.0,
noise_constrained: bool = False,
+ mass_count: bool = False,
) -> SolveTraceResult
```
@@ -223,8 +224,9 @@ calab.solve_trace(
| `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`. |
+| `mass_count` | Count each contiguous event by its deconvolved mass (calibrated to the single-spike mass) and refit alpha, instead of the per-bin binarize→bin-sum. Corrects the upsampling spike overcount / halved alpha; preserves multiplicity for resolvable bursts. Knob-free. Default `False`. See `docs/cadecon-mass-count.md`. |
-Returns a `SolveTraceResult` namedtuple with fields: `s_counts`, `alpha`, `baseline`, `threshold`, `pve`, `iterations`, `converged`.
+Returns a `SolveTraceResult` namedtuple with fields: `s_counts`, `s_rate`, `alpha`, `baseline`, `threshold`, `pve`, `iterations`, `converged`. `s_rate` is a calibrated **continuous firing-rate** estimate (graded, same absolute scale as `s_counts` but not rounded); it is populated only when `mass_count=True` and is an empty array otherwise. See `docs/masscount_R_metrics.pdf`.
Raises `ValueError` if `trace` (or `warm_counts`) contains a non-finite value (`NaN` or `Inf`) — the FFI boundary rejects non-finite input rather than returning garbage.
diff --git a/python/src/calab/_bridge/_apps.py b/python/src/calab/_bridge/_apps.py
index 1167283..a3abac3 100644
--- a/python/src/calab/_bridge/_apps.py
+++ b/python/src/calab/_bridge/_apps.py
@@ -208,6 +208,7 @@ def decon(
upsample_target: int | None = None,
hp_filter_enabled: bool | None = None,
lp_filter_enabled: bool | None = None,
+ mass_count: bool | None = None,
max_iterations: int | None = None,
convergence_tol: float | None = None,
num_subsets: int | None = None,
@@ -247,6 +248,10 @@ def decon(
Enable high-pass filter.
lp_filter_enabled : bool, optional
Enable low-pass filter.
+ mass_count : bool, optional
+ Mass-based count readout: count events by deconvolved mass and refit
+ alpha, correcting the upsampling spike overcount / halved alpha. See
+ ``docs/cadecon-mass-count.md``.
max_iterations : int, optional
Maximum solver iterations (1–200).
convergence_tol : float, optional
@@ -273,6 +278,7 @@ def decon(
upsample_target=upsample_target,
hp_filter_enabled=hp_filter_enabled,
lp_filter_enabled=lp_filter_enabled,
+ mass_count=mass_count,
max_iterations=max_iterations,
convergence_tol=convergence_tol,
num_subsets=num_subsets,
diff --git a/python/src/calab/_bridge/_models.py b/python/src/calab/_bridge/_models.py
index 3481889..15f87eb 100644
--- a/python/src/calab/_bridge/_models.py
+++ b/python/src/calab/_bridge/_models.py
@@ -19,6 +19,11 @@ class DeconConfig(BaseModel):
)
hp_filter_enabled: bool | None = Field(None, description="Enable high-pass filter")
lp_filter_enabled: bool | None = Field(None, description="Enable low-pass filter")
+ mass_count: bool | None = Field(
+ None,
+ description="Mass-based count readout: count events by deconvolved mass and "
+ "refit alpha, correcting the upsampling spike overcount / halved alpha",
+ )
max_iterations: int | None = Field(
None, gt=0, le=200, description="Maximum solver iterations (1-200)"
)
diff --git a/python/src/calab/_compute.py b/python/src/calab/_compute.py
index 713b557..209d501 100644
--- a/python/src/calab/_compute.py
+++ b/python/src/calab/_compute.py
@@ -284,6 +284,10 @@ class SolveTraceResult(NamedTuple):
----------
s_counts : np.ndarray
Spike counts at the original sampling rate, shape ``(n_timepoints,)``, float32.
+ s_rate : np.ndarray
+ Calibrated continuous firing-rate estimate (graded, same absolute scale as
+ ``s_counts`` but not rounded). Populated only when ``mass_count=True``; an
+ empty array otherwise. See ``docs/masscount_R_metrics.pdf``.
alpha : float
Amplitude scaling factor.
baseline : float
@@ -299,6 +303,7 @@ class SolveTraceResult(NamedTuple):
"""
s_counts: np.ndarray
+ s_rate: np.ndarray
alpha: float
baseline: float
threshold: float
@@ -356,6 +361,7 @@ def solve_trace(
warm_counts: np.ndarray | None = None,
lambda_: float = 0.0,
noise_constrained: bool = False,
+ mass_count: bool = False,
) -> SolveTraceResult:
"""Run the InDeCa pipeline on a single trace. Delegates to Rust.
@@ -385,6 +391,12 @@ def solve_trace(
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.
+ mass_count : bool
+ Mass-based count readout: count each contiguous event by its deconvolved
+ mass (calibrated to the single-spike mass) and refit alpha, instead of
+ the per-bin binarize→bin-sum. Corrects the coherent-grid spike overcount
+ that inflates counts and halves alpha at high upsampling. Knob-free.
+ Default False. See ``docs/cadecon-mass-count.md``.
Returns
-------
@@ -395,14 +407,15 @@ def solve_trace(
if warm_counts is not None:
warm = np.ascontiguousarray(warm_counts, dtype=np.float64)
- s_counts, alpha, baseline, threshold, pve, iterations, converged = _indeca_solve_trace(
+ s_counts, s_rate, alpha, baseline, threshold, pve, iterations, converged = _indeca_solve_trace(
trace_1d, tau_rise, tau_decay, fs,
upsample_factor, max_iters, tol,
hp_enabled, lp_enabled, warm, lambda_,
- noise_constrained,
+ noise_constrained, mass_count,
)
return SolveTraceResult(
s_counts=np.asarray(s_counts),
+ s_rate=np.asarray(s_rate),
alpha=float(alpha),
baseline=float(baseline),
threshold=float(threshold),
diff --git a/python/tests/test_solve_trace.py b/python/tests/test_solve_trace.py
index 369ddeb..7a9eb21 100644
--- a/python/tests/test_solve_trace.py
+++ b/python/tests/test_solve_trace.py
@@ -111,6 +111,47 @@ def test_noise_constrained_detects_events_with_noise(self):
assert result.s_counts.sum() >= 1
assert result.pve > 0.5
+ def test_mass_count_accepted(self):
+ # The mass_count 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, mass_count=True)
+ assert isinstance(result, SolveTraceResult)
+ assert result.s_counts.shape == (300,)
+ assert result.s_counts.sum() >= 0
+
+ def test_mass_count_recovers_alpha_at_upsample(self):
+ # At upsampling the default readout inflates the count and halves alpha;
+ # mass_count should recover an alpha much closer to truth on clean, well-
+ # separated spikes.
+ alpha_true = 5.0
+ trace = _make_trace(
+ 0.1, 0.6, 30.0, 1200, [150, 450, 750, 1050], alpha=alpha_true, baseline=2.0
+ )
+ off = solve_trace(trace, 0.1, 0.6, 30.0, upsample_factor=10)
+ on = solve_trace(trace, 0.1, 0.6, 30.0, upsample_factor=10, mass_count=True)
+ # Default path under-estimates alpha; mass_count lands near truth.
+ assert on.alpha > 1.5 * off.alpha
+ assert abs(on.alpha - alpha_true) / alpha_true < 0.2
+
+ def test_s_rate_graded_and_mass_count_only(self):
+ # s_rate is the graded, calibrated continuous rate: populated only under
+ # mass_count, fractional (not integer), same length as the trace, and
+ # integrating to roughly the true spike count.
+ n_true = 5
+ trace = _make_trace(
+ 0.1, 0.6, 30.0, 1500, [150, 450, 750, 1050, 1350], alpha=5.0, baseline=2.0
+ )
+ off = solve_trace(trace, 0.1, 0.6, 30.0, upsample_factor=10)
+ on = solve_trace(trace, 0.1, 0.6, 30.0, upsample_factor=10, mass_count=True)
+ # empty without mass_count, populated with it
+ assert off.s_rate.size == 0
+ assert on.s_rate.shape == (1500,)
+ # graded: has fractional values strictly between 0 and 1
+ assert bool(((on.s_rate > 1e-4) & (on.s_rate < 0.99)).any())
+ # calibrated: integrates to roughly the true count
+ assert abs(on.s_rate.sum() - n_true) < 0.5 * n_true
+
# ---------------------------------------------------------------------------
# estimate_kernel