PIPELINE_PORT_STATUS.md

MEA-NAP Python Port — Status Handoff

Living status doc for the MATLAB → Python port of the MEA-NAP pipeline. Read this first before making changes to src/meanap/pipeline/ or the Pipeline/Test-pipeline GUI wiring — it captures decisions and gotchas that aren't obvious from the code alone.

Where things stand

Step	MATLAB source	Python status
1. Spike detection	Functions/WATERS-master/*, detectSpikes*.m	Done, validated against MATLAB reference output, wired into the GUI (Run + Test pipeline)
2. Neuronal activity (firing rates, burst detection)	Functions/firingRatesBursts.m, Functions/singleChannelBurstDetection.m	Done, validated against MATLAB reference output (100% parity on recording- and node-level fields), wired into runner.py. CSV export now done too (NeuronalActivity_RecordingLevel.csv/NeuronalActivity_NodeLevel.csv, port of saveEphysStats.m) — a 2026-07-08 audit found the ephys dict itself was complete but was only ever written to ephys_results.json, never flattened into the two CSVs MATLAB's own pipeline produces.
3. Functional connectivity (STTC)	Functions/generateAdjMs.m, Functions/STTCandThresholding/*	Core (STTC) done, exact parity. Probabilistic thresholding ported but inherently non-bit-reproducible (see below)
4. Network metrics	Functions/ExtractNetMet.m, Functions/2019_03_03_BCT/*	Deterministic subset done (ND, NS, MEW, Dens, CC_raw, PL_raw, Eglob, Eloc, BC, NE), 100% parity. Modularity-dependent subset also done (Ci/Q/nMod via Louvain + consensus clustering, raw + normalized PC, Z, node cartography 6-role classification, Hub3/Hub4, rich club RC) — 100% parity for everything downstream of a fixed Ci (and, for PC-normalization, a fixed PC_norm too); the stochastic pieces themselves (Ci, PC_norm's null-model randomization) aren't bit-reproducible, same situation as Step 3. Controllability (aveControl/modalControl) also done. Small-worldness also done — SW/SWw and the saved, null-model-normalized CC/PL (formula assembly has exact parity against MATLAB given the same A/R/L, see test_pipeline_small_worldness.py; the null models themselves, randmio_und_v2/latmio_und_v2, aren't bit-reproducible, same situation as everywhere else in this table) — see network_metrics.py. effRank done (network_metrics.effective_rank, port of calEffRank.m). NMF (num_nnmf_components/nComponentsRelNS/nnmf_residuals/nnmf_var_explained) also done — see nmf.py; not just RNG-different from MATLAB but algorithm-different (sklearn NMF solvers vs. MATLAB's nnmf), so treat this one as looser-than-usual parity. Record-level summary-stat scalars now done too (NDmean, NDtop25, NSmean, sigEdgesMean, sigEdgesTop10, PCmean, PCmeanTop10, PCmeanBottom10, percentZscoreGreaterThanZero, percentZscoreLessThanZero) — a 2026-07-08 audit against MATLAB's full default netMetToCal list (see "Auditing for silently-missing metrics" below) found these were in plotting_step4.py's display-name dict (implying they were intended) but never actually computed anywhere in the port.
5. Statistics & machine learning	Functions/doStats.m, featureCorrelation.m, doLDA.m, doClassification.m	Done, and deliberately not a literal port — the MATLAB version is unfinished (its two-way ANOVA is commented out, its LME's random effect has one observation per level, its cross-validation leaks cultures across folds, and there is no regression at all). See "Step 5: statistics and machine learning" below for what was kept, what was changed, and why.

All four pipeline steps are wired into runner.py and reachable from the GUI (start/stop_analysis_step now goes up to 4). Step 5 is deliberately not part of run_pipeline: it analyses a run that has already finished, so it has its own entry points (meanap-stats, the Stats & ML tab) and can be pointed at any output folder or bundle, including one from an earlier session. Output folder structure (the static tree CreateOutputFolders.m builds) is fully ported and used automatically whenever the pipeline runs — diffed directly against OutputData03Mar2026/ (a real MATLAB run) and matches.

End-to-end CSV parity report (2026-07-15/16): python/parityRun/PARITY_REPORT.md. Ran the full Test-pipeline dataset through both ports with matched parameters and compared all four CSVs. Headline: every deterministic stage now reproduces MATLAB exactly — spike detection (bior1p5/thr5 identical counts and times, F1 = 1.0000, after fixing 3 real bugs — see "Spike detection gotchas"), step-2 activity metrics (FRmean diff = 0 end to end), and step-4 network metrics given the same adjacency (worst diff 2.8e-14, controllability included). The only remaining end-to-end divergence is step 3's probabilistic thresholding RNG, which changes aN and hence the network CSVs; effRank (1.2e-05) and aveControlMean (1.4e-04) are near exact because they barely depend on marginal edges. Still not comparable 1:1: per-node network values (r=0.68-0.83) and the stochastic metrics (Q/nMod/SW/PC/Z/NCpn*); num_nnmf_components differs by solver, not RNG. Two fixes landed from this work: the CWT rewrite (also ~2x faster, step 1 272s -> 141s) and NetworkActivity_NodeLevel.csv's Channel column, which was a position index rather than an electrode ID. The report also lists the remaining structural CSV mismatches and 4 further headless-MATLAB traps beyond the 5 in MEApipeline_timingBenchmark.m.

NB: OutputData03Mar2026/ no longer exists on this machine, so test_pipeline_step1/step4.py and the fixtures that reference it can no longer run as written (step 4 skips every check and reports "No checks ran" — a pass/fail that silently means nothing was tested). python/parityRun/ contains a fresh MATLAB reference run plus a regenerated step-4 fixture generator (gen_step4_parity_reference.m) if these need rebasing.

MATLAB is available on this machine (/usr/local/bin/matlab, licensed, confirmed working). This was the key that unlocked exact-parity testing for steps 3-4: MATLAB's own pipeline never persists the intermediate values needed to isolate each step's arithmetic (e.g. the raw unthresholded STTC matrix, or network metrics computed from a fixed adjacency matrix), so python/test_fixtures/gen_sttc_reference.m and gen_step4_reference.m call the relevant MATLAB/BCT functions directly to generate ground truth, saved as .npz fixtures (see "How to verify changes" below to regenerate them). Don't assume this MATLAB install will be present in every future environment — treat the fixtures as the durable artifact, the .m scripts as how to regenerate them if the underlying algorithms ever need re-validating.

Step 5: statistics and machine learning (src/meanap/stats/, 2026-08-30)

The MATLAB pipeline's last stage (MEApipeline.m:1356-1378) asks four questions of a finished run's metric tables. All four are ported; a fifth, which MATLAB does not ask, is added.

Question	MATLAB	Python
Does a metric differ by age or genotype?	doStats.m	stats/comparisons.py
Which metrics measure the same thing?	featureCorrelation.m	stats/correlation.py
Can genotype/age be decoded from the features?	doClassification.m, doLDA.m	stats/decoding.py
How much of the target does each feature explain?	—	stats/regression.py
Which features carry the genotype signal at each age, and does that change?	—	stats/decoding.py (decoding_shapley)
Is it activity, correlation strength, or genuine network topology?	—	stats/decoding.py (family_shapley)
Does the topology difference survive at matched connection density?	—	stats/density_sweep.py

Entry points: meanap-stats <run-or-bundle> (CLI), the Stats & ML tab (gui/panels/stats.py + gui/stats_worker.py), or run_stats() directly. It reads a finished run — an output folder or a .meanap bundle — so it works on runs from earlier sessions and on bundles someone sent, and it writes 5_StatsAndML/<lag>/ with a CSV and a figure for everything it computes.

The one substantive change: recordings are not independent samples

Read this before changing anything here. The same culture is imaged at several DIVs. On the Yin timecourse that is 378 recordings from 121 cultures, median 4 recordings each. MATLAB treats the 378 rows as 378 independent samples everywhere:

• doStats's LME uses (1|recordingName), giving each recording its own random intercept. Each recording appears once, so that random effect has one observation per level and absorbs nothing — it is a no-op dressed as a repeated-measures model.
• doClassification calls crossval(..., 'KFold', 2) with no grouping, so the same culture lands in train and test. The classifier recognises cultures it has already met at another DIV, and the reported accuracy measures culture identity as much as genotype.

Here, culture identity is derived once (dataset.derive_culture_ids — drop the DIV token and any date token from the recording name; falls back to one-culture-per-recording when the names do not parse, so an unrecognised naming scheme makes the analysis conservative, never wrong) and used throughout: as the random-effect grouping in every mixed model, and as the grouping for StratifiedGroupKFold/GroupKFold in every cross-validation and in the permutation null (which permutes labels between cultures, so a culture keeps one genotype as it does in the real data).

test_stats.py's Section A5 has the tripwire for this: a synthetic feature that is constant within a culture and random between them, carrying no group signal. With culture-grouped folds it decodes genotype at chance, as it must. With row-wise folds it would score near-perfectly.

Other deliberate divergences from doStats.m

• Omnibus before pairwise. A three-level genotype gets a joint Wald test across its dummies; the age×genotype interaction gets a joint Wald test of the interaction terms. (MATLAB fits the interaction model, compares it to the main-effects model, then discards the comparison and uses the main-effects model regardless — see the commented-out block at doStats.m:117-125.)
• Multiple comparisons. ~50 metrics × several tests is ~2,900 p-values on the Yin run. Every p-value carries a Benjamini-Hochberg partner corrected within its own named family. MATLAB corrects nothing.
• Hedges' g, not MATLAB's "d-prime", which divides by the mean of the two SDs rather than their pooled value and has no small-sample correction.
• Mixed model instead of RMAOV1. A repeated-measures ANOVA needs every subject at every level and silently drops the rest; the mixed model uses the unbalanced data as it is. (MATLAB works around this by intersecting recording names across DIVs first, which throws away most of the data.)
• Kruskal-Wallis rather than anova1 for the per-age omnibus: several of these metrics are proportions bounded at 0/1 and visibly non-normal. Pairwise follow-ups are Welch t-tests.
• Balanced accuracy, not the misclassification rate, and chance is 1/n_classes, not the majority-class rate MATLAB draws as "random chance".

The new part: variance attribution (stats/regression.py)

The effective dimensionality of the ~50 recording-level metrics on the Yin run is about 4.6 — 37 pairs correlate above |r| = 0.9, and aN and numActiveElec correlate at 0.99997 (they are the same quantity computed by two steps). Under that much collinearity a coefficient or a single-model importance is close to arbitrary. So each feature gets three numbers:

• Marginal — R² of that feature alone (its ceiling);
• Unique — drop in R² when it is removed from the full model (its floor, ≈0 for anything with a near-duplicate);
• Shapley — its average marginal contribution over feature orderings (the LMG decomposition, Monte Carlo over n_orderings random permutations). This is the number to quote: non-negative, and the column sums exactly to the full model's R², so each entry is a genuine share of explained variance.

Worked example from the Yin run: Dens has Marginal = 0.40 (density alone explains 40% of the age variance) but Unique = 0.000008 and Shapley = 0.043. Reading Unique alone would say density is irrelevant to age, which is wrong; the Shapley value says it holds a 4% share of a signal it splits with its correlates.

test_stats.py checks this against a case with an exactly known answer: with orthogonal predictors the Shapley value must equal both the marginal and the unique R² for every feature, and it does to 1e-9.

Attributing the decoding across features, per age (2026-08-30)

variance_decomposition partitions R² for a continuous target. decoding.decoding_shapley is its counterpart for a categorical one: it partitions how far a decoder gets above chance, separately within each age, so the question "which metrics carry the genotype difference, and does that change as the culture matures" has an answer with numbers on it. Two figures come out — 5C6_genotype_shapley_by_age (features × ages heatmap) and 5C7_genotype_shapley_across_ages (trajectories over the total they partition) — plus three CSVs.

Three properties worth not breaking:

• The shares at one age sum exactly to that age's decodability (cross-validated balanced accuracy − chance). value(∅) is defined as 0 rather than fitted, which is what makes the identity hold. Asserted in test_stats.py.
• Shapley values here can be negative, unlike the regression version. Nested least-squares R² can only rise as features are added; cross-validated accuracy cannot make that promise, so a negative share is a real finding — the metric cost the decoder accuracy at that age — not a bug.
• The feature set is chosen once, on the pooled data, without looking at the labels. select_representatives clusters on 1 - |r| and keeps the most central member of each cluster (50 metrics → 15 on the Yin run, Density standing for 13 of them). Target-blind is the load-bearing part: screening features by how well they separate the groups and then reporting how well those features separate the groups would inflate every total here. Choosing once rather than per age is what makes the trajectories comparable.

Cost. Roughly n_ages × n_orderings × max_features subset evaluations, each a small cross-validated fit (~11 ms), with subset scores cached — the short prefixes every ordering shares are paid for once. The defaults (shapley_max_features=15, shapley_orderings=100) take a couple of minutes on a five-age run; --shapley-orderings 0 skips it.

Gotcha: don't assert where a Shapley value peaks. A share is a marginal contribution to accuracy, and accuracy saturates — once a feature nearly separates the classes on its own, its share stops growing however much stronger the underlying effect gets. A test that planted a linearly-growing effect and asserted the share peaked at the oldest age failed for exactly this reason. The robust claim is the crossover between features, which is what the test checks now.

Result on the Yin timecourse (worth knowing as a sanity check on the method): decodability is highest at DIV 14 (0.37 above chance) and roughly halves by DIV 28. Density and mean node degree carry the early signal and fade; mean event rate is ≈0 at DIV 14 and becomes the largest single contributor by DIV 35-42. The genotype signature moves from network structure to activity rate over development.

Activity vs correlation strength vs topology (2026-08-30)

The per-feature attribution says which metric carries the genotype signal. It cannot say whether an apparent difference in network organisation is real, or is what you would expect anyway from a culture that fires more and is more strongly correlated — density and global efficiency correlate at 0.98 on the Yin run, so no per-feature ranking separates those two readings.

decoding.family_shapley does, by making a whole family one player:

Family	What it is	Metrics on the Yin run
activity	firing properties, computed before any pairwise measure exists	11
coupling	how much correlation there is — density, degree, node strength, significant edges	6
topology	how that correlation is arranged — efficiency, clustering, path length, modularity, cartography, controllability	33

The mapping is the explicit dict FEATURE_FAMILIES, overridable per call and written out as feature_family_membership.csv so it can be checked and argued with. Anything unlisted becomes an other family rather than being guessed at or silently folded in. Two placements are genuinely arguable and are commented in the source: aN/numActiveElec (how many cells are active is a fact about firing) and effRank/num_nnmf_components (computed from the activity matrix, not the adjacency matrix).

Two figures: 5C8_feature_family_contributions (grouped bars per age, summing to the total) and 5C9_feature_family_alone_vs_unique.

Why this is more robust than the per-feature view. Collinearity inside a family stops mattering — thirteen metrics that are all density in disguise are one player, so no credit is split arbitrarily between them, and this needs none of the redundancy reduction decoding_shapley requires. Collinearity between families is what the Shapley value is for: the signal density and global efficiency share is split in proportion to how often each is the one that adds it, over all orders of entry.

Exact, not sampled. Three families means eight subsets per age, so the values come from enumeration (_exact_shapley) rather than Monte Carlo. The whole thing takes under a second on the Yin run, and re-running gives bitwise identical numbers — asserted in test_stats.py.

Gotcha: alone and unique do NOT bracket the Shapley value. That bracket holds for a submodular game; cross-validated accuracy is not one. Families can be complementary — on the test fixture, activity's unique (0.65) exceeds its alone (0.48) at DIV 28, because the classifier can only use it well once the other families are present. A test asserting the bracket failed for exactly this reason, and the figure's caption was rewritten to give the unique > alone case its own reading rather than describing the gap as "shared signal" in both directions.

Result on the Yin timecourse. Network topology's alone is 0.39 at DIV 14 but its unique is 0.000, and it goes negative (−0.08) by DIV 35: on this dataset topology carries essentially nothing about genotype that activity and correlation strength do not already carry. Correlation strength leads the Shapley split at DIV 14 (0.165, 44% of the total) and fades. Activity's unique contribution grows monotonically with age (−0.01 at DIV 14 to +0.11 at DIV 42) and is the only family with independent signal by DIV 42.

Density sweep: topology at matched connection density (2026-08-30)

stats/density_sweep.py. The control for the confound the section above documents: every graph metric depends on how many edges the network has, and MEA-NAP's null-model normalisation does not remove that (a degree-preserving null preserves density exactly, so the ratio is density-conditioned, and it saturates as the graph approaches complete).

The sweep thresholds every recording to a common proportion of edges and repeats across a range — 2% to 40% in 2% steps, the grid Schroeter et al. (2015, J Neurosci) used on the closest design to MEA-NAP's (dissociated cultures on MEAs, density climbing over four weeks in vitro). Two figures split by group and by age, one context figure, three CSVs, and a set of cost-integrated features (Ginestet et al. 2011): each metric integrated over the density range, divided by the range so it stays on the metric's own scale, giving one density-controlled number per recording named <metric>_costInt. Those are registered in FEATURE_FAMILIES as topology, so they flow into the family decomposition unchanged.

Off by default (StatsSettings.density_sweep, --density-sweep, or the Density sweep checkbox on the GUI's Stats & ML tab): it is a robustness analysis, and unlike the other four it reads the per-recording adjacency matrices rather than the metric CSVs. The GUI disables the checkbox outright when the chosen run has no ExperimentMatFiles (panels/stats.py's has_adjacency, which handles both folders and bundles) — finding that out after twenty minutes of other analyses would be a poor way to learn it.

Binary at each density, deliberately. Surviving edges still carry weights, and those weights still differ between groups, so a weighted sweep removes only half the confound. Discarding them is what proportional thresholding means in this literature and is the only version that isolates topology. It also makes Dens/ND/NS/MEW constant by construction, which is why they are not swept.

It reads the pre-thresholding matrix. Step 3 stores both adjM{lag}mslag and adjM{lag}mslag_raw; the sweep takes _raw where it exists (CAT-NAP stores one matrix, adj__-prefixed, which is pre-threshold). Imposing 2% density on a significance-thresholded matrix would be thresholding twice, and the sparsest targets would be unreachable for exactly the recordings whose connectivity is weakest — the bias van den Heuvel et al. (2017) warn about.

Gotchas

• lccFraction is not decoration. These networks fragment badly at sparse densities — on the Yin run the largest component holds only ~33% of nodes at 2%, reaching ~0.95 around 26%. charpath (like BCT's) averages over connected pairs only, so a more fragmented network can show a shorter path length. Below the density where lccFraction approaches 1, read Eglob (which scores a disconnected pair as zero) and not PL. The run summary reports density_where_network_intact for this reason.
• The swept range sits entirely below the observed one (2-40% against an observed median of 0.86). That is the point — proportional thresholding imposes a common density regardless — but it does mean the networks being compared are not the ones step 4 measured, and the context figure exists to make that visible rather than implicit.
• Ties are all kept at the cutoff, so a realised density can exceed its target slightly. Breaking ties arbitrarily would make the result depend on node ordering.
• Matching density does not match SIZE. van Wijk et al.'s title is "networks of different size and connectivity density"; thresholding addresses the second half only. It matters: on the Yin run node counts differ 2.4x between genotypes (median 143 KO, 93 Het, 59 WT), and with density alone controlled the Spearman correlation with node count at 20% density is Eglob -0.01, Q -0.04, ElocMean +0.34, PL +0.42, CC +0.43, nMod +0.47, BCmean -0.64 — only Eglob and Q are clean. sweep_n_nodes closes that half: every network is reduced to a random induced subgraph of a common size before thresholding (that order matters — the reverse lands at an uncontrolled density), repeated over sweep_subsamples draws and averaged. Default "auto", which takes the sweep_node_percentile (20th) percentile of the run's own counts.
• Subsampling trades the size confound for a SELECTION one, and the trade has to be watched. Recordings with fewer nodes than the target are dropped, and smallness is not randomly distributed: on the Yin run a target of 80 nodes keeps 76% of KO recordings against 34% of wildtype, which would compare the largest wildtype networks with typical KO ones. retention() reports kept-fraction per group, it is written to density_sweep_retention.csv, and the runner logs a warning when the spread between groups exceeds 20 points. The default percentile is deliberately low for this reason, and _resolve_target_nodes refuses to subsample below 15 nodes at all rather than compare near-empty graphs.
• BCmean is normalised by ((n-1)(n-2)), as step 4 and MATLAB's ExtractNetMet.m both do. The sweep originally omitted that and the unnormalised mean correlated with node count at r = +0.89 — a size measure wearing a topology name. Normalised it is r = -0.64, still the most size-dependent metric in the set, which is worth knowing before quoting it.
• Cost is dominated by the big networks and by modularity. Consensus Louvain at rep_num=50 (what step 4 uses) is ~0.8s per density per recording; the sweep defaults to sweep_modularity_reps=20. On the Yin run (378 recordings, median 131 nodes but up to 652, 20 densities) the whole sweep took roughly 15 minutes across 16 cores. Do not estimate it from a small subset: wall time on a parallel subset is the slowest single recording, not the mean, and dividing that by the subset size then re-multiplying double-counts the parallelism (an error made and corrected during this work). betweenness_wei is the hot spot and is pure Python — only the distance loop is numba-compiled.

Topology under control, fed back into the family decomposition

run.py's _controlled_families runs the family decomposition a second time with the topology family represented only by the sweep's cost-integrated features, and draws 5E4_topology_controlled comparing the two. Activity and correlation strength stay raw in both, deliberately: they are the confounders topology is being tested against, and controlling them away would defeat the comparison. Runs automatically whenever both feature_families and density_sweep are on.

Result on the Yin timecourse. Controlling topology does not rescue it: topology's mean unique contribution is -0.009 controlled against -0.019 raw, and per age it is 0.000 / -0.013 / -0.036 / -0.007 / +0.012 — zero everywhere. Alone actually falls (0.221 -> 0.159). A feature-count control matters here and was run: matching the raw set to the same six quantities (CC_rawMean, PL_raw, Eglob, ElocMean, Q, nMod) moves raw unique from -0.019 to -0.004, so a little of the original gap was the 33-vs-7 asymmetry — but the controlled version is still no better than the matched raw one.

This reconciles with the sweep curves rather than contradicting them: there is a real WT > KO clustering difference at matched density and size (~0.28 SD), and it carries no genotype information beyond what activity and connection density already carry. Topology here is redundant, not absent — which is a stronger claim than "mismeasured", and only sayable because the control was run.

Choosing the subsample target: separate selection from resolution

Raising sweep_n_nodes changes two things at once — topology is better resolved and a different (smaller, non-random) set of recordings survives — so comparing two targets directly confounds them. Run three conditions instead:

	what it is
A	target N over every recording that reaches it
B	a larger target M over every recording that reaches it
C	target N over only B's recordings

A vs C is then the selection effect alone, and C vs B the resolution effect alone. Done on the Yin run at N=22 / M=50 (WT−KO gap, in SDs of each metric's own spread across the sweep):

	CC	Eglob	ElocMean	Q	nMod	BCmean
A (N=22, 345 recs)	0.284	0.012	0.232	0.052	−0.199	0.074
C (N=22, B's 278)	0.262	0.015	0.221	0.166	−0.203	0.139
B (N=50, 278 recs)	0.341	0.046	0.292	0.420	−0.252	0.285

Reading: CC, ElocMean and nMod barely move on selection (A→C ≤0.02) and keep sign and rough magnitude at the larger target, so those differences are real. Q and BCmean move a lot on selection alone (0.052→0.166 and 0.074→0.139), so a Q or BC difference measured at a high target is partly the retention imbalance — at N=50 that was WT 63% vs KO 87%. Resolution (C→B) is the larger term for every metric, which means effect sizes are not comparable between targets and one target's numbers should not be quoted as absolute.

Result on the Yin timecourse, with density and size controlled: a modest but robust WT > KO difference in clustering and local efficiency, and more modules in KO; nothing in Eglob. Far smaller than the uncontrolled comparison showed — the large effects were density and size — but not zero, which amends the family-Shapley conclusion rather than contradicting it.

Measured runtimes (378 recordings, 20 densities, 20 draws, 16 cores): N=50 took 25 min, N=22 over 278 recordings took 9 min. Subsampling multiplies the sweep by the draw count, so it is much more expensive than the unsubsampled sweep (~15 min) despite the smaller graphs. Per-recording cost is close to linear in the target, not cubic as the betweenness term suggests — measured over a five-target ladder at 10 draws on 32 cores: N=20 129s, N=30 178s, N=40 242s, N=60 351s, N=80 452s, 23 min for all five.

Selection and resolution, in code (2026-09-02)

The A/B/C design above is no longer a manual recipe. The asymmetry that makes it cheap is that C is free: a recording's swept result depends only on its matrix, the target and the seed, never on which cohort it sits in, so C is a row subset of A rather than a second sweep. Only B needs one.

• The selection half runs automatically whenever subsampling is on. density_sweep.selection_sensitivity re-reads the sweep just computed over only the recordings clearing a higher node count (sweep_selection_percentile, p40 by default), recomputes every pairwise group gap, and writes density_sweep_selection_check.csv. A gap is flagged when it moves by at least 0.1 SD and either changes sign or moves by half its own size — both conditions, since a large relative move on a gap of 0.02 is noise and a small relative move on a gap of 1.0 is not worth a warning. Flagged gaps are logged and land in the run summary under density_sweep_selection.
• The full decomposition is opt-in: --sweep-nodes 22,50 (or sweep_n_nodes_high) sweeps at both targets and writes density_sweep_target_decomposition.csv with GapA/GapC/GapB and the two terms they split into, Selection (C − A) and Resolution (B − C), plus a Dominant column naming the larger. This doubles the sweep's cost.
• Two figures cover the size half (2026-09-02), which until then had none — the density half had three while size, the larger confound here, was CSV-only. 5E5_size_control is the counterpart to 5E3: node counts per recording against the target they were reduced to, with everything below it marked as dropped, and retention per group beside it, titled in the warning colour when the spread exceeds the same 20 points the runner warns on. 5E6_selection_check draws the check above as one dumbbell per metric per group pair — open marker over every swept recording, filled over the restricted cohort — with flagged comparisons in the warning colour. Both are offered by the viewer and the HTML report, and both redraw from the CSVs alone: _load_sweep recovers the target from NNodesUsed and _load_selection rebuilds the check from its table. The target decomposition is still CSV-only. 5E7-5E9 cover the effects under control; see the section below.
• Every gap is standardised within its own condition, which is what an analysis restricted to that cohort would actually report. So a delta mixes a change in group means with a change in spread — deliberately, because both are things the restriction did to the number you would quote. It also means a delta is never exactly zero even when the means are untouched.
• group_gaps excludes lccFraction: comparing it between groups is a statement about connectedness, not topology.

Result on the Yin timecourse (378 recordings, 10 draws, N = 20/30/40/60/80). Holding the cohort fixed at N=80's 206 recordings, so only resolution varies, the WT−KO gap barely moves across a 4x range in target:

	N=20	N=30	N=40	N=60	N=80
CC	0.888	1.005	1.048	1.029	0.994
ElocMean	1.016	1.138	1.144	1.110	1.073
nMod	−0.820	−0.824	−0.758	−0.716	−0.705
BCmean	0.264	0.265	0.295	0.320	0.311
Q	0.371	0.342	0.482	0.538	0.609
Eglob	0.124	0.127	0.109	0.183	0.162
PL	0.074	0.085	0.141	0.145	0.144

Nothing flips sign and the headline metrics move 12–19%; Q is the exception, nearly doubling. (Standardised per condition against the cost-integrated feature's own spread, so these are on a different scale from the N=22/50 table above — compare within a table, not between.)

The selection term is the larger one. Natural cohort minus fixed cohort at N=20: CC −0.115, ElocMean −0.171, nMod +0.194, BCmean −0.164, Q −0.235 — for CC that exceeds the entire N=20→80 resolution term (0.106). Its sign also runs opposite to the earlier N=22/50 pair, so it depends on which recordings the cut removes and has to be measured per dataset rather than assumed. This is the case for running the check rather than picking a target and hoping.

Subsampling does remove the size dependence it was built for. Residual Spearman against the original node count, fixed cohort, flat across N:

	CC	PL	ElocMean	BCmean	nMod	Eglob	Q
density only	+0.43	+0.42	+0.34	−0.64	+0.47	−0.01	−0.04
+ subsampling	−0.05	−0.03	−0.10	−0.09	+0.12	−0.20	−0.22

BCmean −0.64 → −0.05 is the striking one. But subsampling introduces a ~−0.2 residual in Eglob and Q, the two that were already clean under density control alone — worth knowing before trusting Q under this control.

Effects before and after the controls (2026-09-02)

The sweep figures answer "do the curves separate". The question a reader of the 5A heatmaps actually asks is "does this effect survive the control", and density_sweep.control_effects answers it directly: the same statistic 5A reports, recomputed on the controlled features. Standardised betas from the mixed model and its "SD change across age range" for the pooled terms; Hedges' g from the per-age Welch tests for the per-age genotype contrasts, kept apart by a Scope column because the two are not on one scale and must never share a panel. control_values reports the metrics themselves — mean, SD, SEM, N, median per (control, metric, group, age) — because an effect size in units of each condition's own spread cannot say what either group measured.

Three levels of control, written to control_effects.csv and control_metric_values.csv:

	what it is
none	step 4's own metrics, each network at its own density and size
density	cost-integrated over the density grid, size left alone
density+size	the same, on networks first cut to a common node count

none→density is the density control alone and density→density+size is the node-subsampling effect alone. The middle condition needs the sweep run twice, so it is opt-in: --sweep-density-only / sweep_density_only. Without it the outputs carry two conditions and cannot attribute a change to subsampling; the figure says so in its own caption rather than letting the reader assume otherwise.

Raw counterparts are the unnormalised ones. RAW_COUNTERPART pairs CC with CC_rawMean and PL with PL_raw, not with the null-model-normalised NetMet.CC/NetMet.PL saved under the same short names — the sweep measures raw quantities on binary graphs, and pairing against the normalised ones would compare two different quantities. BCmean has no recording-level counterpart in step 4 at all and so appears only in the controlled conditions; 5E9 labels that panel rather than leaving it blank.

Every condition is put on one cohort (match_cohort, on by default). A subsampled sweep drops every recording too small to reach the target while the raw metrics and an unsubsampled sweep keep all of them, so without this the step between two conditions mixes the control being applied with the cohort changing underneath it. This is not a hypothetical: on the Yin run the unmatched tables read CC's pooled age effect as none +1.09, density −0.53, density+size −1.06, which looks like subsampling doubling the effect. Matched on the 306 recordings all three conditions have, it is +1.24, −1.13, −1.06 — the density control does essentially all of the work and subsampling adds nothing. The unmatched reading would have attributed to node subsampling something that was the retention cut. The cost is that the none column no longer matches the 5A tables, which is the honest trade: those answer "what does this run show", this one answers "what did the control do".

Result on the Yin timecourse (pooled age effect, SD change across the age range, 306 recordings, matched cohort):

	none	density	density+size
CC	+1.24	−1.13	−1.06
ElocMean	+1.16	−1.11	−1.08
Eglob	+1.67	−0.14	+0.18
PL	−1.75	−1.02	−0.66
Q	−0.72	−0.84	−0.70
nMod	−0.93	+0.56	+0.98
BCmean	—	−0.54	−0.51

CC, ElocMean and nMod change sign under control: clustering and local efficiency rise steeply with age at each network's own density and fall once density is matched, which is the confound stated as plainly as it can be. Almost all of that is the density control; subsampling moves nMod (+0.56 → +0.98) and PL (−1.02 → −0.66) and little else. Eglob's large raw age effect (+1.67) essentially vanishes.

Three figures. 5E7_control_effects is the pooled view, one dumbbell per metric per term with filled markers surviving FDR. 5E8_genotype_effect_by_age is the per-age view — one panel per (metric, genotype pair), age on the x axis, a line per control — which is where a difference that only appears late, or only before the control, becomes visible as a shape. 5E9_metric_values_by_control is the metrics themselves, mean ± SEM against age with one line per group, repeated per control; its controlled columns share a y axis so the subsampling step reads directly, while the uncontrolled column is scaled on its own because it is a different measurement (raw nMod is ~2 against ~10-35 controlled, and a shared axis flattens it to a line).

How small can the target be? Draw noise vs attenuation

Two different costs, and only one of them is fixed by more draws. Measured on 30 large Yin recordings, 20 independent subgraphs each at a fixed 20% density, variance split into draw noise and real between-recording signal:

N	CC signal SD	reliability, 1 draw	10 draws	20 draws
20	0.047	0.29	0.80	0.89
30	0.050	0.40	0.87	0.93
40	0.058	0.61	0.94	0.97
60	0.058	0.73	0.96	0.98
80	0.059	0.86	0.98	0.99

• A single draw at a small target is worthless — at N=20 the draw-to-draw SD exceeds the between-recording SD (ratio 1.5 for CC, 1.6 for ElocMean). The sweep_subsamples=20 default is load-bearing, not polish.
• Averaging fixes noise but not attenuation. Real signal SD still rises 0.047 → 0.059 from N=20 to N=80 after draw noise is removed: a 20-node subgraph genuinely makes different networks look more alike. So effect sizes shrink toward zero as the target falls — a null at a small target is weak evidence, a positive one is conservative.
• For per-recording features, N≥40. Rank agreement with the well-resolved N=80 values: at N=20 it is CC 0.83, ElocMean 0.76, Q 0.75, PL 0.67; by N=40 everything except PL is ≥0.90. Group means tolerate N=20; feeding *_costInt into decoding or regression does not. This is why sweep_node_percentile was raised from p10 to p20 on 2026-09-02: p10 landed at 23 nodes on this run, below that line, where p20 lands at 37. The extra retention loss that costs falls unevenly, which is the trade being made: measured over the 370 recordings that reach the sweep at all, wildtype goes 88% → 74% while KO only goes 95% → 93%, widening the between-group spread from 7 to 19 points. That sits just under the runner's own 20-point warning, so p20 does not trip it on this dataset — but it is close enough that a run with more small wildtype recordings would, and the selection check now reports what the tighter cohort does to every gap either way.

Verify: uv run python python/test_density_sweep.py [<extracted run folder>] — 68 checks, of which sections A8–A10 cover the gap tables, the selection check and the decomposition (including that Selection + Resolution == Total, and that C really is A read over B's recordings rather than a second measurement). Section A works on ring lattices and random graphs, where the answer is known by construction: the sweep must separate a lattice from a random graph at the same density, and must give identical results for two graphs that differ only in edge strength.

Where the figures show up (report + viewer, 2026-08-30)

Step 5's figures reach the same two places every other step's do, through one mechanism worth understanding before changing anything here.

One catalogue, three consumers. stats/figures.py owns the list of which figures exist, what each is called, what its caption says, and how to draw it. The step (stats/run.py), the bundle exporter (pipeline/export.py) and the viewer (pipeline/render.py → viewer/server.py) all go through it. Adding a figure means adding it there and nowhere else; adding one to stats/plots.py alone means the viewer will never offer it. test_stats_report.py asserts the catalogue and the files actually written are the same set, in both directions.

Bundles carry the tables, not the pictures. 5_StatsAndML is the first entry in bundle._DATA_ONLY_DIRS: its CSVs travel, its PNGs do not (~100 kB against ~4 MB on the Yin run). That is only safe because every figure is a pure function of the CSVs beside it, so figures.load_results() reads them back and draw_stats_figure() reproduces the figure byte-for-byte — checked in test_stats_report.py, and the check that must keep passing if a plot ever starts reading something the tables do not carry. Two analyses are recomputed rather than stored (the correlation structure and the discriminant projection): both are deterministic and sub-second, and storing them would mean inventing a file format for a stack of matrices and a set of eigenvalues.

Re-bundling a folder that has since been through step 5 tops up the manifest's reconstructable list with 5_stats (bundle._claim_stats), because that list is otherwise static and was written before the folder had any statistics in it.

The report picks the figures up by walking the folder. Its captions come from figures.report_patterns() — the same prose the viewer shows — so the two cannot drift; the folder description and the ~20 CSV descriptions are in report.py alongside the other steps'. Filenames are numbered (5A1_, 5A2_, … 5D3_) because both the report and a file browser sort alphabetically, and without the digits the decoding section opens on its per-age breakdown and the comparisons section on DIV 14 rather than on the model over all of them.

The viewer gets a Statistics tab, hidden when the run has no 5_StatsAndML. It groups the figures by analysis, shows each one's caption underneath, and honours the group/age colour overrides (/api/stats?lag=…&key=…), so a genotype is the same colour there as on every other figure in the session. Figures are cached per address like every other render.

Gotcha: a k-class problem has k-1 discriminant axes. A two-group study — the commonest design — gets one, so plot_lda_projection draws a one-axis strip rather than the scatter. It used to return None there, which meant the catalogue offered a figure that could not be drawn; the test that caught it is the catalogue-vs-files check above.

Cell types: the X-linked mosaic (stats/celltypes.py, 2026-10-14)

The first step-5 analysis whose observation is a cell, not a recording. A heterozygous culture of an X-linked gene (Mecp2) is a mosaic, and staining says which cells express it, so marker- and marker+ cells share one network. Four classes — WT, Het+, Het-, KO — give four contrasts from one mixed model per feature: cell-autonomous (Het- vs Het+), environment on a + cell (Het+ vs WT), environment on a - cell (KO vs Het-), genotype (KO vs WT), plus the change in the cell-autonomous gap per week. Also a WT-vs-KO "KO-likeness" classifier applied to Het cells, a Het-/Het+ decoder on within-recording-normalised features, and per-class trajectories. Figures 5G1-5G8, CSVs mosaic_*.csv and celltype_inventory.csv, in each lag folder; on by default (StatsSettings.celltype_mosaic, --no-mosaic), skipped when the run has no 2_NeuronalActivity/CellTypeCalls.csv or no group mixing + and - cells.

Decisions worth keeping:

• Groups are found from the calls, not their names — mosaic = 10-90% +, references = >=95% + / <=5% +. In the references the marker follows from the genotype, so all their cells are used, labelled or not: on the Yin tracked run only 3 WT cultures have Mecp2 calls, against 37 with cells.
• Two versions, as the user asked: all (every labelled cell-day, recording nested in culture) and tracked (cells from TrackedCellMetrics.csv, with a cell random effect too). Labels always come from CellTypeCalls.csv; tracking supplies only identity (chain#cluster).
• Features are log-scaled when non-negative with skew > 1, then standardised per version, so every estimate is in pooled-SD units. Network node metrics are not normalised within recording for the effects — that would force the environment contrasts to zero by construction — so between-culture contrasts of topology inherit the density/size confound of 5E. The decoder, which asks only about within-recording position, does normalise.
• BLAS is pinned to one thread for the whole analysis (threadpoolctl). With the venv's OpenBLAS each small mixed-model fit spread across 16 cores ran ~7x slower (20 s vs 3 s per tracked feature), and the permutation workers oversubscribed 16x16; the step went from ~4 minutes to most of an hour.

First result on local/yin-tracked-final-oct-8.meanap: 17-18/20 features differ KO vs WT, but the Het- vs Het+ difference is 0 for every network feature and the event rate; only event amplitude / area / duration show small cell-autonomous effects (-0.1 SD, all version only). In activity the Het cultures sit with WT (the KO-vs-Het- environment contrast carries it); in topology both Het classes sit with KO (Het+ vs WT carries it). Marker state is barely decodable within a recording (0.53 balanced accuracy). Tests: python/test_stats_mosaic.py (planted cell-autonomous vs network effects plus a negative control), and the 5G figures go through test_stats_report.py.

Gotchas

• statsmodels' default lbfgs optimiser fails on these mixed models. It drives the culture variance to the boundary and then raises LinAlgError inverting a singular Hessian — on every metric of the Yin run. _MIXEDLM_METHODS tries powell first, which converges on the same data to a variance component of 0.08-0.35. Don't reorder that list without re-checking on real data; the failure is silent in the sense that it looks like "this metric has no effect".
• There is a clustered-OLS fallback (_fit_cluster_ols) for when no optimiser converges. Rows say which was used in the Test column (LME vs OLS-clustered), and the interaction test is a joint Wald rather than a likelihood ratio precisely so it stays valid under the fallback.
• feature_matrix(dropna="rows") is the default and is what MATLAB does. On the Yin run dropna="columns" leaves 2 columns — the calcium imaging activity metrics are missing for enough recordings that requiring complete columns destroys the feature set. Don't switch the default.
• The decoding null is the runtime. 200 permutations × 6 models ≈ 2 minutes on 16 cores; it is parallelised with joblib over permutations, with estimator n_jobs forced to 1 inside the workers (a forest asking for every core from inside every worker is slower than no parallelism).
• CAT-NAP runs name their columns Grp/DIV; MATLAB's doStats expects eGrp/AgeDiv. dataset.META_COLUMNS knows both.

Hundreds of sorted units — 7x end to end (2026-10-24)

Spike-sorted recordings make every node-count-dependent cost bite. Measured on the MaxOne retinal organoid in local/test-spike-sorted-data/ (306 units, 919k spikes, 126 s; steps 1-4, lags 10/15/25 ms, 200 shuffles, full figures, one recording, 16-core box):

	Before	After
Step 3	1494s	121s
Step 4	1749s	332s
Total	3258s	463s

(The "before" run was under cProfile; its step 3 is ~2.2s x 603 STTC matrices unprofiled, so the comparison holds.) Adjacency matrices and both network CSVs are identical before and after; the NMF arrays agree to 1e-14 (see below).

Where it went, and what changed:

Cost (before)	Change	After
STTC pair sweep, 1489s	get_sttc's pair loop is one numba call, rows across cores (prange, rows paired from both ends to balance). The work is O(units x total spikes) per matrix, so compiling alone gained nothing; the cores did.	2.2s -> 0.22s per matrix
NMF, 933s	fits run on threads (sklearn's CD update releases the GIL), each with a seed drawn up front in the old order; the per-rank sweep stops at 95% variance explained or a plateau (nmf_plateau_gain, default 0.005 per 10 ranks; 0 = every rank, as MATLAB). Here it stops at rank 228 of 306.	26s
participation_coef_norm, 379s	its 100 null models are built on threads. Each model's random numbers depend only on the original network's edge count, so they are all drawn first, in the sequential order (_draw_null_plan), then the rewiring (randmio_und_signed, now nogil) runs in parallel. Same models and same generator state as the loop — this is what the 2026-08-26 note above said could not be done without changing results.	33s
randmio_und_v2/latmio_und_v2, <1s	compiled, fed batches from the caller's generator so the draws keep their order; identical output. Negligible here (they were cheap at 306 nodes), done for completeness.	~0s

What is left in step 4 is figures (~250s, of which plot_network is ~190s for 72 network plots at 306 nodes); its computation is 86s. Express mode skips them.

Results. Everything above is exact except one thing: NMF now holds BLAS to one thread for every fit (threadpoolctl). Before, it ran with whatever BLAS share the machine gave it, and a multi-threaded BLAS sums in a thread-count-dependent order, so NMF already differed between machines and worker counts in the last bits. It no longer does; against an old run its arrays move by ~1e-14. num_nnmf_components did not move here, but on another recording such noise could flip it by one, as it could have between machines before.

Threads. parallel.compute_threads() gives a process its share: a map_recordings worker's pinned BLAS share, otherwise every physical core. STTC sets numba's thread count from it; NMF and the null models size their thread pools from it. If numba cannot load a threading layer (a frozen build missing one), get_sttc warns once and runs the same kernel serially; the app smoke test now exercises the parallel kernel, and the PyInstaller spec collects numba.np.ufunc.

Checked by python/test_pipeline_fast_paths.py [3]-[6] (each rewrite against the form it replaced, at one thread and several, including generator state afterwards), plus the step-3/step-4 MATLAB parity tests, CAT-NAP, sorted import and stim-network tests.

Spike detection — 2.5x on step 1 (2026-10-25, 1.15.1)

Step 1 on the ExampleData benchmark (2 recordings, 64 channels, 10 min at 12.5 kHz) went from 39.2s to 16.0s; the whole full run from 105.3s to 80.8s and express from 53.2s to 29.3s. Every output file, spike-check data included, is byte-identical to the run before.

Per recording it had been: 6.3s loading, 9.6s detecting (57s on one thread, so 15 threads bought only 6x), 2.7s summarising the check figures.

Change	Where	Effect
Inflate the v7.3 .mat chunks on threads	io._read_deflate_parallel	6.3s -> 1.0s per recording; identical array
CWT by direct convolution instead of oaconvolve	spike_detection._cwt_bior15	detection 8.8s -> 5.8s per recording
Thresholding keeps c > DTh instead of a scales x samples matrix	detect_spikes_wavelet	~60 ms per channel; identical
Spike counts per bin by bincount	plotting._frequency_curves	1.6s -> ~0 per recording; identical
Pre-draw the example-trace panels, filter their channels on threads	plotting.compute_spike_check_data	~0.7s per recording; identical draws

Loading. MATLAB writes v7.3 variables in tiny gzip chunks (64 x 127 samples: ~59,000 for one ten-minute recording) and h5py inflates them one at a time. zlib releases the GIL, so the raw chunk bytes are read through a file handle per thread and inflated in parallel, straight into the float32 array. That also skips the float64 copy the plain read made (~1.9 GB less peak RAM per recording). Only deflate-only, fully written chunked data takes this path, and only with 3+ threads (below that h5py's own inflate is faster); anything else reads as before.

CWT. The bior1.5 filters are 18-63 taps at 12.5 kHz. At that length overlap-add FFT convolution loses to a direct pass: 1.3x on one thread, and 2.5x with every core busy, since the FFT version's block spectra saturate memory bandwidth (the reason detection scaled so poorly). This is the one change that is not bit-exact: coefficients move by ~1e-16 relative, the same order by which oaconvolve already differed from MATLAB's conv. Spike times on all 128 channels of both recordings are identical, and the MATLAB step-1 parity report is unchanged line for line. Across the seven further recordings below, the only channel it changes is one whose noise floor is itself at rounding level.

More recordings. Step 1 alone, before vs after, on seven more files:

Recording	Rate, channels, format	Before	After	Outputs
HP_tc043_DIV21 (hippocampal)	25 kHz, 60, v7.3	42.2s	14.6s	identical
HCNT26_DIV58_E2	12.5 kHz, 16, v7.3	8.2s	4.7s	identical
20230208_P1_DIV21(000)_A1 (paper example set)	12.5 kHz, 64, v7.3	19.7s	8.1s	identical
MCS N103262_..._div4	25 kHz, 60, v7.3	25.0s	8.0s	1 channel differs (below)
MCS N103262_..._div14	25 kHz, 60, v7.3	26.0s	10.3s	identical
MCS N102512_..._div6	10 kHz, 60, v7	13.5s	10.5s	identical
MCS N102512_..._div14	10 kHz, 60, v7	13.6s	10.8s	identical

The v7 files gain least: scipy inflates a v7 variable as one zlib stream, which cannot be split, so their loading is unchanged.

The one difference is channel 32 of div4, a channel that sits at exactly one value for 94% of the recording. Its filtered trace is zero up to rounding there, so the wavelet's robust noise estimate is ~1e-19 and the detection threshold is rounding error: ~450,000 samples per scale clear it, and which ones do depends on the last bit of each coefficient. Its wavelet spikes went 14,315 -> 14,311, and only the frequency curve in the check data moved with them. Every channel with a real noise floor is identical. The detections on such a channel are rounding noise in either version (and in MATLAB's), which Sigmaj == 0 was meant to catch but only does for an exactly flat channel. Treating a sigma that small as zero would be a behaviour change and is not part of this work.

Left: detection is now ~5.8s per recording, of which per channel the CWT is ~300 ms, trace_noise's two medians ~120 ms and filtfilt ~85 ms.

Checked by python/test_pipeline_fast_paths.py [7]-[9], plus test_pipeline_step1 (MATLAB parity), the spike-detection scale/noise tests, test_raw_io_mcs, test_check_figures_bundle, test_bundle_render, step 2 and the stim tests.

How to verify changes

uv run python python/test_stats.py          # 47 checks, ~2 min — the analyses
uv run python python/test_stats_report.py   # 31 checks, ~1 min — report + viewer
uv run python python/test_stats_mosaic.py   # 22 checks, ~1 min — cell-type mosaic
uv run meanap-stats local/YinThesisRun.meanap --quick -o /tmp/statscheck

Section B of test_stats.py runs the real Yin bundle end to end and skips when local/ is absent. test_stats_report.py builds its own synthetic run, so it needs nothing outside the repo.

Plot parity (MATLAB vs Python)

Goal: the plots each pipeline writes should match in content (values + which plots exist) and be at least as good stylistically — not pixel-identical.

Tool: python/plot_parity_report.py. Walks a MATLAB output tree and a Python output tree, pairs the plots (handling the two pipelines' different file/folder naming via an embedded MATLAB-label→Python-code map + a Lag10ms↔10mslag alias), and writes a self-contained side-by-side HTML report. Re-run after any plotting change:

python python/plot_parity_report.py --matlab OutputData24Dec2025 \
    --python python/parityRun/OutputData_FreshCheck --out plot_parity_report.html

Keep the report at the repo root — it references images by path relative to itself. python/parityRun/run_step1_freshcheck.py + run_steps34_freshcheck.py regenerate a matched Python run (unit=V, lags 10/25/50, recording_workers=1 to avoid the step-3 process-pool crash) for the Python side.

As of 2026-07-16 the Test-pipeline dataset reports 477 paired / 21 MATLAB-only / 12 Python-only (was 378/120/38 before this pass). Fixes landed:

• Half-violin group plots (plot_half_violin_by_x, port of plotHalfViolinByX.m/HalfViolinPlot.m) replaced generic seaborn violins in step-2 and step-4 group comparisons: one-sided KDE fill + jittered scatter + black mean/SEM, subplot-per-group (x=DIV, viridis) / -per-DIV (x=group). Empty metrics (e.g. burst metrics on a no-burst recording) now emit a placeholder axes with a "N/M recordings had no data" note instead of skipping the file.
• Edge-thresholding checks (plotting_step3.py, port of significance_distribution_plots.m) — threshold snapshots come from probabilistic_threshold.adjm_thr(..., collect_check_snapshots=True).
• GraphMetricsByLag + NodeCartographyByLag + DensityLandscape (plotting_step4.py) — the three 4B group plots Python was missing.
• Node cartography scatter now draws the data-driven k-means boundaries (stored on each recording's metrics as cartographyBoundaries) + the NodeCartographyDiagram.jpg bottom panel (bundled in pipeline/assets/), instead of the fixed params thresholds + a side legend.
• Waveforms: align_peaks now uses a separate waveform_width=25 (51-sample waveforms, matching MATLAB) instead of the peak-search win; added a time scale bar. Spike times/counts unchanged (waveforms are plotting-only).
• 4A combined_ plot naming now matches MATLAB (<n>_combined_MEA_NetworkPlot[_<colour name>]).

Remaining MATLAB-only (all in 4A, genuine missing types, not naming): 8_..NullModels (plotNullModelIterations.m) and the two 9_..NodeCartography_modules module-reordered cartography network plots. Boundary/role values differ from MATLAB where the upstream metric is RNG-dependent (step-3 thresholding, k-means) — expected, see the parity table above.

Key files

• src/meanap/stats/ — the step-5 package: dataset.py (loads a run into one feature table, derives culture identity), comparisons.py, correlation.py, decoding.py, regression.py, plots.py, run.py (writes 5_StatsAndML/), cli.py (meanap-stats). See "Step 5" above.
• src/meanap/pipeline/io.py — HDF5/v7.3 .mat I/O: load_raw_recording, load_spike_times_mat, save_spike_times_npz, load_spike_times_npz.
• src/meanap/pipeline/spike_detection.py — the ported detection algorithms (threshold + bior1.5 wavelet CWT). ~540 lines, see "Spike detection gotchas" below before touching this.
• src/meanap/pipeline/output_folders.py — create_output_folders(), a literal port of CreateOutputFolders.m's folder list.
• src/meanap/pipeline/spreadsheet.py — read_recording_csv() / parse_spreadsheet_range(), port of pipelineReadCSV.m. Also parse_ground_electrodes()/ground_spike_times_dict(), port of groundSpikeTimes.m — zeroes out spike times for channels listed in a recording's optional 4th spreadsheet column (Ground, comma-separated channel IDs — RecordingInfo.ground was already parsed from this column, but nothing consumed it until the 2026-07-08 audit found it unwired). Called independently in step2.py/step3.py/step4.py right after each loads its own spike times — MATLAB instead grounds once (in Step 2's formatSpikeTimes.m) and lets the grounded spike times propagate through its single chained .mat file; this port's steps each load spike times fresh from Step 1's raw output, so grounding has to be reapplied at each load site to reach the same outcome. Only matters for spreadsheets that actually have a populated Ground column — the bundled example dataset's exampleData.csv doesn't, so this is a no-op for the Test Pipeline button and wasn't (and can't be, without a fixture that has grounded channels) validated against a real end-to-end run; validated instead via a direct unit test of ground_spike_times_dict() (matches by channel ID against the channels array, zeroes the matched channels' spike-time arrays, leaves everything else untouched).
• src/meanap/pipeline/example_data.py — download_example_data(), port of downloadExampleData.m (Dropbox source only, not the Harvard Dataverse path).
• src/meanap/pipeline/runner.py — run_pipeline(params, log), the new top-level orchestrator that MEApipeline.m corresponds to. Creates the output folders, reads the recording spreadsheet, then runs whichever steps fall in [params.start_analysis_step, params.stop_analysis_step]. All four steps are implemented (each guarded independently by start <= N <= stop, so e.g. start=3, stop=4 correctly skips 1-2 and reads 1-2's saved outputs from disk instead — this was a latent bug for step 2 before this session, which ran whenever stop >= 2 regardless of start). Also implements Params.time_processes (port of MEApipeline.m's Params.timeProcesses — the field and its GUI checkbox already existed in params.py/pipeline.py, but nothing actually timed anything until this session, same "wired but never built" pattern as several other findings in this doc). When set, times each step with time.perf_counter(), logs "Step N duration (seconds): X" per step in MATLAB's own format plus a "Total pipeline duration (seconds): X" line MATLAB has no equivalent of, and writes step_durations.json into the output folder (MATLAB has no single chained output file for this port to append timings to, so this is the durable, programmatically-readable record instead of scraping the log) — built for a MATLAB-vs-Python speed comparison on the Test Pipeline dataset.
• src/meanap/pipeline/step2.py — _run_step2_neuronal_activity(), called from runner.py when stop_analysis_step >= 2. Loads each recording's Step 1 .npz spike times (filtered to params.spikes_method), peeks the raw .mat file's HDF5 shape to get duration_s, calls firing_rates_bursts(), writes per-recording check plots, and dumps every recording's ephys dict to 2_NeuronalActivity/ephys_results.json. (step2_runner_snippet.py alongside it is a stale draft of the same function — safe to delete, superseded by step2.py.)
• src/meanap/pipeline/firing_rates.py — firing_rates_bursts(), the port of firingRatesBursts.m: per-channel firing rates + active-electrode filtering, network burst stats (rate, mean length, ISI within/outside, CV of inter-burst interval, fraction of spikes in bursts), and single-channel burst stats. Returns the ephys dict with MATLAB-matching field names (FRmean, NBurstRate, channelBurstDur, etc.) — see python/test_pipeline_step2.py for the full field list.
• src/meanap/pipeline/burst_detection.py — get_isin_threshold() (Bakkum's automatic ISI_N threshold via histogram peak-finding), burst_detect_isin() (the core ISI_N burst assignment), burst_detect_network() (all-channel network bursts), single_channel_burst_detection() (per-channel bursts). Port of singleChannelBurstDetection.m + the network-burst half of firingRatesBursts.m. See "Burst detection gotchas" below.
• src/meanap/pipeline/plotting_step2.py, parula.py — Step 2 check plots (firing rate heatmap, raster, burst heatmaps); parula.py is a manual port of MATLAB's parula colormap (not in matplotlib) used by the heatmaps.
• src/meanap/pipeline/sttc.py — get_sttc(), the deterministic Spike Time Tiling Coefficient computation, port of get_sttc.m + sttc_m.m (Cutts & Eglen 2014). Vectorized with np.searchsorted rather than a literal two-pointer transliteration — see "STTC gotchas" below for why that's still exactly equivalent. Bit-exact parity (~1e-15 max diff, floating-point noise) confirmed against MATLAB's own get_sttc.m (called directly via python/test_fixtures/gen_sttc_reference.m) on both example recordings at 3 lag values.
• src/meanap/pipeline/probabilistic_threshold.py — adjm_thr() / circular_shift_spikes(), port of adjM_thr_parallel.m's significance thresholding (circular-shift surrogates + upper-tail cutoff). Not bit-reproducible against MATLAB — see the module's own docstring and "STTC gotchas" below.
• src/meanap/pipeline/step3.py — _run_step3_functional_connectivity(). Loads each recording's Step 1 spike times, computes adjm_thr() per lag in params.func_con_lag_val, and saves both the thresholded (adjM{lag}mslag) and raw (adjM{lag}mslag_raw) matrices to ExperimentMatFiles/<recording>_adjM.npz.
• src/meanap/pipeline/network_metrics.py — deterministic BCT-equivalent functions (find_node_deg_edge_weight, strengths_und, density_und, clustering_coef_wu, distance_wei, charpath, efficiency_wei_global, efficiency_wei_local, betweenness_wei), hand-ported from Functions/2019_03_03_BCT/*.m (not via networkx — see "Network metrics gotchas" below for why), plus the modularity-dependent functions: participation_coef (raw, Guimera-Amaral, deterministic given Ci), participation_coef_norm (this is what MATLAB actually saves as NetMet.PC — see "Network metrics gotchas"; additionally stochastic, 100 null-model randomizations via null_models.py on top of Ci), module_degree_zscore, rich_club_wu, classify_node_cartography (6-role classification from fixed Params thresholds), hub_classification (Hub3/Hub4), average_controllability/ modal_controllability (Bassett Lab ave_control.m/modal_control.m), effective_rank, and small_worldness_rl_wu (port of small_worldness_RL_wu.m — deterministic given a real network plus two null models built from it; see null_models.py for those). Read this module's docstring before adding more metrics; it explains exactly which NetMet fields are and aren't in scope and why.
• src/meanap/pipeline/louvain.py — community_louvain(), single-run Louvain modularity optimization, hand-ported from community_louvain.m line-by-line. Sanity-checked against known theoretical values (two disconnected triangles → Q=0.5 exactly), not against MATLAB directly (inherently stochastic — see the module's docstring). Gotcha already found and fixed once: module-label renumbering (unique(Mb)) happens exactly once per hierarchical level, after the local-moving phase fully converges — not after every single node-sweep. Getting this wrong doesn't crash, it just silently produces a worse/wrong partition.
• src/meanap/pipeline/modularity.py — mod_consensus_cluster_iterate(), port of the same-named .m file: runs louvain.py 50× per round, builds a consensus co-classification matrix, re-clusters on the thresholded consensus matrix, repeats until block-diagonal (stable). ~0.3s and 2 outer iterations on the real 64-node example data — fast enough that this isn't a performance concern the way Step 3's 200-shuffle thresholding is.
• src/meanap/pipeline/null_models.py — randmio_und_signed() (Maslov & Sneppen 2002 degree-preserving double-edge-swap rewiring, port of randmio_und_signed.m) and null_model_und_sign() (degree-exact, strength-approximate weighted randomization, port of null_model_und_sign.m's wei_freq<1 periodic-resort branch — the default in any modern MATLAB this codebase targets). Used by participation_coef_norm below. Performance: rng.choice(n, size=4, replace=False) profiled at >85% of randmio_und_signed's runtime for realistic network sizes (array-setup overhead dominates for such a tiny sample) — replaced with a manual batched-candidate-and-filter approach for a ~3.5x speedup; see the inline comment if this ever needs revisiting. Validated via structural invariants only (python/ test_pipeline_null_models.py) — exact degree preservation and total weight preservation are mathematical guarantees of the algorithm (verified empirically), strength preservation is approximate by design (MATLAB's own docs say so), and none of it is bit-reproducible against MATLAB (stochastic — different RNG stream). Also randmio_und_v2()/latmio_und_v2() (ports of the same-named .m files in Functions/CC_PL_SW/) — a different edge-swap algorithm (picks two existing edges by index rather than four random node indices) used to build the random/lattice-like null models that network_metrics.small_worldness_rl_wu normalizes small-worldness against. latmio_und_v2 additionally takes a "distance" matrix (MEA-NAP passes squareform(pdist(adjM)) — Euclidean distance between each node's connectivity profile, not spatial electrode distance) that biases which swaps count as more lattice-like. Both batch their random draws the same way randmio_und_signed does, for the same performance reason. At realistic recording sizes (n≤64) and MATLAB's own iteration counts (10000 for the lattice model, 5000 for the random model), both run in well under a second — sparse active-node subnetworks mean few edges to rewire per attempt. Validated via the same structural-invariants approach, plus (unlike the signed variants above) an exact MATLAB parity check on the deterministic formula that consumes their output — see small_worldness_rl_wu below.
• src/meanap/pipeline/nmf.py — cal_nmf(), port of calNMF.m (num_nnmf_components/nComponentsRelNS/nnmf_residuals/ nnmf_var_explained). Read this module's docstring before touching it — of everything in this port, this is the one place where even the algorithm differs from MATLAB (sklearn.decomposition.NMF's coordinate- descent solver standing in for MATLAB's built-in nnmf, which defaults to Alternating Least Squares), so num_nnmf_components itself — not just the underlying factor matrices — can legitimately differ between the two. Also deliberately diverges from calNMF.m's literal implementation for tractability: bins spike times directly into the target-downsampled- resolution bins instead of building MATLAB's huge intermediate matrix at native sampling rate first — mathematically identical whenever MATLAB's own downSampleSum wouldn't have errored (see the module docstring for the exact-divisibility condition this relies on). Performance: the nnmf_residuals/nnmf_var_explained sweep (one NNMF fit per rank from 1 to the active-electrode count) is unconditional — always runs regardless of Params.includeNMFcomponents — and is the dominant cost, ~55s for a real 64-channel/600s recording on this environment's CPU. Computed once per recording (not once per lag).
• src/meanap/pipeline/step4.py — _run_step4_network_metrics() / compute_network_metrics(). Reproduces ExtractNetMet.m's active-node subsetting (nodeStrength != 0 & activityLevel >= minActivityLevel) then calls the network_metrics.py functions on the subset, and (if params is passed) the modularity-dependent metrics too — Ci/Q/nMod, PC/Z/RC/NdCartDiv/NCpn1-6/Hub3/Hub4. Reads Step 3's thresholded adjacency matrices and Step 1's spike counts; writes 4_NetworkActivity/ netmet_results.json and, per recording/lag, check plots under 4A_IndividualNetworkAnalysis/<group>/<recording>/<lag>mslag/ (see plotting_step4.py). compute_network_metrics()'s result dict includes an adjMsub key (the active-node-subsetted adjacency matrix) purely so plots can be drawn — it's stripped out before the JSON dump, don't expect it to round-trip. Field-naming gotcha: this module computes clustering coefficient and path length twice, for two genuinely different quantities that happen to share MATLAB's CC/PL names in ExtractNetMet.m's docstring — the raw, unnormalized values (CC_raw per-node array + CC_rawMean scalar, PL_raw scalar; independently useful/testable, but NOT what MATLAB saves), and the null-model-normalized scalars from small_worldness_rl_wu (CC/PL — what MATLAB actually saves into NetMet.CC/NetMet.PL, alongside the new SW/SWw). The small-worldness block additionally uses MATLAB's own oddball gate, aN > minNumberOfNodesToCalNetMet (strictly greater — every other block in ExtractNetMet.m uses >=), faithfully replicated rather than "fixed", since the goal is matching MATLAB's actual behavior. Also computes the record-level summary-stat scalars NDmean/NDtop25/ NSmean/sigEdgesMean/sigEdgesTop10 (from ND/NS/the raw adjacency values) and, inside the modularity block, PCmean/PCmeanTop10/ PCmeanBottom10/percentZscoreGreaterThanZero/ percentZscoreLessThanZero (from the normalized PC/Z) — added in the 2026-07-08 audit (see "Auditing for silently-missing metrics" below), these were previously computed nowhere despite being in plotting_step4.py's display-name dict and MATLAB's default netMetToCal. _run_step4_network_metrics() additionally calls nmf.cal_nmf() once per recording (lag-independent, matching ExtractNetMet.m's if e == 1 gate) and merges its result into every lag's metrics dict via metrics.update(nmf_result). The NMF result's rank-indexed arrays (nnmf_residuals, nnmf_var_explained, randResidualPerComponent, and the optional factor matrices) are NOT node-indexed — _NMF_NON_NODE_KEYS explicitly excludes them from the CSV-export loop's generic "any array with size>1 gets spread across NetworkActivity_NodeLevel.csv rows by channel index" logic, since their length can coincidentally match the active-node count and would otherwise get silently mis-attributed to channels.
• src/meanap/pipeline/plotting_step4.py — two plots:
  • plot_connectivity_stats(), port of plotConnectivityProperties.m (adjacency-matrix heatmap + max/mean STTC bars + ND/NS/edge-weight histograms, laid out via GridSpec to mirror MATLAB's 6x6 tiledlayout). Visually cross-checked against MATLAB's own 1_adjM10msConnectivityStats.png from OutputData24Dec2025 — same layout, same shapes; the specific bin heights/edges differ because that comparison used a different MATLAB run's stochastically-thresholded adjM (expected, see Step 3 notes), not a plotting bug.
  • plot_spatial_network(), port of the base 2_MEA_NetworkPlot.png + 3_MEA_NetworkPlotNodedegree<Metric>.png from StandardisedNetworkPlot.m — reuses network_plot.py's plot_network() (built for the Network Viewer GUI tab, already a generic MEA renderer) rather than duplicating it. Needs electrode coordinates — see channel_layout.py below. Visually cross-checked against MATLAB's 2_MEA_NetworkPlot.png from OutputData24Dec2025: node positions and overall edge topology (which nodes are hubs, the shape of the long diagonal high-weight edges) matches closely; exact edge set differs only because of Step 3's inherent stochastic thresholding (different MATLAB run).
  • Same plot_spatial_network() also produces 4_MEA_NetworkPlotNodedegreeParticipationcoefficient.png, colored by the normalized PC (metrics["PC"], i.e. participation_coef_norm's 1st output — matches what MATLAB actually feeds this specific plot, traced through PlotIndvNetMet.m line 119).
  • plot_node_cartography(), port of the top (data) panel of NodeCartography.m — PC vs. Z scatter colored by the 6 cartography roles, with the 5 fixed Params boundary lines. Also colored by the normalized PC now (was raw PC in an earlier version of this port — see "Network metrics gotchas" for why that was a real bug, not just an approximation choice). The bottom panel in MATLAB (a static NodeCartographyDiagram.jpg explaining the roles) is not reproduced — nothing data-driven to check parity against. Visually cross-checked against MATLAB's 9_adjM10msNodeCartography.png from OutputData24Dec2025 — same legend/colors/boundary-line mechanics; the specific point positions differ only because Ci and the PC normalization are both stochastic (expected, documented).
  • Also produces 5_MEA_NetworkPlotNodestrengthLocalefficiency.png (colored by Eloc). Node size here is node strength (NS), not node degree — the one plot in this family that sizes differently, confirmed from PlotIndvNetMet.m directly (nodeSizeMetricsToPlot{end+1} = lagNetMet.NS specifically for the Eloc entry, everything else uses lagNetMet.ND) rather than trusting the (wrong, "node degree") wording in docs/meanap-outputs.rst's own Figure 5 caption — the filename MATLAB itself generates (Nodestrength...) was the tell. Finding this surfaced a real bug in network_plot.py's plot_network(): its legend hardcoded the label "node degree:" and formatted values as 2-digit integers regardless of what z actually was, and applied a max(z, 3) floor intended to keep integer-ND legend divisions sane — which silently shrank every node to a sliver when z was node strength (typically < 1). Fixed by adding a z_name parameter that both relabels the legend and switches off the degree-specific floor/integer-formatting whenever z_name != "node degree" (or z doesn't look integer-valued). Defaults preserve the Network Viewer GUI's existing behavior exactly (it always sizes by ND). If you add another plot that sizes by anything other than ND, pass z_name — don't skip it even if the plot "looks fine" at a glance, since a too-small-node bug is easy to miss visually until you compare side-by-side against MATLAB's own legend. Saved per recording/lag: 1_adjM{lag}msConnectivityStats.png, 2_MEA_NetworkPlot.png (flat, no color metric — matches MATLAB's base variant), 3_MEA_NetworkPlotNodedegreeBetweennesscentrality.png (colored by BC), 4_MEA_NetworkPlotNodedegreeParticipationcoefficient.png (colored by normalized PC), 5_MEA_NetworkPlotNodestrengthLocalefficiency.png (colored by Eloc, sized by NS), 6_circular_NetworkPlotNodedegreeModule.png (circular layout, nodes colored by community assignment Ci from mod_consensus_cluster_iterate, sized by ND — see plotting_step4.py's plot_circular_module_network and the notes in this doc's "⚠️ Network metrics gotchas" section for why Ci is stochastic), 9_adjM{lag}msNodeCartography.png, 10_MEA_NetworkPlotNodedegreeAveragecontrollability.png, 11_MEA_NetworkPlotNodedegreeModalcontrollability.png (controllability plots, including their batch-scaled variants — see "Batch scaling" below). Still not ported: null-model panels (small-worldness — see "Network metrics gotchas").
• src/meanap/pipeline/channel_layout.py — get_coords_from_layout(), port of getCoordsFromLayout.m. Supports MCS60old, MCS60, MCS59 (electrode-grid layouts that drop 1-4 grounded/corner channels via a specific hardcoded reordering — see the module for the transcribed channelsOrdering arrays), Axion64, Axion16 (simpler row/col grid formulas, no channel dropping). Exact parity confirmed against MATLAB for all 5 layouts (matlab -batch calling getCoordsFromLayout.m directly — see python/test_fixtures/gen_coords_reference.m if that needs re-verifying). The example dataset uses Axion64 (Parameters_*.csv: channelLayout=Axion64), which is why all 64 raw channels appear in every step 1-4 output with none dropped — MCS-layout corner/ground-channel dropping was never actually exercised against real pipeline data, only validated against MATLAB's lookup table in isolation. MATLAB's 'Custom' layout (user-drawn random coordinates) isn't ported — raises ValueError if requested.
• src/meanap/gui/main_window.py — _on_run runs the pipeline on a background QThread (see item 1 under "Known limitations" for the cancellation mechanism). _on_test_pipeline downloads the example dataset, points the Data tab at it, sets start_step=1, stop_step=4, then calls _on_run — i.e. it runs the full 4-step pipeline. This mirrors runPipelineApp.m's TestPipelineButton handler, which also falls through into a full pipeline run using the example-data settings.
• src/meanap/gui/panels/pipeline.py — has both "Start at step" and "Stop at step" spin boxes now (kept mutually consistent — moving one past the other drags it along). test_btn ("🧪 Test pipeline") sits left of Run/Stop.
• python/test_pipeline_step1.py — standalone parity script; downloads example data via meanap.pipeline.example_data, runs step 1, and diffs spike times against the MATLAB reference in OutputData03Mar2026/.
• python/test_pipeline_step2.py — standalone parity script for step 2. Unlike the step 1 script, it feeds firing_rates_bursts() the MATLAB reference spike times (loaded straight from OutputData03Mar2026's _spikes.mat, method bior1p5) rather than Python's own step 1 output — this isolates step 2 arithmetic from step 1's known ~82-84% wavelet-approximation gap. Diffs the resulting ephys dict against NeuronalActivity_RecordingLevel.csv and NeuronalActivity_NodeLevel.csv (MATLAB's own flattened summary CSVs, conveniently already keyed by the same field names as ephys). Currently 100% parity (1052/1052 checks) on both example recordings.
• python/test_pipeline_step3.py — two-part test. (1) Exact parity of get_sttc() against the python/test_fixtures/*_sttc_reference.npz fixtures (generated via MATLAB directly, see "Key files" above) — currently exact (~1e-15 max diff). (2) Structural sanity checks on adjm_thr()'s thresholding (symmetric, zero diagonal, thresholding only removes edges) — NOT a MATLAB parity check, since the RNG isn't shared.
• python/test_pipeline_step4.py — feeds network_metrics.py the actual MATLAB-thresholded adjacency matrices from a real run (OutputData03Mar2026/ExperimentMatFiles/*_OutputData03Mar2026.mat, adjMs.adjM{lag}mslag) rather than Python's own Step 3 output, isolating Step 4's arithmetic from Step 3's stochasticity. Diffs against python/test_fixtures/*_step4_reference.npz (MATLAB's own BCT functions, called directly — see gen_step4_reference.m). Currently 100% parity (60/60 checks) across both recordings × 3 lags × 10 metrics.
• python/test_pipeline_cartography.py — feeds network_metrics.py's modularity-dependent functions (raw participation coefficient, z-score, rich club, cartography classification, hub classification) a fixed, MATLAB-generated Ci (python/test_fixtures/*_cartography_reference.npz, via gen_cartography_reference.m) — isolating "does the deterministic math match" from "does the same community assignment get found" (a different, inherently non-bit-reproducible question — see modularity.py's docstring). Cartography classification and hub counting are additionally fed a fixed, MATLAB-generated normalized PC (PCnorm in the fixture) since that's what MATLAB's real pipeline actually feeds them, not the raw PC (see "Network metrics gotchas"). Currently 100% parity (42/42 checks).
• python/test_pipeline_null_models.py — structural tests (not MATLAB parity, since these are stochastic) for null_models.py and participation_coef_norm: exact degree preservation, exact/approximate weight and strength preservation, boundedness, no-NaN. Also a real-data timing smoke test (~15-35s for 100 iterations at real recording sizes).
• python/test_pipeline_small_worldness.py — three-part test for small_worldness_rl_wu + randmio_und_v2/latmio_und_v2. (1) Exact MATLAB parity for the deterministic formula assembly: feeds MATLAB's own A/R/L (python/test_fixtures/small_worldness_reference.npz, via gen_small_worldness_reference.m) into small_worldness_rl_wu and diffs SW/SWw/CC/PL against MATLAB's own small_worldness_RL_wu.m output on the same inputs. Unlike the other gen_*_reference.m scripts, this one doesn't need a real MATLAB pipeline run's saved adjacency matrices — it builds a small fixed-seed random network directly, since the point is isolating formula correctness, not re-deriving a specific recording's numbers. (2) Structural invariants on Python's own randmio_und_v2/ latmio_und_v2 (exact degree/weight preservation, symmetry) — NOT a MATLAB parity check, same rationale as test_pipeline_null_models.py. (3) An end-to-end smoke test chaining Python's own null models into small_worldness_rl_wu. All three currently pass.
• python/test_pipeline_nmf.py — structural/sanity tests for nmf.py (NOT a MATLAB parity check — impossible here even in principle, since the underlying NNMF solver is a different algorithm, not just a different RNG stream, see nmf.py's docstring). Checks: spike-count preservation through randomise_spike_train/_bin_spike_times; on a synthetic low-rank-plus-noise recording, num_nnmf_components lands in a sane range, nnmf_var_explained is bounded in [0, 1] and approaches ~1 at full rank, nnmf_residuals roughly decreases with rank; include_nmf_ components=True returns non-negative factor matrices. All currently pass.
• python/test_fixtures/ — .npz ground-truth fixtures for steps 3-4 plus the .m scripts that generated them (gen_sttc_reference.m, gen_step4_reference.m, gen_cartography_reference.m, gen_small_worldness_reference.m) — re-run those directly in MATLAB if the fixtures ever need regenerating (e.g. after an algorithm change).

HTML output report (report.py)

src/meanap/pipeline/report.py's generate_report(output_root) walks a completed (or partial) MEA-NAP output folder and writes a single self-contained report.html at its root — a folder tree on the left, an image gallery with captions on the right. No server, no external JS/CSS, no new dependencies; opens directly over file://. Wired into the GUI as a "🌐 View report" button next to Run/Stop on the Run tab (main_window.py's _on_view_report) — regenerates the report from whatever's on disk and opens it via webbrowser.open().

• Captions are sourced from docs/meanap-outputs.rst (MEA-NAP's own Sphinx figure-legend reference) wherever it documents a figure, reworded to describe what this Python port's plot actually shows — MATLAB's originals often also render "scaled to whole dataset" and "combined" side-by-side variants this port doesn't produce. Six step-2 burst-heatmap figures (3_BurstRate_heatmap.png through 8_BurstDetectionInfo.png) have no MATLAB documentation anywhere in the repo (confirmed via repo-wide search — MEApipeline.m only lists the filenames, and the MATLAB plotting functions' own docstrings are placeholder boilerplate); their captions in report.py are original, written to match the sibling group-comparison figures' documented semantics for the same metric.
• Deep-linkable: report.html#4_NetworkActivity/.../10mslag auto-expands the sidebar tree and selects that folder on load (openHashPath() in the embedded JS). Every folder node's path (relative to the output root) is baked into the JSON tree at generation time and both drives this and lets you verify a filename maps to the caption you expect by reading _PLOT_PATTERNS'/FOLDER_DESCRIPTIONS' regexes directly rather than reverse-engineering the JS.
• Data files (.npz/.json/.csv/.mat) are listed (with a short caption) rather than embedded as images — clicking one just opens/downloads it via the browser's normal file:// handling.
• Verified by rendering with headless Chrome (google-chrome --headless --screenshot, available on this machine) and visually inspecting — no automated test script for this one, since there's no MATLAB ground truth to diff an HTML report against. If the tree/gallery structure changes, re-screenshot a few folders by hand rather than trusting it blind.

Spike sorting (Params.spike_source = "sort", 2026-09-16)

Step 1 can sort spikes into units instead of detecting them per electrode; each unit is then a node through steps 2–4. Everything about it — the survey of sorters, why low-density MEAs are a different problem from Neuropixels, the synthetic-ground-truth benchmark that made Tridesclous2 the default, the file format (channels = parent electrode per unit, unit_* arrays, method sorted), and the follow-ups — is in python/SPIKE_SORTING_PLAN.md. User docs: docs/python/spike-sorting.md. Tests: python/test_spike_sorting.py; benchmark: python/benchmark_spike_sorting.py.

Gotchas that are easy to re-break:

• active_spike_method(params), not params.spikes_method, is what steps 2/3/4 (and stim, render) read: it is sorted in sort mode.
• Never spikeinterface[full] (CUDA torch, and still no hdbscan).
• Preprocessing is saved in-process: a spawned worker gets a pickled copy of an in-memory NumpyRecording.
• Workers are spawned, set through si.set_global_job_kwargs, because the sorters' internal stages do not all take job kwargs and a fork copies the pipeline process.
• Curate on the refractory-violation fraction; SI's isi_violations_ratio is Poisson-normalised and reads 4–60 on bursting cultures.

Spike detection gotchas (don't re-discover these)

• The wavelet detectors never had a refractory period — in MATLAB (detectSpikesCWT.m passes refPeriod only to detectSpikesThreshold) and, until 2026-09-17, in this port. The CWT fires at more than one scale on one spike and alignPeaks then puts both on the same trough. Measured on HP_tc043_DIV21 (dense hippocampal culture): every bior1.5 pair closer than 0.5 ms was one trough detected twice (442/442, 435/435, 370/370 on three busy electrodes; 5–6.5 % of their spikes), while pairs 0.5–1 ms apart were 95 % distinct troughs and everything above 1 ms was 100 % distinct — population spikes, which refPeriod (2 ms) would discard as the threshold detectors do (bior1.5 on electrode 57: 7925 → 7484 at 0.5 ms, → 4468 at 2 ms, below thr4's 5499). Hence Params.wavelet_ref_period_ms, default 0.5 ms, applied after alignment; None restores the old behaviour and is what test_pipeline_step1.py uses for MATLAB parity. Analysis script: local/spikesort_test/real/bior_doubles.py (not committed).

• bior1.5 CWT sign: PyWavelets' wavefun() gives the analysis wavelet with the opposite sign to MATLAB's Wavelet Toolbox CWT. _cwt_bior15() must return -coeffs at the end, or F1 collapses from 0.82 → 0.10. Obsolete (2026-07-15) — that sign flip was compensating for the old FFT-based CWT, which has been replaced by a literal port of MATLAB's own algorithm (see below). The -sqrt(a) factor in MATLAB's formula now carries the sign; there is no separate negation, and PyWavelets' psi_d has the same sign as MATLAB's (they agree to 1.6e-15).

• Threshold naming: MATLAB does strcat('thr', num2str(t)). The Python _thr_name() helper must special-case whole numbers (thr4 not thr4.0) — getting this wrong silently produces F1 = 0.000 because saved/loaded keys never match.

• Threshold sign convention (fixed this session — was a live bug): the detection formula is threshold = median − multiplier × MAD, so multiplier must be positive to push the threshold below the median for negative-going spikes. Params.thresholds used to default to [-3.5, -4.5, -5.5] (copied from an assumption that never got exercised end-to-end); running step 1 with those values would have pushed the threshold the wrong way and detected near-nothing. Fixed default to [3.0, 4.0, 5.0], matching MATLAB's own Params.thresholds = {'3','4','5'} convention. If you see thresholds go negative anywhere again, that's a bug.

• HDF5 loading: all .mat files here are v7.3 (HDF5), so scipy.io.loadmat can't read them — everything goes through h5py. Raw data is stored (n_channels, n_samples) and needs .T. Spike times are stored as a (64, 1) array of object references into per-channel Groups keyed by method name (bior1p5, thr4, thr5), values in seconds.

• Parity numbers: thr4/thr5 F1 = 1.000 (exact match), bior1.5 F1 ≈ 0.82–0.84 (PyWavelets approximates the wavelet via cascade; MATLAB uses its native CWT — expected, not a bug to chase further unless parity requirements tighten). Wrong diagnosis, and now fixed (2026-07-15): bior1.5 F1 = 1.0000, exact. The premise was false in two ways. (1) PyWavelets' bior1.5 is not an approximation of MATLAB's: the filters are byte-identical, wavefun's grid is identical, and psi_d agrees to 1.6e-15. (2) MATLAB's "native CWT" for this call is not the modern Morse-wavelet cwt — cwt(x, scales, 'bior1.5') dispatches to the legacy algorithm, which is itself cascade-based and fully reproducible:

  [psi, xval] = intwave(wname, 10);  step = xval(2)-xval(1);
  j = 1 + floor((0:a*(xmax-xmin)) / (a*step));
  f = fliplr(psi(j));
  coefs = -sqrt(a) * wkeep1(diff(wconv1(signal, f)), len);

  (verified bit-identical to R2024b's cwt, max diff 0). _cwt_bior15 is now a literal port of that — integrated wavelet, convolve, differentiate — and matches MATLAB's cwt to ~1e-15 relative on real traces. Two further bugs were fixed alongside it, both of which had to be right before F1 reached 1.0:

  • _determine_scales compared 0-based zero-crossing indices against MATLAB's >500/<500 bounds, which are 1-based. A crossing at 1-based index 501 fell into neither bin, corrupting the narrowest scale's width entry and yielding scales [3 3 4 6 7] instead of MATLAB's [2 3 4 6 7] — so the port simply never looked at the narrowest spikes (precision ~1.0, recall 0.93: a pure subset of MATLAB's spikes).
  • Sigmaj centred the subsampled coefficients on their own mean; MATLAB centres them on the mean of the full row (mean(c(i,:))). Beware np.round vs MATLAB round here too (half-to-even vs half-away-from-zero) — the scale interpolation lands on exact .5 values.

• Known 2-sample offset (deliberate, not a bug): Python's spike times are a constant 2 samples (0.16 ms) earlier than MATLAB's, for every method. Cause: MATLAB reports 1-based frame indices as spikeFrames/fs (so a spike on the first sample gets 1/fs, not 0), and alignPeaks.m computes newSpikeTime = spikeTimes(i)+pos-win with a 1-based pos, landing on the true peak +1 even for a perfectly centred peak. Python reports the true peak at idx/fs. Decision (2026-07-15, with the user): keep Python correct rather than replicate MATLAB's off-by-one — same call as the setUpSpreadSheet.m coordinate bug above. This affects no downstream metric (a uniform shift leaves firing rates and STTC unchanged); it only means spike times need a +2-sample shift to diff against MATLAB's directly.

• Open, minor: thr4 finds 2 extra spikes per recording vs MATLAB (24943 vs 24941; F1 still 1.0000 to 4 dp). Predates the CWT work and is unrelated to it (threshold detection doesn't touch the CWT) — likely the refractory-period loop. thr5 and bior1p5 are exact.

Burst detection gotchas (don't re-discover these)

• array_fracInBursts is not a median of per-channel fractions (fixed this session — was a live bug). MATLAB computes it as sum(total_num_sp_in_bst(bursting_electrodes)) / sum(sum(raster)) — total spikes-in-bursts summed across only the bursting electrodes, divided by total spikes across all electrodes (singleChannelBurstDetection.m line 152). burst_detection.py's single_channel_burst_detection() used to compute np.nanmedian(all_fracs_in_burst) instead, which returns NaN when zero channels burst — MATLAB returns 0.0 there (0 spikes-in-bursts over a nonzero total). Fixed by accumulating total_sp_in_bst_sum / total_all_spikes directly instead of medianing per-channel ratios. The per-channel all_fracsInBursts (channelFracSpikesInBursts in ephys) was already correct — only the aggregate was wrong.
• Weak test coverage for the actual burst-matching logic: on both example recordings, firing rates are low enough (~0.16-0.20 Hz mean) that MATLAB itself detects zero bursts, network or single-channel, for either spike-detection method. test_pipeline_step2.py therefore only proves Python agrees with MATLAB that there are no bursts here — it does not exercise whether burst_detect_isin() / get_isin_threshold() produce the same burst boundaries as MATLAB when bursts do exist. If a future dataset with genuine bursting activity becomes available, re-run the parity test against it before trusting burst-timing outputs.

STTC gotchas (don't re-discover these)

• _run_p is vectorized via searchsorted, not a literal port of the two-pointer loop in sttc_m.m. The MATLAB loop advances a monotonic pointer j into the second spike train and never backtracks; this is mathematically equivalent to "does the nearest spike in train 2 (by insertion point) fall within dt", since both trains are sorted ascending — so checking the immediate left/right neighbours via np.searchsorted per spike gives identical results without the pointer bookkeeping. Confirmed empirically to bit-exact match MATLAB's C-mex sttc.c output, not just argued from theory — if you ever suspect this equivalence breaks (e.g. for unusual/degenerate spike time inputs), regenerate python/test_fixtures/*_sttc_reference.npz and check before trusting a change here.
• Probabilistic thresholding (adjMci) is fundamentally not bit-reproducible against MATLAB. adjM_thr_parallel.m never seeds MATLAB's RNG, and even two MATLAB runs of the same recording can give different thresholded matrices (worse if the Parallel Computing Toolbox's parfor is active — each worker gets its own substream). Don't spend time trying to match adjMci bit-for-bit; the raw adjM (via get_sttc) is where exact parity lives, and it's what step 4 downstream should prefer once cross-step chaining matters more than matching MATLAB's specific thresholded output for a specific past run.
• Performance: get_sttc() on the full 64-channel example recording takes ~0.07s (vectorized). adjm_thr() with rep_num=200 (MATLAB's default in the reference run) costs ~200x that per lag — budget for ~15-80s per (recording, lag) depending on activity level; step3.py does this for every recording × lag combination in params.func_con_lag_val, so a full run can take minutes, similar to step 1's wavelet CWT.

Network metrics gotchas (don't re-discover these)

• MATLAB's saved NetMet.CC / NetMet.PL are NOT the plain clustering coefficient / path length. They come from small_worldness_RL_wu.m, which computes the real network's C/PL and then normalizes them against randomized (randmio_und_v2, 5000 iterations) and lattice (latmio_und_v2, 10000 iterations) null models: CC = C/Cl, PL(saved) = PL/PLr. Those null models involve randomized edge rewiring (stochastic, like step 3's thresholding) and are NOT ported. This module's CC/PL are the raw, unnormalized values (what small_worldness_RL_wu calls C and PL internally, before the null-model division) — useful and fully deterministic, but a different number than what MATLAB's NetMet.CC / NetMet.PL fields contain. Don't compare them directly without remembering this distinction.
• Hand-ported from BCT, not networkx, even though networkx is a project dependency. Tried & measured: BCT's weighted formulas (especially clustering_coef_wu's Onnela geometric-mean triangle intensity, and efficiency_wei's modified local efficiency) don't have exact networkx equivalents with matching normalization — porting the BCT .m source directly (available at Functions/2019_03_03_BCT/) and validating against real MATLAB output was more reliable than trying to reverse-engineer an equivalence. distance_wei (Dijkstra) and betweenness_wei (Brandes') are vectorized-per-source-node but otherwise structurally identical to their .m counterparts — see the inline comments in network_metrics.py if you need to verify the port line-by-line against the MATLAB source.
• Ground truth requires a fixed adjacency matrix, which MATLAB's real pipeline never persists on its own (only Ephys, Info, Params, spikeTimes, and the thresholded adjMs are saved to ExperimentMatFiles/*.mat — never a NetMet struct for the example dataset, since step 4 was never run on it). The fixtures instead take the real adjMs.adjM{lag}mslag matrices already sitting in OutputData03Mar2026/ExperimentMatFiles/ and run MATLAB's own BCT functions on them directly, replicating ExtractNetMet.m's active-node subsetting logic by hand in gen_step4_reference.m. If step4.py's subsetting logic ever changes, that MATLAB script's subsetting logic needs to change in lockstep or the fixture stops being a valid ground truth.
• Q/Ci/nMod (modularity), PC (raw + normalized)/Z, RC, node cartography (NdCartDiv/NCpn1-6), Hub3/Hub4 are now ported — see louvain.py (single-run Louvain), modularity.py (consensus clustering → Ci/Q/nMod), null_models.py (randmio_und_signed, null_model_und_sign — the randomization participation_coef_norm needs), and network_metrics.py's participation_coef, participation_coef_norm, module_degree_zscore, rich_club_wu, classify_node_cartography, hub_classification. Everything downstream of a fixed Ci (or, for the normalized PC, a fixed Ci and a fixed PC_norm) has 100% parity against MATLAB (python/test_pipeline_cartography.py); the stochastic pieces themselves (Ci, and separately the 100 randomizations inside participation_coef_norm) aren't bit-reproducible — same situation as Step 3's thresholding. participation_coef_norm costs ~15-35s per (recording, lag) at real data sizes (100 randomizations × Maslov-Sneppen rewiring) — see null_models.py's "Performance" note above if this ever needs to be faster.
• Auto-set cartography boundaries (TrialLandscapeDensity.m) is now ported — network_metrics.trial_landscape_density + step4._apply_cartography_boundaries (Phase B reducer, gated on params.auto_set_cartography_boundaries, per-lag by default to match MEApipeline.m). Before this, step4.py classified nodes with the fixed default boundaries (peri_part_coef=0.625 …); against real data almost every node's normalized PC sits below 0.625, so ~62-95% of nodes piled into cartography role 1 (NCpn1) and roles 2-6 were near-empty — a systematic, non-RNG divergence from MATLAB's NCpn1-6 CSV columns (found via end-to-end CSV verification vs OutputData24Dec2025). The port pools PC/Z across recordings and re-derives all five boundaries by optimal 1-D k-means (_optimal_1d_split_boundaries, a deterministic O(k·n²) DP — MATLAB's kmeans is random-seeded and not itself bit-reproducible, so exact boundary parity is unavailable, same class as Step 3). Validated against MATLAB's stored per-lag boundaries (several match to ~1e-16 where MATLAB's kmeans also hit the optimum) and by test_pipeline_landscape.py (48 checks, incl. a regression guard that data-driven boundaries un-pile role 1 across all 6 fixture recordings/lags). Deliberate divergence from MATLAB: TrialLandscapeDensity.m line 61 reloads ExpName{1} every loop iteration, so MATLAB actually pools only the first recording's PC/Z (duplicated) — an upstream bug. The port pools all recordings (the documented intent); this is the largest remaining source of NCpn1-3 difference vs MATLAB and is expected/bounded, not a defect.
• Found and fixed a real bug from an earlier session (not a MATLAB bug — this one was mine): node cartography classification and Hub3/Hub4 were wired to the raw PC instead of the normalized PC in an earlier version of step4.py and the corresponding test fixture (gen_cartography_reference.m). Traced ExtractNetMet.m → PlotIndvNetMet.m line 119 to confirm MATLAB's real pipeline feeds cartography/hub classification the normalized PC (the same one that colors 4_MEA_NetworkPlotNodedegreeParticipationcoefficient.png), not the raw Guimera-Amaral formula. Both step4.py and the test fixture now use the normalized PC for these two — confirmed 100% parity re-holds (python/test_pipeline_cartography.py, 42/42). If a metric name in MEA-NAP's saved output has both a "PC"-like short name and a longer formula-named function producing a different number (raw vs. normalized, as here — or the CC/PL/SW situation above), always trace which one the actual plotting/downstream code consumes before assuming the obvious match.
• Cmcblty (communicability) is correctly NOT ported — it's not a gap. MATLAB's own code path that would compute it (fcn_find_hubs_wu.m, called from a block in ExtractNetMet.m) is commented out (wrapped in %{ %}) in the current MATLAB source. It's dead code upstream too.
• SW/SWw (small-worldness) still NOT ported — needs the randmio_und_v2/latmio_und_v2 null models discussed above. Lower priority than everything else in this list since it was the most stochastic and most computationally expensive (10000/5000 rewiring iterations) piece to begin with.

⚠️ Electrode coordinate bug — found in MATLAB, not just the Python port

Functions/setUpSpreadSheet.m has a real indexing bug that misplaces electrodes in every spatial plot (heatmaps, network plots) for layouts where the raw recording's channel order differs from getCoordsFromLayout.m's own internal generation order — which is the case for Axion64 (and likely the MCS-family layouts too). This was found and confirmed while investigating a user report that the step-4 network plots looked rotated relative to step 2's firing-rate heatmap.

Root cause (setUpSpreadSheet.m, ~line 68):

[channels, coords] = getCoordsFromLayout(Params.channelLayoutPerRecording{nRecording});
recordingChannels = recordingChannelData.channels;   % raw file order, e.g. 11,12,13,...,88
subsetIndex = find(ismember(channels, recordingChannels));
Params.channels{nRecording} = recordingChannels;      % raw file order
Params.coords{nRecording} = coords(subsetIndex, :);   % getCoordsFromLayout's OWN order (11,21,31,...)

find(ismember(channels, recordingChannels)) does not reorder coords to align with recordingChannels— since every channel is present in both arrays, it just returns 1:64 in channels's own (unreordered) order. The result: Params.coords{n}[i] and Params.channels{n}[i] are for two different physical electrodes whenever the two source orderings differ. Every spatial plot that reads Params.coords{ExN} (heatmaps, network plots, node cartography — electrodeHeatMaps.m, StandardisedNetworkPlot.m, and others all consume this same variable with no correction) inherits the misplacement.

getCoordsFromLayout.m's Axion64 formula itself is correct — confirmed against Axion Biosystems' own AxIS Navigator software manual, which states "[t]he electrode ID consists of two digits, the column and row" and that column/row map to the x-/y-axis respectively. That matches channels(i) = colIdx*10 + rowIdx, coords(i,:) = [spacing(colIdx), spacing(rowIdx)] exactly. The bug is entirely in the reindexing step, not the coordinate formula.

Proof: extracted the real Params.coords{1} from a genuine MATLAB run (OutputData24Dec2025/ExperimentMatFiles/*.mat) and reproduced the exact find(ismember(...)) logic in Python — it reproduces MATLAB's actual saved coordinates exactly. Cross-checked against MATLAB's own Ephys.FR (channel 23 = true max firing rate for that recording) against pixel-measured positions in the real 2_Heatmap.png — the bug's predicted (wrong) position matches where the brightest dot actually renders; the semantically-correct position does not.

Decision (2026-07-02, discussed with the user): the Python port (channel_layout.py, plotting_step4.py's plot_spatial_network, plotting_step2.py's plot_heatmap) uses the correct electrode-to- coordinate mapping, not a bug-for-bug replication of MATLAB's current output. This means Python's spatial plots will not pixel-match MATLAB's current Axion64 output, but will be internally self-consistent (heatmap and network plots agree with each other) and consistent with Axion's documented convention. plotting_step2.py previously had its own ad-hoc, undocumented coordinate function (get_coords_from_chs, since removed) that turned out to be the exact transpose of the correct mapping — likely someone's attempt to compensate for a symptom of this same bug by eye, applied to only one of the two plot types, which is what produced the apparent heatmap-vs-network-plot rotation mismatch the user noticed. Both now go through channel_layout.py consistently.

Not fixed: the user asked not to touch setUpSpreadSheet.m for now. If this gets revisited, the fix is to properly reorder coords to match recordingChannels, e.g. [~, subsetIndex] = ismember(recordingChannels, channels); Params.coords{nRecording} = coords(subsetIndex, :); (note ismember's second output gives the actual per-element position, unlike find(ismember(...)) which only gives a same-order subset).

Open question, follow-up needed: the user mentioned their stimulation- data pipeline runs in MATLAB look correctly oriented already. Stimulation plotting functions (electrodeHeatMaps.m et al.) consume the same Params.coords{ExN} — no separate code path was found (confirmed via repo-wide grep for stim.*coord). The likely reconciliation: stim-data raw files may happen to have their channels field stored in an order that already coincides with getCoordsFromLayout's internal generation order, which would make the bug not manifest for those files by coincidence, not because of different code. Worth checking directly against a real stimulation dataset before trusting this hypothesis.

Resuming from a previous run (Params.prior_analysis, 2026-07-31)

MATLAB's Params.priorAnalysis + startAnalysisStep workflow now works in the port. pipeline/resume.py owns it. Semantics, mirroring MATLAB: output still goes to this run's own (fresh, dated) folder — the prior run is only ever read.

Only three artefacts cross step boundaries in this port, so resuming is a search path over a few lookups rather than MATLAB's chained-.mat machinery:

Produced by	Consumed by	File
Step 1	Steps 2, 3, 4	1_SpikeDetection/1A_SpikeDetectedData/<rec>_spikes.npz
Step 3	Step 4	ExperimentMatFiles/<rec>_adjM.npz
CAT-NAP step 2	CAT-NAP step 4	ExperimentMatFiles/<rec>_catnap.npz

InputLocator resolves both against, in order: this run's output folder → Params.spike_detected_data (spike files only, MATLAB's spikeDetectedData) → Params.prior_analysis_path. It's a frozen dataclass of plain paths so it pickles into the spawned Step 3/4 workers; each step rebuilds its own via build_input_locator(params, output_root) rather than threading it through the task tuples.

Two things this needed on top of the lookup itself:

• Fail fast. runner.py calls missing_step_inputs() before any compute and raises if no recording has what the starting step needs (it warns and continues if only some do). Previously a typo'd path just produced "SKIP: spike data not found" per recording and a "successful" run with empty CSVs.
• Duration no longer requires the raw recording. Steps 2/3/4 each used to re-derive duration_s by reopening the raw .mat, which made resuming useless without the raw data mounted. Step 1 now stores duration_s in the spike .npz and everything reads it back via io.resolve_duration_s(), falling back to the raw file only for .npz files written before this change. That also collapsed three divergent implementations into one — step 2's had a hardcoded if n_samples == 64 transposed-check (Axion64-specific, so a channels-major MCS60 file would have given duration_s = 60/fs and firing rates wrong by orders of magnitude) while steps 3/4 compared against the real channel count. Step 2 also used to invent a duration (max spike time, else literally 60 s) when the raw file was unreadable; it now skips the recording like steps 3/4 do, since every firing rate scales with this number.

Verified end to end on the example data: steps 1-3 into folder A, then step 3 alone and step 4 alone resumed from A with raw_data pointed at a non-existent folder — same seed reproduced A's adjacency bit-for-bit, and step 4 produced its CSVs from A's spikes + adjacency.

CAT-NAP resume (2026-08-06)

The calcium path had none of this: runner.py returned straight out of the suite2p_mode branch before the step-range logic, so start_analysis_step, stop_analysis_step and prior_analysis were silently ignored — and nothing intermediate was written, so there was nothing to resume from even in principle. (Only denoising was reused, and via the raw suite2p folder: Fdenoised.npy & friends are written next to the inputs, not into the output folder.) MATLAB has always supported this: its step 2 appends adjMs to ExperimentMatFiles/<rec>_<folder>.mat, and step 4's priorAnalysis == 1 && startAnalysisStep == 4 branch loads it back.

Ported in catnap/store.py: phase 1 writes one <rec>_catnap.npz per recording (adjacency, coords, channels, per-node activity counts, duration, the coordinate normalisation, and the calc_twop_activity_stats dict), and start_analysis_step >= 4 reads it via InputLocator.catnap_file() instead of recomputing. Step 4 is the only resumable boundary here, because CAT-NAP has no step 1 or 3 — adjacency is built in step 2, as in MATLAB.

Deliberate choices worth knowing:

• The fluorescence matrices are not stored. MATLAB saves F/denoisedF/ spks because its step 4 re-derives the activity matrix from them; this port stores the per-node quantity step 4 actually consumes (spike_counts) and keeps the file small — a resumed run rewrites it into the new output folder so that folder is resumable in turn. Phase 3 still re-reads the suite2p folder for the trace figures and the mean-projection backdrop, and already degrades gracefully when the raw data isn't mounted.
• Cell-type groups are re-read, not stored, so a resumed run picks up an edited spreadsheet rather than freezing the first run's grouping.
• Different filename from the ephys _adjM.npz. The two hold different things; pointing a CAT-NAP resume at an ephys output folder now reports a missing input instead of failing deep inside a load.
• The phase-2 cartography barrier is batch-wide. Resuming a subset of the batch re-derives the boundaries from that subset alone, so the roles won't match the original run. Documented rather than worked around — pooling is the point of that barrier.

Two bugs fell out of the same code and were fixed alongside:

1. The metric loop was driven by Params.funcConLagval. suite2pToAdjm derives a single lag round(1000 / fs) for the F/spks/denoised F activity types and ignores funcConLagval entirely (MATLAB even overwrites Params.FuncConLagval with it). Iterating the Params list meant the key adjM{lag}mslag matched nothing and those three activity types produced no network metrics at all. The loop now follows the adjacency that was actually built (store.sorted_adjm_items). Only twop_activity="peaks" was ever exercised end to end, which is why this survived.
2. One RNG stream for the whole batch. Unlike steps 3/4, CAT-NAP drew every stochastic stage from a single make_rng(seed, "catnap"), so a recording's metrics depended on how much randomness the recordings before it consumed — and a resumed run, having skipped the thresholding draws, would have produced different numbers from the run it resumed. Now per-recording: ("catnap", rec) for adjacency, ("step4", rec) for network metrics (same label as the ephys path — same work, same recording), ("catnap_subnetwork", rec) for the subnetwork analysis.

stop_analysis_step is still not honoured on this path — the calcium flow is one step, not four — but the runner now says so in the log instead of ignoring it, as it does for a start_analysis_step of 2 or 3.

Tests: python/test_catnap_resume.py (46 checks — round-trip including the None-vs-NaN distinction in the calcium imaging stats, locator resolution, fail-fast, resuming with no suite2p folder on disk at all, and a determinism check that was confirmed to fail against the old batch-wide RNG before being kept).

Random seed (Params.random_seed, 2026-07-31)

pipeline/rng.py's make_rng(seed, *labels). None (the default) keeps the old behaviour — fresh OS entropy per run, matching MATLAB, which seeds nothing. Set an integer and every stochastic stage becomes reproducible within the Python port (this can't and doesn't make anything match MATLAB — see item 4 under "Known limitations").

Two non-obvious constraints drove the design:

1. Results must not depend on scheduling. Steps 3/4 run recordings across a process pool auto-sized from the machine's cores/RAM, so one shared stream would make the numbers depend on the host. Every generator is instead derived from (seed, "step3"|"step4", rec.filename) — a recording gets the same stream whether it ran first, last, or alone.
2. Labels are hashed with BLAKE2b, not hash(). Python randomizes string hashing per process, so hash("rec_A") differs between the parent and a spawned worker — exactly the failure this is meant to prevent. Verified under two PYTHONHASHSEED values.

Wired at: step3.py (thresholding surrogates), step4.py (modularity, PC-norm and small-worldness null models, NMF), plotting.py (which example channels the step-1 check figure shows), stim_step.py figures, and the CAT-NAP path.

One extra fix was needed to make "same seed → same numbers" literally true. average_controllability/modal_controllability are deterministic metrics, but both call svds(A, 1), and ARPACK draws its start vector from NumPy's global RNG — so they wobbled in the last ~2 ulps between runs on identical input, and no random_seed could have fixed it (the draw doesn't come from any generator this pipeline owns). _dominant_singular_value() now pins v0 to a fixed vector. MATLAB parity is unaffected: python/parityRun/compare_step4_same_adjacency.py still reports worst-diff 2.842e-14 with aveControl at ~2.9e-15, unchanged.

With that in place, two runs at random_seed=42 produce byte-identical NetworkActivity_RecordingLevel.csv and NetworkActivity_NodeLevel.csv (all 48 columns, Q/SW/PC/num_nnmf_components included); random_seed=7 differs, as it should. No BLAS thread pinning was needed.

Verbose level (Params.verbose_level, 2026-09-04)

Was stored, written to params.json and never read — the GUI offered three levels and every run printed the same thing. pipeline/verbosity.py now makes it mean something.

The design constraint is that the levels are additive: Verbose adds to what Normal prints, Debug adds to both, and nothing is ever removed. A level that could suppress a line would make a Debug log unsafe to reason from (is that warning missing, or was it never printed?) and would make raising the level a decision rather than a free one. Concretely: none of the ~150 existing log(...) call sites had to be audited or changed.

RunLog is callable, so log("…") still means "print at every level"; the new channels are log.detail(…) (Verbose+), log.debug(…) (Debug only) and log.timed(label) (a context manager that reports elapsed time at Verbose+). run_pipeline promotes whatever callable it was given once, at the top, so every step and every worker downstream inherits the level.

Two things worth knowing:

1. Workers get the level through Params, not through the log. Steps 3 and 4 compute in spawned processes and buffer their own lines into a list that the parent replays. Each promotes its own buffer with as_run_log(logs.append, params.verbose_level) — the level travels in the pickled Params that was going there anyway.
2. as_run_log never re-levels an existing RunLog. A function further down the stack guessing a level from its own arguments is how the run's level and a step's level would drift apart.

MATLAB's Params.verboseLevel is 'Normal' | 'High' | 'Silent'; normalise_level accepts those names (High → Verbose, Silent → Normal, there being no quieter level here) so a params file from that pipeline keeps its intent instead of silently landing on the default. Anything unrecognised — a hand-edited file, a None — logs normally rather than raising.

What each level adds, and where:

Level	Adds	Where
Normal	(unchanged)	everywhere
Verbose	the batch and its groups/DIVs; every setting that differs from its default; per-method spike counts and mean rate; active electrodes, mean FR, burst rate; edges surviving thresholding per lag; active nodes, density, modules/Q, small-worldness, global efficiency per lag; per-stage timings	runner._log_run_header, _describe_spikes, step2._describe_activity, step3._describe_edges, step4._describe_network, catnap._describe_adjms
Debug	file paths, sampling rates and durations, detection/threshold settings, worker counts, Python/platform/library versions, and every setting rather than the changed ones	the same places, via log.debug

Checked by python/test_verbose_level.py (~5 s), whose load-bearing assertion is that the netmet JSON is byte-identical across all three levels: turning Debug on to diagnose a run must not invalidate it.

Spreadsheet range convention

Params.spreadsheet_range is a string like "A2:A100000". The parser (parse_spreadsheet_range) extracts the two numbers and treats them as 1-indexed file line numbers, header = line 1 — i.e. "A2:A3" reads the first two data rows, matching MATLAB's csvRange = [2, 3] / DataLines semantics exactly. Column letters are accepted but ignored (kept for readability/familiarity, not functionally meaningful). Default is "A2:A100000" — not "A1:Z1000", which would include the header as a data row.

Batch scaling (scaled-to-whole-dataset plots)

MATLAB renders a "scaled to entire data batch" variant of several plots so activity/metric levels are visually comparable across recordings, not just within one. It does this with a two-pass structure: a first pass over all recordings computes batch-wide maxima/bounds, a second pass plots each recording against those shared scales. This port now mirrors that:

• Step 2 raster (plotting_step2.py's plot_raster, port of rasterPlot.m): now a two-panel figure — top scaled to the recording's own raster_plot_upper_percentile, bottom scaled to spike_freq_max, the batch-wide max per-channel firing rate (MATLAB's maxValStruct.FR). step2.py computes spike_freq_max after its compute loop and threads it into a second plotting loop. (The old single-panel raster + avg-FR histogram was replaced — rasterPlot.m has no histogram panel.)
• Step 2 electrode + burst heatmaps (plotting_step2.py's plot_heatmap, port of electrodeHeatMaps.m / plotNodeHeatmap.m): now a two-panel figure when a batch max is supplied — left scaled to the recording (99th-percentile color axis), right scaled to the entire dataset (color axis = that metric's batch-wide max). step2.py computes a batch_max dict over all six of MATLAB's valsTogetMax metrics (FR, channelBurstRate, channelBurstDur, channelFracSpikesInBursts, channelISIwithinBurst, channeISIoutsideBurst) in the same compute pass and threads it into the plotting loop. Covers 2_Heatmap.png (FR) + the five burst heatmaps (3_BurstRate … 7_ISIoutside). Verified on A2/A3: the right ("entire dataset") panels share the same color ceiling (~2.75 Hz FR) across both recordings while the left panels use each recording's own — A3 looks vividly active scaled to itself but correctly reads as much quieter than A2 on the shared batch scale. (MATLAB additionally log-scales the BurstDur and ISIoutsideBurst heatmaps via useLogScale; this port keeps them linear for now — see limitation note. Burst heatmaps are all ~zero on the example data anyway, since it has no detected bursts.)
• Step 4 network plots (plotting_step4.py's plot_spatial_network → network_plot.py's plot_network): each spatial plot now also writes a _scaled sibling and a _combined side-by-side figure (plot_spatial_network_combined, port of MATLAB's N_combined_MEA_NetworkPlot — individual scale left, batch scale right). step4.py was restructured into compute-then-plot: the compute loop collects plot_jobs, then _batch_metric_bounds() pools each node-level metric (ND/NS/BC/PC/Eloc) across all recordings and all lags (port of findMinMaxNetMetTable.m, which reads the whole NodeLevel.csv column), then _plot_recording_lag() draws both variants. The scaled variant passes three overrides into plot_network: z_scale_override (batch max of the size metric → node sizes comparable across recordings), z2_bounds_override (batch min/max of the color metric), and edge_bounds_override fixed to (0.1, 1.0) — MATLAB hardcodes minMax.EW = [0.1, 1]. These mirror MATLAB's useMinMaxBoundsForPlots / Params.metricsMinMax mechanism exactly. Verified: on the two example recordings (A2, A3), the _scaled plots share identical legend scales (node degree max 17, BC 0–0.263, fixed edge weights) while the individual plots differ (A2's own ND max is 14) — i.e. the batch scale is genuinely shared. Because the example data is low-connectivity (edges ~0.04–0.10, well under the fixed [0.1, 1] edge scale), the scaled variants' edges render faint — this is faithful to MATLAB, not a bug.
• Not batch-scaled yet: connectivity-stats / cartography (cartography has fixed [0,1]×Z axes in MATLAB anyway, so there's nothing to batch-scale). Also the log-scaling MATLAB applies to the BurstDur and ISIoutsideBurst heatmaps (useLogScale = [0,1,0,0,1] in MEApipeline.m) is not ported — the Python heatmaps are linear. Low priority: it's a cosmetic axis transform, and it can't be exercised/verified on the example data (no bursts → all-zero burst metrics). Add it to plot_heatmap (a log_scale flag transforming values + color axis to log10) if a bursting dataset makes it matter.

Known limitations / next steps

Last audited 2026-07-08 — if you fix or port something in this list, please either delete the item or strike it through (like the coordinate-lookup one below) rather than leaving it to silently go stale like the previous version of this list did.

1. Pipeline runs synchronously on the Qt UI thread. Done — the pipeline now runs on a background QThread (src/meanap/gui/pipeline_worker.py, PipelineWorker). run_pipeline gained a should_cancel: CancelCheck parameter (see src/meanap/pipeline/cancellation.py) that it polls at every step boundary and once per recording (and per lag inside the expensive step-3/4 inner loops); when it returns True the run unwinds by raising PipelineCancelled, which the worker turns into a cancelled signal. The GUI's Stop button (main_window.py's _on_stop) now calls worker.request_cancel() for a real cooperative cancel — it halts at the next checkpoint (i.e. after the in-flight recording/lag finishes, not mid-computation), not just resets the buttons. Log lines and the final outcome come back via Qt signals (log_message/finished_ok/cancelled/ failed) delivered on the UI thread, so no more QApplication.processEvents() pumping. closeEvent requests cancel + wait(5000) so closing the window mid-run doesn't destroy a live QThread. Granularity caveat: because cancellation is cooperative (polled at loop boundaries), a click during a single long operation — one recording's wavelet CWT, or one lag's 200-shuffle thresholding — won't take effect until that operation returns.
2. Controllability not ported. Done — average_controllability()/ modal_controllability() (ports of the Bassett Lab ave_control.m/ modal_control.m) are in network_metrics.py, computed in step4.py (aveControl/modalControl + mean/top-25%/threshold summary stats), and plotted as 10_…Averagecontrollability.png / 11_…Modalcontrollability.png (see item 3 below). Small-worldness not ported. Done — small_worldness_rl_wu() (port of small_worldness_RL_wu.m) is in network_metrics.py, fed by two new null-model rewiring functions in null_models.py (randmio_und_v2/latmio_und_v2, ports of the same-named .m files in Functions/CC_PL_SW/ — distinct from randmio_und_signed/ null_model_und_sign above, which serve participation_coef_norm instead). Wired into step4.py's compute_network_metrics as SW/SWw/ CC/PL (the last two are what MATLAB actually saves into NetMet.CC/NetMet.PL — the raw, unnormalized values this module already computed are now under CC_raw/PL_raw/CC_rawMean instead, to avoid colliding with the real field names; see network_metrics.py's docstring for why both are kept). MATLAB's own gate for this block is aN > minNumberOfNodesToCalNetMet (strictly greater — unlike every other block's aN >= ...), faithfully replicated in step4.py, not a typo. Formula-assembly parity against MATLAB (fixed A/R/L) is exact — see python/test_pipeline_small_worldness.py and python/test_fixtures/gen_small_worldness_reference.m. The null models themselves aren't bit-reproducible against MATLAB (independent RNG streams), validated instead via structural invariants (exact degree/ weight preservation), same situation as null_models.py's other functions. NMF-based num_nnmf_components not ported. Done — nmf.py's cal_nmf(), port of calNMF.m. Unlike everything else in this port, this one is not just RNG-stream-different from MATLAB but algorithm-different: MATLAB's nnmf defaults to Alternating Least Squares, nmf.py uses sklearn.decomposition.NMF (coordinate descent, the closest available equivalent) — different solvers can converge to different local optima and even pick a different num_nnmf_components, since that value depends on exactly where each solver's residual crosses a phase-randomized reference's. Also computes the (mathematically identical, see nmf.py's docstring) downsampled spike matrix directly at the target resolution rather than materializing MATLAB's huge intermediate native-fs matrix first. Performance: nnmf_residuals/nnmf_var_explained are unconditional (saved regardless of Params.includeNMFcomponents) and require one NNMF fit per rank 1..(active-electrode count), which is the dominant cost — ~55s for a real 64-channel/600s recording (measured on this environment's CPU), independent of lag count since it's computed once per recording, not once per lag. Spatial/temporal autocorrelation not ported. Still true, and staying that way for now: SA_lambda/SA_inf/TA_regional/TA_global aren't in MATLAB's own default netMetToCal list either (AdvancedSettings.m calls them out as "other optional ones"), and ExtractNetMet.m's own temporal-autocorrelation code path is an explicit % TODO stub — there's no complete MATLAB reference behavior to port yet. Everything else in ExtractNetMet.m's own docstring list of metrics is done with either exact or "deterministic given a fixed stochastic input" parity — see the "Where things stand" table and network_metrics.py's docstring for the precise breakdown. Cmcblty (communicability) needs no work: MATLAB's own code path for it is commented out upstream.
3. Step 4 plots: 9 of MATLAB's 4A_IndividualNetworkAnalysis plots are ported (connectivity stats, base/BC/PC/Eloc-colored network plots, circular module-colored network plot, node cartography, average/modal controllability — see plotting_step4.py). Batch-scaled ("scaled to entire data batch") variants are now ported for the four spatial network plots (2/3/4/5) plus the two controllability plots (10/11) — each writes a sibling N_scaled_MEA_NetworkPlot….png whose node size / node color / edge-weight scales come from batch-wide pooled bounds instead of the recording's own range (see the "Batch scaling" note below). The N_combined_MEA_NetworkPlot….png side-by-side figure is now ported too (plot_spatial_network_combined — left panel scaled to the recording, right panel scaled to the batch, each with its own inline legend/colorbar; rendered as one two-axis figure rather than MATLAB's copyobj merge of two separate figures). Still not ported: the null-model-iterations panel (plotNullModelIterations.m — a diagnostic plot of the met/met2 convergence traces latmio_und_v2/randmio_und_v2 collect every 1000 rewiring iterations; item 2's null-model calculation is done, this is purely the unported diagnostic plot on top of it).
4. Step 3's probabilistic thresholding is not bit-reproducible against MATLAB by design (see "STTC gotchas" above) — this is not a bug to fix, it's inherent to comparing two independently-seeded RNG streams. The deterministic STTC computation itself has exact parity. Same situation for Step 4's modularity (Ci) and the null-model randomization inside the normalized participation coefficient. Run-to-run reproducibility within the Python port is now available via Params.random_seed (see the "Random seed" section above) — that is a different guarantee from parity with MATLAB, and it does not move any of these closer to MATLAB's numbers.
5. No Parameters_*.csv/.mat or channel_layout.png generation. MATLAB writes these as part of its run; Python writes .npz/.json outputs per step instead (see "Key files" above for exactly what each step writes) and the ExperimentMatFiles/*_adjM.npz per-recording adjacency file from step 3, but nothing analogous to MATLAB's single combined <recording>_<output>.mat (Info/Params/spikeTimes/Ephys/adjMs all in one file).
6. wid_ms / n_scales / abs_thresholds not exposed in the GUI or Params→SpikeDetectionParams wiring in runner.py. The runner uses SpikeDetectionParams's own defaults (wid_ms=(0.4, 0.8), n_scales=5) rather than reading them from Params, because Params has no fields for them yet. abs_thresholds is stored in Params but never read by the runner — only relative thresholds are wired up.
7. filter_high_pass default is 8000.0 Hz in Params, whereas the parity-validated example run used 6150.0 Hz — both are legitimate, user-tunable values; just don't be surprised if the GUI-run parity looks slightly different from the test_pipeline_step1.py numbers.
8. The example dataset barely exercises step 3/4 downstream logic. Firing rates are low enough (~0.16-0.20 Hz mean) that density stays under ~10% and many node pairs are disconnected (NE=0 for those nodes, PL/Eglob dominated by very sparse paths) — plausible and validated against MATLAB, but not a stress test. If a busier reference dataset shows up, it's worth re-running all the parity tests against it.
9. No full group-level (Step 2B/4B 4_RecordingsByGroup/5_GraphMetricsByLag/ etc.) analysis in the port. Most outputs are per-recording, per-lag. The one exception is batch scaling (see the "Batch scaling" section above): steps 2 and 4 now do a first cross-recording pass to compute batch-wide plot bounds, then a second plotting pass — but this only feeds the scaled plot variants, it does not produce MATLAB's group-comparison figures/CSVs.
10. No electrode coordinate lookup table Done — channel_layout.py ports getCoordsFromLayout.m (exact parity, all 5 layouts).
11. No HTML/output browser Done — report.py, wired into the GUI as "🌐 View report" on the Run tab. See its own section above.

Auditing for silently-missing metrics

A real gap in this port (found 2026-07-08, now fixed) went undetected for a while because it wasn't a crash or a wrong number — it was fields that were simply never computed, with no test asserting they existed. Controllability and small-worldness were caught because someone happened to remember seeing plots for them; the record-level summary-stat scalars (NDmean, NDtop25, NSmean, sigEdgesMean, sigEdgesTop10, PCmean, PCmeanTop10, PCmeanBottom10, percentZscoreGreaterThanZero, percentZscoreLessThanZero) were only found by deliberately auditing against MATLAB's ground truth for "what should exist" rather than trusting this doc's own "still missing" list (which had gone stale — it only tracked things someone had noticed were missing).

The audit that found them: Functions/ExtractNetMet.m's docstring (lines 26-50) lists every NetMet field name, and AdvancedSettings.m's Params.netMetToCal/Params.networkLevelNetMetToPlot/ Params.unitLevelNetMetToPlot defaults are the authoritative list of what a default pipeline run actually produces. Cross-referencing that full list against grep -n '"X"\]\|result\["X"\]' src/meanap/pipeline/step4.py network_metrics.py for every field name is the reliable check — checking this doc's own "known limitations" prose against the code is not, since prose gets stale in exactly the way that just happened. A second useful signal: plotting_step4.py's NETMET_REC_METRICS/NETMET_NODE_METRICS display-name dicts were built ahead of the metrics they label (i.e. someone added the label anticipating the field, then the field itself never got written) — grepping for display-dict keys with no corresponding result["X"] = ... in step4.py catches exactly this pattern. If you're auditing again later, also grep ExtractNetMet.m for every distinct function call (grep -oE '[a-zA-Z_][a-zA-Z0-9_]*\(' Functions/ ExtractNetMet.m | sort -u) and check each one has a Python equivalent somewhere in network_metrics.py/nmf.py/null_models.py — that's how randmio_und_v2/latmio_und_v2/calNMF/calEffRank were confirmed to be the only remaining unaccounted-for function calls at the time of this audit (excluding MATLAB builtins, the checkIfRecomputeMetric/ prevNetMet incremental-caching machinery Python doesn't replicate, the commented-out fcn_find_hubs_wu block, and binary-adjacency-matrix-only functions like betweenness_bin/distance_bin, since this port only supports Params.adjMtype == 'weighted').

Extended to Steps 1-3 the same day. Same methodology, different authoritative lists: firingRatesBursts.m's own field-by-field assignments to Ephys.* (there's no netMetToCal-equivalent list for Step 2 — the closest thing is saveEphysStats.m's hardcoded NetMetricsE/NetMetricsC cell arrays, which turned out to be the more useful reference since they're also the CSV column whitelist) and batchDetectSpikes.m's varsList (what actually gets saved to each recording's _spikes.mat). Found two more real gaps, both now fixed:

1. Step 2 had no CSV export at all — firing_rates_bursts()'s ephys dict was complete and correct (already validated at 100% parity), but step2.py only ever wrote it to ephys_results.json; the NeuronalActivity_RecordingLevel.csv/NeuronalActivity_NodeLevel.csv files MATLAB's saveEphysStats.m produces (called unconditionally in MEApipeline.m, not gated behind an optional flag) never got written by the Python port. Fixed — see _save_ephys_stats_csv() in step2.py.
2. Electrode grounding was parsed but never applied — RecordingInfo.ground (the spreadsheet's optional 4th Ground column) was already being read by read_recording_csv(), but nothing consumed it anywhere in the pipeline. Fixed — see spreadsheet.py's parse_ground_electrodes()/ground_spike_times_dict() note above. This one's a good example of why grepping call sites matters more than grepping field/parameter existence: Params.ground-equivalent fields existing in params.py/RecordingInfo proved nothing about whether they were wired to anything.

Also checked and confirmed not gaps (dead code in MATLAB itself, so correctly unported): Ephys.channelHasBurst (computed by firingRatesBursts.m but not in saveEphysStats.m's CSV whitelist, same "computed but never surfaced" pattern as Step 4's BCmeantop5); spikeAmps/getSpikeAmp.m (only consumed by the stimulation-analysis code path, Params.stimulationMode == 1, out of this port's scope); Params.ProbThreshPlotChecks's diagnostic replicate-stability plot (adjM_thr_checkreps.m — optional, off by default, and Params .prob_thresh_plot_checks already exists in params.py unwired to match that same off-by-default state, not silently dropped).

MATLAB vs Python speed comparison

Current: Python ~7x faster (re-measured 2026-10-25, MEA-NAP 1.15.0)

Same dataset and settings as the original comparison below (2 Axion64 recordings, lags [10, 25, 50]ms, 200 shuffles, every figure drawn). MATLAB ran python/parityRun/MEApipeline_timingBenchmark.m (R2025b, 16-worker parpool), then Python ran python/benchmark_express.py's full-run settings immediately afterwards on the same machine:

Step	MATLAB	Python	Ratio
1. Spike detection	445s	39.6s	11x
2. Neuronal activity	62s	7.6s	8x
3. Functional connectivity	9s	1.8s	5x
4. Network activity	377s	76.6s	5x
Total	893s (~15 min)	125.6s (~2 min)	~7x

An earlier Python run that day on the same box took 105.6s (step 4 56.5s), which would make the ratio 8.5x. Step 4 is the step that varies with load on this shared machine, so "7-8x" is the honest range. MATLAB's own total is in line with July's 931s.

Express mode, which MATLAB has no equivalent of, finishes the same analysis in 53.3s (step 4 12.7s, about 17x MATLAB's full run), and every value it computes is byte-identical to the full run's.

What changed since the 2.66x below: step 1 got faster as the spike-detection channel threading matured, and steps 3/4 gained the sorted-scale work ("Hundreds of sorted units" above): STTC across cores, threaded NMF with an early stop, and threaded exact null models. Against the code just before that work (356fc0f), the same benchmark went 236.1s -> 105.6s full and 168.7s -> 53.3s express, with all CSVs, adjacency matrices and spike files byte-identical to the older code in both modes.

Figures are now most of what is left in a full run's step 4 (about 44s of 56.5s), and step 1 is most of an express run.

Since then (1.15.1) step 1 is 2.5x faster ("Spike detection — 2.5x on step 1" above): the same full run takes 80.8s on this machine (express 29.3s). That was not measured back to back with MATLAB, but against the 893s above it would be ~11x.

Original comparison (2026-07-09)

Full steps-1-4 run on the Test Pipeline dataset (ExampleData/, recordings NGN2_20230208_P1_DIV14_A2/_A3, Axion64 layout, ProbThreshRepNum=200, lags [10, 25, 50]ms), timed via MATLAB's Params.timeProcesses and the Python port's newly-implemented equivalent (Params.time_processes, see runner.py's entry above), both run alone (no concurrent load) on the same machine:

Step	MATLAB	Python	Ratio
1. Spike detection	580s	568.8s	Python ~2% faster
2. Neuronal activity	53s	6.2s	Python ~8.5x faster
3. Functional connectivity	15s	80.2s	Python ~5.3x slower
4. Network activity	283s	511.7s	Python ~1.8x slower
Total	931s (~15.5 min)	1166.9s (~19.5 min)	Python ~25% slower overall

Not a uniform "Python is slower" story — it's faster than MATLAB at the two purely-numerical steps, slower specifically where MATLAB parallelizes:

• Step 3: MATLAB's log shows it spinning up an 8-worker parallel pool for the 200-shuffle probabilistic thresholding; adjm_thr() (this port) runs single-threaded. Likely the dominant cause of the 5.3x gap — not yet investigated further (e.g. whether MATLAB is also using a compiled/mex STTC core on top of the parallelism).
• Step 4: NMF alone costs ~55s per recording (~110s of the 511.7s total, measured earlier this session on real data) — sklearn's coordinate-descent solver run serially for every rank from 1 to the active-electrode count, see nmf.py's docstring. MATLAB likely also parallelizes some of its BCT calls here (unconfirmed). The other ~400s is the rest of step 4's compute
  • plotting, not yet broken down further.
• Step 2: Python's vectorized numpy firing-rate/burst computation is dramatically faster than MATLAB's implementation, which does per-channel LOESS regression fitting for the burst ISI threshold (getISInTh.m → fLOESS.m — also the source of the "Matrix is close to singular" warnings MATLAB prints throughout Step 2 on this low-firing-rate dataset).

Obvious next lever if closing the gap matters: parallelize steps 3/4 (e.g. multiprocessing/joblib across recordings or lags) — that's where essentially all of the total deficit lives.

Optimization pass — 3.3x faster end-to-end (2026-07-09)

A follow-up performance pass took the full run from 1166.9s to 350.0s — Python is now ~2.66x faster than MATLAB (931s), not 25% slower. Same machine, same config (2 recordings, [10, 25, 50]ms, ProbThreshRepNum=200):

Step	Before	After	Speedup	What changed
1. Spike detection	568.8s	103.6s	5.5x	Channel-level threading (detect_spikes_recording(max_workers=)): per-channel wavelet CWT + bandpass release the GIL, so threads parallelize over one shared ~3.8 GB dat array with no extra RAM. Bit-identical output.
2. Neuronal activity	6.2s	6.4s	—	Unchanged (already fast).
3. Functional connectivity	80.2s	42.7s	1.9x	Recording-level process map.
4. Network activity	511.7s	197.3s	2.6x	numba-JIT randmio_und_signed (~28x on the PC-norm null models); collection-based plotting (~5x plot phase, single-threaded); recording-level process map.
Total	1166.9s	350.0s	3.3x

Key pieces (all default-on, with serial fallbacks; see src/meanap/pipeline/parallel.py):

• RAM/CPU-aware worker sizing. Step 1 is RAM-bound (~3.8 GB/recording) → threads over channels on one shared array (no per-worker copy). Steps 3/4 are CPU-bound/low-RAM → processes over recordings. Counts auto-size to physical cores and free RAM (psutil) with headroom; override via params.spike_detection_channel_workers / params.recording_workers (None = auto). Per-worker BLAS threads pinned to avoid oversubscription.
• numba is now a dependency (forces numpy<2.5); randmio_und_signed has an @njit core + pure-Python fallback (runs without numba, just slower). All step-3/4 parity fixtures pass under numpy 2.4.6.
• Plotting (network_plot.py): edges were individual ax.plot() Line2D artists and nodes individual add_patch circles (hundreds per figure × ~380 figures) → one LineCollection + PatchCollection(match_original=True). Cut the plot phase ~5x and removed a get_text_width_height_descent layout pathology (458s→1.1s in profile — the per-edge artists were making tight_layout/bbox re-measure text catastrophically). Output pixel-validated unchanged.
• Step 4 restructured into A (parallel compute) → B (serial batch-bounds reduce) → C (parallel plot). Deterministic metrics bit-identical serial-vs-parallel. Phase B is also where the cross-recording node-cartography boundary clustering runs (_apply_cartography_boundaries → network_metrics.trial_landscape_density, port of TrialLandscapeDensity.m's default kmeans branch — see the node-cartography note below).

Caveats: measured on a shared 16-core box where step-4 process-parallelism scales poorly (memory-bandwidth/BLAS contention → only ~1.1-1.4x across recordings), so most of step 4's win is the numba + plotting changes, which are single-threaded and portable. On a dedicated ~8-core machine the recording map should contribute more. Remaining portable candidates: numba randmio_und_v2/latmio_und_v2 (small-worldness, same @njit pattern); cutting NMF's fit count (init="nndsvda" was tried — ~8% slower since the per-rank sweep pays NNDSVD's SVD-init cost on every fit, reverted).

How this was measured — headless (Params.guiMode=0) execution of MATLAB MEA-NAP is not a well-trodden path; getting a clean end-to-end run required finding and working around 5 bugs/config traps that the interactive GUI (runPipelineApp.m/getParamsFromApp.m) silently avoids. See MEApipeline_timingBenchmark.m's header comment (repo root, not part of the Python port — a throwaway benchmarking copy of MEApipeline.m) for the full list and exact fixes; briefly:

1. batchDetectSpikes() called with 5 args in MEApipeline.m's non-GUI branch, but its arguments block requires 6.
2. Params.electrodesToGroundPerRecording is left as a non-cell [] by setUpSpreadSheet.m when there's no spreadsheet Ground column, but Step 2/3/4 unconditionally brace-index it.
3. AdvancedSettings.m's default Params.netMetToCal includes 'PC_raw', whose computation is dead code — crashes Step 4's eval-assignment loop.
4. Several Params fields (minActivityLevel, nodeScalingMethod, networkPlotEdgeThresholdMethod, etc.) are only ever set by the GUI, with no AdvancedSettings.m fallback — silently undefined in headless runs.
5. Params.nodeLayout='Original' (not 'MEA') is StandardisedNetworkPlot .m's actual sentinel for "use real electrode coordinates"; several other plotting/scaling params (nodeScalingMethod='Linear' not 'degree', maxNodeSize=1 not 0.06, networkPlotEdgeThresholdMethod='Absolute Value' not 'percentile', etc.) also don't match the Python port's same-sounding defaults — the two ports' vocabularies/scales for these particular fields diverged. The values finally used were extracted from Parameters_OutputData26Jun2025.csv, a completed GUI-driven MATLAB run on this exact dataset — i.e. "what MEA-NAP actually computes in practice," not a guess.

Items 1-4 are genuine MATLAB-side bugs (would affect anyone scripting MEA-NAP outside the GUI, e.g. for HPC batch jobs); item 5 is this benchmark's own config error, not a MATLAB bug — noted here so a future re-run doesn't have to rediscover it.

Dense-network hot loops — 6.1x on CAT-NAP (2026-08-26)

A 378-recording CAT-NAP batch (one 1000ms lag, express + Dropbox) took a collaborator 99h48m, ~951s per recording. Profiling one 270-cell recording locally reproduced it almost exactly (960s), so this was compute, not the network link.

The cause is network density, not calcium imaging as such. A 1000ms coincidence window leaves the thresholded STTC matrix nearly complete — that batch had a median density of 0.86, with 55 of 378 recordings at >=0.99 and 17 at exactly 1.000. Several BCT routines are fine at a 64-electrode MEA's ~2000 edges and quadratic-or-worse at the ~26000 edges a 270-cell recording reaches at that density:

Fix	What was wrong	Gain
null_models.null_model_und_sign	the weight-assignment loop rebuilt four Python lists of length m per pass, and it makes ~m/10 passes	3.3x
null_models._stable_ranks	that loop argsorted every remaining edge to read wei_period (10) random ranks; argpartition places just those in O(m)	5.0x cumulative
sttc.get_sttc	run_T reads a single spike train, but sttc_pair recomputed it per pair — O(n^2) calls where O(n) suffice, 201 times per lag	3.0x
network_metrics.distance_wei	pure-Python Dijkstra; efficiency_wei_local runs one per node over that node's neighbourhood, so at density ~1 it is O(n^4)	18x on local efficiency

All four are exactly equivalent, not approximations, and are verified as such: python/test_pipeline_fast_paths.py checks each against the form it replaced (including the tie case that _stable_ranks falls back on, which real data never produces), and an end-to-end run over 5 real recordings reproduced every adjM and all 58 metric fields per recording bit-for-bit against the previous commit.

Measured end to end on those 5 recordings: 1062s -> 173s. The largest (N=270, one lag) went 1052s -> 166s. Modelled over the collaborator's actual N/density distribution, her batch drops from ~46h to ~8h of compute on the box these numbers were taken on.

What is left. participation_coef_norm is still ~87% of the remainder (~38s of that is randmio_und_signed, the rest the weight-assignment loop), because it runs 100 null models per recording per lag. Cutting that further means either parallelising the 100 draws or taking fewer of them, and both change results — the draws come from one sequential RNG stream, so any parallel split changes the values, and n_iter is a straight accuracy/time trade. Neither has been done.

Parallel network metrics, streaming-safe (2026-08-26)

Parallelising across recordings is free numerically, since each is already seeded from its own filename — so that is what landed (catnap 1.6.0).

The obstacle was the remote path. map_recordings wants the whole task list up front, and the CAT-NAP stream cannot give it one: recordings arrive from a source that fetches, analyses and releases each before the next, which is what keeps peak local storage at one or two recordings rather than the whole dataset. What makes this work anyway is that network metrics need only the adjacency matrix, never the suite2p folder — so the metrics can go to a pool while the fetch/denoise/STTC half stays serial and keeps its storage guarantee. Workers never touch raw data, so the parallelism costs no disk.

pipeline/parallel.py gained StreamingPool for this: submit as the stream yields, collect in completion order, and block the producer once saturated (without that bound the parent queues every remaining recording's adjacency matrices — gigabytes on a large batch). Results are re-keyed into spreadsheet order before the cartography barrier, batch bounds or CSVs see them, so nothing downstream can notice the reordering. Cancellation is checked between completions as well as between submissions, so a stop during the drain drops queued work instead of waiting out the batch.

Verified in python/test_pipeline_streaming_pool.py: 1, 3 and 4 workers give identical metrics, reversed submission order gives identical metrics, and a real 13-recording batch produced the same netmet_results.json digest at 1 and 8 workers.

The __main__ trap, and why probing does not find it. spawn re-imports the parent's __main__ in every worker, so a caller with no if __name__ == "__main__": guard re-runs its whole script per worker — for a pipeline script, re-running the pipeline into the same output folder. This bit test_viewer_run_styling.py / test_viewer_trace_download.py as soon as CAT-NAP started spawning, and it is a trap map_recordings had all along; the ephys tests just never had enough recordings to open a pool.

Python's recursion check does not reliably stop it. A child whose nested pool attempt raises, but which catches it, goes on to complete its task — so the parent sees a healthy pool and never learns that every worker re-ran the script. A probe task therefore reports success in exactly the case it is meant to catch (measured: four MARKER: top-level ran lines, spawn_usable True). Only reading the entry point works, so spawn_usable() parses __main__ for a module-level guard and keeps unguarded callers serial with an explanatory warning. map_recordings uses it too. StreamingPool._degrade remains as the backstop for pools that die for other reasons (an OOM kill looks identical from the parent), finishing outstanding tasks in-process rather than failing the run.

The ceiling is the largest recording. A batch cannot finish before its slowest single member, and on the collaborator's dataset the ten biggest recordings are 61% of the compute — the local 13-recording cache is worse still (one 137s recording and twelve trivial ones), which is why it only moves 178s -> 148s. Getting past that ceiling requires the within-recording parallelism above, which changes results.

On a batch of equally large networks — which is what the collaborator's biggest recordings look like — the pool scales close to linearly. Eight copies of the N=270, 25964-edge recording, one lag, on a 16-physical-core box:

Workers	Wall	Speedup	Efficiency
1	1101.6s	1.00x	—
4	310.1s	3.55x	89%
8	166.2s	6.63x	83%
16	168.1s	6.55x	(only 8 tasks — nothing left to parallelise)

Same digest at every worker count.

The density is worth raising with whoever runs these. At 86-100% density the network is close to complete, which makes modularity, participation coefficient and small-worldness near-degenerate — and it is also precisely what makes these routines quadratic. A shorter lag or a stricter threshold would address the science and the runtime together.

Where the ephys side's time actually goes (2026-08-28)

The dense-network fixes above are in shared modules (sttc.py, null_models.py, network_metrics.py) and CAT-NAP calls step4.compute_network_metrics itself, so the ephys path got all of them for free. Re-measuring the documented benchmark dataset (2 recordings, 3 lags, 200 shuffles, steps 2-4 from prior spike detection):

Step	2026-07-09 pass	2026-08-28
3. Functional connectivity	42.7s	6.0s
4. Network activity	197.3s	135.9s

The get_sttc hoist is worth 4.2x at 64 channels on real data (26ms -> 6ms per lag, bit-identical), which is ~24s of step 3's drop; it runs 201 times per lag per recording.

What the density fixes could not do much for is the metrics themselves: a real ephys network here is 64 nodes at density 0.10, the opposite end from the 0.86-1.00 that made those routines quadratic. Every network metric for both recordings and all three lags now costs 6.6s. The remaining named candidates (randmio_und_v2/latmio_und_v2 under numba) are 0.19s per recording per 3 lags — worth having at CAT-NAP scale, not here.

The ephys hotspot is the two dimensionality metrics. Same dataset, step 4 only, express, one worker:

Step 4 compute	Time
with NMF + effective rank	69.5s
without	6.6s

i.e. 90% of it, for two fields, one of which (num_nnmf_components) is not even algorithm-identical to MATLAB's (see nmf.py). Per recording: NMF 27.5s, effective rank 4.4s, every network metric 3.7s.

Switching the dimensionality metrics off (compute_nmf / compute_eff_rank, 2026-08-28)

ExtractNetMet.m gates both on Params.netMetToCal; this port declares net_met_to_cal and has never read it. Gating on that list now would be a trap: every params.json ever written carries a copy of its old default, which lists neither metric, so reproducing an old run would silently switch NMF off. Hence two fields of their own, both defaulting to what the pipeline has always done. Measured effect on the benchmark dataset: step 4's compute 69.5s -> 6.6s.

Off means the fields are absent, matching MATLAB (a metric missing from netMetToCal is never assigned) and leaving NaN to mean what it always did — computed, and undefined. The recording-level CSV loses exactly those columns; _plot_violin already skips a metric that is not in the frame.

One behavioural change comes with this. NMF drew from the same generator as the network metrics, so switching it off would have moved Ci/Q/PC/ SW — a metric you remove must not move the metrics you keep. It now draws from make_rng(seed, "step4-dimensionality", filename), the same separation "step4-plots" already used, and CAT-NAP has structurally (its lag-independent metrics run in the producer, its network metrics in a worker with its own stream). The consequence: seeded ephys runs from before this commit reproduce different Ci/Q/PC/SW values, because those metrics no longer inherit a stream that NMF had drawn from first. Nothing deterministic changes, and no MATLAB parity fixture is affected (they cover the deterministic metrics only).

CAT-NAP is untouched: it keeps twop_nmf (default off) and always computes effective rank. The Connectivity tab's new "Network metrics" box says so and goes quiet in CAT-NAP mode, rather than offering a second switch for one metric.

Checked by python/test_step4_optional_metrics.py: the call is genuinely not made, the fields are absent rather than NaN, and every other metric is bit-for-bit what the full run produced.

Express mode — same dataset, figures skipped (2026-08-06)

python/benchmark_express.py re-runs the same dataset and config as the table above (ExampleData/, NGN2_20230208_P1_DIV14_A2/_A3, Axion64, lags [10, 25, 50]ms, ProbThreshRepNum=200) twice — once normally, once with Params.express_mode — so the two sets of numbers are directly comparable. Deliberately not the smaller local/testBurstDetection recording: a 16-channel single recording flatters express mode, because its figure count barely shrinks while its compute does.

Step	Full	Express	Saved
1. Spike detection	69.8s	69.9s	— (untouched)
2. Neuronal activity	7.7s	0.0s	7.7s
3. Functional connectivity	14.6s	14.6s	— (untouched)
4. Network activity	159.0s	113.4s	45.6s
Total	251.1s	198.0s	53.1s (21%)

	figures	folder	bundle
Full	483	56.4 MB	—
Express	6	3.2 MB	2.2 MB

Size is the real win — 17x off the folder, 25x off it as a single shareable file. Time is a secondary benefit and bounded: steps 1 and 3 draw almost nothing, so 84.4s of the 251.1s can never be saved no matter how many figures are skipped. Step 4 keeps 113.4s of compute (null models, NMF, modularity) after its plotting is removed. Expect ~20%, not the ~36% a small single-recording dataset suggests.

The 6 remaining figures are the spike-detection checks (3 per recording) — the only family that can't be rebuilt, since it needs the raw voltage.

Verified alongside the timings, not assumed:

• both runs write byte-identical CSVs at the recording and node level for steps 2 and 4 — express changes what is drawn, never what is computed;
• all 66 per-recording 4A figures (2 recordings × 3 lags × 11 figures) re-render from the bundle pixel-identical to the full run's, in 13.6s;
• the bundle is a complete resume artifact: step 4 re-run from the .meanap alone, with raw_data pointed at a non-existent folder, reproduced the full run's metrics exactly. The ephys spike .npz is small enough (~0.3 MB per recording here) to travel with it.

Measured on a 32-core box, so these absolute numbers are not comparable to the 350.0s in the table above (shared 16-core machine) — the same full run takes 251.1s here. The full-vs-express comparison is the meaningful part, since both halves ran back to back on the same machine.

How to verify changes

# Standalone parity checks against MATLAB reference (no GUI, no example-data path config needed)
uv run python python/test_pipeline_step1.py
uv run python python/test_pipeline_step2.py
uv run python python/test_pipeline_step3.py
uv run python python/test_pipeline_step4.py
uv run python python/test_pipeline_cartography.py
uv run python python/test_pipeline_null_models.py
uv run python python/test_pipeline_channel_layout.py

# Full runner smoke test (creates real output folder tree + runs all 4 steps)
uv run python -c "
from pathlib import Path
from meanap.params import Params
from meanap.pipeline.runner import run_pipeline

p = Params()
p.home_dir = str(Path.cwd())
p.raw_data = str(Path.cwd() / 'ExampleData')
p.spreadsheet_file_name = str(Path.cwd() / 'ExampleData' / 'exampleData.csv')
p.spreadsheet_range = 'A2:A3'
p.output_data_folder = '/tmp/some_test_dir'
p.output_data_folder_name = 'OutputDataTest'
p.start_analysis_step = 1
p.stop_analysis_step = 4
print(run_pipeline(p, log=print))
"

# GUI smoke test
uv run meanap-gui
# Data tab: set "MEA-NAP folder" to the repo root, then Run tab → 🧪 Test pipeline

Note: full spike detection on the two example recordings (64 channels each, wavelet CWT) takes several minutes on a single core — it's CPU-bound, not stuck. Step 3 with the default prob_thresh_rep_num=200 adds another minute-ish per lag value on top of that (see "STTC gotchas" above). Step 4's normalized participation coefficient adds ~15-35s per (recording, lag) on top of that too (100 degree-preserving randomizations — see null_models.py "Performance" note) — a full 2-recording × 3-lag step 4 run alone takes ~2-4 minutes; don't assume a hang.

If MATLAB is available on the machine (which matlab), the step 3/4 fixtures can be regenerated with:

matlab -batch "run('python/test_fixtures/gen_sttc_reference.m')"
matlab -batch "run('python/test_fixtures/gen_step4_reference.m')"
matlab -batch "run('python/test_fixtures/gen_cartography_reference.m')"
matlab -batch "run('python/test_fixtures/gen_coords_reference.m')"
# then convert the resulting .mat files to .npz — see the comment at the
# top of gen_sttc_reference.m for the exact conversion snippet.