Fix CIFModel.from_linear_terms to use Poisson log link

.github/workflows/ci.yml

@@ -233,3 +233,31 @@ jobs:
           python -m pip install -e .[dev]
       - name: Verify audited MATLAB-facing runtime surface
         run: pytest -q tests/test_matlab_symbol_surface.py tests/test_class_fidelity_audit.py tests/test_parity_report.py
+
+  regenerate-figures:
+    runs-on: ubuntu-latest
+    env:
+      NSTAT_DATA_DIR: ${{ github.workspace }}/.nstat_data_cache
+
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+      - name: Cache example data
+        uses: actions/cache@v4
+        with:
+          path: ${{ env.NSTAT_DATA_DIR }}
+          key: nstat-example-data-v1
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          python -m pip install -e .[dev]
+      - name: Regenerate all paper figures
+        run: python examples/paper/regenerate_all_figures.py
+      - name: Upload regenerated figures
+        uses: actions/upload-artifact@v4
+        with:
+          name: paper-figures
+          path: docs/figures/
+          retention-days: 30

examples/paper/example01_mepsc_poisson.py

@@ -140,12 +140,36 @@ def run_example01(*, export_figures: bool = False, export_dir: Path | None = Non
     print(f"  AIC: {resultConst.AIC}")
     print(f"  BIC: {resultConst.BIC}")
 
-    # --- Figure 1: Constant Mg2+ diagnostics (Matlab-matching plotResults) ---
-    # Matlab calls resultConst.plotResults which creates a 2x4 grid:
-    #   [KSPlot (2-wide)]  [InvGausTrans]  [SeqCorr]
-    #   [plotCoeffs (2-wide)]               [plotResidual (2-wide)]
-    fig1 = plt.figure(figsize=(14, 9))
-    resultConst.plotResults(handle=fig1)
+    # --- Figure 1: Constant Mg2+ diagnostics (Matlab-matching 2x2 layout) ---
+    # Matlab uses subplot(2,2,...) with: raster, InvGausTrans, KSPlot, lambda
+    fig1, axes1 = plt.subplots(2, 2, figsize=(14, 9))
+
+    # (2,2,1): Neural raster
+    ax = axes1[0, 0]
+    spikeCollConst.plot(handle=ax)
+    ax.set_title("Neural Raster with constant Mg$^{2+}$ Concentration",
+                 fontweight="bold", fontsize=12)
+    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
+    ax.set_ylabel("mEPSCs", fontsize=12, fontweight="bold")
+    ax.set_yticks([0, 1])
+
+    # (2,2,2): Inverse Gaussian transform (ACF)
+    resultConst.plotInvGausTrans(fit_num=None, handle=axes1[0, 1])
+
+    # (2,2,3): KS plot
+    resultConst.KSPlot(fit_num=None, handle=axes1[1, 0])
+
+    # (2,2,4): Lambda estimate
+    ax = axes1[1, 1]
+    lam = resultConst.lambda_signal
+    ax.plot(np.asarray(lam.time, dtype=float),
+            np.asarray(lam.data[:, 0], dtype=float),
+            "b", linewidth=2)
+    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
+    ax.set_ylabel(lam.ylabel if hasattr(lam, "ylabel") else "spikes/sec",
+                  fontsize=12, fontweight="bold")
+    ax.legend(["$\\lambda_{const}$"], loc="upper right")
+
     fig1.suptitle("Example 01 — Figure 1: Constant Mg$^{2+}$ Summary",
                   fontsize=14, fontweight="bold")
     fig1.tight_layout()
@@ -229,10 +253,39 @@ def run_example01(*, export_figures: bool = False, export_dir: Path | None = Non
     print(f"  AIC: {resultWashout.AIC}")
     print(f"  BIC: {resultWashout.BIC}")
 
-    # --- Figure 3: Piecewise model diagnostics (Matlab-matching plotResults) ---
-    # Matlab calls resultWashout.plotResults which creates the same 2x4 grid.
-    fig3 = plt.figure(figsize=(14, 9))
-    resultWashout.plotResults(handle=fig3)
+    # --- Figure 3: Piecewise model diagnostics (Matlab-matching 2x2 layout) ---
+    # Matlab uses subplot(2,2,...) with: raster+epoch lines, InvGausTrans, KSPlot, lambda comparison
+    fig3, axes3 = plt.subplots(2, 2, figsize=(14, 9))
+
+    # (2,2,1): Neural raster with epoch boundary lines
+    ax = axes3[0, 0]
+    spikeCollWashout.plot(handle=ax)
+    ax.set_title("Neural Raster with decreasing Mg$^{2+}$ Concentration",
+                 fontweight="bold", fontsize=12)
+    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
+    ax.set_yticklabels([])
+    ax.axvline(495, color="r", linewidth=4)
+    ax.axvline(765, color="r", linewidth=4)
+
+    # (2,2,2): Inverse Gaussian transform (ACF) — all fits
+    resultWashout.plotInvGausTrans(fit_num=None, handle=axes3[0, 1])
+
+    # (2,2,3): KS plot — all fits
+    resultWashout.KSPlot(fit_num=None, handle=axes3[1, 0])
+
+    # (2,2,4): Lambda comparison (two models overlaid)
+    ax = axes3[1, 1]
+    lam = resultWashout.lambda_signal
+    t = np.asarray(lam.time, dtype=float)
+    ax.plot(t, np.asarray(lam.data[:, 0], dtype=float), "b", linewidth=2)
+    if lam.data.shape[1] > 1:
+        ax.plot(t, np.asarray(lam.data[:, 1], dtype=float), "g", linewidth=2)
+    ax.set_ylim(0, 5)
+    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
+    ax.set_ylabel(lam.ylabel if hasattr(lam, "ylabel") else "spikes/sec",
+                  fontsize=12, fontweight="bold")
+    ax.legend(["$\\lambda_{const}$", "$\\lambda_{const-epoch}$"], loc="upper right")
+
     fig3.suptitle("Example 01 — Figure 3: Piecewise Baseline Comparison",
                   fontsize=14, fontweight="bold")
     fig3.tight_layout()

examples/paper/example03_psth_and_ssglm.py

@@ -494,45 +494,30 @@ def run_part_b(data_dir, export_dir=None):
     estStimEffect = np.exp(basisMat @ xK) / delta  # (T, K)
 
     # ------------------------------------------------------------------
-    # Figure 5: True/PSTH/SSGLM stimulus effect surfaces (3D mesh)
+    # Figure 5: True/PSTH/SSGLM stimulus effect surfaces
+    # Match MATLAB: mesh(trial, time, data) with view([90 -90]) → top-down
+    # Using pcolormesh for clean 2D rendering matching MATLAB's top-down mesh
     # ------------------------------------------------------------------
-    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
-
-    fig5 = plt.figure(figsize=(10, 12))
+    fig5, axes5 = plt.subplots(3, 1, figsize=(14, 9))
     trial_axis = np.arange(1, numRealizations + 1)
     T_act = min(actStimEffect.shape[0], len(basis_time))
-    T_mesh, K_mesh = np.meshgrid(
-        basis_time[:T_act], trial_axis, indexing="ij"
-    )
 
-    ax = fig5.add_subplot(3, 1, 1, projection="3d")
-    ax.plot_surface(K_mesh, T_mesh, actStimEffect[:T_act, :],
-                    cmap="viridis", edgecolor="none")
-    ax.view_init(elev=-90, azim=90)
-    ax.set_xticks([])
-    ax.set_yticks([])
-    ax.set_title("True Stimulus Effect", fontweight="bold", fontsize=14)
-
-    ax = fig5.add_subplot(3, 1, 2, projection="3d")
-    ax.plot_surface(K_mesh, T_mesh, psthSurface2D[:T_act, :],
-                    cmap="viridis", edgecolor="none")
-    ax.view_init(elev=-90, azim=90)
-    ax.set_xticks([])
-    ax.set_yticks([])
-    ax.set_title("PSTH Estimated Stimulus Effect", fontweight="bold",
-                 fontsize=14)
-
-    ax = fig5.add_subplot(3, 1, 3, projection="3d")
-    ax.plot_surface(K_mesh, T_mesh, estStimEffect[:T_act, :],
-                    cmap="viridis", edgecolor="none")
-    ax.view_init(elev=-90, azim=90)
-    ax.set_xticks([])
-    ax.set_yticks([])
-    ax.set_title("SSGLM Estimated Stimulus Effect", fontweight="bold",
-                 fontsize=14)
+    surfaces = [
+        (actStimEffect[:T_act, :], "True Stimulus Effect"),
+        (psthSurface2D[:T_act, :], "PSTH Estimated Stimulus Effect"),
+        (estStimEffect[:T_act, :], "SSGLM Estimated Stimulus Effect"),
+    ]
+    for ax, (data, title_str) in zip(axes5, surfaces):
+        # data is (T, K) — time on rows, trials on columns
+        # MATLAB: imagesc shows trial on x, time on y
+        ax.pcolormesh(trial_axis, basis_time[:T_act], data, cmap="viridis",
+                      shading="auto")
+        ax.set_ylabel("time [s]")
+        ax.set_title(title_str, fontweight="bold", fontsize=14)
+    axes5[-1].set_xlabel("Trial [k]")
 
     fig5.tight_layout()
-    print("  Figure 5: Stimulus effect surfaces (3D mesh)")
+    print("  Figure 5: Stimulus effect surfaces (top-down mesh)")
 
     # ------------------------------------------------------------------
     # 6. Learning-trial analysis: spike rate CIs

examples/paper/example04_place_cells_continuous_stimulus.py

@@ -269,19 +269,19 @@ def run_example04(*, export_figures: bool = False, export_dir: Path | None = Non
     fig2, axes2 = plt.subplots(1, 3, figsize=(14, 5))
 
     axes2[0].boxplot([dKS1[np.isfinite(dKS1)], dKS2[np.isfinite(dKS2)]],
-                     tick_labels=["Animal 1", "Animal 2"])
+                     labels=["Animal 1", "Animal 2"])
     axes2[0].set_ylabel(r"$\Delta$KS (Gaussian - Zernike)")
     axes2[0].set_title("KS Statistic Difference")
     axes2[0].axhline(0, color="gray", linestyle="--", linewidth=0.5)
 
     axes2[1].boxplot([dAIC1[np.isfinite(dAIC1)], dAIC2[np.isfinite(dAIC2)]],
-                     tick_labels=["Animal 1", "Animal 2"])
+                     labels=["Animal 1", "Animal 2"])
     axes2[1].set_ylabel(r"$\Delta$AIC (Zernike - Gaussian)")
     axes2[1].set_title("AIC Difference")
     axes2[1].axhline(0, color="gray", linestyle="--", linewidth=0.5)
 
     axes2[2].boxplot([dBIC1[np.isfinite(dBIC1)], dBIC2[np.isfinite(dBIC2)]],
-                     tick_labels=["Animal 1", "Animal 2"])
+                     labels=["Animal 1", "Animal 2"])
     axes2[2].set_ylabel(r"$\Delta$BIC (Zernike - Gaussian)")
     axes2[2].set_title("BIC Difference")
     axes2[2].axhline(0, color="gray", linestyle="--", linewidth=0.5)

