Re-measure MATLAB vs Python speed (~7x) and fix benchmark_express for express bundles

python/PIPELINE_PORT_STATUS.md

@@ -2158,7 +2158,45 @@ path, `Params.stimulationMode == 1`, out of this port's scope);
 .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 (2026-07-09)
+## 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.
+
+### 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`,

python/benchmark_express.py

@@ -4,7 +4,7 @@
 
     uv run python python/benchmark_express.py --out /path/for/outputs
 
-Runs the Test Pipeline dataset (``local/testBurstDetection``) through steps 1-4
+Runs the benchmark dataset (``ExampleData/``, recordings A2/A3) through steps 1-4
 twice — once producing every figure, once producing only the data plus the
 spike-detection checks — and reports the per-step durations and the size of
 what each leaves behind. Both runs use the same seed, so the numbers they
@@ -14,9 +14,10 @@
 from __future__ import annotations
 
 import argparse
-import json
+import re
 import sys
 import time
+import zipfile
 from pathlib import Path
 
 REPO_ROOT = Path(__file__).parent.parent
@@ -56,32 +57,44 @@ def _params(out_dir: Path, name: str, *, express: bool) -> Params:
     )
 
 
-def _dir_size(path: Path) -> int:
-    return sum(p.stat().st_size for p in path.rglob("*") if p.is_file())
+#: The runner's timing lines. Read from the log rather than step_durations.json
+#: because express mode deletes the output folder (json included) once its
+#: .meanap bundle is written.
+_DURATION = re.compile(r"^(Step (\d)|Total pipeline) duration \(seconds\): ([\d.]+)")
 
 
-def _count(path: Path, pattern: str) -> int:
-    return sum(1 for _ in path.rglob(pattern))
+def _files(root: Path) -> list[tuple[str, int]]:
+    """(relative path, size) of every file in a run folder or .meanap bundle."""
+    if root.is_dir():
+        return [(p.relative_to(root).as_posix(), p.stat().st_size)
+                for p in root.rglob("*") if p.is_file()]
+    with zipfile.ZipFile(root) as zf:
+        return [(i.filename, i.file_size) for i in zf.infolist() if not i.is_dir()]
 
 
 def _run(out_dir: Path, name: str, *, express: bool) -> dict:
     print(f"\n{'=' * 70}\n{name} (express={express})\n{'=' * 70}", flush=True)
+    durations = {}
+
+    def log(msg: str) -> None:
+        print(msg, flush=True)
+        m = _DURATION.match(str(msg))
+        if m:
+            durations[f"step{m.group(2)}" if m.group(2) else "total"] = float(m.group(3))
+
     t0 = time.perf_counter()
-    root = run_pipeline(_params(out_dir, name, express=express),
-                        log=lambda m: print(m, flush=True))
+    # A full run returns its folder; an express run returns the bundle that
+    # replaced it.
+    root = run_pipeline(_params(out_dir, name, express=express), log=log)
     wall = time.perf_counter() - t0
 
-    durations = {}
-    durations_path = root / "step_durations.json"
-    if durations_path.exists():
-        durations = json.loads(durations_path.read_text())
-
-    bundle = root.with_suffix(".meanap")
+    files = _files(root)
     return {
         "name": name, "root": root, "wall": wall, "durations": durations,
-        "bytes": _dir_size(root),
-        "figures": _count(root, "*.png"),
-        "bundle_bytes": bundle.stat().st_size if bundle.exists() else None,
+        "files": files,
+        "bytes": root.stat().st_size if root.is_file() else sum(s for _, s in files),
+        "figures": sum(1 for f, _ in files if f.endswith(".png")),
+        "bundle_bytes": root.stat().st_size if root.is_file() else None,
     }
 
 
@@ -113,29 +126,32 @@ def main() -> int:
         print(row)
 
     print()
-    print(f"{'':<14}{'figures':>12}{'folder MB':>12}{'bundle MB':>12}")
+    print(f"{'':<14}{'figures':>12}{'output MB':>12}")
     for r in results:
-        b = "-" if r["bundle_bytes"] is None else f"{r['bundle_bytes'] / 1e6:.1f}"
-        print(f"{r['name']:<14}{r['figures']:>12}{r['bytes'] / 1e6:>12.1f}{b:>12}")
+        kind = "bundle" if r["bundle_bytes"] is not None else "folder"
+        print(f"{r['name']:<14}{r['figures']:>12}{r['bytes'] / 1e6:>12.1f}  ({kind})")
 
     full, exp = results
     if full["durations"].get("total") and exp["durations"].get("total"):
         saved = full["durations"]["total"] - exp["durations"]["total"]
         pct = saved / full["durations"]["total"] * 100
         print(f"\ntime saved: {saved:.1f}s ({pct:.0f}%)")
-    print(f"folder size: {full['bytes'] / 1e6:.1f} MB → {exp['bytes'] / 1e6:.1f} MB "
-          f"({full['bytes'] / max(exp['bytes'], 1):.1f}x smaller)")
-    if exp["bundle_bytes"]:
-        print(f"bundle:      {exp['bundle_bytes'] / 1e6:.1f} MB "
-              f"({full['bytes'] / exp['bundle_bytes']:.0f}x smaller than the full folder)")
+    print(f"output size: {full['bytes'] / 1e6:.1f} MB folder → {exp['bytes'] / 1e6:.1f} MB "
+          f"{'bundle' if exp['bundle_bytes'] is not None else 'folder'} "
+          f"({full['bytes'] / max(exp['bytes'], 1):.0f}x smaller)")
 
-    # Per-folder breakdown of where the bytes go.
+    # Per-folder breakdown of where the bytes go (uncompressed, for a bundle).
     print()
     for r in results:
         print(f"{r['name']}:")
-        for sub in sorted(p for p in r["root"].iterdir() if p.is_dir()):
-            print(f"    {sub.name:<28}{_dir_size(sub) / 1e6:>8.1f} MB"
-                  f"{_count(sub, '*.png'):>7} png")
+        tops: dict[str, list[int]] = {}
+        for f, size in r["files"]:
+            if "/" in f:
+                t = tops.setdefault(f.split("/", 1)[0], [0, 0])
+                t[0] += size
+                t[1] += f.endswith(".png")
+        for top, (size, pngs) in sorted(tops.items()):
+            print(f"    {top:<28}{size / 1e6:>8.1f} MB{pngs:>7} png")
     return 0
 
 

python/parityRun/MEApipeline_timingBenchmark.m

@@ -25,7 +25,10 @@ function MEApipeline_timingBenchmark(InputParamsFilePath)
 %   Step 2 duration (seconds): 53
 %   Step 3 duration (seconds): 15
 %   Step 4 duration (seconds): 283
-%   (total 931s) — see python/PIPELINE_PORT_STATUS.md's "MATLAB vs Python
+%   (total 931s). Re-run 2026-10-25 (R2025b, 16-worker parpool):
+%   445 / 62 / 9 / 377, total 893s. Run as
+%   matlab -batch "addpath('<repo>/matlab'); addpath('<repo>/python/parityRun'); MEApipeline_timingBenchmark"
+%   — see python/PIPELINE_PORT_STATUS.md's "MATLAB vs Python
 %   speed comparison" section for the full writeup incl. the Python side.
 %
 % Getting here required finding and working around 5 bugs/config traps in