cajigaslab · iahncajigas · Mar 4, 2026 · Mar 4, 2026 · Mar 4, 2026 · Mar 4, 2026
diff --git a/.github/workflows/performance-parity.yml b/.github/workflows/performance-parity.yml
@@ -59,10 +59,11 @@ jobs:
             --python-report output/performance/python_performance_report.json \
             --matlab-report tests/performance/fixtures/matlab/performance_baseline_470fde8.json \
             --policy parity/performance_gate_policy.yml \
-            --previous-python-report tests/performance/fixtures/python/performance_baseline_20260303.json \
+            --previous-python-report tests/performance/fixtures/python/performance_baseline_linux_20260304.json \
             --report-out output/performance/performance_parity_report.json \
             --csv-out output/performance/performance_parity_report.csv \
-            --fail-on-regression
+            --fail-on-regression \
+            --require-regression-env-match
 
       - name: Run pytest-benchmark smoke suite
         env:
@@ -80,5 +81,5 @@ jobs:
             output/performance/*.json
             output/performance/*.csv
             tests/performance/fixtures/matlab/performance_baseline_470fde8.json
-            tests/performance/fixtures/python/performance_baseline_20260303.json
+            tests/performance/fixtures/python/performance_baseline_linux_20260304.json
           if-no-files-found: warn
diff --git a/README.md b/README.md
@@ -187,10 +187,11 @@ python tools/performance/compare_matlab_python_performance.py \
   --python-report output/performance/python_performance_report.json \
   --matlab-report tests/performance/fixtures/matlab/performance_baseline_470fde8.json \
   --policy parity/performance_gate_policy.yml \
-  --previous-python-report tests/performance/fixtures/python/performance_baseline_20260303.json \
+  --previous-python-report tests/performance/fixtures/python/performance_baseline_linux_20260304.json \
   --report-out parity/performance_parity_report.json \
   --csv-out parity/performance_parity_report.csv \
-  --fail-on-regression
+  --fail-on-regression \
+  --require-regression-env-match
 ```
 
 Generate MATLAB baseline report (controlled environment):

diff --git a/parity/CYCLE_VALIDATION_CHECKLIST.md b/parity/CYCLE_VALIDATION_CHECKLIST.md
@@ -0,0 +1,47 @@
+# Cycle Validation Checklist (2026-03-04)
+
+Commands used each cycle:
+- `pytest -q`
+- `python tools/parity/build_numeric_drift_report.py --fixtures-manifest tests/parity/fixtures/matlab_gold/manifest.yml --thresholds parity/numeric_drift_thresholds.yml --report-out parity/numeric_drift_report.json --fail-on-violation`
+- `python tools/parity/check_functional_parity_progress.py --report parity/function_example_alignment_report.json --policy parity/functional_gate_policy.yml`
+- `python tools/parity/check_example_output_spec.py --report parity/function_example_alignment_report.json --spec parity/example_output_spec.yml`
+- `python tools/reports/generate_validation_pdf.py --repo-root "$PWD" --matlab-help-root /Users/iahncajigas/Library/CloudStorage/Dropbox/Research/Matlab/nSTAT_currentRelease_Local/helpfiles --notebook-group all --timeout 900 --skip-command-tests --parity-mode gate --enforce-unique-images --min-unique-images-per-topic 1 --max-cross-topic-reuse-ratio 1.0`
+- `python tools/reports/generate_validation_pdf.py --repo-root "$PWD" --matlab-help-root /Users/iahncajigas/Library/CloudStorage/Dropbox/Research/Matlab/nSTAT_currentRelease_Local/helpfiles --notebook-group all --timeout 900 --skip-command-tests --parity-mode image --skip-parity-check`
+- `python tools/reports/build_image_parity_pdfs.py --report-json <latest-json> --python-out output/pdf/image_mode_parity/python_pages.pdf --matlab-out output/pdf/image_mode_parity/matlab_pages.pdf --pairs-json output/pdf/image_mode_parity/pairs.json`
+- `python tools/reports/check_pdf_image_parity.py --python-pdf output/pdf/image_mode_parity/python_pages.pdf --matlab-pdf output/pdf/image_mode_parity/matlab_pages.pdf --out-dir output/pdf/image_mode_parity --dpi 150 --ssim-threshold 0.70 --max-failing-pages 0`
+- `python tools/performance/run_python_benchmarks.py --tiers S --repeats 5 --warmup 1 --out-json output/performance/python_performance_report.json --out-csv output/performance/python_performance_report.csv`
+- `python tools/performance/compare_matlab_python_performance.py --python-report output/performance/python_performance_report.json --matlab-report tests/performance/fixtures/matlab/performance_baseline_470fde8.json --policy parity/performance_gate_policy.yml --previous-python-report tests/performance/fixtures/python/performance_baseline_linux_20260304.json --report-out output/performance/performance_parity_report.json --csv-out output/performance/performance_parity_report.csv --fail-on-regression --require-regression-env-match`
+- Local macOS reruns use `tests/performance/fixtures/python/performance_baseline_20260303.json` with the same command to satisfy strict env matching.
+
+## Cycle 1
+- Log: `output/cycle/cycle1.log`
+- `pytest`: PASS
+- numeric drift (0 failed topics): PASS
+- functional parity (no gaps/partials): PASS
+- example output spec: PASS
+- gate-mode validation PDF (0 parity failures, 0 uniqueness violations): PASS
+- image-mode parity (0 failing pages): PASS
+- performance-parity (0 regression failures): PASS
+- Fixes applied in cycle: comparator option to require regression env match + regression test coverage.
+
+## Cycle 2
+- Log: `output/cycle/cycle2.log`
+- `pytest`: PASS
+- numeric drift (0 failed topics): PASS
+- functional parity (no gaps/partials): PASS
+- example output spec: PASS
+- gate-mode validation PDF (0 parity failures, 0 uniqueness violations): PASS
+- image-mode parity (0 failing pages): PASS
+- performance-parity (0 regression failures): PASS
+- Fixes applied in cycle: Linux baseline + strict regression env matching in workflow/tests, decoding `computeSpikeRateCIs` vectorization, and added deterministic performance workloads for `nspikeTrain.getSigRep` and `Analysis.fitGLM`.
+
+## Cycle 3
+- Log: `output/cycle/cycle3.log`
+- `pytest`: PASS
+- numeric drift (0 failed topics): PASS
+- functional parity (no gaps/partials): PASS
+- example output spec: PASS
+- gate-mode validation PDF (0 parity failures, 0 uniqueness violations): PASS
+- image-mode parity (0 failing pages): PASS
+- performance-parity (0 regression failures): PASS
+- Fixes applied in cycle: none required; full acceptance suite rerun clean after Cycle 2 changes.
diff --git a/src/nstat/compat/matlab/__init__.py b/src/nstat/compat/matlab/__init__.py
@@ -3383,36 +3383,37 @@ def _compute_spike_rate_cis_matlab(
                 chol_m = DecodingAlgorithms._chol_like_matlab(Wku_temp)
                 if chol_m.shape != (K, K):
                     raise ValueError("Wku covariance slice must be KxK")
-            for c in range(int(Mc)):
-                z = rng.normal(0.0, 1.0, size=(K,))
-                xK_draw[r, :, c] = xK_arr[r, :] + (chol_m @ z)
+            # Preserve MATLAB-parity draw ordering by sampling (Mc, K), where each row is one Monte-Carlo draw.
+            z_draw = rng.normal(0.0, 1.0, size=(int(Mc), K))
+            xK_draw[r, :, :] = xK_arr[r, :][:, None] + (chol_m @ z_draw.T)
 
-        lambda_delta = np.zeros((dN_arr.shape[1], K, int(Mc)), dtype=float)
         spike_rate = np.zeros((int(Mc), K), dtype=float)
-        for c in range(int(Mc)):
+        mask = (time >= float(t0)) & (time <= float(tf))
+        interval = max(float(tf - t0), np.finfo(float).eps)
+        integrate_fn = getattr(np, "trapezoid", None)
+        if integrate_fn is None:
+            integrate_fn = getattr(np, "trapz", None)  # pragma: no cover - NumPy<2 fallback
+
+        if window_vals.size > 0 and np.any(np.abs(gamma_vec) > 0.0):
+            hist_term = np.zeros((dN_arr.shape[1], K), dtype=float)
             for k in range(K):
-                stim_k = basis_mat @ xK_draw[:, k, c]
-                if window_vals.size > 0 and np.any(np.abs(gamma_vec) > 0.0):
-                    hk = Hk[k]
-                    cols = min(hk.shape[1], gamma_vec.size)
-                    hist_lin = hk[:, :cols] @ gamma_vec[:cols]
-                else:
-                    hist_lin = np.zeros(stim_k.shape[0], dtype=float)
-                eta = stim_k + hist_lin
-                if fit_type == "poisson":
-                    lam = np.exp(eta)
-                else:
-                    exp_eta = np.exp(eta)
-                    lam = exp_eta / (1.0 + exp_eta)
-                lambda_delta[:, k, c] = lam
-            rates = lambda_delta[:, :, c] / float(delta)
-            mask = (time >= float(t0)) & (time <= float(tf))
+                hk = Hk[k]
+                cols = min(hk.shape[1], gamma_vec.size)
+                hist_term[:, k] = hk[:, :cols] @ gamma_vec[:cols]
+        else:
+            hist_term = np.zeros((dN_arr.shape[1], K), dtype=float)
+
+        for c in range(int(Mc)):
+            stim_ck = basis_mat @ xK_draw[:, :, c]
+            eta = stim_ck + hist_term
+            if fit_type == "poisson":
+                rates = np.exp(eta) / float(delta)
+            else:
+                exp_eta = np.exp(eta)
+                rates = (exp_eta / (1.0 + exp_eta)) / float(delta)
             if np.sum(mask) < 2:
                 integral_vals = np.zeros(K, dtype=float)
             else:
-                integrate_fn = getattr(np, "trapezoid", None)
-                if integrate_fn is None:
-                    integrate_fn = getattr(np, "trapz", None)  # pragma: no cover - NumPy<2 fallback
                 if integrate_fn is None:  # pragma: no cover - extreme fallback
                     dt_vec = np.diff(time[mask]).reshape(-1, 1)
                     y0 = rates[mask, :][:-1, :]
@@ -3423,7 +3424,7 @@ def _compute_spike_rate_cis_matlab(
                         integrate_fn(rates[mask, :], x=time[mask], axis=0),
                         dtype=float,
                     )
-            spike_rate[c, :] = integral_vals / max(float(tf - t0), np.finfo(float).eps)
+            spike_rate[c, :] = integral_vals / interval
 
         CIs = np.zeros((K, 2), dtype=float)
         for k in range(K):
@@ -3451,10 +3452,9 @@ def _compute_spike_rate_cis_matlab(
         ci_obj.setColor("b")
         spike_rate_sig.setConfInterval(ci_obj)
 
-        prob_mat = np.zeros((K, K), dtype=float)
-        for k in range(K):
-            for m in range(k + 1, K):
-                prob_mat[k, m] = float(np.sum(spike_rate[:, m] > spike_rate[:, k])) / float(Mc)
+        # prob_mat(k,m) = P(rate_m > rate_k), with MATLAB-style upper-triangle usage.
+        prob_full = np.mean(spike_rate[:, None, :] > spike_rate[:, :, None], axis=0)
+        prob_mat = np.triu(np.asarray(prob_full, dtype=float), k=1)
         sig_mat = (prob_mat > (1.0 - float(alphaVal))).astype(float)
         return spike_rate_sig, prob_mat, sig_mat
 

diff --git a/src/nstat/performance_workloads.py b/src/nstat/performance_workloads.py
@@ -7,7 +7,7 @@
 
 import numpy as np
 
-from nstat.compat.matlab import CIF, Covariate, DecodingAlgorithms, History, nstColl
+from nstat.compat.matlab import Analysis, CIF, Covariate, DecodingAlgorithms, History, nspikeTrain, nstColl
 
 
 TIER_ORDER = ("S", "M", "L")
@@ -17,6 +17,8 @@
     "history_design_matrix",
     "simulate_cif_thinning",
     "decoding_spike_rate_cis",
+    "nspiketrain_get_sigrep",
+    "analysis_fit_glm_pipeline",
 )
 
 
@@ -36,6 +38,10 @@ class CaseConfig:
     n_bins: int = 120
     mc_draws: int = 30
     decode_delta_s: float = 0.01
+    sigrep_bin_s: float = 0.001
+    glm_n_samples: int = 1000
+    glm_n_features: int = 6
+    glm_dt_s: float = 0.001
 
 
 def get_case_config(case: str, tier: str) -> CaseConfig:
@@ -74,6 +80,18 @@ def get_case_config(case: str, tier: str) -> CaseConfig:
             "M": dict(num_basis=6, num_trials=8, n_bins=200, mc_draws=50, decode_delta_s=0.01),
             "L": dict(num_basis=8, num_trials=12, n_bins=320, mc_draws=80, decode_delta_s=0.01),
         }
+    elif case == "nspiketrain_get_sigrep":
+        vals = {
+            "S": dict(n_spikes=800, duration_s=2.0, sigrep_bin_s=0.002),
+            "M": dict(n_spikes=3000, duration_s=3.0, sigrep_bin_s=0.001),
+            "L": dict(n_spikes=9000, duration_s=5.0, sigrep_bin_s=0.001),
+        }
+    elif case == "analysis_fit_glm_pipeline":
+        vals = {
+            "S": dict(glm_n_samples=900, glm_n_features=6, glm_dt_s=0.001),
+            "M": dict(glm_n_samples=1800, glm_n_features=8, glm_dt_s=0.001),
+            "L": dict(glm_n_samples=3200, glm_n_features=10, glm_dt_s=0.001),
+        }
     else:
         raise ValueError(f"Unknown case: {case}")
 
@@ -101,6 +119,23 @@ def _deterministic_decode_inputs(cfg: CaseConfig) -> tuple[np.ndarray, np.ndarra
     return xk, wku, d_n
 
 
+def _deterministic_glm_inputs(cfg: CaseConfig) -> tuple[np.ndarray, np.ndarray]:
+    n = int(cfg.glm_n_samples)
+    p = int(cfg.glm_n_features)
+    t = np.linspace(0.0, 1.0, n, dtype=float)
+    X = np.zeros((n, p), dtype=float)
+    for j in range(p):
+        f = float(j + 1)
+        X[:, j] = np.sin(2.0 * np.pi * f * t) + 0.35 * np.cos(2.0 * np.pi * (f + 0.5) * t)
+
+    beta = np.linspace(-0.25, 0.30, p, dtype=float)
+    eta = -2.0 + X @ beta
+    mu = np.exp(np.clip(eta, -25.0, 25.0)) * float(cfg.glm_dt_s)
+    phase = np.sin(np.arange(n, dtype=float) * 0.071) + 1.0
+    y = np.floor(mu + 0.35 * phase).astype(float)
+    return X, y
+
+
 def run_python_workload(case: str, tier: str, seed: int = 20260303) -> dict[str, float]:
     """Execute one deterministic Python workload and return summary metrics."""
 
@@ -183,4 +218,26 @@ def run_python_workload(case: str, tier: str, seed: int = 20260303) -> dict[str,
             "rate_mean": float(np.mean(rate)),
         }
 
+    if case == "nspiketrain_get_sigrep":
+        spikes = _deterministic_spike_times(cfg.n_spikes, cfg.duration_s)
+        train = nspikeTrain(spikes, t_start=0.0, t_end=float(cfg.duration_s), name="perf_unit")
+        sig_binary = np.asarray(train.getSigRep(binSize_s=cfg.sigrep_bin_s, mode="binary"), dtype=float)
+        sig_count = np.asarray(train.getSigRep(binSize_s=cfg.sigrep_bin_s, mode="count"), dtype=float)
+        return {
+            "n_bins": float(sig_binary.size),
+            "binary_sum": float(np.sum(sig_binary)),
+            "count_sum": float(np.sum(sig_count)),
+        }
+
+    if case == "analysis_fit_glm_pipeline":
+        X, y = _deterministic_glm_inputs(cfg)
+        fit = Analysis.fitGLM(X=X, y=y, fitType="poisson", dt=float(cfg.glm_dt_s))
+        pred = np.asarray(fit.predict(X), dtype=float)
+        return {
+            "coeff_norm": float(np.linalg.norm(fit.coefficients)),
+            "intercept": float(fit.intercept),
+            "log_likelihood": float(fit.log_likelihood),
+            "pred_mean": float(np.mean(pred)),
+        }
+
     raise ValueError(f"Unhandled workload case: {case}")
diff --git a/tests/performance/fixtures/python/performance_baseline_linux_20260304.csv b/tests/performance/fixtures/python/performance_baseline_linux_20260304.csv
@@ -0,0 +1,16 @@
+case,tier,repeats,median_runtime_ms,mean_runtime_ms,std_runtime_ms,median_peak_memory_mb,summary
+unit_impulse_basis,S,7,1.7549330000008467,1.767814000002334,0.03712841168356124,0.39142704010009766,"{""cols"": 50.0, ""rows"": 501.0, ""total_mass"": 500.0}"
+unit_impulse_basis,M,7,5.230034999996747,5.244222714285992,0.04359156024491504,3.076793670654297,"{""cols"": 100.0, ""rows"": 2001.0, ""total_mass"": 2000.0}"
+unit_impulse_basis,L,7,10.598025999996707,10.61451614285959,0.03639945260812847,18.373775482177734,"{""cols"": 200.0, ""rows"": 6001.0, ""total_mass"": 6000.0}"
+covariate_resample,S,7,0.3187809999900537,0.31151314285742565,0.027315777992366848,0.06198883056640625,"{""cols"": 1.0, ""rows"": 1001.0, ""signal_energy"": 0.5195204795204795}"
+covariate_resample,M,7,0.3733120000077861,0.3833688571432958,0.027679721024068783,0.1355915069580078,"{""cols"": 1.0, ""rows"": 3001.0, ""signal_energy"": 0.5198042747802833}"
+covariate_resample,L,7,0.4388640000030364,0.43310742857321266,0.016053381349413916,0.2376575469970703,"{""cols"": 1.0, ""rows"": 6001.0, ""signal_energy"": 0.5199200133311115}"
+history_design_matrix,S,7,0.34009099999821046,0.3473752857193598,0.02510644529088935,0.0890655517578125,"{""cols"": 4.0, ""rows"": 1000.0, ""total_count"": 9737.0}"
+history_design_matrix,M,7,0.7935309999993478,0.8029685714266829,0.02089668436143026,0.43688201904296875,"{""cols"": 4.0, ""rows"": 5000.0, ""total_count"": 243740.0}"
+history_design_matrix,L,7,1.6487759999961327,1.6487702857140059,0.03152991976246771,0.8869400024414062,"{""cols"": 4.0, ""rows"": 10000.0, ""total_count"": 1462420.0}"
+simulate_cif_thinning,S,7,11.90902599999788,11.841307428570401,0.2313883554103292,0.07111740112304688,"{""mean_spikes_per_unit"": 13.8, ""num_units"": 5.0, ""total_spikes"": 69.0}"
+simulate_cif_thinning,M,7,46.7842029999872,46.623392285716086,0.34392725136599767,0.13444900512695312,"{""mean_spikes_per_unit"": 23.5, ""num_units"": 10.0, ""total_spikes"": 235.0}"
+simulate_cif_thinning,L,7,138.86579399999732,138.25536671428398,2.017449271316567,0.20075607299804688,"{""mean_spikes_per_unit"": 36.65, ""num_units"": 20.0, ""total_spikes"": 733.0}"
+decoding_spike_rate_cis,S,7,29.966428000008705,29.846530285713666,0.3328735707401364,0.2328357696533203,"{""num_trials"": 6.0, ""prob_mean"": 0.1509259259259259, ""rate_mean"": 50.4457886636761, ""sig_count"": 0.0}"
+decoding_spike_rate_cis,M,7,62.59277799999552,63.25780771428567,1.3782261025910456,0.7640476226806641,"{""num_trials"": 8.0, ""prob_mean"": 0.18562499999999998, ""rate_mean"": 50.12398439148756, ""sig_count"": 0.0}"
+decoding_spike_rate_cis,L,7,138.46276700000715,138.051728285717,1.750384804578685,2.7336788177490234,"{""num_trials"": 12.0, ""prob_mean"": 0.21328124999999998, ""rate_mean"": 50.073736692667104, ""sig_count"": 0.0}"