examples/paper/example05_decoding_ppaf_pphf.py

@@ -455,7 +455,7 @@ def _plot_part_b(result):
     fig4, axes4 = plt.subplots(1, 4, figsize=(14, 4))
     for d, (ax, lab) in enumerate(zip(axes4, labels)):
         data = [result["rmse_free"][:, d], result["rmse_goal"][:, d]]
-        bp = ax.boxplot(data, tick_labels=["Free", "Goal"])
+        bp = ax.boxplot(data, labels=["Free", "Goal"])
         ax.set_title(f"RMSE: {lab}")
         ax.set_ylabel("RMSE")
 

examples/paper/regenerate_all_figures.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python3
+"""Regenerate all paper example figures.
+
+Runs each example script with --export-figures and saves PNGs to
+docs/figures/exampleNN/.  Called by CI on every push and can also be
+run locally:
+
+    python examples/paper/regenerate_all_figures.py
+"""
+from __future__ import annotations
+
+import importlib
+import sys
+import traceback
+from pathlib import Path
+
+import matplotlib
+matplotlib.use("Agg")
+
+THIS_DIR = Path(__file__).resolve().parent
+REPO_ROOT = THIS_DIR.parents[1]
+FIGURES_ROOT = REPO_ROOT / "docs" / "figures"
+
+if str(REPO_ROOT) not in sys.path:
+    sys.path.insert(0, str(REPO_ROOT))
+
+EXAMPLES = [
+    ("example01_mepsc_poisson", "example01", "run_example01"),
+    ("example02_whisker_stimulus_thalamus", "example02", "run_example02"),
+    ("example03_psth_and_ssglm", "example03", "run_example03"),
+    ("example04_place_cells_continuous_stimulus", "example04", "run_example04"),
+    ("example05_decoding_ppaf_pphf", "example05", "run_example05"),
+]
+
+
+def main() -> int:
+    # Ensure example data is available
+    from nstat.data_manager import ensure_example_data
+    ensure_example_data(download=True)
+
+    failed = 0
+    for mod_name, dir_name, run_fn_name in EXAMPLES:
+        export_dir = FIGURES_ROOT / dir_name
+        print(f"\n{'='*60}")
+        print(f"  {mod_name}")
+        print(f"{'='*60}")
+        try:
+            mod = importlib.import_module(
+                f"examples.paper.{mod_name}"
+            )
+            run_fn = getattr(mod, run_fn_name)
+            run_fn(export_figures=True, export_dir=export_dir)
+            print(f"  OK: figures saved to {export_dir}")
+        except Exception:
+            traceback.print_exc()
+            failed += 1
+
+    total = len(EXAMPLES)
+    print(f"\n=== Done. {total - failed}/{total} examples succeeded ===")
+    return 1 if failed else 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

nstat/cif.py

@@ -319,9 +319,7 @@ def from_linear_terms(
         name: str = "lambda",
     ) -> "CIFModel":
         eta = intercept + np.asarray(design_matrix, dtype=float) @ np.asarray(coefficients, dtype=float)
-        p = np.exp(np.clip(eta, -20.0, 20.0))
-        p = p / (1.0 + p)
-        rate = p / max(float(dt), 1e-12)
+        rate = np.exp(np.clip(eta, -20.0, 20.0))  # Poisson log link: lambda = exp(eta)
         return cls(np.asarray(time, dtype=float).reshape(-1), rate, name)