From 8a22a4fecaefdb3107dfe556ae59d4b04a8a60d9 Mon Sep 17 00:00:00 2001 From: Iahn Cajigas Date: Sun, 8 Mar 2026 00:52:08 -0500 Subject: [PATCH] Tighten source-driven analysis and notebook audit parity --- notebooks/ExplicitStimulusWhiskerData.ipynb | 18 ++--- notebooks/HippocampalPlaceCellExample.ipynb | 14 ++-- notebooks/HybridFilterExample.ipynb | 20 +++--- notebooks/PPSimExample.ipynb | 38 +++++----- notebooks/StimulusDecode2D.ipynb | 18 ++--- notebooks/TrialExamples.ipynb | 41 ++++------- notebooks/ValidationDataSet.ipynb | 48 ++++++++----- nstat/analysis.py | 38 ++++++++-- nstat/matlab_reference.py | 33 +++++++++ nstat/parity_report.py | 18 ++++- parity/class_fidelity.yml | 8 ++- parity/notebook_fidelity.yml | 21 +++--- parity/report.md | 24 +++++-- parity/simulink_fidelity.yml | 68 +++++++++++++++++- .../matlab_gold/analysis_exactness.mat | Bin 0 -> 924 bytes .../fixtures/matlab_gold/cif_exactness.mat | Bin 614 -> 614 bytes .../matlab_gold/nspiketrain_exactness.mat | Bin 1451 -> 1451 bytes .../matlab_gold/point_process_exactness.mat | Bin 366 -> 366 bytes .../matlab_gold/signalobj_exactness.mat | Bin 1310 -> 1310 bytes tests/test_matlab_gold_fixtures.py | 36 +++++++++- tests/test_matlab_reference.py | 23 +++++- tests/test_notebook_fidelity_audit.py | 4 +- tests/test_notebook_parity_notes.py | 4 +- tests/test_parity_report.py | 8 +-- tests/test_simulink_fidelity_audit.py | 17 +++++ .../build_foundational_help_notebooks.py | 19 ++--- .../build_helpfile_fidelity_notebooks.py | 24 +++++-- tools/notebooks/parity_notes.yml | 8 +-- .../matlab/export_matlab_gold_fixtures.m | 27 +++++++ 29 files changed, 415 insertions(+), 162 deletions(-) create mode 100644 tests/parity/fixtures/matlab_gold/analysis_exactness.mat diff --git a/notebooks/ExplicitStimulusWhiskerData.ipynb b/notebooks/ExplicitStimulusWhiskerData.ipynb index 325ba9d7..19c4d8fd 100644 --- a/notebooks/ExplicitStimulusWhiskerData.ipynb +++ b/notebooks/ExplicitStimulusWhiskerData.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "8b8aa493", + "id": "a51db33a", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "dbf8e486", + "id": "9f51d683", "metadata": {}, "outputs": [], "source": [ @@ -83,7 +83,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e296bd4d", + "id": "4ee3a32e", "metadata": {}, "outputs": [], "source": [ @@ -106,7 +106,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6a18af75", + "id": "35afed3f", "metadata": {}, "outputs": [], "source": [ @@ -134,7 +134,7 @@ { "cell_type": "code", "execution_count": null, - "id": "fd207a34", + "id": "214e2ec3", "metadata": {}, "outputs": [], "source": [ @@ -149,7 +149,7 @@ { "cell_type": "code", "execution_count": null, - "id": "2afb535d", + "id": "782d0257", "metadata": {}, "outputs": [], "source": [ @@ -170,7 +170,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7a48f375", + "id": "494c31d8", "metadata": {}, "outputs": [], "source": [ @@ -198,7 +198,7 @@ { "cell_type": "code", "execution_count": null, - "id": "043ef33a", + "id": "b49274b0", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +230,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d66ac872", + "id": "c1129e2d", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/HippocampalPlaceCellExample.ipynb b/notebooks/HippocampalPlaceCellExample.ipynb index 9789a3ac..e5682b19 100644 --- a/notebooks/HippocampalPlaceCellExample.ipynb +++ b/notebooks/HippocampalPlaceCellExample.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "39256c2a", + "id": "04918a27", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "2b552e70", + "id": "cd6ae412", "metadata": {}, "outputs": [], "source": [ @@ -85,7 +85,7 @@ { "cell_type": "code", "execution_count": null, - "id": "02b3faec", + "id": "bc94da8b", "metadata": {}, "outputs": [], "source": [ @@ -105,7 +105,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3e5df761", + "id": "7383932a", "metadata": {}, "outputs": [], "source": [ @@ -125,7 +125,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22988fc3", + "id": "3d6944b1", "metadata": {}, "outputs": [], "source": [ @@ -152,7 +152,7 @@ { "cell_type": "code", "execution_count": null, - "id": "2bb2d885", + "id": "f52ed4fd", "metadata": {}, "outputs": [], "source": [ @@ -179,7 +179,7 @@ { "cell_type": "code", "execution_count": null, - "id": "34d08f71", + "id": "ca3ffc88", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/HybridFilterExample.ipynb b/notebooks/HybridFilterExample.ipynb index f3293504..d4e8eb2d 100644 --- a/notebooks/HybridFilterExample.ipynb +++ b/notebooks/HybridFilterExample.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "deae92f2", + "id": "799ac706", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `HybridFilterExample.mlx`\n", "- Fidelity status: `high_fidelity`\n", - "- Remaining justified differences: The notebook now reproduces the hybrid-filter simulation, single-run decoding, and averaged summary figures with real outputs; the Python port still uses the current hybrid-filter implementation instead of every MATLAB-specific reporting branch." + "- Remaining justified differences: The notebook now reproduces the hybrid-filter simulation, single-run decoding, and averaged summary figures with real outputs; the Python port still uses the current hybrid-filter implementation instead of every MATLAB-specific reporting branch.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "ad07c855", + "id": "ac9d53a9", "metadata": {}, "outputs": [], "source": [ @@ -63,7 +63,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8ad6c4c4", + "id": "bf21531d", "metadata": {}, "outputs": [], "source": [ @@ -87,7 +87,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d3676bdb", + "id": "e005f11e", "metadata": {}, "outputs": [], "source": [ @@ -98,7 +98,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c0ac4ff8", + "id": "f52cdb38", "metadata": {}, "outputs": [], "source": [ @@ -109,7 +109,7 @@ { "cell_type": "code", "execution_count": null, - "id": "20954cd6", + "id": "3c6317b7", "metadata": {}, "outputs": [], "source": [ @@ -161,7 +161,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e59c36e8", + "id": "c0a4b72b", "metadata": {}, "outputs": [], "source": [ @@ -172,7 +172,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1da0a680", + "id": "e28a4b2e", "metadata": {}, "outputs": [], "source": [ @@ -248,4 +248,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/PPSimExample.ipynb b/notebooks/PPSimExample.ipynb index 3ac2cdb2..1199c4b9 100644 --- a/notebooks/PPSimExample.ipynb +++ b/notebooks/PPSimExample.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "7eb3e8aa", + "id": "75db499a", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9b87979a", + "id": "344b225c", "metadata": {}, "outputs": [], "source": [ @@ -71,7 +71,7 @@ { "cell_type": "code", "execution_count": null, - "id": "11e65873", + "id": "e586ec40", "metadata": {}, "outputs": [], "source": [ @@ -82,7 +82,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5a9c6419", + "id": "e7759808", "metadata": {}, "outputs": [], "source": [ @@ -93,7 +93,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5637fb8b", + "id": "1066f0cb", "metadata": {}, "outputs": [], "source": [ @@ -105,7 +105,7 @@ { "cell_type": "code", "execution_count": null, - "id": "79f48bf9", + "id": "c10cdb5a", "metadata": {}, "outputs": [], "source": [ @@ -116,7 +116,7 @@ { "cell_type": "code", "execution_count": null, - "id": "739175fd", + "id": "1bc35d15", "metadata": {}, "outputs": [], "source": [ @@ -127,7 +127,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9a5d334e", + "id": "c6fa075a", "metadata": {}, "outputs": [], "source": [ @@ -143,7 +143,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1cf6bf85", + "id": "2c6f4f13", "metadata": {}, "outputs": [], "source": [ @@ -157,7 +157,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d1587532", + "id": "cb839486", "metadata": {}, "outputs": [], "source": [ @@ -173,7 +173,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4476c77c", + "id": "f728a04e", "metadata": {}, "outputs": [], "source": [ @@ -185,7 +185,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b208e63b", + "id": "c20b4fcd", "metadata": {}, "outputs": [], "source": [ @@ -196,7 +196,7 @@ { "cell_type": "code", "execution_count": null, - "id": "cc4bacb4", + "id": "2058e202", "metadata": {}, "outputs": [], "source": [ @@ -208,7 +208,7 @@ { "cell_type": "code", "execution_count": null, - "id": "81244323", + "id": "76209533", "metadata": {}, "outputs": [], "source": [ @@ -220,7 +220,7 @@ { "cell_type": "code", "execution_count": null, - "id": "da3c3380", + "id": "11b2922d", "metadata": {}, "outputs": [], "source": [ @@ -232,7 +232,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9768bcd9", + "id": "2666d94d", "metadata": {}, "outputs": [], "source": [ @@ -244,7 +244,7 @@ { "cell_type": "code", "execution_count": null, - "id": "91a6eb4a", + "id": "246bce8c", "metadata": {}, "outputs": [], "source": [ @@ -258,7 +258,7 @@ { "cell_type": "code", "execution_count": null, - "id": "626975dc", + "id": "481ca8c7", "metadata": {}, "outputs": [], "source": [ @@ -272,7 +272,7 @@ { "cell_type": "code", "execution_count": null, - "id": "54d342c8", + "id": "dc217235", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/StimulusDecode2D.ipynb b/notebooks/StimulusDecode2D.ipynb index 469531c6..ecff1c14 100644 --- a/notebooks/StimulusDecode2D.ipynb +++ b/notebooks/StimulusDecode2D.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "832f42b3", + "id": "249987ba", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `StimulusDecode2D.mlx`\n", - "- Fidelity status: `high_fidelity`\n", - "- Remaining justified differences: The notebook now reproduces the 2-D stimulus-decoding workflow with simulated receptive fields and decoded trajectories; the current Python decoder uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear filter." + "- Fidelity status: `partial`\n", + "- Remaining justified differences: The notebook reproduces the MATLAB section order, figure inventory, simulated receptive fields, and decoded-trajectory presentation, but the current Python decoder still uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear filter.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "610a2d16", + "id": "3989e7ba", "metadata": {}, "outputs": [], "source": [ @@ -122,7 +122,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7b19bc85", + "id": "e5b0164f", "metadata": {}, "outputs": [], "source": [ @@ -136,7 +136,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c4494ced", + "id": "44c9df66", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +178,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7e531f49", + "id": "94f8bcd5", "metadata": {}, "outputs": [], "source": [ @@ -196,7 +196,7 @@ { "cell_type": "code", "execution_count": null, - "id": "435982bd", + "id": "71c1f658", "metadata": {}, "outputs": [], "source": [ @@ -252,4 +252,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/TrialExamples.ipynb b/notebooks/TrialExamples.ipynb index 566b9389..21ee977e 100644 --- a/notebooks/TrialExamples.ipynb +++ b/notebooks/TrialExamples.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "f4474479", + "id": "728c6277", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `TrialExamples.mlx`\n", "- Fidelity status: `high_fidelity`\n", - "- Remaining justified differences: The notebook now mirrors the MATLAB Trial workflow with executable object construction, masking, history extraction, and plotting; the closing analysis section uses one representative Python `Analysis` run instead of linking out to separate MATLAB help pages.\n" + "- Remaining justified differences: The notebook now mirrors the MATLAB Trial workflow with executable object construction, masking, history extraction, and plotting; only minor Python plotting defaults differ from the published MATLAB help output.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "a61a3753", + "id": "62d9faec", "metadata": {}, "outputs": [], "source": [ @@ -35,7 +35,7 @@ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", - "from nstat import Analysis, ConfigColl, CovColl, Covariate, Events, History, Trial, TrialConfig, nspikeTrain, nstColl\n", + "from nstat import CovColl, Covariate, Events, History, Trial, nspikeTrain, nstColl\n", "from nstat.notebook_figures import FigureTracker\n", "\n", "np.random.seed(7)\n", @@ -124,7 +124,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1d69a420", + "id": "509b204a", "metadata": {}, "outputs": [], "source": [ @@ -140,7 +140,7 @@ { "cell_type": "code", "execution_count": null, - "id": "57a90e50", + "id": "374bcec9", "metadata": {}, "outputs": [], "source": [ @@ -153,7 +153,7 @@ { "cell_type": "code", "execution_count": null, - "id": "46786f1c", + "id": "afc89310", "metadata": {}, "outputs": [], "source": [ @@ -165,7 +165,7 @@ { "cell_type": "code", "execution_count": null, - "id": "521b6cf2", + "id": "0d3c2fa7", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +178,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1209c432", + "id": "7b522976", "metadata": {}, "outputs": [], "source": [ @@ -191,7 +191,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6a2f219e", + "id": "207e8239", "metadata": {}, "outputs": [], "source": [ @@ -203,7 +203,7 @@ { "cell_type": "code", "execution_count": null, - "id": "51befd6a", + "id": "af1d8f85", "metadata": {}, "outputs": [], "source": [ @@ -219,34 +219,23 @@ { "cell_type": "code", "execution_count": null, - "id": "8072393e", + "id": "efa7bea7", "metadata": {}, "outputs": [], "source": [ "# SECTION 8: Example 2: Analyzing Trial Data\n", - "cfg = TrialConfig([[\"Position\", \"x\"], [\"Force\", \"f_x\"]], sampleRate=ctx[\"sample_rate\"], history=[0.0, 0.05, 0.1], name=\"Position+Force+History\")\n", - "cfgColl = ConfigColl([cfg])\n", - "fit = Analysis.RunAnalysisForNeuron(trial1, 1, cfgColl)\n", - "fit_stats = fit.computeKSStats()\n", - "print(\n", - " {\n", - " \"config_name\": fit.configNames[0],\n", - " \"aic\": round(float(fit.AIC[0]), 3),\n", - " \"bic\": round(float(fit.BIC[0]), 3),\n", - " \"ks_stat\": round(float(fit_stats[\"ks_stat\"]), 4),\n", - " }\n", - ")\n" + "print(\"Examples of neural spike analysis using AnalysisExamples2 (Neural Spike Analysis Toolbox) or AnalysisExamples (standard methods).\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "e781b407", + "id": "4d808c62", "metadata": {}, "outputs": [], "source": [ "# SECTION 9: Related analysis workflows\n", - "print(\"For larger model-comparison walkthroughs, continue with AnalysisExamples and AnalysisExamples2.\")\n", + "print({\"recommended_next_notebooks\": [\"AnalysisExamples2\", \"AnalysisExamples\"]})\n", "__tracker.finalize()\n" ] } diff --git a/notebooks/ValidationDataSet.ipynb b/notebooks/ValidationDataSet.ipynb index 3337a575..e6292b4e 100644 --- a/notebooks/ValidationDataSet.ipynb +++ b/notebooks/ValidationDataSet.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "a034081e", + "id": "5d534eeb", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `ValidationDataSet.mlx`\n", "- Fidelity status: `high_fidelity`\n", - "- Remaining justified differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; the Python port uses shorter deterministic simulations than MATLAB so the notebook remains stable in CI." + "- Remaining justified differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; local execution uses the MATLAB-scale simulation sizes, while CI switches to a documented shorter deterministic fast path for stability.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "835e27fe", + "id": "0177a98f", "metadata": {}, "outputs": [], "source": [ @@ -30,6 +30,8 @@ "if str(SRC_PATH) not in sys.path:\n", " sys.path.insert(0, str(SRC_PATH))\n", "\n", + "import os\n", + "\n", "import matplotlib\n", "matplotlib.use(\"Agg\")\n", "import matplotlib.pyplot as plt\n", @@ -58,7 +60,12 @@ " return time, data\n", "\n", "\n", - "def _simulate_constant_case(seed=0, *, p=0.01, n_samples=20001, delta=0.001):\n", + "CI_FAST_PATH = os.environ.get(\"CI\", \"\").strip().lower() in {\"1\", \"true\", \"yes\"}\n", + "\n", + "\n", + "def _simulate_constant_case(seed=0, *, p=0.01, n_samples=None, delta=0.001):\n", + " if n_samples is None:\n", + " n_samples = 20001 if CI_FAST_PATH else 100001\n", " rng = np.random.default_rng(seed)\n", " total_time = n_samples * delta\n", " time = np.linspace(0.0, total_time, n_samples)\n", @@ -85,7 +92,11 @@ " }\n", "\n", "\n", - "def _simulate_piecewise_case(seed=1, *, p1=0.001, p2=0.01, n1=20000, n2=20000, delta=0.001):\n", + "def _simulate_piecewise_case(seed=1, *, p1=0.001, p2=0.01, n1=None, n2=None, delta=0.001):\n", + " if n1 is None:\n", + " n1 = 20000 if CI_FAST_PATH else 100000\n", + " if n2 is None:\n", + " n2 = 20000 if CI_FAST_PATH else 100000\n", " rng = np.random.default_rng(seed)\n", " t1 = np.linspace(0.0, n1 * delta, n1 + 1)\n", " t2 = np.linspace(n1 * delta, (n1 + n2) * delta, n2 + 1)[1:]\n", @@ -143,17 +154,18 @@ { "cell_type": "code", "execution_count": null, - "id": "31374daa", + "id": "13e27c11", "metadata": {}, "outputs": [], "source": [ "# SECTION 0: Software Validation Data Set\n", - "# This notebook follows the MATLAB validation helpfile with deterministic simulations for CI-stable execution.\n", + "# This notebook follows the MATLAB validation helpfile; CI uses a documented short fast path while local runs use MATLAB-scale sample counts.\n", "plt.close(\"all\")\n", "constant_case = _simulate_constant_case()\n", "piecewise_case = _simulate_piecewise_case()\n", "print(\n", " {\n", + " \"ci_fast_path\": CI_FAST_PATH,\n", " \"constant_lambda_hz\": round(float(constant_case[\"lambda_hz\"]), 4),\n", " \"piecewise_lambda1_hz\": round(float(piecewise_case[\"lambda1_hz\"]), 4),\n", " \"piecewise_lambda2_hz\": round(float(piecewise_case[\"lambda2_hz\"]), 4),\n", @@ -164,7 +176,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f56503ae", + "id": "d94296d5", "metadata": {}, "outputs": [], "source": [ @@ -175,7 +187,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3c7479dd", + "id": "b3d520d3", "metadata": {}, "outputs": [], "source": [ @@ -187,7 +199,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5c357cd0", + "id": "8b792ab5", "metadata": {}, "outputs": [], "source": [ @@ -201,7 +213,7 @@ { "cell_type": "code", "execution_count": null, - "id": "15efa1ec", + "id": "681433fa", "metadata": {}, "outputs": [], "source": [ @@ -223,7 +235,7 @@ { "cell_type": "code", "execution_count": null, - "id": "69d4e988", + "id": "1ec098fd", "metadata": {}, "outputs": [], "source": [ @@ -245,7 +257,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ff7b1cfd", + "id": "daf22292", "metadata": {}, "outputs": [], "source": [ @@ -258,7 +270,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d2e52928", + "id": "1286f142", "metadata": {}, "outputs": [], "source": [ @@ -288,7 +300,7 @@ { "cell_type": "code", "execution_count": null, - "id": "52b4802a", + "id": "b35f530e", "metadata": {}, "outputs": [], "source": [ @@ -299,7 +311,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4119e5a0", + "id": "75542bfe", "metadata": {}, "outputs": [], "source": [ @@ -334,7 +346,7 @@ { "cell_type": "code", "execution_count": null, - "id": "a0a82984", + "id": "ad6eec57", "metadata": {}, "outputs": [], "source": [ @@ -378,4 +390,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/nstat/analysis.py b/nstat/analysis.py index 25ec22d8..472d58d1 100644 --- a/nstat/analysis.py +++ b/nstat/analysis.py @@ -82,6 +82,32 @@ def _fit_lambda_matrix_to_covariate(lambda_time: np.ndarray, lambda_columns: lis ) +def _glm_deviance(y: np.ndarray, mean_counts: np.ndarray, distribution: str) -> float: + observed = np.asarray(y, dtype=float).reshape(-1) + expected = np.clip(np.asarray(mean_counts, dtype=float).reshape(-1), 1e-12, None) + if observed.shape != expected.shape: + raise ValueError("observed and expected counts must have matching shapes") + + dist = str(distribution).lower() + if dist == "poisson": + ratio = np.ones_like(observed) + positive = observed > 0.0 + ratio[positive] = observed[positive] / expected[positive] + return float(2.0 * np.sum(observed * np.log(ratio) - (observed - expected))) + + if dist == "binomial": + prob = np.clip(expected, 1e-12, 1.0 - 1e-12) + return float( + 2.0 + * np.sum( + observed * np.log(np.clip(observed / prob, 1e-12, None)) + + (1.0 - observed) * np.log(np.clip((1.0 - observed) / (1.0 - prob), 1e-12, None)) + ) + ) + + raise ValueError(f"Unsupported GLM distribution for deviance: {distribution}") + + def _benjamini_hochberg(p_values: np.ndarray, alpha: float) -> np.ndarray: p = np.asarray(p_values, dtype=float).reshape(-1) if p.size == 0: @@ -157,25 +183,25 @@ def GLMFit( y = np.concatenate(stacked_y) if stacked_y else np.array([], dtype=float) lambda_time_full = np.concatenate(lambda_time_segments) if lambda_time_segments else np.array([], dtype=float) sample_rate = float(tObj.sampleRate) - dt = 1.0 / max(sample_rate, 1e-12) if algorithm == "BNLRCG": glm_res = fit_binomial_glm(X, y, include_intercept=False, l2=l2, max_iter=max_iter) - probability = glm_res.predict_probability(X) - rate_hz = probability * sample_rate + lambda_delta = np.clip(glm_res.predict_probability(X), 1e-12, 1.0 - 1e-9) + rate_hz = lambda_delta * sample_rate distribution = "binomial" b = np.asarray(glm_res.coefficients, dtype=float).reshape(-1) + logLL = float(np.sum(y * np.log(lambda_delta) + (1.0 - y) * np.log(1.0 - lambda_delta))) + dev = _glm_deviance(y, lambda_delta, distribution) else: glm_res = fit_poisson_glm(X, y, include_intercept=False, l2=l2, max_iter=max_iter) lambda_delta = glm_res.predict_rate(X) rate_hz = lambda_delta * sample_rate distribution = "poisson" b = np.asarray(glm_res.coefficients, dtype=float).reshape(-1) + logLL = float(glm_res.log_likelihood) + dev = _glm_deviance(y, lambda_delta, distribution) n_params = int(b.size) - prob = np.clip(rate_hz * dt, 1e-12, 1.0 - 1e-9) - logLL = float(np.sum(y * np.log(prob) + (1.0 - y) * np.log(1.0 - prob))) - dev = float(-2.0 * logLL) AIC = float(2.0 * n_params + dev) BIC = float(np.log(max(y.shape[0], 1)) * n_params + dev) diff --git a/nstat/matlab_reference.py b/nstat/matlab_reference.py index 3d38cfe1..dd2864b3 100644 --- a/nstat/matlab_reference.py +++ b/nstat/matlab_reference.py @@ -119,8 +119,41 @@ def run_simulated_network_reference(*, matlab_repo: str | Path | None = None, se } +def run_analysis_reference(*, matlab_repo: str | Path | None = None) -> dict[str, np.ndarray]: + repo = Path(matlab_repo) if matlab_repo is not None else _default_matlab_repo() + if not repo.exists(): + raise FileNotFoundError(f"MATLAB reference repo not found at {repo}") + engine = start_matlab_engine() + _add_repo_to_path(engine, repo) + engine.eval( + """ + t = (0:0.1:1.0)'; + stim = Covariate(t, sin(2*pi*t), 'Stimulus', 'time', 's', '', {'stim'}); + spikeTrain = nspikeTrain([0.1 0.4 0.7], '1', 0.1, 0.0, 1.0, 'time', 's', '', '', -1); + trial = Trial(nstColl({spikeTrain}), CovColl({stim})); + cfg = TrialConfig({{'Stimulus', 'stim'}}, 10, [], []); + cfg.setName('stim'); + fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl({cfg})); + analysisAIC = fit.AIC(1); + analysisBIC = fit.BIC(1); + analysisLogLL = fit.logLL(1); + analysisCoeffs = fit.getCoeffs(1)'; + analysisLambdaHead = fit.lambda.data(1:5, 1)'; + """, + nargout=0, + ) + return { + "aic": _to_numpy(engine.workspace["analysisAIC"]).reshape(-1), + "bic": _to_numpy(engine.workspace["analysisBIC"]).reshape(-1), + "logll": _to_numpy(engine.workspace["analysisLogLL"]).reshape(-1), + "coeffs": _to_numpy(engine.workspace["analysisCoeffs"]).reshape(-1), + "lambda_head": _to_numpy(engine.workspace["analysisLambdaHead"]).reshape(-1), + } + + __all__ = [ "matlab_engine_available", + "run_analysis_reference", "run_point_process_reference", "run_simulated_network_reference", "start_matlab_engine", diff --git a/nstat/parity_report.py b/nstat/parity_report.py index 3541d2c4..d349bb05 100644 --- a/nstat/parity_report.py +++ b/nstat/parity_report.py @@ -90,6 +90,12 @@ def render_parity_report(repo_root: Path | None = None) -> str: notebook_partial = iter_outstanding_notebook_fidelity(notebook_fidelity) simulink_counts = summarize_simulink_strategies(simulink_fidelity) simulink_outstanding = iter_outstanding_simulink_items(simulink_fidelity) + simulink_reference_only = [ + row + for row in simulink_fidelity.get("items", []) + if str(row.get("current_python_status", "")).strip() == "reference_only" + or str(row.get("python_strategy", "")).strip() == "reference_only" + ] lines = [ "# nSTAT Python Parity Report", "", @@ -181,6 +187,10 @@ def render_parity_report(repo_root: Path | None = None) -> str: lines.append( f"- Simulink fidelity: {len(simulink_outstanding)} Simulink-backed assets still rely on partial, fallback, or unsupported Python execution paths." ) + elif simulink_reference_only: + lines.append( + f"- Simulink fidelity: native Python coverage exists for the required published workflows, and {len(simulink_reference_only)} inventoried MATLAB assets remain reference-only." + ) else: lines.append("- Simulink fidelity: all inventoried Simulink-backed workflows have an explicit Python execution strategy.") lines.extend(["", "## Remaining Mapping Deltas", ""]) @@ -228,13 +238,17 @@ def render_parity_report(repo_root: Path | None = None) -> str: lines.append(f"- `{label}` -> `{python_target}` [{row['status']}]: {detail}") lines.extend(["", "## Simulink Fidelity Deltas", ""]) - if not simulink_outstanding: - lines.append("No partial, fallback, or unsupported Simulink execution paths remain in the audit.") + if not simulink_outstanding and not simulink_reference_only: + lines.append("No partial, reference-only, fallback, or unsupported Simulink execution paths remain in the audit.") else: for row in simulink_outstanding: lines.append( f"- `{row['model_name']}` -> `{row['model_path']}` [{row['python_strategy']}/{row['current_python_status']}]: {row['chosen_interoperability_strategy']}" ) + for row in simulink_reference_only: + lines.append( + f"- `{row['model_name']}` -> `{row['model_path']}` [{row['python_strategy']}/{row['current_python_status']}]: {row['chosen_interoperability_strategy']}" + ) lines.extend(["", "## Justified Non-Applicable Items", ""]) non_applicable = _iter_non_applicable_rows(payload) diff --git a/parity/class_fidelity.yml b/parity/class_fidelity.yml index a4ff1489..37c663eb 100644 --- a/parity/class_fidelity.yml +++ b/parity/class_fidelity.yml @@ -1,5 +1,5 @@ version: 1 -generated_on: 2026-03-07 +generated_on: 2026-03-08 source_repositories: matlab: https://github.com/cajigaslab/nSTAT python: https://github.com/cajigaslab/nSTAT-python @@ -250,8 +250,12 @@ items: known_remaining_differences: - Advanced MATLAB algorithm-selection, cross-validation, and some report-layout branches are still lighter than MATLAB. + - The canonical Poisson GLM path is now fixture-matched for coefficients, lambda + traces, AIC, and BIC on the baseline workflow, but the stored `logLL` convention + is still not numerically exact to MATLAB. required_remediation: - - Add dataset-backed numerical parity fixtures for canonical analysis workflows. + - Extend the committed MATLAB-derived fixture coverage beyond the canonical single-neuron + GLM workflow to multi-neuron, validation-window, and alternative algorithm branches. - Port remaining algorithm-selection and validation-option branches from MATLAB. plotting_report_parity: KS, inverse-Gaussian, coefficient, residual, and summary plots now execute on canonical Analysis output; advanced algorithm-selection, diff --git a/parity/notebook_fidelity.yml b/parity/notebook_fidelity.yml index e44ee7f7..a6bc3a87 100644 --- a/parity/notebook_fidelity.yml +++ b/parity/notebook_fidelity.yml @@ -1,5 +1,5 @@ version: 1 -generated_on: '2026-03-07' +generated_on: '2026-03-08' source_repositories: matlab: https://github.com/cajigaslab/nSTAT python: https://github.com/cajigaslab/nSTAT-python @@ -31,9 +31,8 @@ items: python_notebook: notebooks/TrialExamples.ipynb fidelity_status: high_fidelity remaining_differences: The notebook now mirrors the MATLAB Trial workflow with executable - object construction, masking, history extraction, and plotting; the closing analysis - section uses one representative Python `Analysis` run instead of linking out to - separate MATLAB help pages. + object construction, masking, history extraction, and plotting; only minor Python + plotting defaults differ from the published MATLAB help output. python_sections: 9 python_expected_figures: 6 python_uses_figure_tracker: true @@ -240,8 +239,8 @@ items: fidelity_status: high_fidelity remaining_differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; - the Python port uses shorter deterministic simulations than MATLAB so the notebook - remains stable in CI. + local execution uses the MATLAB-scale simulation sizes, while CI switches to a + documented shorter deterministic fast path for stability. python_sections: 11 python_expected_figures: 10 python_uses_figure_tracker: true @@ -258,11 +257,11 @@ items: - topic: StimulusDecode2D source_matlab: StimulusDecode2D.mlx python_notebook: notebooks/StimulusDecode2D.ipynb - fidelity_status: high_fidelity - remaining_differences: The notebook now reproduces the 2-D stimulus-decoding workflow - with simulated receptive fields and decoded trajectories; the current Python decoder - uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear - filter. + fidelity_status: partial + remaining_differences: The notebook reproduces the MATLAB section order, figure + inventory, simulated receptive fields, and decoded-trajectory presentation, but + the current Python decoder still uses regression-based state recovery instead + of MATLAB's symbolic-CIF nonlinear filter. python_sections: 4 python_expected_figures: 6 python_uses_figure_tracker: true diff --git a/parity/report.md b/parity/report.md index f5ccd2d4..30ff20cf 100644 --- a/parity/report.md +++ b/parity/report.md @@ -34,8 +34,8 @@ Generated from `parity/manifest.yml`, `parity/class_fidelity.yml`, and `tools/no | Status | Count | |---|---:| | `exact` | 0 | -| `high_fidelity` | 13 | -| `partial` | 0 | +| `high_fidelity` | 12 | +| `partial` | 1 | ## Simulink Fidelity Summary @@ -46,16 +46,17 @@ Generated from `parity/manifest.yml`, `parity/class_fidelity.yml`, and `tools/no | `packaged_runtime` | 0 | | `matlab_engine_fallback` | 0 | | `unsupported` | 0 | -| `reference_only` | 4 | +| `reference_only` | 10 | ## Coverage Notes - Public API: no missing MATLAB public APIs remain; only the MATLAB help-browser utility is explicitly non-applicable. - Help/notebook parity: all inventoried MATLAB help workflows are mapped to Python notebooks or equivalents. -- Notebook fidelity: all tracked MATLAB-helpfile notebook ports are marked high fidelity or exact. +- Notebook fidelity: workflow coverage is complete, but 1 MATLAB-helpfile notebook ports are still marked partial in `tools/notebooks/parity_notes.yml`. +- Notebook fidelity audit: structural section/figure comparisons plus placeholder/tracker-only cell detection are recorded in `parity/notebook_fidelity.yml`. - Paper examples and docs gallery: all canonical paper examples and committed gallery directories are mapped. - Class fidelity: the class audit reports no partial, wrapper-only, or missing items. -- Simulink fidelity: all inventoried Simulink-backed workflows have an explicit Python execution strategy. +- Simulink fidelity: native Python coverage exists for the required published workflows, and 10 inventoried MATLAB assets remain reference-only. ## Remaining Mapping Deltas @@ -63,7 +64,7 @@ No partial or missing items remain in the mapping inventory. ## Remaining Notebook-Fidelity Deltas -No partial notebook-fidelity items remain in `tools/notebooks/parity_notes.yml`. +- `StimulusDecode2D` -> `notebooks/StimulusDecode2D.ipynb` [partial]: The notebook reproduces the MATLAB section order, figure inventory, simulated receptive fields, and decoded-trajectory presentation, but the current Python decoder still uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear filter. ## Remaining Class-Fidelity Deltas @@ -71,7 +72,16 @@ No partial, wrapper-only, or missing class-fidelity items remain. ## Simulink Fidelity Deltas -No partial, fallback, or unsupported Simulink execution paths remain in the audit. +- `PointProcessSimulationCont` -> `PointProcessSimulationCont.slx` [reference_only/reference_only]: Keep as reference while the Python port uses the native discrete simulation path. +- `PointProcessSimulationLegacy2010b` -> `PointProcessSimulation.mdl.r2010b` [reference_only/reference_only]: Treat as a compatibility/reference asset because the native Python port targets the current `PointProcessSimulation.slx` behavior rather than every historic MATLAB model format. +- `PointProcessSimulationLegacy2011a` -> `PointProcessSimulation.mdl.r2011a` [reference_only/reference_only]: Treat as a compatibility/reference asset alongside the current `.slx` model. +- `PointProcessSimulationLegacy2011b` -> `PointProcessSimulation.mdl.r2011b` [reference_only/reference_only]: Treat as a compatibility/reference asset alongside the current `.slx` model. +- `PointProcessSimulationLegacy2013a` -> `PointProcessSimulation.mdl.r2013a` [reference_only/reference_only]: Treat as a compatibility/reference asset alongside the current `.slx` model. +- `PointProcessSimulationLegacySLX2013a` -> `PointProcessSimulation.slx.r2013a` [reference_only/reference_only]: Treat as a versioned MATLAB reference asset rather than a distinct Python execution target. +- `PointProcessSimulationThinningLegacy2011a` -> `PointProcessSimulationThinning.mdl.r2011a` [reference_only/reference_only]: Keep as a MATLAB reference asset while the Python port validates thinning behavior through `CIF.simulateCIFByThinning` and MATLAB Engine comparisons. +- `PointProcessSimulationCache` -> `PointProcessSimulation.slxc` [reference_only/reference_only]: Treat as a MATLAB build artifact, not as a Python execution target. +- `HelpPointProcessSimulationCache` -> `helpfiles/PointProcessSimulation.slxc` [reference_only/reference_only]: Treat as a published-help artifact only. +- `SimulatedNetwork2Cache` -> `helpfiles/SimulatedNetwork2.slxc` [reference_only/reference_only]: Treat as a MATLAB build artifact, not a Python target. ## Justified Non-Applicable Items diff --git a/parity/simulink_fidelity.yml b/parity/simulink_fidelity.yml index 4e793344..40739395 100644 --- a/parity/simulink_fidelity.yml +++ b/parity/simulink_fidelity.yml @@ -1,5 +1,5 @@ version: 1 -generated_on: 2026-03-07 +generated_on: 2026-03-08 source_repositories: matlab: https://github.com/cajigaslab/nSTAT python: https://github.com/cajigaslab/nSTAT-python @@ -35,6 +35,72 @@ items: - No Python-executable equivalent currently mirrors the continuous Simulink model directly. validation_plan: - Audit whether any MATLAB helpfile or paper workflow depends on the continuous model before promoting it beyond reference-only. + - model_name: PointProcessSimulationLegacy2010b + model_path: PointProcessSimulation.mdl.r2010b + purpose: Legacy pre-SLX discrete point-process model retained for MATLAB compatibility and regression reference. + matlab_usage: used_for_reference + python_strategy: reference_only + current_python_status: reference_only + chosen_interoperability_strategy: Treat as a compatibility/reference asset because the native Python port targets the current `PointProcessSimulation.slx` behavior rather than every historic MATLAB model format. + fidelity_risks: + - Legacy `.mdl` files may embed solver or block-version behavior that differs from the current `.slx` model. + validation_plan: + - Keep inventoried and compare only if a MATLAB helpfile or published baseline is shown to depend on the legacy model specifically. + - model_name: PointProcessSimulationLegacy2011a + model_path: PointProcessSimulation.mdl.r2011a + purpose: Legacy pre-SLX discrete point-process model retained for MATLAB compatibility and regression reference. + matlab_usage: used_for_reference + python_strategy: reference_only + current_python_status: reference_only + chosen_interoperability_strategy: Treat as a compatibility/reference asset alongside the current `.slx` model. + fidelity_risks: + - Block behavior may differ subtly from the current `.slx` model even when the published workflow is unchanged. + validation_plan: + - Keep inventoried and promote only if a workflow is shown to require this exact legacy model. + - model_name: PointProcessSimulationLegacy2011b + model_path: PointProcessSimulation.mdl.r2011b + purpose: Legacy pre-SLX discrete point-process model retained for MATLAB compatibility and regression reference. + matlab_usage: used_for_reference + python_strategy: reference_only + current_python_status: reference_only + chosen_interoperability_strategy: Treat as a compatibility/reference asset alongside the current `.slx` model. + fidelity_risks: + - Block behavior may differ subtly from the current `.slx` model even when the published workflow is unchanged. + validation_plan: + - Keep inventoried and promote only if a workflow is shown to require this exact legacy model. + - model_name: PointProcessSimulationLegacy2013a + model_path: PointProcessSimulation.mdl.r2013a + purpose: Legacy pre-SLX discrete point-process model retained for MATLAB compatibility and regression reference. + matlab_usage: used_for_reference + python_strategy: reference_only + current_python_status: reference_only + chosen_interoperability_strategy: Treat as a compatibility/reference asset alongside the current `.slx` model. + fidelity_risks: + - Block behavior may differ subtly from the current `.slx` model even when the published workflow is unchanged. + validation_plan: + - Keep inventoried and promote only if a workflow is shown to require this exact legacy model. + - model_name: PointProcessSimulationLegacySLX2013a + model_path: PointProcessSimulation.slx.r2013a + purpose: Version-specific discrete point-process SLX asset retained for MATLAB compatibility and regression reference. + matlab_usage: used_for_reference + python_strategy: reference_only + current_python_status: reference_only + chosen_interoperability_strategy: Treat as a versioned MATLAB reference asset rather than a distinct Python execution target. + fidelity_risks: + - Version-specific SLX files are not stable Python execution targets and may diverge from the current published model. + validation_plan: + - Keep inventoried and compare only if a published workflow is tied to this asset specifically. + - model_name: PointProcessSimulationThinningLegacy2011a + model_path: PointProcessSimulationThinning.mdl.r2011a + purpose: Legacy thinning-based point-process simulation model retained for MATLAB compatibility/reference. + matlab_usage: used_for_reference + python_strategy: reference_only + current_python_status: reference_only + chosen_interoperability_strategy: Keep as a MATLAB reference asset while the Python port validates thinning behavior through `CIF.simulateCIFByThinning` and MATLAB Engine comparisons. + fidelity_risks: + - The legacy thinning model may encode branch behavior that is not separately exercised by current Python CI. + validation_plan: + - Keep inventoried and add direct MATLAB-reference comparisons if a helpfile or paper workflow depends on thinning-specific semantics. - model_name: PointProcessSimulationCache model_path: PointProcessSimulation.slxc purpose: Simulink compiled cache artifact for `PointProcessSimulation`. diff --git a/tests/parity/fixtures/matlab_gold/analysis_exactness.mat b/tests/parity/fixtures/matlab_gold/analysis_exactness.mat new file mode 100644 index 0000000000000000000000000000000000000000..8e1e89f1f4ee0f98fed609a2c1fa88460c36d40e GIT binary patch literal 924 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X01nwjV*I2WZRmZYXAM7+9GYSQ%O<7#SFuDG&=7V1UunmmkOu1>%Z1kCPJ;UN9t?DLiYq zq;P=iX)&{*vczVCDJjpE6*1iu(^+ zk}a$fc(qN4DXPngiZWzoRu8$?;BVIY^t@3uLyRoz5|B1LA;Y$@Sg3e%kkDgqqvGI$ zuPhoZniv_=&ayPl{U-2h$=*&wj+nMj`L~=yY&WoTOWSYrXKqdp_&NPs`uTL{t3@oPEBuFXYrS z)Hb#@kS-CwT&?~LbGcYYE(=)A3&SNe^$UX)8%{wGC LSQz&F<~R!g)PPiD literal 0 HcmV?d00001 diff --git a/tests/parity/fixtures/matlab_gold/cif_exactness.mat b/tests/parity/fixtures/matlab_gold/cif_exactness.mat index a9b37b028b8ed3a7bbc412c0a62c9c2b4a003baf..ddd9e0bff43c9ca423e1ba1e0c2dfee0ca19a63e 100644 GIT binary patch delta 40 wcmaFH@{DDIiFj$Af^TAxf`WyDfq|8Yft8_=f{}rd*~CEQi3w~QOFl6I0QV3K1ONa4 delta 40 wcmaFH@{DDIiFjg(f^TAxf`Yk%k)f4=k(Hs5f{}rd*~CEQi3w~QOFl6I0QEf!_5c6? diff --git a/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat b/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat index bdcfad9b3ffa8a060e78055d108198063f3e5780..b43b32f5477ce7e6c354160510d55f8bf9c9c81e 100644 GIT binary patch delta 40 vcmZ3@y_$Q1iFj$Af^TAxf`WyDfq|8Yft8`5f{}rd*~CEQi3w~QOUzjT<){ke delta 40 vcmZ3@y_$Q1iFjg(f^TAxf`Yk%k)f55sg;45f{}rd*~CEQi3w~QOUzjTdXu3 diff --git a/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat b/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat index b15a0a7406133ab5bf57743ede528f6897a29466..f3e667f3f35c79937df7006827445f4ea76bfd96 100644 GIT binary patch delta 40 wcmbQoHIHk8iFj$Af^TAxf`WyDfq|8Yft8`5f{}rd*~CEQi3w~QOIEP}0OCsvIsgCw delta 40 wcmbQoHIHk8iFjg(f^TAxf`Yk%k)f55sg;4Lf{}rd*~CEQi3w~QOIEP}0O1G=Gynhq diff --git a/tests/test_matlab_gold_fixtures.py b/tests/test_matlab_gold_fixtures.py index 48c2ad98..8c87b029 100644 --- a/tests/test_matlab_gold_fixtures.py +++ b/tests/test_matlab_gold_fixtures.py @@ -5,7 +5,7 @@ import numpy as np from scipy.io import loadmat -from nstat import CIF, Covariate, SignalObj, nspikeTrain +from nstat import Analysis, CIF, ConfigColl, CovColl, Covariate, SignalObj, Trial, TrialConfig, nspikeTrain, nstColl REPO_ROOT = Path(__file__).resolve().parents[1] @@ -24,6 +24,18 @@ def _vector(payload: dict[str, np.ndarray], key: str) -> np.ndarray: return np.asarray(payload[key], dtype=float).reshape(-1) +def _string(payload: dict[str, np.ndarray], key: str) -> str: + value = payload[key] + if isinstance(value, bytes): + return value.decode("utf-8") + if isinstance(value, str): + return value + arr = np.asarray(value) + if arr.shape == (): + return str(arr.item()) + return str(arr.reshape(-1)[0]) + + def test_signalobj_matches_matlab_gold_fixture() -> None: payload = _load_fixture("signalobj_exactness.mat") signal = SignalObj(_vector(payload, "time"), np.asarray(payload["data"], dtype=float), "sig", "time", "s", "u", ["x1", "x2"]) @@ -111,6 +123,28 @@ def test_cif_eval_surface_matches_matlab_gold_fixture() -> None: np.testing.assert_allclose(cif.evalJacobianLog(stim_val), np.asarray(payload["jacobian_log"], dtype=float), rtol=1e-8, atol=1e-10) +def test_analysis_fit_surface_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("analysis_exactness.mat") + time = _vector(payload, "time") + stim_data = _vector(payload, "stim_data") + spike_times = _vector(payload, "spike_times") + sample_rate = _scalar(payload, "sample_rate") + + stim = Covariate(time, stim_data, "Stimulus", "time", "s", "", ["stim"]) + spike_train = nspikeTrain(spike_times, "1", 0.1, 0.0, 1.0, "time", "s", "", "", -1) + trial = Trial(nstColl([spike_train]), CovColl([stim])) + cfg = TrialConfig([["Stimulus", "stim"]], sample_rate, [], [], name="stim") + fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl([cfg])) + + np.testing.assert_allclose(fit.getCoeffs(1), _vector(payload, "coeffs"), rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(fit.lambdaSignal.time, _vector(payload, "lambda_time"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(fit.lambdaSignal.data[:, 0], _vector(payload, "lambda_data"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(fit.AIC[0]), _scalar(payload, "AIC"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(fit.BIC[0]), _scalar(payload, "BIC"), rtol=1e-8, atol=1e-10) + assert np.isfinite(float(fit.logLL[0])) + assert fit.fitType[0] == _string(payload, "distribution") + + def test_point_process_lambda_trace_matches_matlab_gold_fixture() -> None: payload = _load_fixture("point_process_exactness.mat") diff --git a/tests/test_matlab_reference.py b/tests/test_matlab_reference.py index 83d8cc4f..086ad06e 100644 --- a/tests/test_matlab_reference.py +++ b/tests/test_matlab_reference.py @@ -5,9 +5,10 @@ import numpy as np import pytest -from nstat import CIF, Covariate, simulate_two_neuron_network +from nstat import Analysis, CIF, ConfigColl, CovColl, Covariate, Trial, TrialConfig, nspikeTrain, nstColl, simulate_two_neuron_network from nstat.matlab_reference import ( matlab_engine_available, + run_analysis_reference, run_point_process_reference, run_simulated_network_reference, ) @@ -73,3 +74,23 @@ def test_native_network_simulation_preserves_matlab_connectivity_layout_when_eng np.testing.assert_allclose(native.actual_network, matlab_ref["actual_network"]) np.testing.assert_allclose(native.lambda_delta[:5], matlab_ref["prob_head"], rtol=1e-6, atol=1e-8) assert np.all((matlab_ref["state_head"] == 0.0) | (matlab_ref["state_head"] == 1.0)) + + +@pytest.mark.skipif(not MATLAB_REPO_ROOT.exists(), reason="MATLAB reference repo not available") +def test_native_analysis_fit_matches_matlab_reference_when_engine_is_available() -> None: + if not matlab_engine_available(): + pytest.skip("MATLAB Engine for Python is not installed") + + time = np.arange(0.0, 1.0 + 0.1, 0.1) + stim = Covariate(time, np.sin(2 * np.pi * time), "Stimulus", "time", "s", "", ["stim"]) + spike_train = nspikeTrain([0.1, 0.4, 0.7], "1", 0.1, 0.0, 1.0, "time", "s", "", "", -1) + trial = Trial(nstColl([spike_train]), CovColl([stim])) + cfg = TrialConfig([["Stimulus", "stim"]], 10.0, [], [], name="stim") + fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl([cfg])) + matlab_ref = run_analysis_reference(matlab_repo=MATLAB_REPO_ROOT) + + np.testing.assert_allclose(fit.getCoeffs(1), matlab_ref["coeffs"], rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(fit.lambdaSignal.data[:5, 0], matlab_ref["lambda_head"], rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(np.asarray(fit.AIC, dtype=float)[:1], matlab_ref["aic"], rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(np.asarray(fit.BIC, dtype=float)[:1], matlab_ref["bic"], rtol=1e-6, atol=1e-8) + assert np.isfinite(np.asarray(fit.logLL, dtype=float)[0]) diff --git a/tests/test_notebook_fidelity_audit.py b/tests/test_notebook_fidelity_audit.py index f36b68a3..40aac0d8 100644 --- a/tests/test_notebook_fidelity_audit.py +++ b/tests/test_notebook_fidelity_audit.py @@ -45,10 +45,10 @@ def test_notebook_fidelity_audit_marks_upgraded_ports_as_high_fidelity() -> None } <= high_fidelity_topics -def test_notebook_fidelity_audit_has_no_partial_or_placeholder_notebooks() -> None: +def test_notebook_fidelity_audit_tracks_only_known_partial_notebooks() -> None: audit = yaml.safe_load(AUDIT_PATH.read_text(encoding="utf-8")) or {} partial_topics = {row["topic"] for row in audit.get("items", []) if row["fidelity_status"] in {"partial", "placeholder", "missing"}} - assert not partial_topics + assert partial_topics == {"StimulusDecode2D"} def test_high_fidelity_notebooks_have_no_placeholder_or_tracker_only_cells() -> None: diff --git a/tests/test_notebook_parity_notes.py b/tests/test_notebook_parity_notes.py index a105ad01..f49ca12f 100644 --- a/tests/test_notebook_parity_notes.py +++ b/tests/test_notebook_parity_notes.py @@ -38,9 +38,9 @@ def test_target_notebooks_start_with_machine_readable_parity_note() -> None: assert row["remaining_differences"] in source -def test_notebook_parity_notes_have_no_partial_statuses() -> None: +def test_notebook_parity_notes_track_only_known_partial_statuses() -> None: partial = [row["topic"] for row in _load_notes() if row["fidelity_status"] == "partial"] - assert not partial + assert partial == ["StimulusDecode2D"] def test_high_fidelity_parity_notes_do_not_admit_placeholder_or_tracker_only_status() -> None: diff --git a/tests/test_parity_report.py b/tests/test_parity_report.py index d1938ce4..cbbb1bf8 100644 --- a/tests/test_parity_report.py +++ b/tests/test_parity_report.py @@ -23,13 +23,13 @@ def test_parity_report_highlights_current_constraints() -> None: assert "Notebook Fidelity Summary" in text assert "Simulink Fidelity Summary" in text assert "Remaining Notebook-Fidelity Deltas" in text - assert "all tracked MATLAB-helpfile notebook ports are marked high fidelity or exact" in text - assert "No partial notebook-fidelity items remain in `tools/notebooks/parity_notes.yml`." in text + assert "workflow coverage is complete, but 1 MATLAB-helpfile notebook ports are still marked partial" in text + assert "`StimulusDecode2D` -> `notebooks/StimulusDecode2D.ipynb` [partial]" in text assert "No partial or missing items remain in the mapping inventory." in text assert "Remaining Class-Fidelity Deltas" in text assert "the class audit reports no partial, wrapper-only, or missing items" in text assert "No partial, wrapper-only, or missing class-fidelity items remain." in text assert "Simulink Fidelity Deltas" in text - assert "all inventoried Simulink-backed workflows have an explicit Python execution strategy" in text - assert "No partial, fallback, or unsupported Simulink execution paths remain in the audit." in text + assert "native Python coverage exists for the required published workflows" in text + assert "reference-only" in text assert "nstatOpenHelpPage" in text diff --git a/tests/test_simulink_fidelity_audit.py b/tests/test_simulink_fidelity_audit.py index 9843d88b..a26d47f6 100644 --- a/tests/test_simulink_fidelity_audit.py +++ b/tests/test_simulink_fidelity_audit.py @@ -17,6 +17,17 @@ "unsupported", "reference_only", } +REQUIRED_MODEL_PATHS = { + "PointProcessSimulation.slx", + "PointProcessSimulationCont.slx", + "PointProcessSimulation.mdl.r2010b", + "PointProcessSimulation.mdl.r2011a", + "PointProcessSimulation.mdl.r2011b", + "PointProcessSimulation.mdl.r2013a", + "PointProcessSimulation.slx.r2013a", + "PointProcessSimulationThinning.mdl.r2011a", + "helpfiles/SimulatedNetwork2.mdl", +} def _load_audit() -> dict: @@ -40,6 +51,12 @@ def test_simulink_fidelity_audit_records_required_execution_fields() -> None: assert row["validation_plan"] +def test_simulink_fidelity_audit_covers_required_model_inventory() -> None: + payload = _load_audit() + model_paths = {row["model_path"] for row in payload["items"]} + assert REQUIRED_MODEL_PATHS <= model_paths + + def test_simulink_fidelity_audit_paths_exist_when_matlab_repo_is_available() -> None: if not MATLAB_REPO_ROOT.exists(): pytest.skip(f"MATLAB reference repo not available at {MATLAB_REPO_ROOT}") diff --git a/tools/notebooks/build_foundational_help_notebooks.py b/tools/notebooks/build_foundational_help_notebooks.py index b84ccfb4..77f3c316 100644 --- a/tools/notebooks/build_foundational_help_notebooks.py +++ b/tools/notebooks/build_foundational_help_notebooks.py @@ -50,7 +50,7 @@ def _write_notebook( ## MATLAB Parity Note - Source MATLAB helpfile: `TrialExamples.mlx` - Fidelity status: `high_fidelity` -- Remaining justified differences: The notebook now mirrors the MATLAB Trial workflow with executable object construction, masking, history extraction, and plotting; the closing analysis section uses one representative Python `Analysis` run instead of linking out to separate MATLAB help pages. +- Remaining justified differences: The notebook now mirrors the MATLAB Trial workflow with executable object construction, masking, history extraction, and plotting; only minor Python plotting defaults differ from the published MATLAB help output. """ @@ -72,7 +72,7 @@ def _write_notebook( import matplotlib.pyplot as plt import numpy as np - from nstat import Analysis, ConfigColl, CovColl, Covariate, Events, History, Trial, TrialConfig, nspikeTrain, nstColl + from nstat import CovColl, Covariate, Events, History, Trial, nspikeTrain, nstColl from nstat.notebook_figures import FigureTracker np.random.seed(7) @@ -205,22 +205,11 @@ def _build_trial(): """, """ # SECTION 8: Example 2: Analyzing Trial Data - cfg = TrialConfig([["Position", "x"], ["Force", "f_x"]], sampleRate=ctx["sample_rate"], history=[0.0, 0.05, 0.1], name="Position+Force+History") - cfgColl = ConfigColl([cfg]) - fit = Analysis.RunAnalysisForNeuron(trial1, 1, cfgColl) - fit_stats = fit.computeKSStats() - print( - { - "config_name": fit.configNames[0], - "aic": round(float(fit.AIC[0]), 3), - "bic": round(float(fit.BIC[0]), 3), - "ks_stat": round(float(fit_stats["ks_stat"]), 4), - } - ) + print("Examples of neural spike analysis using AnalysisExamples2 (Neural Spike Analysis Toolbox) or AnalysisExamples (standard methods).") """, """ # SECTION 9: Related analysis workflows - print("For larger model-comparison walkthroughs, continue with AnalysisExamples and AnalysisExamples2.") + print({"recommended_next_notebooks": ["AnalysisExamples2", "AnalysisExamples"]}) __tracker.finalize() """, ] diff --git a/tools/notebooks/build_helpfile_fidelity_notebooks.py b/tools/notebooks/build_helpfile_fidelity_notebooks.py index 0bc53d8a..a4be05ce 100644 --- a/tools/notebooks/build_helpfile_fidelity_notebooks.py +++ b/tools/notebooks/build_helpfile_fidelity_notebooks.py @@ -272,7 +272,7 @@ def _plot_ks(ax, ideal, empirical, ci, *, label, color): ## MATLAB Parity Note - Source MATLAB helpfile: `ValidationDataSet.mlx` - Fidelity status: `high_fidelity` -- Remaining justified differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; the Python port uses shorter deterministic simulations than MATLAB so the notebook remains stable in CI. +- Remaining justified differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; local execution uses the MATLAB-scale simulation sizes, while CI switches to a documented shorter deterministic fast path for stability. """ @@ -289,6 +289,8 @@ def _plot_ks(ax, ideal, empirical, ci, *, label, color): if str(SRC_PATH) not in sys.path: sys.path.insert(0, str(SRC_PATH)) + import os + import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt @@ -317,7 +319,12 @@ def _lambda_columns(fit_result): return time, data - def _simulate_constant_case(seed=0, *, p=0.01, n_samples=20001, delta=0.001): + CI_FAST_PATH = os.environ.get("CI", "").strip().lower() in {"1", "true", "yes"} + + + def _simulate_constant_case(seed=0, *, p=0.01, n_samples=None, delta=0.001): + if n_samples is None: + n_samples = 20001 if CI_FAST_PATH else 100001 rng = np.random.default_rng(seed) total_time = n_samples * delta time = np.linspace(0.0, total_time, n_samples) @@ -344,7 +351,11 @@ def _simulate_constant_case(seed=0, *, p=0.01, n_samples=20001, delta=0.001): } - def _simulate_piecewise_case(seed=1, *, p1=0.001, p2=0.01, n1=20000, n2=20000, delta=0.001): + def _simulate_piecewise_case(seed=1, *, p1=0.001, p2=0.01, n1=None, n2=None, delta=0.001): + if n1 is None: + n1 = 20000 if CI_FAST_PATH else 100000 + if n2 is None: + n2 = 20000 if CI_FAST_PATH else 100000 rng = np.random.default_rng(seed) t1 = np.linspace(0.0, n1 * delta, n1 + 1) t2 = np.linspace(n1 * delta, (n1 + n2) * delta, n2 + 1)[1:] @@ -401,12 +412,13 @@ def _plot_isi_hist(ax, train, lambda_hz, *, title): """, """ # SECTION 0: Software Validation Data Set - # This notebook follows the MATLAB validation helpfile with deterministic simulations for CI-stable execution. + # This notebook follows the MATLAB validation helpfile; CI uses a documented short fast path while local runs use MATLAB-scale sample counts. plt.close("all") constant_case = _simulate_constant_case() piecewise_case = _simulate_piecewise_case() print( { + "ci_fast_path": CI_FAST_PATH, "constant_lambda_hz": round(float(constant_case["lambda_hz"]), 4), "piecewise_lambda1_hz": round(float(piecewise_case["lambda1_hz"]), 4), "piecewise_lambda2_hz": round(float(piecewise_case["lambda2_hz"]), 4), @@ -956,8 +968,8 @@ def _plot_raster(ax, time_s, spikes, *, max_cells=18): ## MATLAB Parity Note - Source MATLAB helpfile: `StimulusDecode2D.mlx` -- Fidelity status: `high_fidelity` -- Remaining justified differences: The notebook now reproduces the 2-D stimulus-decoding workflow with simulated receptive fields and decoded trajectories; the current Python decoder uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear filter. +- Fidelity status: `partial` +- Remaining justified differences: The notebook reproduces the MATLAB section order, figure inventory, simulated receptive fields, and decoded-trajectory presentation, but the current Python decoder still uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear filter. """ diff --git a/tools/notebooks/parity_notes.yml b/tools/notebooks/parity_notes.yml index 06a43428..fa1b4c8d 100644 --- a/tools/notebooks/parity_notes.yml +++ b/tools/notebooks/parity_notes.yml @@ -9,7 +9,7 @@ notes: file: notebooks/TrialExamples.ipynb source_matlab: TrialExamples.mlx fidelity_status: high_fidelity - remaining_differences: The notebook now mirrors the MATLAB Trial workflow with executable object construction, masking, history extraction, and plotting; the closing analysis section uses one representative Python `Analysis` run instead of linking out to separate MATLAB help pages. + remaining_differences: The notebook now mirrors the MATLAB Trial workflow with executable object construction, masking, history extraction, and plotting; only minor Python plotting defaults differ from the published MATLAB help output. - topic: AnalysisExamples file: notebooks/AnalysisExamples.ipynb source_matlab: AnalysisExamples.mlx @@ -59,9 +59,9 @@ notes: file: notebooks/ValidationDataSet.ipynb source_matlab: ValidationDataSet.mlx fidelity_status: high_fidelity - remaining_differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; the Python port uses shorter deterministic simulations than MATLAB so the notebook remains stable in CI. + remaining_differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; local execution uses the MATLAB-scale simulation sizes, while CI switches to a documented shorter deterministic fast path for stability. - topic: StimulusDecode2D file: notebooks/StimulusDecode2D.ipynb source_matlab: StimulusDecode2D.mlx - fidelity_status: high_fidelity - remaining_differences: The notebook now reproduces the 2-D stimulus-decoding workflow with simulated receptive fields and decoded trajectories; the current Python decoder uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear filter. + fidelity_status: partial + remaining_differences: The notebook reproduces the MATLAB section order, figure inventory, simulated receptive fields, and decoded-trajectory presentation, but the current Python decoder still uses regression-based state recovery instead of MATLAB's symbolic-CIF nonlinear filter. diff --git a/tools/parity/matlab/export_matlab_gold_fixtures.m b/tools/parity/matlab/export_matlab_gold_fixtures.m index 2c4fe953..0b4d1fa6 100644 --- a/tools/parity/matlab/export_matlab_gold_fixtures.m +++ b/tools/parity/matlab/export_matlab_gold_fixtures.m @@ -21,6 +21,7 @@ function export_matlab_gold_fixtures(repoRoot, matlabRepoRoot) export_signalobj_fixture(fixtureRoot); export_nspiketrain_fixture(fixtureRoot); export_cif_fixture(fixtureRoot); +export_analysis_fixture(fixtureRoot); export_point_process_fixture(fixtureRoot); end @@ -114,6 +115,32 @@ function export_cif_fixture(fixtureRoot) save(fullfile(fixtureRoot, 'cif_exactness.mat'), '-struct', 'payload'); end +function export_analysis_fixture(fixtureRoot) +t = (0:0.1:1.0)'; +stimData = sin(2*pi*t); +stim = Covariate(t, stimData, 'Stimulus', 'time', 's', '', {'stim'}); +spikeTrain = nspikeTrain([0.1 0.4 0.7], '1', 0.1, 0.0, 1.0, 'time', 's', '', '', -1); +trial = Trial(nstColl({spikeTrain}), CovColl({stim})); +cfg = TrialConfig({{'Stimulus', 'stim'}}, 10, [], []); +cfg.setName('stim'); +fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl({cfg})); + +payload = struct(); +payload.time = t; +payload.stim_data = stimData; +payload.spike_times = spikeTrain.spikeTimes; +payload.sample_rate = trial.sampleRate; +payload.coeffs = fit.getCoeffs(1); +payload.lambda_time = fit.lambda.time; +payload.lambda_data = fit.lambda.data(:,1); +payload.AIC = fit.AIC(1); +payload.BIC = fit.BIC(1); +payload.logLL = fit.logLL(1); +payload.distribution = fit.fitType{1}; + +save(fullfile(fixtureRoot, 'analysis_exactness.mat'), '-struct', 'payload'); +end + function export_point_process_fixture(fixtureRoot) rng(5); Ts = .001;