diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9b7a81e4..fd9c33f3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -196,7 +196,7 @@ jobs: python -m pip install --upgrade pip python -m pip install -e .[dev] - name: Run API surface checks - run: pytest -q tests/test_api_surface.py + run: pytest -q tests/test_api_surface.py tests/test_cleanroom_boundary.py symbol-surface-audit: runs-on: ubuntu-latest diff --git a/notebooks/ExplicitStimulusWhiskerData.ipynb b/notebooks/ExplicitStimulusWhiskerData.ipynb index 19c4d8fd..436fab61 100644 --- a/notebooks/ExplicitStimulusWhiskerData.ipynb +++ b/notebooks/ExplicitStimulusWhiskerData.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "a51db33a", + "id": "1039ac25", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9f51d683", + "id": "ea575816", "metadata": {}, "outputs": [], "source": [ @@ -83,7 +83,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4ee3a32e", + "id": "c32ee0d2", "metadata": {}, "outputs": [], "source": [ @@ -106,7 +106,7 @@ { "cell_type": "code", "execution_count": null, - "id": "35afed3f", + "id": "3e7f4cb5", "metadata": {}, "outputs": [], "source": [ @@ -134,7 +134,7 @@ { "cell_type": "code", "execution_count": null, - "id": "214e2ec3", + "id": "79a0c40f", "metadata": {}, "outputs": [], "source": [ @@ -149,7 +149,7 @@ { "cell_type": "code", "execution_count": null, - "id": "782d0257", + "id": "88ef749f", "metadata": {}, "outputs": [], "source": [ @@ -170,7 +170,7 @@ { "cell_type": "code", "execution_count": null, - "id": "494c31d8", + "id": "8fb09e81", "metadata": {}, "outputs": [], "source": [ @@ -198,7 +198,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b49274b0", + "id": "1a7627c6", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +230,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c1129e2d", + "id": "3097f64f", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/HippocampalPlaceCellExample.ipynb b/notebooks/HippocampalPlaceCellExample.ipynb index e5682b19..5b7a7816 100644 --- a/notebooks/HippocampalPlaceCellExample.ipynb +++ b/notebooks/HippocampalPlaceCellExample.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "04918a27", + "id": "9cc6c119", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "cd6ae412", + "id": "3743d52f", "metadata": {}, "outputs": [], "source": [ @@ -85,7 +85,7 @@ { "cell_type": "code", "execution_count": null, - "id": "bc94da8b", + "id": "8804c316", "metadata": {}, "outputs": [], "source": [ @@ -105,7 +105,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7383932a", + "id": "f088de8e", "metadata": {}, "outputs": [], "source": [ @@ -125,7 +125,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3d6944b1", + "id": "f42a4e6c", "metadata": {}, "outputs": [], "source": [ @@ -152,7 +152,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f52ed4fd", + "id": "c75cfbaa", "metadata": {}, "outputs": [], "source": [ @@ -179,7 +179,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ca3ffc88", + "id": "366f0e8c", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/HybridFilterExample.ipynb b/notebooks/HybridFilterExample.ipynb index d4e8eb2d..72346edf 100644 --- a/notebooks/HybridFilterExample.ipynb +++ b/notebooks/HybridFilterExample.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "799ac706", + "id": "51ab8762", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ac9d53a9", + "id": "b9c53e65", "metadata": {}, "outputs": [], "source": [ @@ -63,7 +63,7 @@ { "cell_type": "code", "execution_count": null, - "id": "bf21531d", + "id": "2aeada10", "metadata": {}, "outputs": [], "source": [ @@ -87,7 +87,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e005f11e", + "id": "1ca4154b", "metadata": {}, "outputs": [], "source": [ @@ -98,7 +98,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f52cdb38", + "id": "bc311d6b", "metadata": {}, "outputs": [], "source": [ @@ -109,7 +109,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3c6317b7", + "id": "9a25e97a", "metadata": {}, "outputs": [], "source": [ @@ -161,7 +161,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c0a4b72b", + "id": "78a4e6b5", "metadata": {}, "outputs": [], "source": [ @@ -172,7 +172,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e28a4b2e", + "id": "90fd8d80", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/PPSimExample.ipynb b/notebooks/PPSimExample.ipynb index 1199c4b9..820a2183 100644 --- a/notebooks/PPSimExample.ipynb +++ b/notebooks/PPSimExample.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "75db499a", + "id": "d32212af", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "344b225c", + "id": "10c55dde", "metadata": {}, "outputs": [], "source": [ @@ -71,7 +71,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e586ec40", + "id": "82e2e6d2", "metadata": {}, "outputs": [], "source": [ @@ -82,7 +82,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e7759808", + "id": "bb8a8dca", "metadata": {}, "outputs": [], "source": [ @@ -93,7 +93,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1066f0cb", + "id": "7aa83848", "metadata": {}, "outputs": [], "source": [ @@ -105,7 +105,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c10cdb5a", + "id": "e8660e03", "metadata": {}, "outputs": [], "source": [ @@ -116,7 +116,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1bc35d15", + "id": "7b3aa59d", "metadata": {}, "outputs": [], "source": [ @@ -127,7 +127,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c6fa075a", + "id": "a2f76513", "metadata": {}, "outputs": [], "source": [ @@ -143,7 +143,7 @@ { "cell_type": "code", "execution_count": null, - "id": "2c6f4f13", + "id": "5dd18704", "metadata": {}, "outputs": [], "source": [ @@ -157,7 +157,7 @@ { "cell_type": "code", "execution_count": null, - "id": "cb839486", + "id": "d1867d92", "metadata": {}, "outputs": [], "source": [ @@ -173,7 +173,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f728a04e", + "id": "393f2a17", "metadata": {}, "outputs": [], "source": [ @@ -185,7 +185,7 @@ { "cell_type": "code", "execution_count": null, - "id": "c20b4fcd", + "id": "59d8878a", "metadata": {}, "outputs": [], "source": [ @@ -196,7 +196,7 @@ { "cell_type": "code", "execution_count": null, - "id": "2058e202", + "id": "e13231cb", "metadata": {}, "outputs": [], "source": [ @@ -208,7 +208,7 @@ { "cell_type": "code", "execution_count": null, - "id": "76209533", + "id": "9298a78b", "metadata": {}, "outputs": [], "source": [ @@ -220,7 +220,7 @@ { "cell_type": "code", "execution_count": null, - "id": "11b2922d", + "id": "1bebdc33", "metadata": {}, "outputs": [], "source": [ @@ -232,7 +232,7 @@ { "cell_type": "code", "execution_count": null, - "id": "2666d94d", + "id": "37097000", "metadata": {}, "outputs": [], "source": [ @@ -244,7 +244,7 @@ { "cell_type": "code", "execution_count": null, - "id": "246bce8c", + "id": "a6cdb847", "metadata": {}, "outputs": [], "source": [ @@ -258,7 +258,7 @@ { "cell_type": "code", "execution_count": null, - "id": "481ca8c7", + "id": "78d2200b", "metadata": {}, "outputs": [], "source": [ @@ -266,20 +266,20 @@ "summary = FitResSummary(results)\n", "fig = _figure(\"Summary.plotSummary\", figsize=(10.0, 4.5))\n", "summary.plotSummary(handle=fig)\n", - "print({\"fit_names\": summary.fitNames, \"mean_AIC\": np.asarray(summary.AIC, dtype=float).round(3).tolist()})\n" + "print({\"fit_names\": summary.fitNames, \"mean_AIC\": np.asarray(summary.meanAIC, dtype=float).round(3).tolist()})\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "dc217235", + "id": "b920b8cb", "metadata": {}, "outputs": [], "source": [ "# SECTION 17: Summarize model selection\n", "fig = _figure(\"bar(summary.AIC)\", figsize=(8.0, 4.5))\n", "ax = fig.subplots(1, 1)\n", - "ax.bar(np.arange(len(summary.fitNames)), np.asarray(summary.AIC, dtype=float), color=[\"0.6\", \"tab:blue\", \"tab:green\"])\n", + "ax.bar(np.arange(len(summary.fitNames)), np.asarray(summary.meanAIC, dtype=float), color=[\"0.6\", \"tab:blue\", \"tab:green\"])\n", "ax.set_xticks(np.arange(len(summary.fitNames)), summary.fitNames, rotation=20)\n", "ax.set_ylabel(\"mean AIC\")\n", "ax.set_title(\"Model comparison across realizations\")\n", diff --git a/notebooks/StimulusDecode2D.ipynb b/notebooks/StimulusDecode2D.ipynb index a5b6b0b7..e90ff61a 100644 --- a/notebooks/StimulusDecode2D.ipynb +++ b/notebooks/StimulusDecode2D.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "0e2c1f81", + "id": "4f37c1b6", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `StimulusDecode2D.mlx`\n", "- Fidelity status: `high_fidelity`\n", - "- Remaining justified differences: The notebook now follows the MATLAB nonlinear-CIF decoding workflow and uses `DecodingAlgorithms.PPDecodeFilter` before the same documented linear fallback branch as MATLAB. Exact decoded traces and figure styling can still vary modestly because Python's symbolic/numeric stack and random streams are not byte-identical to MATLAB." + "- Remaining justified differences: The notebook now follows the MATLAB nonlinear-CIF decoding workflow and uses `DecodingAlgorithms.PPDecodeFilter` before the same documented linear fallback branch as MATLAB. Exact decoded traces and figure styling can still vary modestly because Python's symbolic/numeric stack and random streams are not byte-identical to MATLAB.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "ef0ce881", + "id": "f94aa3da", "metadata": {}, "outputs": [], "source": [ @@ -163,7 +163,7 @@ { "cell_type": "code", "execution_count": null, - "id": "318e16a5", + "id": "4c32de3a", "metadata": {}, "outputs": [], "source": [ @@ -184,7 +184,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ea4c43c8", + "id": "25960ea2", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +230,7 @@ { "cell_type": "code", "execution_count": null, - "id": "19e02f21", + "id": "b0784e40", "metadata": {}, "outputs": [], "source": [ @@ -248,7 +248,7 @@ { "cell_type": "code", "execution_count": null, - "id": "f1d2fde2", + "id": "714372bb", "metadata": {}, "outputs": [], "source": [ @@ -303,4 +303,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/TrialExamples.ipynb b/notebooks/TrialExamples.ipynb index 21ee977e..ffe918b2 100644 --- a/notebooks/TrialExamples.ipynb +++ b/notebooks/TrialExamples.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "728c6277", + "id": "914238fa", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "62d9faec", + "id": "8ffbef9e", "metadata": {}, "outputs": [], "source": [ @@ -124,7 +124,7 @@ { "cell_type": "code", "execution_count": null, - "id": "509b204a", + "id": "3ada46f2", "metadata": {}, "outputs": [], "source": [ @@ -140,7 +140,7 @@ { "cell_type": "code", "execution_count": null, - "id": "374bcec9", + "id": "c5ad1d59", "metadata": {}, "outputs": [], "source": [ @@ -153,7 +153,7 @@ { "cell_type": "code", "execution_count": null, - "id": "afc89310", + "id": "87bae7f2", "metadata": {}, "outputs": [], "source": [ @@ -165,7 +165,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0d3c2fa7", + "id": "2bb06d98", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +178,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7b522976", + "id": "9ecba539", "metadata": {}, "outputs": [], "source": [ @@ -191,7 +191,7 @@ { "cell_type": "code", "execution_count": null, - "id": "207e8239", + "id": "828b38f9", "metadata": {}, "outputs": [], "source": [ @@ -203,7 +203,7 @@ { "cell_type": "code", "execution_count": null, - "id": "af1d8f85", + "id": "24ada846", "metadata": {}, "outputs": [], "source": [ @@ -219,7 +219,7 @@ { "cell_type": "code", "execution_count": null, - "id": "efa7bea7", + "id": "406f0e5e", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +230,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4d808c62", + "id": "2a99d085", "metadata": {}, "outputs": [], "source": [ diff --git a/notebooks/ValidationDataSet.ipynb b/notebooks/ValidationDataSet.ipynb index e6292b4e..120e63a0 100644 --- a/notebooks/ValidationDataSet.ipynb +++ b/notebooks/ValidationDataSet.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "5d534eeb", + "id": "d15a297a", "metadata": {}, "source": [ "\n", @@ -15,7 +15,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0177a98f", + "id": "f2698163", "metadata": {}, "outputs": [], "source": [ @@ -154,7 +154,7 @@ { "cell_type": "code", "execution_count": null, - "id": "13e27c11", + "id": "702ec99c", "metadata": {}, "outputs": [], "source": [ @@ -176,7 +176,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d94296d5", + "id": "a65040f3", "metadata": {}, "outputs": [], "source": [ @@ -187,7 +187,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b3d520d3", + "id": "1f04cbfd", "metadata": {}, "outputs": [], "source": [ @@ -199,7 +199,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8b792ab5", + "id": "1a464c1b", "metadata": {}, "outputs": [], "source": [ @@ -213,7 +213,7 @@ { "cell_type": "code", "execution_count": null, - "id": "681433fa", + "id": "fb4f5ef8", "metadata": {}, "outputs": [], "source": [ @@ -235,7 +235,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1ec098fd", + "id": "c286b8c0", "metadata": {}, "outputs": [], "source": [ @@ -257,7 +257,7 @@ { "cell_type": "code", "execution_count": null, - "id": "daf22292", + "id": "879b1951", "metadata": {}, "outputs": [], "source": [ @@ -270,7 +270,7 @@ { "cell_type": "code", "execution_count": null, - "id": "1286f142", + "id": "58f09d75", "metadata": {}, "outputs": [], "source": [ @@ -300,7 +300,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b35f530e", + "id": "6928e1f4", "metadata": {}, "outputs": [], "source": [ @@ -311,7 +311,7 @@ { "cell_type": "code", "execution_count": null, - "id": "75542bfe", + "id": "3522b3b9", "metadata": {}, "outputs": [], "source": [ @@ -346,7 +346,7 @@ { "cell_type": "code", "execution_count": null, - "id": "ad6eec57", + "id": "f576be75", "metadata": {}, "outputs": [], "source": [ @@ -355,10 +355,10 @@ "fig = _prepare_figure(\"Summary.plotSummary\", figsize=(8.5, 4.5))\n", "axs = fig.subplots(1, 2)\n", "xloc = np.arange(len(summary.fitNames))\n", - "axs[0].bar(xloc, summary.AIC, color=[\"tab:blue\", \"tab:green\"])\n", + "axs[0].bar(xloc, summary.meanAIC, color=[\"tab:blue\", \"tab:green\"])\n", "axs[0].set_xticks(xloc, summary.fitNames)\n", "axs[0].set_title(\"Mean AIC across neurons\")\n", - "axs[1].bar(xloc, summary.BIC, color=[\"tab:blue\", \"tab:green\"])\n", + "axs[1].bar(xloc, summary.meanBIC, color=[\"tab:blue\", \"tab:green\"])\n", "axs[1].set_xticks(xloc, summary.fitNames)\n", "axs[1].set_title(\"Mean BIC across neurons\")\n", "\n", diff --git a/nstat/__init__.py b/nstat/__init__.py index 35b3e7b5..9695c97b 100644 --- a/nstat/__init__.py +++ b/nstat/__init__.py @@ -16,7 +16,6 @@ from .fit import FitResSummary, FitResult, FitSummary from .glm import PoissonGLMResult, fit_poisson_glm from .history import History, HistoryBasis -from .matlab_reference import matlab_engine_available, run_point_process_reference, run_simulated_network_reference from .paper_examples_full import run_full_paper_examples from .signal import Covariate, Signal from .simulation import simulate_poisson_from_rate @@ -85,7 +84,6 @@ def __getattr__(name: str): "FitSummary", "History", "HistoryBasis", - "matlab_engine_available", "NetworkSimulationResult", "ParityValidationError", "PointProcessSimulation", @@ -106,8 +104,6 @@ def __getattr__(name: str): "nstat_install", "psth", "run_full_paper_examples", - "run_point_process_reference", - "run_simulated_network_reference", "simulate_point_process", "simulate_poisson_from_rate", "simulate_two_neuron_network", diff --git a/nstat/analysis.py b/nstat/analysis.py index 472d58d1..55484313 100644 --- a/nstat/analysis.py +++ b/nstat/analysis.py @@ -6,10 +6,10 @@ from scipy.stats import chi2, norm from .SignalObj import SignalObj -from .fit import FitResult, _SingleFit +from .fit import FitResult, _SingleFit, _matlab_compute_ks_arrays from .glm import fit_binomial_glm, fit_poisson_glm from .signal import Covariate -from .trial import ConfigCollection, Trial +from .trial import ConfigCollection, SpikeTrainCollection, Trial def psth(spike_trains: Sequence[object], bin_edges: np.ndarray) -> tuple[np.ndarray, np.ndarray]: @@ -129,6 +129,20 @@ class Analysis: colors = ["b", "g", "r", "c", "m", "y", "k"] + @staticmethod + def _collapse_spike_input(nspikeObj): + if isinstance(nspikeObj, SpikeTrainCollection): + return nspikeObj.toSpikeTrain() + if isinstance(nspikeObj, Sequence) and not hasattr(nspikeObj, "spikeTimes"): + if len(nspikeObj) == 0: + raise ValueError("Spike input sequence must not be empty") + if any(not hasattr(item, "spikeTimes") for item in nspikeObj): + raise ValueError("Python Analysis expects sequences of MATLAB-style nspikeTrain objects") + if len(nspikeObj) == 1: + return nspikeObj[0] + return SpikeTrainCollection(list(nspikeObj)).toSpikeTrain() + return nspikeObj + @staticmethod def psth(spike_trains: Sequence[object], bin_edges: np.ndarray) -> tuple[np.ndarray, np.ndarray]: return psth(spike_trains, bin_edges) @@ -190,7 +204,6 @@ def GLMFit( 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) @@ -198,9 +211,13 @@ def GLMFit( 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) + # MATLAB stores logLL using the legacy per-bin convention + # `sum(y.*log(data*delta) + (1-y).*(1-data*delta))` for both GLM branches. + matlab_bin_mass = np.maximum(rate_hz / max(sample_rate, 1e-12), 1e-12) + logLL = float(np.sum(y * np.log(matlab_bin_mass) + (1.0 - y) * (1.0 - matlab_bin_mass))) + n_params = int(b.size) AIC = float(2.0 * n_params + dev) BIC = float(np.log(max(y.shape[0], 1)) * n_params + dev) @@ -371,50 +388,9 @@ def RunAnalysisForAllNeurons(tObj: Trial, configs: ConfigCollection, makePlot=1, return fits[0] if len(fits) == 1 else fits @staticmethod - def computeKSStats(nspikeObj, lambdaInput: Covariate, DTCorrection: int = 1): - del DTCorrection - if isinstance(nspikeObj, Sequence) and not hasattr(nspikeObj, "spikeTimes"): - if len(nspikeObj) != 1: - raise ValueError("Python computeKSStats currently expects a single spike train") - nspikeObj = nspikeObj[0] - - time = np.asarray(lambdaInput.time, dtype=float).reshape(-1) - data = np.asarray(lambdaInput.data, dtype=float) - if data.ndim == 1: - data = data[:, None] - dt = float(np.median(np.diff(time))) if time.size > 1 else 1.0 / max(lambdaInput.sampleRate, 1.0) - edges = np.concatenate([time, [time[-1] + dt]]) - counts = np.asarray(nspikeObj.to_binned_counts(edges), dtype=float).reshape(-1) - - z_cols: list[np.ndarray] = [] - u_cols: list[np.ndarray] = [] - x_cols: list[np.ndarray] = [] - ks_cols: list[np.ndarray] = [] - stats: list[float] = [] - for col in range(data.shape[1]): - lam_per_bin = np.clip(data[:, col].reshape(-1) * dt, 1e-12, None) - z = _time_rescaled_z(counts, lam_per_bin) - u = 1.0 - np.exp(-z) - ks_sorted = np.sort(u) - n_events = ks_sorted.shape[0] - if n_events: - x_axis = ((np.arange(1, n_events + 1, dtype=float) - 0.5) / n_events) - ks_stat = float(np.max(np.abs(ks_sorted - x_axis))) - else: - x_axis = np.asarray([], dtype=float) - ks_stat = 1.0 - z_cols.append(z) - u_cols.append(u) - x_cols.append(x_axis) - ks_cols.append(ks_sorted) - stats.append(ks_stat) - - Z = np.column_stack(z_cols) if z_cols and z_cols[0].size else np.zeros((0, data.shape[1]), dtype=float) - U = np.column_stack(u_cols) if u_cols and u_cols[0].size else np.zeros((0, data.shape[1]), dtype=float) - xAxis = np.column_stack(x_cols) if x_cols and x_cols[0].size else np.zeros((0, data.shape[1]), dtype=float) - KSSorted = np.column_stack(ks_cols) if ks_cols and ks_cols[0].size else np.zeros((0, data.shape[1]), dtype=float) - ks_stat = np.asarray(stats, dtype=float) - return Z, U, xAxis, KSSorted, ks_stat if ks_stat.size > 1 else float(ks_stat[0]) + def computeKSStats(nspikeObj, lambdaInput: Covariate, DTCorrection: int = 1, *, random_values=None): + nspikeObj = Analysis._collapse_spike_input(nspikeObj) + return _matlab_compute_ks_arrays(nspikeObj, lambdaInput, dt_correction=DTCorrection, random_values=random_values) @staticmethod def computeInvGausTrans(Z): @@ -446,21 +422,36 @@ def computeInvGausTrans(Z): @staticmethod def computeFitResidual(nspikeObj, lambdaInput: Covariate, windowSize: float = 0.01): - time = np.asarray(lambdaInput.time, dtype=float).reshape(-1) - data = np.asarray(lambdaInput.data, dtype=float) - if data.ndim == 1: - data = data[:, None] - dt = float(np.median(np.diff(time))) if time.size > 1 else 1.0 / max(lambdaInput.sampleRate, 1.0) - window = max(float(windowSize), dt) - edges = np.arange(float(time[0]), float(time[-1]) + window, window, dtype=float) - if edges[-1] < time[-1]: - edges = np.append(edges, time[-1] + window) - counts = np.asarray(nspikeObj.to_binned_counts(edges), dtype=float).reshape(-1) - out = np.zeros((counts.shape[0], data.shape[1]), dtype=float) - for col in range(data.shape[1]): - rate = np.interp(edges[:-1], time, data[:, col].reshape(-1)) - out[:, col] = counts - rate * window - return Covariate(edges[:-1], out, "M(t_k)", lambdaInput.xlabelval, lambdaInput.xunits, lambdaInput.yunits, list(lambdaInput.dataLabels)) + nspikeObj = Analysis._collapse_spike_input(nspikeObj) + + nCopy = nspikeObj.nstCopy() + nCopy.resample(lambdaInput.sampleRate) + nCopy.setMinTime(lambdaInput.minTime) + nCopy.setMaxTime(lambdaInput.maxTime) + + sumSpikes = nCopy.getSigRep(windowSize) + windowTimes = np.linspace(float(nCopy.minTime), float(nCopy.maxTime), sumSpikes.time.size, dtype=float) + lambdaInt = lambdaInput.integral() + lambdaIntVals = ( + lambdaInt.getValueAt(windowTimes[1:]).reshape(-1, lambdaInt.dimension) + - lambdaInt.getValueAt(windowTimes[:-1]).reshape(-1, lambdaInt.dimension) + ) + spike_window_data = np.asarray(sumSpikes.data, dtype=float) + if lambdaIntVals.shape[0] == spike_window_data.shape[0]: + sumSpikesOverWindow = spike_window_data + else: + sumSpikesOverWindow = spike_window_data[1:, :] + mdata = np.asarray(sumSpikesOverWindow, dtype=float) - np.asarray(lambdaIntVals, dtype=float) + out = np.vstack([np.zeros((1, mdata.shape[1]), dtype=float), mdata]) + return Covariate( + windowTimes, + out, + "M(t_k)", + lambdaInt.xlabelval, + lambdaInt.xunits, + lambdaInt.yunits, + list(lambdaInput.dataLabels), + ) @staticmethod def KSPlot(fitResults: FitResult, DTCorrection: int = 1, makePlot: int = 1): diff --git a/nstat/cif.py b/nstat/cif.py index afbd10a2..16331ae7 100644 --- a/nstat/cif.py +++ b/nstat/cif.py @@ -88,6 +88,56 @@ def _sigmoid(values: np.ndarray) -> np.ndarray: return 1.0 / (1.0 + np.exp(-np.clip(values, -20.0, 20.0))) +def _prepare_uniform_matrix( + values, + num_realizations: int, + num_draws: int, + *, + name: str, +) -> np.ndarray | None: + if values is None: + return None + arr = np.asarray(values, dtype=float) + if arr.ndim == 1: + if arr.size < num_draws: + raise ValueError(f"{name} must contain at least {num_draws} draws") + return np.tile(arr.reshape(1, -1), (int(num_realizations), 1)) + if arr.ndim != 2: + raise ValueError(f"{name} must be one- or two-dimensional") + if arr.shape[0] == int(num_realizations) and arr.shape[1] >= int(num_draws): + return arr + if arr.shape[1] == int(num_realizations) and arr.shape[0] >= int(num_draws): + return arr.T + raise ValueError( + f"{name} must have shape ({num_realizations}, >= {num_draws}) or " + f"(>= {num_draws}, {num_realizations})" + ) + + +def _prepare_optional_uniform_rows(values, num_realizations: int, *, name: str) -> np.ndarray | None: + if values is None: + return None + arr = np.asarray(values, dtype=float) + if arr.ndim == 1: + return np.tile(arr.reshape(1, -1), (int(num_realizations), 1)) + if arr.ndim != 2: + raise ValueError(f"{name} must be one- or two-dimensional") + if arr.shape[0] == int(num_realizations): + return arr + if arr.shape[1] == int(num_realizations): + return arr.T + raise ValueError(f"{name} must align with numRealizations={num_realizations}") + + +def _as_single_lambda_trace(lambda_covariate: Covariate) -> np.ndarray: + lambda_data = np.asarray(lambda_covariate.data, dtype=float) + if lambda_data.ndim == 1: + return lambda_data.reshape(-1) + if lambda_data.ndim == 2 and lambda_data.shape[1] == 1: + return lambda_data[:, 0].reshape(-1) + raise ValueError("simulateCIFByThinningFromLambda requires a single lambda trace") + + def _ordered_independent_names( xnames: Sequence[str] | None, stim_names: Sequence[str] | None, @@ -610,19 +660,120 @@ def simulateCIFByThinningFromLambda( maxTimeRes: float | None = None, *, seed: int | None = None, + random_values: np.ndarray | None = None, + thinning_values: np.ndarray | None = None, + return_details: bool = False, ) -> SpikeTrainCollection: - model = CIFModel(lambda_covariate.time, np.asarray(lambda_covariate.data, dtype=float).reshape(-1), getattr(lambda_covariate, "name", "lambda")) - coll = model.simulate(num_realizations=numRealizations, seed=seed) - if maxTimeRes is not None: - rounded = [] - for idx in range(1, coll.numSpikeTrains + 1): - train = coll.getNST(idx).nstCopy() - spikes = np.unique(np.ceil(train.spikeTimes / float(maxTimeRes)) * float(maxTimeRes)) - rounded.append(nspikeTrain(spikes, name=train.name, minTime=lambda_covariate.minTime, maxTime=lambda_covariate.maxTime, makePlots=-1)) - coll = SpikeTrainCollection(rounded) + if int(numRealizations) < 1: + raise ValueError("numRealizations must be >= 1") + + lambda_trace = _as_single_lambda_trace(lambda_covariate) + Tmax = float(lambda_covariate.maxTime) + lambda_bound = float(np.max(lambda_trace)) if lambda_trace.size else 0.0 + proposal_count = int(np.ceil(lambda_bound * (1.5 * Tmax))) if lambda_bound > 0.0 and Tmax > 0.0 else 0 + + arrival_uniforms = _prepare_uniform_matrix( + random_values, + int(numRealizations), + proposal_count, + name="random_values", + ) + thinning_uniforms = _prepare_optional_uniform_rows( + thinning_values, + int(numRealizations), + name="thinning_values", + ) + rng = np.random.default_rng(seed) + + trains: list[nspikeTrain] = [] + details = { + "lambda": lambda_trace.copy(), + "time": np.asarray(lambda_covariate.time, dtype=float).reshape(-1).copy(), + "lambda_bound": float(lambda_bound), + "proposal_count": int(proposal_count), + "arrival_uniforms": [], + "interarrival_times": [], + "candidate_spike_times": [], + "lambda_ratio": [], + "thinning_uniforms": [], + "accepted_spike_times": [], + } + + for realization in range(int(numRealizations)): + if proposal_count <= 0: + arrival = np.zeros(0, dtype=float) + interarrival = np.zeros(0, dtype=float) + candidate_spike_times = np.zeros(0, dtype=float) + lambda_ratio = np.zeros(0, dtype=float) + thinning = np.zeros(0, dtype=float) + accepted_spike_times = np.zeros(0, dtype=float) + else: + arrival = ( + np.asarray(arrival_uniforms[realization, :proposal_count], dtype=float).reshape(-1) + if arrival_uniforms is not None + else rng.random(proposal_count, dtype=float) + ) + arrival = np.clip(arrival, np.finfo(float).tiny, 1.0) + interarrival = -np.log(arrival) / lambda_bound + candidate_spike_times = np.cumsum(interarrival) + candidate_spike_times = candidate_spike_times[candidate_spike_times <= Tmax] + + if candidate_spike_times.size: + lambda_ratio = np.asarray(lambda_covariate.getValueAt(candidate_spike_times), dtype=float).reshape(-1) + lambda_ratio = np.clip(lambda_ratio / lambda_bound, 0.0, 1.0) + if thinning_uniforms is not None: + if thinning_uniforms.shape[1] < lambda_ratio.size: + raise ValueError( + "thinning_values must contain at least as many draws as surviving proposal times" + ) + thinning = np.asarray( + thinning_uniforms[realization, : lambda_ratio.size], + dtype=float, + ).reshape(-1) + else: + thinning = rng.random(lambda_ratio.size, dtype=float) + thinning = np.clip(thinning, 0.0, 1.0) + accepted_spike_times = candidate_spike_times[lambda_ratio >= thinning] + else: + lambda_ratio = np.zeros(0, dtype=float) + thinning = np.zeros(0, dtype=float) + accepted_spike_times = np.zeros(0, dtype=float) + + if maxTimeRes is not None and accepted_spike_times.size: + accepted_spike_times = np.unique( + np.ceil(accepted_spike_times / float(maxTimeRes)) * float(maxTimeRes) + ) + + train = nspikeTrain( + accepted_spike_times, + name="1", + minTime=lambda_covariate.minTime, + maxTime=lambda_covariate.maxTime, + makePlots=-1, + ) + train.setName("1") + trains.append(train) + + details["arrival_uniforms"].append(arrival.copy()) + details["interarrival_times"].append(interarrival.copy()) + details["candidate_spike_times"].append(candidate_spike_times.copy()) + details["lambda_ratio"].append(lambda_ratio.copy()) + details["thinning_uniforms"].append(thinning.copy()) + details["accepted_spike_times"].append(accepted_spike_times.copy()) + + coll = SpikeTrainCollection(trains) coll.setMinTime(lambda_covariate.minTime) coll.setMaxTime(lambda_covariate.maxTime) - return coll + if not return_details: + return coll + + if int(numRealizations) == 1: + detail_payload = { + key: (value[0] if isinstance(value, list) else value) + for key, value in details.items() + } + return coll, detail_payload + return coll, details @staticmethod def simulateCIFByThinning( @@ -664,6 +815,8 @@ def simulateCIF( *, seed: int | None = None, return_lambda: bool = False, + random_values: np.ndarray | None = None, + return_details: bool = False, ): if int(numRealizations) < 1: raise ValueError("numRealizations must be >= 1") @@ -698,8 +851,18 @@ def simulateCIF( raise ValueError("simType must be either poisson or binomial") lambda_data = np.zeros((time.size, int(numRealizations)), dtype=float) + lambda_delta_data = np.zeros((time.size, int(numRealizations)), dtype=float) + hist_effect_data = np.zeros((time.size, int(numRealizations)), dtype=float) + eta_data = np.zeros((time.size, int(numRealizations)), dtype=float) + if random_values is None: + rng = np.random.default_rng(seed) + draws = rng.random((time.size, int(numRealizations))) + else: + draws = np.asarray(random_values, dtype=float) + if draws.shape != (time.size, int(numRealizations)): + raise ValueError("random_values must have shape (len(time), numRealizations)") trains: list[nspikeTrain] = [] - rng = np.random.default_rng(seed) + spike_matrix = np.zeros((time.size, int(numRealizations)), dtype=float) mu_val = float(np.asarray(mu, dtype=float).reshape(-1)[0]) for realization in range(int(numRealizations)): @@ -710,15 +873,19 @@ def simulateCIF( if idx - lag >= 0: hist_effect += float(coeff) * float(spikes[idx - lag]) eta = mu_val + float(stim_drive[idx]) + float(ens_drive[idx]) + hist_effect + hist_effect_data[idx, realization] = hist_effect + eta_data[idx, realization] = eta if fit_type == "binomial": lambda_delta = float(_sigmoid(np.asarray([eta], dtype=float))[0]) rate_hz = lambda_delta / max(dt, 1e-12) - spikes[idx] = float(rng.random() < lambda_delta) + spikes[idx] = float(draws[idx, realization] < lambda_delta) else: rate_hz = float(np.exp(np.clip(eta, -20.0, 20.0))) lambda_delta = 1.0 - np.exp(-rate_hz * dt) - spikes[idx] = float(rng.random() < np.clip(lambda_delta, 0.0, 1.0)) + spikes[idx] = float(draws[idx, realization] < np.clip(lambda_delta, 0.0, 1.0)) lambda_data[idx, realization] = rate_hz + lambda_delta_data[idx, realization] = lambda_delta + spike_matrix[idx, realization] = spikes[idx] spike_times = time[spikes > 0.5] train = nspikeTrain(spike_times, name=str(realization + 1), minTime=float(time[0]), maxTime=float(time[-1]), makePlots=-1) trains.append(train) @@ -727,6 +894,20 @@ def simulateCIF( spikeTrainColl.setMinTime(float(time[0])) spikeTrainColl.setMaxTime(float(time[-1])) lambda_cov = Covariate(time, lambda_data, "\\lambda(t|H_t)", "time", "s", "Hz") + if return_details: + details = { + "time": time, + "stimulus_drive": stim_drive.copy(), + "ensemble_drive": ens_drive.copy(), + "history_effect": hist_effect_data, + "eta": eta_data, + "lambda_delta": lambda_delta_data, + "uniform_values": draws, + "spike_indicator": spike_matrix, + } + if return_lambda: + return spikeTrainColl, lambda_cov, details + return spikeTrainColl, details return (spikeTrainColl, lambda_cov) if return_lambda else spikeTrainColl @staticmethod diff --git a/nstat/confidence_interval.py b/nstat/confidence_interval.py index 4ec4c0ab..e8e13b50 100644 --- a/nstat/confidence_interval.py +++ b/nstat/confidence_interval.py @@ -1,26 +1,70 @@ from __future__ import annotations -from dataclasses import dataclass - import numpy as np -@dataclass class ConfidenceInterval: - time: np.ndarray - bounds: np.ndarray - color: str = "b" - - def __init__(self, time, bounds, color: str = "b") -> None: + def __init__(self, time, bounds, *args, color: str | None = None, value: float = 0.95) -> None: t = np.asarray(time, dtype=float).reshape(-1) b = np.asarray(bounds, dtype=float) if b.ndim != 2 or b.shape[1] != 2: raise ValueError("bounds must have shape (n, 2)") if b.shape[0] != t.shape[0]: raise ValueError("bounds rows must match time length") + use_color_alias = bool( + len(args) == 1 + and color is None + and isinstance(args[0], str) + and len(args[0]) <= 2 + ) + self.name = "" if use_color_alias or len(args) < 1 else str(args[0]) + self.xlabelval = "time" if use_color_alias or len(args) < 2 else str(args[1]) + self.xunits = "s" if use_color_alias or len(args) < 3 else str(args[2]) + self.yunits = "" if use_color_alias or len(args) < 4 else str(args[3]) + self.dataLabels = ( + ["lower", "upper"] + if use_color_alias or len(args) < 5 or args[4] is None + else self._normalize_text_sequence(args[4]) + ) + self.plotProps = ( + [] + if use_color_alias or len(args) < 6 or args[5] is None + else self._normalize_plot_props(args[5]) + ) self.time = t self.bounds = b - self.color = color + self.data = self.bounds + self.dimension = int(self.bounds.shape[1]) + self.minTime = float(self.time[0]) if self.time.size else float("nan") + self.maxTime = float(self.time[-1]) if self.time.size else float("nan") + self.sampleRate = ( + float(1.0 / np.mean(np.diff(self.time))) + if self.time.size > 1 and np.all(np.diff(self.time) != 0) + else float("nan") + ) + self.dataMask = np.ones(self.dimension, dtype=int) + self.color = str(args[0]) if use_color_alias else ("b" if color is None else str(color)) + self.value = float(value) + if len(self.plotProps) == 1 and self.dimension > 1: + self.plotProps = self.plotProps * self.dimension + + @staticmethod + def _normalize_text_sequence(values) -> list[str]: + if isinstance(values, (str, bytes)): + return [str(values)] + arr = np.asarray(values, dtype=object) + if arr.shape == (): + return [str(arr.item())] + return [str(item) for item in arr.reshape(-1)] + + @staticmethod + def _normalize_plot_props(values) -> list: + if isinstance(values, (str, bytes)): + return [str(values)] + arr = np.asarray(values, dtype=object) + if arr.shape == (): + return [arr.item()] + return [item.item() if isinstance(item, np.generic) else item for item in arr.reshape(-1)] @property def lower(self) -> np.ndarray: @@ -33,6 +77,29 @@ def upper(self) -> np.ndarray: def setColor(self, color: str) -> None: self.color = str(color) + def setValue(self, value: float) -> None: + self.value = float(value) + + def dataToMatrix(self) -> np.ndarray: + return np.asarray(self.bounds, dtype=float) + + def dataToStructure(self) -> dict: + return { + "time": self.time.tolist(), + "signals": {"values": self.bounds.tolist(), "dimensions": self.dimension}, + "name": self.name, + "dimension": self.dimension, + "minTime": self.minTime, + "maxTime": self.maxTime, + "xlabelval": self.xlabelval, + "xunits": self.xunits, + "yunits": self.yunits, + "dataLabels": list(self.dataLabels), + "dataMask": self.dataMask.tolist(), + "sampleRate": self.sampleRate, + "plotProps": list(self.plotProps), + } + def _coerce_signal_values(self, other) -> np.ndarray: if hasattr(other, "time") and hasattr(other, "data"): other_time = np.asarray(other.time, dtype=float).reshape(-1) @@ -78,8 +145,43 @@ def __rsub__(self, other): def __neg__(self): return ConfidenceInterval(self.time, np.column_stack([-self.upper, -self.lower]), self.color) - def plot(self, color: str | None = None, ax=None): + def toStructure(self) -> dict: + return self.dataToStructure() + + @staticmethod + def fromStructure(structure: dict) -> "ConfidenceInterval": + signals = structure.get("signals", {}) + values = signals.get("values", structure.get("data")) + ci = ConfidenceInterval( + structure["time"], + values, + structure.get("name", ""), + structure.get("xlabelval", "time"), + structure.get("xunits", "s"), + structure.get("yunits", ""), + structure.get("dataLabels"), + structure.get("plotProps"), + ) + if "dataMask" in structure: + ci.dataMask = np.asarray(structure["dataMask"], dtype=int).reshape(-1) + if "sampleRate" in structure: + ci.sampleRate = float(structure["sampleRate"]) + if "minTime" in structure: + ci.minTime = float(structure["minTime"]) + if "maxTime" in structure: + ci.maxTime = float(structure["maxTime"]) + return ci + + def plot(self, color: str | None = None, alphaVal: float = 0.2, drawPatches: int = 0, ax=None): import matplotlib.pyplot as plt axis = plt.gca() if ax is None else ax - return axis.fill_between(self.time, self.lower, self.upper, color=color or self.color, alpha=0.2) + plot_color = color or self.color + if drawPatches: + return axis.fill_between(self.time, self.lower, self.upper, color=plot_color, alpha=alphaVal) + lines = axis.plot(self.time, self.bounds) + for line in lines: + line.set_alpha(alphaVal) + if plot_color is not None and not isinstance(plot_color, (str, bytes)): + line.set_color(plot_color) + return lines diff --git a/nstat/decoding_algorithms.py b/nstat/decoding_algorithms.py index 3e48eaf1..3433391a 100644 --- a/nstat/decoding_algorithms.py +++ b/nstat/decoding_algorithms.py @@ -98,7 +98,11 @@ def _normalize_mu(mu, num_cells: int) -> np.ndarray: def _normalize_beta(beta, num_states: int, num_cells: int) -> np.ndarray: arr = np.asarray(beta, dtype=float) - if arr.ndim == 1: + if arr.ndim == 0: + if num_states != 1 or num_cells != 1: + raise ValueError("scalar beta is only valid for the 1-state, 1-cell MATLAB surface") + arr = arr.reshape(1, 1) + elif arr.ndim == 1: if num_cells == 1 and arr.size == num_states: arr = arr.reshape(num_states, 1) elif arr.size == num_cells and num_states == 1: @@ -774,19 +778,31 @@ def PPDecodeFilter(A, Q, Px0, dN, lambdaCIFColl, binwidth=0.001, x0=None, Pi0=No @staticmethod def PP_fixedIntervalSmoother(A, Q, dN, lags, mu, beta, fitType="poisson", delta=0.001, gamma=None, windowTimes=None, x0=None, Pi0=None): - x_p, W_p, x_u, W_u, *_ = DecodingAlgorithms._ppdecode_filter_linear( - A, - Q, - dN, - mu, - beta, - fitType, - delta, - gamma, - windowTimes, - x0, - Pi0, - ) + obs = _as_observation_matrix(dN) + num_cells, num_steps = obs.shape + num_states = _infer_state_dim(A, beta, num_cells) + mu_vec = _normalize_mu(mu, num_cells) + beta_mat = _normalize_beta(beta, num_states, num_cells) + + x0_vec = np.zeros(num_states, dtype=float) if _is_empty_value(x0) else np.asarray(x0, dtype=float).reshape(-1) + if x0_vec.size != num_states: + raise ValueError("x0 must match the decoding state dimension") + Pi0_mat = np.zeros((num_states, num_states), dtype=float) if _is_empty_value(Pi0) else _as_state_matrix(Pi0, num_states) + + if _is_empty_value(windowTimes): + H_tensor = np.zeros((num_steps, 0, num_cells), dtype=float) + gamma_mat = np.zeros((0, num_cells), dtype=float) + else: + H_tensor = _compute_history_terms(obs, float(delta), windowTimes) + gamma_mat = _normalize_gamma(gamma, H_tensor.shape[1], num_cells) + + x_p = np.zeros((num_states, num_steps + 1), dtype=float) + x_u = np.zeros((num_states, num_steps), dtype=float) + W_p = np.zeros((num_states, num_states, num_steps + 1), dtype=float) + W_u = np.zeros((num_states, num_states, num_steps), dtype=float) + x_p[:, 0] = x0_vec + W_p[:, :, 0] = Pi0_mat + lag_count = max(int(lags), 1) num_states, num_steps = x_u.shape x_uLag = np.zeros_like(x_u) @@ -795,6 +811,26 @@ def PP_fixedIntervalSmoother(A, Q, dN, lags, mu, beta, fitType="poisson", delta= W_pLag = np.zeros_like(W_p) for n in range(num_steps): + x_u[:, n], W_u[:, :, n], _ = DecodingAlgorithms.PPDecode_updateLinear( + x_p[:, n], + W_p[:, :, n], + obs, + mu_vec, + beta_mat, + fitType, + gamma_mat, + H_tensor, + n + 1, + None, + ) + A_t = _select_time_matrix(A, n, num_states) + Q_t = _select_time_matrix(Q, n, num_states) + x_p[:, n + 1], W_p[:, :, n + 1] = DecodingAlgorithms.PPDecode_predict( + x_u[:, n], + W_u[:, :, n], + A_t, + Q_t, + ) if n < lag_count: continue @@ -851,88 +887,159 @@ def PPHybridFilterLinear( if len(Q_models) != n_models: raise ValueError("A and Q must define the same number of hybrid models") - num_cells = obs.shape[0] + num_cells, num_steps = obs.shape state_dims = [_infer_state_dim(A_models[index], beta, num_cells) for index in range(n_models)] + max_dim = max(state_dims) mu_models = _normalize_mu_models(mu, n_models, num_cells) beta_models = _normalize_beta_models(beta, n_models, num_cells, state_dims) - x0_models = _normalize_model_sequence(x0, n_models, lambda index: np.zeros(state_dims[index], dtype=float)) - Pi0_models = _normalize_model_sequence(Pi0, n_models, lambda index: np.zeros((state_dims[index], state_dims[index]), dtype=float)) + x0_models_raw = _normalize_model_sequence(x0, n_models, lambda index: np.zeros(state_dims[index], dtype=float)) + Pi0_models_raw = _normalize_model_sequence(Pi0, n_models, lambda index: np.zeros((state_dims[index], state_dims[index]), dtype=float)) + x0_models = [np.asarray(x0_models_raw[index], dtype=float).reshape(-1) for index in range(n_models)] + Pi0_models = [_as_state_matrix(Pi0_models_raw[index], state_dims[index]) for index in range(n_models)] transition = np.asarray(p_ij, dtype=float) if transition.shape != (n_models, n_models): raise ValueError("p_ij must be an nModels x nModels transition matrix") - model_probs = _normalize_probabilities(Mu0) - if model_probs.size != n_models: - raise ValueError("Mu0 must contain one probability per hybrid model") + row_sums = np.sum(transition, axis=1) + if not np.allclose(row_sums, np.ones(n_models), atol=1e-8): + raise ValueError("State Transition probability matrix must sum to 1 along each row") + + if _is_empty_value(Mu0): + model_probs0 = np.full(n_models, 1.0 / float(n_models), dtype=float) + else: + model_probs0 = _normalize_probabilities(Mu0) + if model_probs0.size != n_models: + raise ValueError("Mu0 must contain one probability per hybrid model") if _is_empty_value(windowTimes): - H_tensor = np.zeros((obs.shape[1], 0, num_cells), dtype=float) + H_tensor = np.zeros((num_steps, 0, num_cells), dtype=float) gamma_mat = np.zeros((0, num_cells), dtype=float) else: H_tensor = _compute_history_terms(obs, float(binwidth), windowTimes) gamma_mat = _normalize_gamma(gamma, H_tensor.shape[1], num_cells) - model_results = [ - DecodingAlgorithms._ppdecode_filter_linear( - A_models[index], - Q_models[index], - obs, - mu_models[index], - beta_models[index], - fitType, - binwidth, - gamma_mat, - windowTimes, - x0_models[index], - Pi0_models[index], - ) - for index in range(n_models) - ] - - max_dim = max(state_dims) - num_steps = obs.shape[1] X = np.zeros((max_dim, num_steps), dtype=float) W = np.zeros((max_dim, max_dim, num_steps), dtype=float) + X_s = [np.zeros((max_dim, num_steps), dtype=float) for _ in range(n_models)] + W_s = [np.zeros((max_dim, max_dim, num_steps), dtype=float) for _ in range(n_models)] + X_u = [np.zeros((state_dims[index], num_steps), dtype=float) for index in range(n_models)] + W_u = [np.zeros((state_dims[index], state_dims[index], num_steps), dtype=float) for index in range(n_models)] + X_p = [np.zeros((state_dims[index], num_steps), dtype=float) for index in range(n_models)] + W_p = [np.zeros((state_dims[index], state_dims[index], num_steps), dtype=float) for index in range(n_models)] MU_u = np.zeros((n_models, num_steps), dtype=float) pNGivenS = np.zeros((n_models, num_steps), dtype=float) - X_s = [result[2] for result in model_results] - W_s = [result[3] for result in model_results] S_est = np.zeros(num_steps, dtype=int) + fit_type = str(fitType) + for time_index in range(num_steps): - predicted_probs = transition.T @ model_probs + if time_index == 0: + MU_p = transition.T @ model_probs0 + prev_probs = model_probs0 + else: + MU_p = transition.T @ MU_u[:, time_index - 1] + prev_probs = MU_u[:, time_index - 1] + + p_ij_s = transition * prev_probs[:, None] + column_norm = np.sum(p_ij_s, axis=0, keepdims=True) + column_norm[column_norm == 0.0] = 1.0 + p_ij_s = p_ij_s / column_norm + + for target_model in range(n_models): + mixed_state = np.zeros(max_dim, dtype=float) + for source_model in range(n_models): + dim_i = state_dims[source_model] + source_state = x0_models[source_model] if time_index == 0 else X_u[source_model][:, time_index - 1] + mixed_state[:dim_i] += source_state * p_ij_s[source_model, target_model] + X_s[target_model][:, time_index] = mixed_state + + mixed_cov = np.zeros((max_dim, max_dim), dtype=float) + for source_model in range(n_models): + dim_i = state_dims[source_model] + source_state = x0_models[source_model] if time_index == 0 else X_u[source_model][:, time_index - 1] + source_cov = Pi0_models[source_model] if time_index == 0 else W_u[source_model][:, :, time_index - 1] + diff = source_state - mixed_state[:dim_i] + mixed_cov[:dim_i, :dim_i] += ( + source_cov + np.outer(diff, diff) + ) * p_ij_s[source_model, target_model] + W_s[target_model][:, :, time_index] = _symmetrize(mixed_cov) + likelihoods = np.zeros(n_models, dtype=float) for model_index in range(n_models): - x_state = model_results[model_index][2][:, time_index] - lambda_delta = _lambda_delta_from_state( - x_state, + dim = state_dims[model_index] + A_t = _select_time_matrix(A_models[model_index], time_index, dim) + Q_t = _select_time_matrix(Q_models[model_index], time_index, dim) + pred_x, pred_W = DecodingAlgorithms.PPDecode_predict( + X_s[model_index][:dim, time_index], + W_s[model_index][:dim, :dim, time_index], + A_t, + Q_t, + ) + upd_x, upd_W, lambda_delta = DecodingAlgorithms.PPDecode_updateLinear( + pred_x, + pred_W, + obs, mu_models[model_index], beta_models[model_index], - str(fitType), + fit_type, gamma_mat, H_tensor, time_index + 1, + None, ) - likelihoods[model_index] = _likelihood_from_lambda(obs[:, time_index], lambda_delta, str(fitType)) - - weighted = likelihoods * predicted_probs - model_probs = _normalize_probabilities(weighted) - MU_u[:, time_index] = model_probs - pNGivenS[:, time_index] = _normalize_probabilities(likelihoods) + X_p[model_index][:, time_index] = pred_x + W_p[model_index][:, :, time_index] = pred_W + X_u[model_index][:, time_index] = upd_x + W_u[model_index][:, :, time_index] = upd_W + + det_ratio = np.sqrt(max(np.linalg.det(upd_W), 0.0)) / max(np.sqrt(max(np.linalg.det(pred_W), 0.0)), 1e-15) + log_term = np.sum(obs[:, time_index] * np.log(np.clip(lambda_delta.reshape(-1), 1e-12, np.inf)) - lambda_delta.reshape(-1)) + likelihoods[model_index] = float(det_ratio * np.exp(np.clip(log_term, -200.0, 50.0))) + + finite_likelihoods = likelihoods.copy() + finite_likelihoods[~np.isfinite(finite_likelihoods)] = 0.0 + pNGivenS[:, time_index] = finite_likelihoods + norm = np.sum(pNGivenS[:, time_index]) + if norm != 0.0 and np.isfinite(norm): + pNGivenS[:, time_index] /= norm + elif time_index > 0: + pNGivenS[:, time_index] = pNGivenS[:, time_index - 1] + else: + pNGivenS[:, time_index] = np.full(n_models, 0.5 if n_models == 2 else 1.0 / float(n_models), dtype=float) + + posterior = MU_p * pNGivenS[:, time_index] + posterior_norm = np.sum(posterior) + if posterior_norm != 0.0 and np.isfinite(posterior_norm): + MU_u[:, time_index] = posterior / posterior_norm + elif time_index > 0: + MU_u[:, time_index] = MU_u[:, time_index - 1] + else: + MU_u[:, time_index] = model_probs0 - best_model = int(np.argmax(model_probs)) + best_model = int(np.argmax(MU_u[:, time_index])) S_est[time_index] = best_model + 1 if MinClassificationError: chosen = best_model - X[: state_dims[chosen], time_index] = model_results[chosen][2][:, time_index] - W[: state_dims[chosen], : state_dims[chosen], time_index] = model_results[chosen][3][:, :, time_index] + dim = state_dims[chosen] + X[:dim, time_index] = X_u[chosen][:, time_index] + W[:dim, :dim, time_index] = W_u[chosen][:, :, time_index] continue + mixed_global_state = np.zeros(max_dim, dtype=float) for model_index in range(n_models): dim = state_dims[model_index] - X[:dim, time_index] += model_probs[model_index] * model_results[model_index][2][:, time_index] - W[:dim, :dim, time_index] += model_probs[model_index] * model_results[model_index][3][:, :, time_index] + mixed_global_state[:dim] += MU_u[model_index, time_index] * X_u[model_index][:, time_index] + X[:, time_index] = mixed_global_state + + mixed_global_cov = np.zeros((max_dim, max_dim), dtype=float) + for model_index in range(n_models): + dim = state_dims[model_index] + diff = X_u[model_index][:, time_index] - mixed_global_state[:dim] + mixed_global_cov[:dim, :dim] += MU_u[model_index, time_index] * ( + W_u[model_index][:, :, time_index] + np.outer(diff, diff) + ) + W[:, :, time_index] = _symmetrize(mixed_global_cov) return S_est, X, W, MU_u, X_s, W_s, pNGivenS diff --git a/nstat/fit.py b/nstat/fit.py index c3b6a9f4..24eef97f 100644 --- a/nstat/fit.py +++ b/nstat/fit.py @@ -8,7 +8,7 @@ matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np -from scipy.stats import kstest +from scipy.stats import norm from .core import Covariate, nspikeTrain @@ -72,15 +72,254 @@ def _time_rescaled_uniforms(y: np.ndarray, lam_per_bin: np.ndarray) -> np.ndarra return np.asarray(uniforms, dtype=float) +def _ksdiscrete( + pk: np.ndarray, + st: np.ndarray, + spikeflag: str, + *, + random_values: np.ndarray | None = None, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, np.ndarray]: + """Port of MATLAB Analysis.m local ksdiscrete helper.""" + + pk_arr = np.asarray(pk, dtype=float).reshape(-1) + if np.any(pk_arr < 0.0) or np.any(pk_arr > 1.0): + raise ValueError("all values for pk must be within [0,1]") + + if spikeflag == "spiketrain": + st_arr = np.asarray(st, dtype=float).reshape(-1) + if pk_arr.shape[0] != st_arr.shape[0]: + raise ValueError("pk and spike train must be same length") + spike_indices = np.flatnonzero(st_arr == 1.0) + 1 + elif spikeflag == "spikeind": + st_arr = np.asarray(st, dtype=float).reshape(-1) + spike_indices = np.unique(np.asarray(st_arr, dtype=int)) + else: + raise ValueError("spikeflag must be 'spiketrain' or 'spikeind'") + + if spike_indices.size == 0: + rst = pk_arr.copy() + return rst, np.sort(rst), np.asarray([], dtype=float), np.nan, np.sort(rst) + + if spike_indices[0] < 1: + raise ValueError("There is at least one spike with index less than 0") + if spike_indices[-1] > pk_arr.shape[0]: + raise ValueError("There is at least one spike with an index greater than the length of pk") + + with np.errstate(divide="ignore", invalid="ignore"): + qk = -np.log(1.0 - pk_arr) + n_spikes = int(spike_indices.size) + rst = np.zeros(max(n_spikes - 1, 0), dtype=float) + rstold = np.zeros_like(rst) + + if random_values is None: + draws = np.random.random_sample(rst.shape[0]) + else: + draws = np.asarray(random_values, dtype=float).reshape(-1) + if draws.shape[0] != rst.shape[0]: + raise ValueError("random_values must match the number of inter-spike intervals") + + for r in range(rst.shape[0]): + ind1 = int(spike_indices[r]) + ind2 = int(spike_indices[r + 1]) + total = float(np.sum(qk[ind1: ind2 - 1])) + qk_ind2 = float(qk[ind2 - 1]) + delta = -(1.0 / qk_ind2) * np.log(1.0 - float(draws[r]) * (1.0 - np.exp(-qk_ind2))) + if delta != 0.0: + total += qk_ind2 * delta + rst[r] = total + rstold[r] = float(np.sum(qk[ind1:ind2])) + + rstsort = np.sort(rst) + inrst = 1.0 / float(max(n_spikes - 1, 1)) + xrst = np.arange(0.5 * inrst, 1.0 - 0.5 * inrst + 0.5 * inrst, inrst, dtype=float) + cb = 1.36 * np.sqrt(inrst) + rstoldsort = np.sort(rstold) + return rst, rstsort, xrst, cb, rstoldsort + + +def _matlab_compute_ks_arrays( + spike_obj: nspikeTrain, + lambda_input: Covariate, + *, + dt_correction: int = 1, + random_values: np.ndarray | None = None, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, float | np.ndarray]: + """Port of MATLAB Analysis.computeKSStats for a single spike train.""" + + n_copy = spike_obj.nstCopy() + lambda_signal = Covariate( + np.asarray(lambda_input.time, dtype=float).reshape(-1), + np.asarray(lambda_input.data, dtype=float), + lambda_input.name, + lambda_input.xlabelval, + lambda_input.xunits, + lambda_input.yunits, + list(lambda_input.dataLabels) if getattr(lambda_input, "dataLabels", None) else [], + ) + + n_copy.resample(lambda_signal.sampleRate) + n_copy.setMinTime(lambda_signal.minTime) + n_copy.setMaxTime(lambda_signal.maxTime) + + rep_bin = n_copy.isSigRepBinary() + if not rep_bin: + lambda_signal = lambda_signal.resample(2.0 * lambda_signal.sampleRate) + n_copy.resample(lambda_signal.sampleRate) + + lambda_data = np.asarray(lambda_signal.data, dtype=float) + if lambda_data.ndim == 1: + lambda_data = lambda_data[:, None] + + n_dims = int(lambda_data.shape[1]) + if random_values is None: + random_cols = [None] * n_dims + else: + rv = np.asarray(random_values, dtype=float) + if rv.ndim == 1: + random_cols = [rv.reshape(-1)] * n_dims + elif rv.ndim == 2 and rv.shape[1] == n_dims: + random_cols = [rv[:, idx].reshape(-1) for idx in range(n_dims)] + else: + raise ValueError("random_values must be 1D or have one column per lambda dimension") + + z_cols: list[np.ndarray] = [] + u_cols: list[np.ndarray] = [] + x_cols: list[np.ndarray] = [] + ks_cols: list[np.ndarray] = [] + stats: list[float] = [] + + if int(dt_correction) == 1 and rep_bin: + pk = np.maximum(lambda_data * (1.0 / max(float(lambda_signal.sampleRate), 1e-12)), 1e-10) + spike_signal = np.asarray(n_copy.getSigRep().data, dtype=float).reshape(-1) + min_dim = min(pk.shape[0], spike_signal.shape[0]) + pk = pk[:min_dim, :] + spike_signal = spike_signal[:min_dim] + + int_cols: list[np.ndarray] = [] + for idx in range(n_dims): + pk_col = np.clip(pk[:, idx], 0.0, 1.0) + z_col, _, _, _, _ = _ksdiscrete( + pk_col, + spike_signal, + "spiketrain", + random_values=random_cols[idx], + ) + int_cols.append(np.asarray(z_col, dtype=float).reshape(-1)) + + if int_cols: + Z = _pad_rows(int_cols, fill_value=np.nan).T + else: + Z = np.zeros((0, n_dims), dtype=float) + else: + lambda_pos = np.maximum(lambda_data, 0.0) + lambda_cov = Covariate( + lambda_signal.time, + lambda_pos, + lambda_signal.name, + lambda_signal.xlabelval, + lambda_signal.xunits, + lambda_signal.yunits, + list(lambda_signal.dataLabels), + ) + lambda_int = lambda_cov.integral() + if n_copy.isSigRepBinary(): + spike_times = np.concatenate([[0.0], np.asarray(n_copy.getSpikeTimes(), dtype=float).reshape(-1)]) + else: + nst_signal = n_copy.getSigRep() + spike_times = np.concatenate([[0.0], np.asarray(nst_signal.time[np.asarray(nst_signal.data).reshape(-1) != 0], dtype=float).reshape(-1)]) + + if spike_times.size: + temp_vals = np.asarray(lambda_int.getValueAt(spike_times), dtype=float) + Z = temp_vals[1:, :] - temp_vals[:-1, :] + else: + Z = np.zeros((1, n_dims), dtype=float) + + U = 1.0 - np.exp(-Z) + if U.ndim == 1: + U = U[:, None] + + for idx in range(U.shape[1]): + ks_sorted = np.sort(np.asarray(U[:, idx], dtype=float).reshape(-1)) + n_events = int(ks_sorted.shape[0]) + if n_events: + x_axis = ((np.arange(1, n_events + 1, dtype=float) - 0.5) / float(n_events)) + ks_stat = float(np.max(np.abs(ks_sorted - x_axis))) + else: + x_axis = np.asarray([], dtype=float) + ks_stat = 1.0 + z_cols.append(np.asarray(Z[:, idx], dtype=float).reshape(-1)) + u_cols.append(np.asarray(U[:, idx], dtype=float).reshape(-1)) + x_cols.append(x_axis) + ks_cols.append(ks_sorted) + stats.append(ks_stat) + + Z_out = _pad_rows(z_cols, fill_value=np.nan).T if z_cols else np.zeros((0, n_dims), dtype=float) + U_out = _pad_rows(u_cols, fill_value=np.nan).T if u_cols else np.zeros((0, n_dims), dtype=float) + x_out = _pad_rows(x_cols, fill_value=np.nan).T if x_cols else np.zeros((0, n_dims), dtype=float) + ks_out = _pad_rows(ks_cols, fill_value=np.nan).T if ks_cols else np.zeros((0, n_dims), dtype=float) + ks_stat = np.asarray(stats, dtype=float) + if ks_stat.size == 1: + return Z_out, U_out, x_out, ks_out, float(ks_stat[0]) + return Z_out, U_out, x_out, ks_out, ks_stat + + def _ks_curve(uniforms: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: u = np.sort(np.asarray(uniforms, dtype=float).reshape(-1)) if u.size == 0: return np.asarray([], dtype=float), np.asarray([], dtype=float), np.asarray([], dtype=float) - ideal = np.linspace(1.0 / u.size, 1.0, u.size, dtype=float) + ideal = (np.arange(1, u.size + 1, dtype=float) - 0.5) / float(u.size) ci = np.full(u.size, 1.36 / np.sqrt(float(u.size)), dtype=float) return ideal, u, ci +def _matlab_kstest2( + x1: np.ndarray, + x2: np.ndarray, + *, + alpha: float = 0.05, + tail: str = "unequal", +) -> tuple[bool, float, float]: + """Port of MATLAB's asymptotic kstest2 implementation.""" + + sample1 = np.asarray(x1, dtype=float).reshape(-1) + sample2 = np.asarray(x2, dtype=float).reshape(-1) + sample1 = sample1[~np.isnan(sample1)] + sample2 = sample2[~np.isnan(sample2)] + if sample1.size == 0 or sample2.size == 0: + raise ValueError("kstest2 requires non-empty samples") + + tail_key = str(tail).lower() + if tail_key not in {"unequal", "smaller", "larger"}: + raise ValueError("tail must be 'unequal', 'smaller', or 'larger'") + + bin_edges = np.concatenate(([-np.inf], np.sort(np.concatenate((sample1, sample2))), [np.inf])) + cdf1 = np.histogram(sample1, bins=bin_edges)[0].cumsum(dtype=float) / float(sample1.size) + cdf2 = np.histogram(sample2, bins=bin_edges)[0].cumsum(dtype=float) / float(sample2.size) + + if tail_key == "unequal": + delta_cdf = np.abs(cdf1 - cdf2) + elif tail_key == "smaller": + delta_cdf = cdf2 - cdf1 + else: + delta_cdf = cdf1 - cdf2 + + ks_statistic = float(np.max(delta_cdf)) + n1 = float(sample1.size) + n2 = float(sample2.size) + n = (n1 * n2) / (n1 + n2) + lam = max((np.sqrt(n) + 0.12 + 0.11 / np.sqrt(n)) * ks_statistic, 0.0) + + if tail_key != "unequal": + p_value = float(np.exp(-2.0 * lam * lam)) + else: + j = np.arange(1.0, 102.0, dtype=float) + p_value = float(2.0 * np.sum(((-1.0) ** (j - 1.0)) * np.exp(-2.0 * lam * lam * j * j))) + p_value = min(max(p_value, 0.0), 1.0) + + h = bool(alpha >= p_value) + return h, p_value, ks_statistic + + def _extract_stat_component(stat: Any, candidates: Sequence[str]) -> Any: if stat is None: return None @@ -593,14 +832,40 @@ def _compute_diagnostics(self, fit_num: int = 1) -> dict[str, np.ndarray | float lam_per_bin = rate_hz * dt residual = counts - lam_per_bin residual_std = residual / np.sqrt(np.maximum(lam_per_bin, 1e-12)) - uniforms = _time_rescaled_uniforms(counts, lam_per_bin) - ideal, empirical, ci = _ks_curve(uniforms) - ks_stat = float(np.max(np.abs(empirical - ideal))) if ideal.size else 0.0 - ks_pvalue = float(kstest(uniforms, "uniform").pvalue) if uniforms.size else np.nan - within = float(np.mean(np.abs(empirical - ideal) <= ci)) if ideal.size else np.nan - lags, acf = _autocorrelation(uniforms, max_lag=25) - acf_ci = 1.96 / np.sqrt(float(uniforms.size)) if uniforms.size else np.nan - gauss = np.clip(uniforms, 1e-6, 1.0 - 1e-6) + + lambda_labels = list(self.lambda_signal.dataLabels) if getattr(self.lambda_signal, "dataLabels", None) else [] + if lambda_labels: + idx = min(max(fit_num - 1, 0), len(lambda_labels) - 1) + selected_labels = [str(lambda_labels[idx])] + else: + selected_labels = [f"\\lambda_{{{fit_num}}}"] + + lambda_signal = Covariate( + time, + rate_hz, + "\\lambda(t)", + self.lambda_signal.xlabelval, + self.lambda_signal.xunits, + self.lambda_signal.yunits, + selected_labels, + ) + Z, U, xAxis, KSSorted, _ = _matlab_compute_ks_arrays(self._primary_spike_train(), lambda_signal, dt_correction=1) + z = np.asarray(Z[:, 0], dtype=float).reshape(-1) if np.asarray(Z).size else np.asarray([], dtype=float) + uniforms = np.asarray(U[:, 0], dtype=float).reshape(-1) if np.asarray(U).size else np.asarray([], dtype=float) + ideal = np.asarray(xAxis[:, 0], dtype=float).reshape(-1) if np.asarray(xAxis).size else np.asarray([], dtype=float) + empirical = np.asarray(KSSorted[:, 0], dtype=float).reshape(-1) if np.asarray(KSSorted).size else np.asarray([], dtype=float) + ci = np.full(ideal.size, 1.36 / np.sqrt(float(ideal.size)), dtype=float) if ideal.size else np.asarray([], dtype=float) + ks_curve_stat = float(np.max(np.abs(empirical - ideal))) if ideal.size else 1.0 + if ideal.size: + different, ks_pvalue, ks_stat = _matlab_kstest2(ideal, empirical) + within = float(not different) + else: + ks_stat = 1.0 + ks_pvalue = np.nan + within = np.nan + gaussianized = norm.ppf(np.clip(uniforms, 1e-6, 1.0 - 1e-6)) + lags, acf = _autocorrelation(gaussianized, max_lag=25) + acf_ci = 1.96 / np.sqrt(float(gaussianized.size)) if gaussianized.size else np.nan coeffs = self.getCoeffs(fit_num) labels = self.covLabels[fit_num - 1] if fit_num - 1 < len(self.covLabels) else [] if coeffs.size == len(labels): @@ -616,33 +881,33 @@ def _compute_diagnostics(self, fit_num: int = 1) -> dict[str, np.ndarray | float "lambda_per_bin": lam_per_bin, "residual": residual, "residual_std": residual_std, + "z": z, "uniforms": uniforms, "ks_ideal": ideal, "ks_empirical": empirical, "ks_ci": ci, "ks_stat": ks_stat, + "ks_curve_stat": ks_curve_stat, "ks_pvalue": ks_pvalue, "within_conf_int": within, "acf_lags": lags, "acf_values": acf, "acf_ci": acf_ci, - "gaussianized": gauss, + "gaussianized": gaussianized, "coefficients": coeffs, "coeff_labels": np.asarray(coeff_labels, dtype=object), } self._diagnostic_cache[fit_num] = diagnostics - self.KSStats[fit_num - 1, 0] = ks_stat + self.setKSStats(z, uniforms, ideal, empirical, np.asarray([ks_stat], dtype=float)) self.KSPvalues[fit_num - 1] = ks_pvalue self.withinConfInt[fit_num - 1] = within - self.U = uniforms - self.Z = gauss - self.X = time + self.X = gaussianized self.Residual = { "time": time, "residual": residual, "standardized": residual_std, } - self.invGausStats = {"rhoSig": acf.tolist(), "confBoundSig": [acf_ci]} + self.invGausStats = {"X": gaussianized, "rhoSig": acf.tolist(), "confBoundSig": [acf_ci]} return diagnostics def computeKSStats(self, fit_num: int = 1) -> dict[str, float]: @@ -654,19 +919,57 @@ def computeKSStats(self, fit_num: int = 1) -> dict[str, float]: } def computeInvGausTrans(self, fit_num: int = 1) -> np.ndarray: - return np.asarray(self._compute_diagnostics(fit_num)["uniforms"], dtype=float) + return np.asarray(self._compute_diagnostics(fit_num)["gaussianized"], dtype=float) def computeFitResidual(self, fit_num: int = 1) -> Covariate: - diag = self._compute_diagnostics(fit_num) - return Covariate( - np.asarray(diag["time"], dtype=float), - np.asarray(diag["residual"], dtype=float), - "fit residual", - "time", - "s", - "counts/bin", - ["residual"], + time, rate_hz = self._lambda_series(fit_num) + if time.size == 0: + residual = Covariate([], [], "M(t_k)", "time", "s", "counts/bin", ["residual"]) + self.setFitResidual(residual) + return residual + + window_size = float(np.median(np.diff(time))) if time.size > 1 else 1.0 + spike_train = self._primary_spike_train().nstCopy() + spike_train.resample(1.0 / max(window_size, 1e-12)) + spike_train.setMinTime(float(time[0])) + spike_train.setMaxTime(float(time[-1])) + sum_spikes = spike_train.getSigRep(window_size, float(time[0]), float(time[-1])) + window_times = np.linspace(float(time[0]), float(time[-1]), sum_spikes.time.size, dtype=float) + + lambda_signal = Covariate( + time, + rate_hz, + "\\lambda(t)", + self.lambda_signal.xlabelval, + self.lambda_signal.xunits, + self.lambda_signal.yunits, + self.lambda_signal.dataLabels if getattr(self.lambda_signal, "dataLabels", None) else ["\\lambda"], + ) + lambda_int = lambda_signal.integral() + lambda_int_vals = ( + lambda_int.getValueAt(window_times[1:]).reshape(-1, lambda_int.dimension) + - lambda_int.getValueAt(window_times[:-1]).reshape(-1, lambda_int.dimension) + ) + + spike_window_data = np.asarray(sum_spikes.data, dtype=float) + if lambda_int_vals.shape[0] == spike_window_data.shape[0]: + sum_spikes_over_window = spike_window_data + else: + sum_spikes_over_window = spike_window_data[1:, :] + + mdata = np.asarray(sum_spikes_over_window, dtype=float) - np.asarray(lambda_int_vals, dtype=float) + residual_data = np.vstack([np.zeros((1, mdata.shape[1]), dtype=float), mdata]) + residual = Covariate( + window_times, + residual_data, + "M(t_k)", + lambda_int.xlabelval, + lambda_int.xunits, + lambda_int.yunits, + list(lambda_signal.dataLabels), ) + self.setFitResidual(residual) + return residual def evalLambda(self, fit_num: int = 1, newData=None) -> np.ndarray: coeffs = self.getCoeffs(fit_num) @@ -723,9 +1026,9 @@ def KSPlot(self, fit_num: int = 1, handle=None): return ax def plotResidual(self, fit_num: int = 1, handle=None): - diag = self._compute_diagnostics(fit_num) ax = handle if handle is not None else plt.subplots(1, 1, figsize=(6.0, 3.5))[1] - ax.plot(np.asarray(diag["time"], dtype=float), np.asarray(diag["residual"], dtype=float), color="tab:purple", linewidth=1.0) + residual = self.computeFitResidual(fit_num) + ax.plot(np.asarray(residual.time, dtype=float), np.asarray(residual.data[:, 0], dtype=float), color="tab:purple", linewidth=1.0) ax.axhline(0.0, color="0.4", linewidth=1.0, linestyle="--") ax.set_xlabel("time [s]") ax.set_ylabel("count residual") @@ -735,12 +1038,12 @@ def plotResidual(self, fit_num: int = 1, handle=None): def plotInvGausTrans(self, fit_num: int = 1, handle=None): diag = self._compute_diagnostics(fit_num) ax = handle if handle is not None else plt.subplots(1, 1, figsize=(6.0, 3.5))[1] - u = np.asarray(diag["uniforms"], dtype=float) - if u.size: - ax.plot(np.arange(1, u.size + 1), u, color="tab:green", linewidth=1.0) - ax.axhline(0.5, color="0.4", linewidth=1.0, linestyle="--") + x = np.asarray(diag["gaussianized"], dtype=float) + if x.size: + ax.plot(np.arange(1, x.size + 1), x, color="tab:green", linewidth=1.0) + ax.axhline(0.0, color="0.4", linewidth=1.0, linestyle="--") ax.set_xlabel("event index") - ax.set_ylabel("time-rescaled transform") + ax.set_ylabel("\\Phi^{-1}(u_i)") ax.set_title("Inverse-Gaussian/Uniform Transform") return ax @@ -805,10 +1108,25 @@ def plotHistCoeffs(self, fit_num: int = 1, sortByEpoch: int = 0, plotSignificanc def setKSStats(self, Z, U, xAxis, KSSorted, ks_stat): self.Z = np.asarray(Z, dtype=float) self.U = np.asarray(U, dtype=float) - self.X = np.asarray(xAxis, dtype=float) + self.KSXAxis = np.asarray(xAxis, dtype=float) self.KSSorted = np.asarray(KSSorted, dtype=float) - value = np.asarray(ks_stat, dtype=float).reshape(-1) - self.KSStats[: value.size, 0] = value + x_axis = np.asarray(xAxis, dtype=float) + ks_sorted = np.asarray(KSSorted, dtype=float) + if x_axis.ndim == 1: + x_axis = x_axis[:, None] + if ks_sorted.ndim == 1: + ks_sorted = ks_sorted[:, None] + + if x_axis.size and ks_sorted.size: + n_cols = min(x_axis.shape[1], ks_sorted.shape[1], self.numResults) + for idx in range(n_cols): + different, p_value, stat = _matlab_kstest2(x_axis[:, idx], ks_sorted[:, idx]) + self.KSStats[idx, 0] = stat + self.KSPvalues[idx] = p_value + self.withinConfInt[idx] = float(not different) + else: + value = np.asarray(ks_stat, dtype=float).reshape(-1) + self.KSStats[: value.size, 0] = value return self def setInvGausStats(self, X, rhoSig, confBoundSig): @@ -923,15 +1241,18 @@ def __init__(self, fit_results: FitResult | Iterable[FitResult]) -> None: self.fitNames = self.fitResCell[max(range(self.numNeurons), key=lambda idx: self.fitResCell[idx].numResults)].configNames self.neuronNumbers = [fr.neuronNumber for fr in self.fitResCell] - aic = _pad_rows([np.asarray(fr.AIC, dtype=float).reshape(-1) for fr in self.fitResCell]) - bic = _pad_rows([np.asarray(fr.BIC, dtype=float).reshape(-1) for fr in self.fitResCell]) - logll = _pad_rows([np.asarray(fr.logLL, dtype=float).reshape(-1) for fr in self.fitResCell]) - ks = _pad_rows([np.asarray(fr.KSStats, dtype=float).reshape(-1) for fr in self.fitResCell], fill_value=np.nan) - - self.AIC = np.nanmean(aic, axis=0) - self.BIC = np.nanmean(bic, axis=0) - self.logLL = np.nanmean(logll, axis=0) - self.KSStats = np.column_stack([np.nanmean(ks, axis=0), np.nanstd(ks, axis=0)]) + self.dev = _pad_rows([np.asarray(fr.dev, dtype=float).reshape(-1) for fr in self.fitResCell]) + self.AIC = _pad_rows([np.asarray(fr.AIC, dtype=float).reshape(-1) for fr in self.fitResCell]) + self.BIC = _pad_rows([np.asarray(fr.BIC, dtype=float).reshape(-1) for fr in self.fitResCell]) + self.logLL = _pad_rows([np.asarray(fr.logLL, dtype=float).reshape(-1) for fr in self.fitResCell]) + self.KSStats = _pad_rows([np.asarray(fr.KSStats, dtype=float).reshape(-1) for fr in self.fitResCell], fill_value=np.nan) + self.KSPvalues = _pad_rows([np.asarray(fr.KSPvalues, dtype=float).reshape(-1) for fr in self.fitResCell], fill_value=np.nan) + self.withinConfInt = _pad_rows([np.asarray(fr.withinConfInt, dtype=float).reshape(-1) for fr in self.fitResCell], fill_value=np.nan) + self.meanAIC = np.nanmean(self.AIC, axis=0) + self.meanBIC = np.nanmean(self.BIC, axis=0) + self.meanlogLL = np.nanmean(self.logLL, axis=0) + self.meanKSStats = np.nanmean(self.KSStats, axis=0) + self.stdKSStats = np.nanstd(self.KSStats, axis=0) self.uniqueCovLabels: list[str] = [] self.coeffMin = np.nan self.coeffMax = np.nan @@ -939,16 +1260,22 @@ def __init__(self, fit_results: FitResult | Iterable[FitResult]) -> None: self.mapCovLabelsToUniqueLabels() def getDiffAIC(self, idx: int = 1) -> np.ndarray: - base = self.AIC[idx - 1] - return self.AIC - base + if self.numResults > 1: + keep = [col for col in range(self.AIC.shape[1]) if col != (idx - 1)] + return self.AIC[:, keep] - self.AIC[:, [idx - 1]] + return self.AIC.copy() def getDiffBIC(self, idx: int = 1) -> np.ndarray: - base = self.BIC[idx - 1] - return self.BIC - base + if self.numResults > 1: + keep = [col for col in range(self.BIC.shape[1]) if col != (idx - 1)] + return self.BIC[:, keep] - self.BIC[:, [idx - 1]] + return self.BIC.copy() def getDifflogLL(self, idx: int = 1) -> np.ndarray: - base = self.logLL[idx - 1] - return self.logLL - base + if self.numResults > 1: + keep = [col for col in range(self.logLL.shape[1]) if col != (idx - 1)] + return self.logLL[:, keep] - self.logLL[:, [idx - 1]] + return self.logLL.copy() def mapCovLabelsToUniqueLabels(self): self.uniqueCovLabels = _ordered_unique( @@ -1035,21 +1362,21 @@ def plotIC(self, handle=None): def plotAIC(self, handle=None): ax = handle if handle is not None else plt.subplots(1, 1, figsize=(5.0, 3.5))[1] - ax.boxplot(_pad_rows([np.asarray(fit.AIC, dtype=float) for fit in self.fitResCell]).T, labels=self.fitNames) + ax.boxplot(self.AIC, labels=self.fitNames) ax.set_ylabel("AIC") ax.set_title("AIC Across Neurons") return ax def plotBIC(self, handle=None): ax = handle if handle is not None else plt.subplots(1, 1, figsize=(5.0, 3.5))[1] - ax.boxplot(_pad_rows([np.asarray(fit.BIC, dtype=float) for fit in self.fitResCell]).T, labels=self.fitNames) + ax.boxplot(self.BIC, labels=self.fitNames) ax.set_ylabel("BIC") ax.set_title("BIC Across Neurons") return ax def plotlogLL(self, handle=None): ax = handle if handle is not None else plt.subplots(1, 1, figsize=(5.0, 3.5))[1] - ax.boxplot(_pad_rows([np.asarray(fit.logLL, dtype=float) for fit in self.fitResCell]).T, labels=self.fitNames) + ax.boxplot(self.logLL, labels=self.fitNames) ax.set_ylabel("log likelihood") ax.set_title("log likelihood Across Neurons") return ax @@ -1075,7 +1402,7 @@ def plotSummary(self, handle=None): labels = list(self.fitNames) for ax, values, title in zip( axes, - (self.AIC, self.BIC, self.logLL), + (self.meanAIC, self.meanBIC, self.meanlogLL), ("AIC", "BIC", "log likelihood"), strict=False, ): @@ -1087,12 +1414,19 @@ def plotSummary(self, handle=None): return fig def boxPlot(self, X, diffIndex: int = 1, h=None, dataLabels=None, **kwargs): - del diffIndex, kwargs + del kwargs ax = h if h is not None else plt.subplots(1, 1, figsize=(6.0, 3.5))[1] values = np.asarray(X, dtype=float) - labels = list(dataLabels) if dataLabels is not None else list(self.fitNames[: values.shape[1] if values.ndim == 2 else 1]) if values.ndim == 1: values = values[:, None] + if dataLabels is not None: + labels = list(dataLabels) + elif values.shape[1] == len(self.fitNames): + labels = list(self.fitNames) + elif values.shape[1] == max(len(self.fitNames) - 1, 1): + labels = [name for idx, name in enumerate(self.fitNames, start=1) if idx != diffIndex] + else: + labels = list(self.fitNames[: values.shape[1]]) ax.boxplot(values, labels=labels) return ax @@ -1102,6 +1436,13 @@ def toStructure(self) -> dict[str, Any]: "numNeurons": self.numNeurons, "numResults": self.numResults, "fitNames": list(self.fitNames), + "dev": self.dev.tolist(), + "AIC": self.AIC.tolist(), + "BIC": self.BIC.tolist(), + "logLL": self.logLL.tolist(), + "KSStats": self.KSStats.tolist(), + "KSPvalues": self.KSPvalues.tolist(), + "withinConfInt": self.withinConfInt.tolist(), } @staticmethod diff --git a/nstat/matlab_reference.py b/nstat/matlab_reference.py deleted file mode 100644 index 576606b9..00000000 --- a/nstat/matlab_reference.py +++ /dev/null @@ -1,191 +0,0 @@ -from __future__ import annotations - -from functools import lru_cache -from pathlib import Path -from typing import Any - -import numpy as np - - -def _repo_root() -> Path: - return Path(__file__).resolve().parents[1] - - -def _default_matlab_repo() -> Path: - return _repo_root().parent / "nSTAT" - - -def matlab_engine_available() -> bool: - try: - import matlab.engine # type: ignore - except Exception: - return False - return True - - -def _to_numpy(value: Any) -> np.ndarray: - if isinstance(value, np.ndarray): - return value - try: - return np.asarray(value, dtype=float) - except Exception: - if hasattr(value, "_data") and hasattr(value, "size"): - return np.asarray(value._data, dtype=float).reshape(value.size, order="F") - raise - - -@lru_cache(maxsize=1) -def start_matlab_engine(): - if not matlab_engine_available(): - raise RuntimeError("MATLAB Engine for Python is not available in this environment") - import matlab.engine # type: ignore - - return matlab.engine.start_matlab() - - -def _add_repo_to_path(engine, matlab_repo: Path) -> None: - engine.addpath(str(matlab_repo), nargout=0) - engine.addpath(str(matlab_repo / "helpfiles"), nargout=0) - - -def run_point_process_reference(*, matlab_repo: str | Path | None = None, seed: int = 5) -> 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( - f""" - rng({int(seed)}); - Ts=.001; tMin=0; tMax=50; t=tMin:Ts:tMax; - mu=-3; - H=tf([-1 -2 -4],[1],Ts,'Variable','z^-1'); - S=tf([1],1,Ts,'Variable','z^-1'); - E=tf([0],1,Ts,'Variable','z^-1'); - u=sin(2*pi*1*t)'; - e=zeros(length(t),1); - stim=Covariate(t',u,'Stimulus','time','s','Voltage',{'sin'}); - ens=Covariate(t',e,'Ensemble','time','s','Spikes',{'n1'}); - [sC, lambda] = CIF.simulateCIF(mu,H,S,E,stim,ens,5,'binomial'); - ppSpikeCounts = zeros(1, sC.numSpikeTrains); - for i=1:sC.numSpikeTrains - ppSpikeCounts(i) = length(sC.getNST(i).spikeTimes); - end - ppLambdaHead = lambda.data(1:10,1)'; - """, - nargout=0, - ) - return { - "spike_counts": _to_numpy(engine.workspace["ppSpikeCounts"]).reshape(-1), - "lambda_head": _to_numpy(engine.workspace["ppLambdaHead"]).reshape(-1), - } - - -def run_simulated_network_reference(*, matlab_repo: str | Path | None = None, seed: int = 4) -> 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( - f""" - rng({int(seed)}); - Ts=.001; tMin=0; tMax=50; t=tMin:Ts:tMax; - mu{1}=-3; mu{2}=-3; - H{1}=tf([-4 -2 -1],[1],Ts,'Variable','z^-1'); - H{2}=tf([-4 -2 -1],[1],Ts,'Variable','z^-1'); - S{1}=tf([1],1,Ts,'Variable','z^-1'); - S{2}=tf([-1],1,Ts,'Variable','z^-1'); - E{1}=tf([1],1,Ts,'Variable','z^-1'); - E{2}=tf([-4],1,Ts,'Variable','z^-1'); - u = sin(2*pi*1*t)'; - stim=Covariate(t',u,'Stimulus','time','s','Voltage',{'sin'}); - assignin('base','S1',S{1}); assignin('base','H1',H{1}); assignin('base','E1',E{1}); assignin('base','mu1',mu{1}); - assignin('base','S2',S{2}); assignin('base','H2',H{2}); assignin('base','E2',E{2}); assignin('base','mu2',mu{2}); - options = simget; - [tout,~,yout] = sim('SimulatedNetwork2',[stim.minTime stim.maxTime],options,stim.dataToStructure); - [h1Num, ~] = tfdata(H{1}, 'v'); - [h2Num, ~] = tfdata(H{2}, 'v'); - [s1Num, ~] = tfdata(S{1}, 'v'); - [s2Num, ~] = tfdata(S{2}, 'v'); - [e1Num, ~] = tfdata(E{1}, 'v'); - [e2Num, ~] = tfdata(E{2}, 'v'); - stateMat = yout(:,1:2); - probMat = zeros(size(stateMat)); - for n = 1:size(stateMat, 1) - hist1 = 0; hist2 = 0; - for lag = 1:length(h1Num) - if n-lag >= 1 - hist1 = hist1 + h1Num(lag) * stateMat(n-lag,1); - hist2 = hist2 + h2Num(lag) * stateMat(n-lag,2); - end - end - ens1 = 0; ens2 = 0; - if n > 1 - ens1 = e1Num(1) * stateMat(n-1,2); - ens2 = e2Num(1) * stateMat(n-1,1); - end - eta1 = mu{1} + hist1 + s1Num(1) * u(n) + ens1; - eta2 = mu{2} + hist2 + s2Num(1) * u(n) + ens2; - probMat(n,1) = exp(eta1) / (1 + exp(eta1)); - probMat(n,2) = exp(eta2) / (1 + exp(eta2)); - end - netSpikeCounts = [sum(yout(:,1)>.5), sum(yout(:,2)>.5)]; - netProbHead = probMat(1:5,:); - netStateHead = yout(1:5,1:2); - netActual = [0 1; -4 0]; - """, - nargout=0, - ) - return { - "spike_counts": _to_numpy(engine.workspace["netSpikeCounts"]).reshape(-1), - "prob_head": _to_numpy(engine.workspace["netProbHead"]), - "state_head": _to_numpy(engine.workspace["netStateHead"]), - "actual_network": _to_numpy(engine.workspace["netActual"]), - } - - -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})); - summary = FitResSummary({fit}); - analysisAIC = fit.AIC(1); - analysisBIC = fit.BIC(1); - analysisLogLL = fit.logLL(1); - analysisCoeffs = fit.getCoeffs(1)'; - analysisLambdaHead = fit.lambda.data(1:5, 1)'; - analysisSummaryAIC = summary.AIC(1); - analysisSummaryBIC = summary.BIC(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), - "summary_aic": _to_numpy(engine.workspace["analysisSummaryAIC"]).reshape(-1), - "summary_bic": _to_numpy(engine.workspace["analysisSummaryBIC"]).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 07e08fe7..fa51512e 100644 --- a/nstat/parity_report.py +++ b/nstat/parity_report.py @@ -173,6 +173,9 @@ def render_parity_report(repo_root: Path | None = None) -> str: "- Public API: no missing MATLAB public APIs remain; only the MATLAB help-browser utility is explicitly non-applicable." ) lines.append("- Help/notebook parity: all inventoried MATLAB help workflows are mapped to Python notebooks or equivalents.") + lines.append( + "- Clean-room boundary: the installed Python package, pure-Python tests, notebooks, examples, and default CI no longer invoke MATLAB; cross-language execution is confined to the MATLAB-side `tests/python_port_fidelity` harness." + ) if notebook_partial: lines.append( f"- Notebook fidelity: workflow coverage is complete, but {len(notebook_partial)} MATLAB-helpfile notebook ports are still marked partial in `tools/notebooks/parity_notes.yml`." @@ -209,7 +212,7 @@ def render_parity_report(repo_root: Path | None = None) -> str: ) 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." + f"- Simulink fidelity: native Python coverage exists for the required published workflows, deterministic injected-input fixtures now back the required native paths, 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.") diff --git a/nstat/release_check.py b/nstat/release_check.py deleted file mode 100644 index 800a11ae..00000000 --- a/nstat/release_check.py +++ /dev/null @@ -1,98 +0,0 @@ -from __future__ import annotations - -import argparse -import shutil -import subprocess -import sys -from pathlib import Path - - -def _repo_root() -> Path: - return Path(__file__).resolve().parents[1] - - -def default_matlab_repo_root(repo_root: Path | None = None) -> Path: - base = _repo_root() if repo_root is None else repo_root.resolve() - return base.parent / "nSTAT" - - -def _matlab_quote(value: str) -> str: - return value.replace("'", "''") - - -def build_release_gate_commands( - repo_root: Path | None = None, - *, - matlab_repo_root: Path | None = None, - skip_matlab: bool = False, -) -> list[list[str]]: - base = _repo_root() if repo_root is None else repo_root.resolve() - matlab_root = default_matlab_repo_root(base) if matlab_repo_root is None else matlab_repo_root.resolve() - - commands: list[list[str]] = [ - [sys.executable, str(base / "tools" / "parity" / "export_matlab_gold_fixtures.py"), "--repo-root", str(base), "--matlab-repo", str(matlab_root)], - [ - sys.executable, - "-m", - "pytest", - "-q", - "tests/test_zernike.py", - "tests/test_matlab_gold_fixtures.py", - "tests/test_signalobj_fidelity.py", - "tests/test_nspiketrain_fidelity.py", - "tests/test_workflow_fidelity.py", - "tests/test_class_fidelity_audit.py", - "tests/test_matlab_symbol_surface.py", - "tests/test_matlab_reference.py", - "tests/test_simulink_fidelity_audit.py", - "tests/test_parity_report.py", - ], - [sys.executable, str(base / "tools" / "notebooks" / "run_notebooks.py"), "--group", "parity_core", "--timeout", "900"], - [sys.executable, str(base / "tools" / "parity" / "build_report.py")], - ] - if not skip_matlab: - commands.append( - [ - "matlab", - "-batch", - ( - f"cd('{_matlab_quote(str(matlab_root))}'); " - "addpath(pwd); " - "addpath(fullfile(pwd,'helpfiles')); " - "results = runtests('tests/python_port_fidelity'); " - "assertSuccess(results); exit" - ), - ] - ) - return commands - - -def run_release_gate( - repo_root: Path | None = None, - *, - matlab_repo_root: Path | None = None, - skip_matlab: bool = False, -) -> None: - base = _repo_root() if repo_root is None else repo_root.resolve() - if not skip_matlab and shutil.which("matlab") is None: - raise FileNotFoundError("`matlab` is not on PATH; rerun with --skip-matlab or install MATLAB CLI access") - for command in build_release_gate_commands(base, matlab_repo_root=matlab_repo_root, skip_matlab=skip_matlab): - subprocess.run(command, cwd=base, check=True) - - -def main(argv: list[str] | None = None) -> int: - parser = argparse.ArgumentParser(description="Run the Python plus MATLAB fidelity release gate.") - parser.add_argument("--repo-root", type=Path, default=None) - parser.add_argument("--matlab-repo", type=Path, default=None) - parser.add_argument("--skip-matlab", action="store_true") - args = parser.parse_args(argv) - - run_release_gate(args.repo_root, matlab_repo_root=args.matlab_repo, skip_matlab=args.skip_matlab) - return 0 - - -__all__ = ["build_release_gate_commands", "default_matlab_repo_root", "main", "run_release_gate"] - - -if __name__ == "__main__": - raise SystemExit(main(sys.argv[1:])) diff --git a/nstat/simulators.py b/nstat/simulators.py index 07611e2f..68a70f1c 100644 --- a/nstat/simulators.py +++ b/nstat/simulators.py @@ -12,6 +12,9 @@ class PointProcessSimulation: time: np.ndarray rate_hz: np.ndarray spikes: SpikeTrain + lambda_delta: np.ndarray | None = None + spike_indicator: np.ndarray | None = None + uniform_values: np.ndarray | None = None @dataclass @@ -25,9 +28,20 @@ class NetworkSimulationResult: stimulus_kernel: np.ndarray ensemble_kernel: np.ndarray baseline_mu: np.ndarray - - -def simulate_point_process(time: np.ndarray, rate_hz: np.ndarray, *, seed: int | None = None) -> PointProcessSimulation: + eta: np.ndarray | None = None + history_effect: np.ndarray | None = None + ensemble_effect: np.ndarray | None = None + spike_indicator: np.ndarray | None = None + uniform_values: np.ndarray | None = None + + +def simulate_point_process( + time: np.ndarray, + rate_hz: np.ndarray, + *, + seed: int | None = None, + uniform_values: np.ndarray | None = None, +) -> PointProcessSimulation: t = np.asarray(time, dtype=float).reshape(-1) r = np.asarray(rate_hz, dtype=float).reshape(-1) if t.shape[0] != r.shape[0]: @@ -40,9 +54,22 @@ def simulate_point_process(time: np.ndarray, rate_hz: np.ndarray, *, seed: int | p = 1.0 - np.exp(-np.clip(r, 0.0, np.inf) * dt) p = np.clip(p, 0.0, 1.0) - rng = np.random.default_rng(seed) - keep = rng.random(t.shape[0]) < p - return PointProcessSimulation(t, r, SpikeTrain(t[keep])) + if uniform_values is None: + rng = np.random.default_rng(seed) + draws = rng.random(t.shape[0]) + else: + draws = np.asarray(uniform_values, dtype=float).reshape(-1) + if draws.shape[0] != t.shape[0]: + raise ValueError("uniform_values must match the length of time") + keep = draws < p + return PointProcessSimulation( + t, + r, + SpikeTrain(t[keep]), + lambda_delta=p, + spike_indicator=keep.astype(float), + uniform_values=draws, + ) def simulate_two_neuron_network( @@ -54,6 +81,7 @@ def simulate_two_neuron_network( ensemble_kernel: tuple[float, float] = (1.0, -4.0), stimulus_frequency_hz: float = 1.0, seed: int | None = 13, + uniform_values: np.ndarray | None = None, ) -> NetworkSimulationResult: """Standalone Python replacement for the MATLAB/Simulink 2-neuron NetworkTutorial.""" if duration_s <= 0 or dt <= 0: @@ -73,9 +101,18 @@ def simulate_two_neuron_network( dtype=float, ) - rng = np.random.default_rng(seed) spikes = np.zeros((time.shape[0], 2), dtype=float) lambda_delta = np.zeros_like(spikes) + eta_trace = np.zeros_like(spikes) + history_effect = np.zeros_like(spikes) + ensemble_effect = np.zeros_like(spikes) + if uniform_values is None: + rng = np.random.default_rng(seed) + draws = rng.random(spikes.shape) + else: + draws = np.asarray(uniform_values, dtype=float) + if draws.shape != spikes.shape: + raise ValueError("uniform_values must have shape (len(time), 2)") for i in range(time.shape[0]): hist_self = np.zeros(2, dtype=float) for lag, coeff in enumerate(history_kernel_arr, start=1): @@ -87,9 +124,12 @@ def simulate_two_neuron_network( ens_effect[0] = ensemble_kernel_arr[0] * float(spikes[i - 1, 1]) ens_effect[1] = ensemble_kernel_arr[1] * float(spikes[i - 1, 0]) eta = baseline_mu_arr + hist_self + (stimulus_kernel_arr * float(drive[i])) + ens_effect + history_effect[i] = hist_self + ensemble_effect[i] = ens_effect + eta_trace[i] = eta lambda_delta[i] = 1.0 / (1.0 + np.exp(-np.clip(eta, -20.0, 20.0))) - spikes[i, 0] = 1.0 if rng.random() < lambda_delta[i, 0] else 0.0 - spikes[i, 1] = 1.0 if rng.random() < lambda_delta[i, 1] else 0.0 + spikes[i, 0] = 1.0 if draws[i, 0] < lambda_delta[i, 0] else 0.0 + spikes[i, 1] = 1.0 if draws[i, 1] < lambda_delta[i, 1] else 0.0 t1 = time[spikes[:, 0] > 0.5] t2 = time[spikes[:, 1] > 0.5] @@ -104,6 +144,11 @@ def simulate_two_neuron_network( stimulus_kernel=stimulus_kernel_arr, ensemble_kernel=ensemble_kernel_arr, baseline_mu=baseline_mu_arr, + eta=eta_trace, + history_effect=history_effect, + ensemble_effect=ensemble_effect, + spike_indicator=spikes, + uniform_values=draws, ) diff --git a/nstat/trial.py b/nstat/trial.py index ee82f954..3bc0183c 100644 --- a/nstat/trial.py +++ b/nstat/trial.py @@ -551,6 +551,67 @@ def getNSTIndicesFromName(self, name: Sequence[str] | str): return matches if len(matches) > 1 else matches[0] return [self.getNSTIndicesFromName(item) for item in name] + def toSpikeTrain( + self, + selectorArray: Sequence[int] | Sequence[str] | str | None = None, + minTime: float | None = None, + maxTime: float | None = None, + windowTimes: Sequence[float] | None = None, + ) -> nspikeTrain: + if self.numSpikeTrains == 0: + raise ValueError("nstColl.toSpikeTrain requires at least one spike train") + + if maxTime is None: + maxTime = self.maxTime + if minTime is None: + minTime = self.minTime + + if selectorArray is None: + selector = self.getIndFromMask() if self.isNeuronMaskSet() else list(range(1, self.numSpikeTrains + 1)) + elif isinstance(selectorArray, str) or _is_string_sequence(selectorArray): + resolved = self.getNSTIndicesFromName(selectorArray) + if isinstance(resolved, list): + selector = [int(item) if not isinstance(item, list) else int(item[0]) for item in resolved] + else: + selector = [int(resolved)] + else: + selector = [int(item) for item in selectorArray] + + if not selector: + raise ValueError("selectorArray resolved to no spike trains") + + delta = 1.0 / max(float(self.sampleRate), 1e-12) + spike_times: list[float] = [] + offset = 0.0 + selected_trains = [self.getNST(index) for index in selector] + name = selected_trains[0].name + + if windowTimes is None or len(windowTimes) == 0: + for idx, train in enumerate(selected_trains): + if idx == 0: + spike_times.extend(np.asarray(train.spikeTimes, dtype=float).reshape(-1).tolist()) + else: + prev_train = selected_trains[idx - 1] + offset += float(prev_train.maxTime) + float(delta) + if np.asarray(train.spikeTimes).size: + spike_times.extend((np.asarray(train.spikeTimes, dtype=float).reshape(-1) + offset).tolist()) + else: + window_arr = np.asarray(windowTimes, dtype=float).reshape(-1) + if len(selector) != window_arr.size - 1: + raise ValueError("Window Times must be 1 row longer than selectorArray") + for idx, train in enumerate(selected_trains): + local_min = float(window_arr[idx]) + delta_tw = float(window_arr[idx + 1] - local_min) + if np.asarray(train.spikeTimes).size: + spike_times.extend((np.asarray(train.spikeTimes, dtype=float).reshape(-1) * delta_tw + local_min).tolist()) + + collapsed = nspikeTrain(spike_times, name, delta, minTime, float(maxTime) * len(selector), "time", "s", "", "", -1) + collapsed.setName(name) + collapsed.setMinTime(float(minTime)) + collapsed.setMaxTime(float(maxTime) * len(selector)) + collapsed.resample(1.0 / max(delta, 1e-12)) + return collapsed + def setMinTime(self, value: float | None = None) -> None: if value is None: value = self.minTime @@ -858,12 +919,13 @@ def toStructure(self) -> dict[str, Any]: @staticmethod def fromStructure(structure: dict[str, Any]) -> "TrialConfig": + # MATLAB's `TrialConfig.fromStructure` omits `ensCovMask` and shifts + # the remaining trailing arguments left by one position. return TrialConfig( structure.get("covMask"), structure.get("sampleRate"), structure.get("history"), structure.get("ensCovHist"), - structure.get("ensCovMask"), structure.get("covLag"), structure.get("name", ""), ) @@ -946,8 +1008,6 @@ def setConfigNames(self, names, index: Sequence[int] | None = None) -> None: while len(self.configNames) < self.numConfigs: self.configNames.append("") self.configNames[target] = names if names else f"Fit {target + 1}" - if isinstance(self.configArray[target], TrialConfig): - self.configArray[target].setName(self.configNames[target]) return if isinstance(names, Sequence) and not isinstance(names, (str, bytes)): if len(index) != len(names): @@ -982,10 +1042,7 @@ def fromStructure(structure: dict[str, Any]) -> "ConfigCollection": configs.append(TrialConfig.fromStructure(row)) else: configs.append(row) - coll = ConfigCollection(configs) - if "configNames" in structure: - coll.configNames = list(structure["configNames"]) - return coll + return ConfigCollection(configs) class Trial: @@ -1294,6 +1351,11 @@ def getSpikeVector(self, *args, neuron_index: int = 1) -> np.ndarray: if not args: return self.nspikeColl.dataToMatrix() first = args[0] + if isinstance(first, (int, np.integer)): + selector = [int(first)] + if len(args) == 1: + return self.nspikeColl.dataToMatrix(selector) + return self.nspikeColl.dataToMatrix(selector, *args[1:]) if isinstance(first, Sequence) and not isinstance(first, (str, bytes, np.ndarray)): bin_edges = np.asarray(first, dtype=float).reshape(-1) return self.nspikeColl.getNST(neuron_index).to_binned_counts(bin_edges) diff --git a/parity/README.md b/parity/README.md index e6b6417f..a4fb6dea 100644 --- a/parity/README.md +++ b/parity/README.md @@ -14,16 +14,16 @@ Generated report: python tools/parity/build_report.py ``` -Refresh MATLAB-derived fixtures: +Refresh MATLAB-derived fixtures from the sibling MATLAB repo: ```bash -python tools/parity/export_matlab_gold_fixtures.py +matlab -batch "cd('../nSTAT'); addpath(fullfile(pwd,'tools','python')); export_python_port_fixtures; exit" ``` -Run the combined Python plus MATLAB release gate: +Run the pure-Python release gate: ```bash -nstat-release-check --matlab-repo ../nSTAT +python tools/release/run_fidelity_gate.py ``` Run the MATLAB-side `pyenv` fidelity suite from the sibling MATLAB repo: diff --git a/parity/class_fidelity.yml b/parity/class_fidelity.yml index 6c55c44d..c18368e6 100644 --- a/parity/class_fidelity.yml +++ b/parity/class_fidelity.yml @@ -128,8 +128,8 @@ items: minTime, maxTime, sampleRate, neuronMask, and neighbors. method_parity: MATLAB-facing collection methods are now first-class, including addToColl, addSingleSpikeToColl, merge, getNST, name/index lookup, masking, neighborhood - management, dataToMatrix, ensemble-covariate helpers, restoreToOriginal, psth, - and psthGLM. + management, dataToMatrix, fixture-backed toSpikeTrain collapsing, ensemble-covariate + helpers, restoreToOriginal, psth, and psthGLM. defaults_parity: Defaults for masks, sample rate, and min/max time now track MATLAB collection semantics closely. indexing_parity: MATLAB-facing one-based getNST is preserved. @@ -141,9 +141,9 @@ items: - Some plotting/statistics helpers and lower-level utility methods from MATLAB are still absent. required_remediation: - - Extend the committed MATLAB-derived fixtures beyond collection naming and - `dataToMatrix` outputs to cover neighbor masks, ensemble covariates, and PSTH - outputs. + - Extend the committed MATLAB-derived fixtures beyond collection naming, + `dataToMatrix`, and `toSpikeTrain` outputs to cover neighbor masks, ensemble + covariates, and PSTH outputs. - Port any remaining collection utilities that surface in MATLAB helpfiles. plotting_report_parity: Raster and PSTH plotting works for core workflows; some collection summary visuals remain unported. @@ -192,8 +192,9 @@ items: handling. property_parity: Core configuration fields and normalized metadata are now exposed in the canonical implementation rather than a dataclass shim. - method_parity: MATLAB-facing methods now include naming, structure round-trip, and - setConfig application against Trial state. + method_parity: MATLAB-facing methods now include naming, fixture-backed structure + round-trip, and setConfig application against Trial state, including the legacy + MATLAB fromStructure argument-shift quirk. defaults_parity: Defaults for empty masks/configs and name handling are close to MATLAB. indexing_parity: N/A for this class. @@ -202,10 +203,11 @@ items: output_type_parity: Returns and mutates canonical TrialConfig/Trial objects as expected. symbol_presence_verified: yes known_remaining_differences: - - Some MATLAB normalization and validation branches remain looser in Python. + - Some malformed-configuration normalization and validation branches remain looser + in Python. required_remediation: - - Add malformed-config fixtures from MATLAB to tighten validation and default coercion - behavior. + - Add malformed-config fixtures from MATLAB before promoting TrialConfig from fixture-backed + high_fidelity to exact. plotting_report_parity: N/A - matlab_name: ConfigColl kind: class @@ -213,13 +215,15 @@ items: python_public_name: nstat.ConfigColl python_impl_path: nstat/trial.py status: high_fidelity - constructor_parity: Supports MATLAB-style collections of TrialConfig objects, string-named - configs, and empty config placeholders. + constructor_parity: Fixture-backed canonical behavior now matches MATLAB for collections + of TrialConfig objects, while Python still preserves extra string/empty convenience + paths beyond the MATLAB round-trip surface. property_parity: numConfigs, configNames, and configArray are exposed with MATLAB-style semantics. method_parity: addConfig, getConfig, setConfig, getConfigNames, setConfigNames, - getSubsetConfigs, and structure round-trip are now implemented in the canonical - collection. + getSubsetConfigs, and the TrialConfig-only structure round-trip now follow the + MATLAB collection behavior, including rebuilt Fit N names after the legacy TrialConfig + round-trip bug. defaults_parity: Empty-config and naming defaults now align closely with MATLAB behavior. indexing_parity: One-based getConfig behavior is preserved. @@ -228,9 +232,11 @@ items: output_type_parity: Returns TrialConfig instances. symbol_presence_verified: yes known_remaining_differences: - - Some MATLAB-specific collection manipulation helpers remain unported. + - Python still preserves extra string/empty convenience branches that are not part + of the canonical MATLAB TrialConfig collection round-trip. required_remediation: - - Add fixture-backed tests for edge-case config coercion and selection semantics. + - Decide whether to retain or remove the extra Python convenience branches before + labeling ConfigColl exact. plotting_report_parity: N/A - matlab_name: Analysis kind: class @@ -248,8 +254,8 @@ items: surface for GLMFit, KS/residual/inverse-Gaussian plotting, history-lag search, ensemble-history coefficients, neighbor analysis, Granger-style comparisons, and spike-triggered averaging. - defaults_parity: Default fitting behavior and Poisson-GLM selection are much closer - to the MATLAB workflow defaults. + defaults_parity: Default fitting behavior and the stored Poisson-GLM AIC/BIC/logLL + convention are now fixture-backed on the canonical single-neuron workflow. indexing_parity: MATLAB-facing one-based neuron numbering remains available through the public entry points. error_warning_parity: Core validation is present, though algorithm-selection and @@ -260,9 +266,13 @@ 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. + - The canonical single-neuron GLM path is now fixture-backed for coefficients, + lambda traces, AIC, BIC, stored logLL, KS statistic, residuals, and the discrete-time + KS arrays under injected within-bin draws. Remaining gaps are now concentrated + in broader algorithm-selection, validation-window, and multi-neuron branches + rather than the canonical baseline diagnostics, and the helper surface now also + accepts MATLAB-style multi-trial spike inputs by collapsing them through fixture-backed + `nstColl.toSpikeTrain` semantics. required_remediation: - Extend the committed MATLAB-derived fixture coverage beyond the canonical single-neuron GLM workflow to multi-neuron, validation-window, and alternative algorithm branches. @@ -282,9 +292,9 @@ items: aliases, config metadata, coefficient arrays, history metadata, AIC/BIC/logLL, validation placeholders, and plotParams scaffolding. method_parity: getCoeffs/getHistCoeffs, subset/merge helpers, label remapping, plot-parameter - computation, validation surface, parameter lookup, KS/inverse-Gaussian diagnostics, - residual computation, coefficient/history plotting, and structure round-trip now - operate on the richer MATLAB-style result surface. + computation, validation surface, parameter lookup, fixture-backed KS/residual + diagnostics, coefficient/history plotting, and structure round-trip now operate + on the richer MATLAB-style result surface. defaults_parity: Default result metadata and placeholder fields are much closer to MATLAB than the earlier lightweight container. indexing_parity: N/A for this class. @@ -294,11 +304,14 @@ items: and list/array fields. symbol_presence_verified: yes known_remaining_differences: - - Plotting/report methods now execute, but their numerical detail and layout remain - lighter than MATLAB. + - Plotting/report methods now execute, Z/U/X semantics now follow MATLAB more closely, + and the canonical baseline fit is fixture-backed for AIC/BIC/logLL, KS statistic, + residual traces, deterministic discrete-time KS arrays, and the stored MATLAB-style + KS p-value. Remaining differences are concentrated in richer report layouts, + validation payloads, and multi-fit branches. required_remediation: - - Add MATLAB-derived golden fixtures for coefficient metadata and validation/report - payloads. + - Add MATLAB-derived golden fixtures for validation/report payloads and the remaining + multi-fit branches. - Tighten report-layout and validation rendering against MATLAB screenshots/fixtures. plotting_report_parity: Result plotting/report methods now exist on the canonical object and cover the MATLAB-facing workflow surface, though visual detail still @@ -310,23 +323,26 @@ items: python_impl_path: nstat/fit.py status: high_fidelity constructor_parity: Summary objects now aggregate MATLAB-style FitResult collections - directly. + directly, including MATLAB-style matrix-valued summary fields. property_parity: Core summary fields exist, including fitResCell, numNeurons, numResults, - fitNames, neuronNumbers, AIC, BIC, logLL, and KSStats. + fitNames, neuronNumbers, dev, AIC, BIC, logLL, KSStats, KSPvalues, + and withinConfInt as MATLAB-style neuron-by-fit matrices. method_parity: MATLAB-style difference helpers, coefficient aggregation, significance - summaries, IC plots, residual summary, box-plot surface, summary structure round-trip, - and plotSummary now operate on canonical FitResult collections. + summaries, IC plots, residual summary, box-plot surface, summary + structure round-trip, and plotSummary now operate on canonical FitResult + collections, and the multi-neuron matrix/diff semantics are fixture-backed. defaults_parity: Summary initialization is close for the implemented metadata surface. indexing_parity: N/A for this class. error_warning_parity: Still lighter than MATLAB for mismatched summary inputs. output_type_parity: Returns canonical FitResSummary/FitSummary objects. symbol_presence_verified: yes known_remaining_differences: - - Summary plotting now exists, but richer MATLAB report/table exports remain visually - lighter than MATLAB. + - Summary plotting now exists and the neuron-by-fit AIC/BIC/logLL and diff + aggregation are fixture-backed, but richer MATLAB report/table exports + remain visually lighter than MATLAB. required_remediation: - - Extend the committed golden fixtures beyond single-fit AIC/BIC aggregation to - multi-neuron summary aggregation and remaining report outputs. + - Extend the committed golden fixtures beyond matrix/diff aggregation to + the remaining MATLAB report/table exports and coefficient-view layouts. plotting_report_parity: Summary plotting and report aggregation now cover the MATLAB-facing workflow surface, though fixture-backed visual parity is still pending. - matlab_name: CIF @@ -354,14 +370,17 @@ items: symbol_presence_verified: yes known_remaining_differences: - Analytic and nonlinear polynomial CIF surfaces are now fixture-backed against - MATLAB, but recursive Simulink-backed stochastic trajectories are still validated - as high-fidelity summaries rather than exact sample-by-sample reproductions. + MATLAB, deterministic recursive point-process traces are now fixture-backed + for history, eta, lambda-delta, and spike-indicator evolution under injected + uniform draws, and the continuous-time `simulateCIFByThinningFromLambda` path + is now fixture-backed for proposal generation, thinning ratios, and rounded + accepted spike times. Seeded Simulink-backed stochastic trajectories still + remain high-fidelity rather than exact sample-by-sample reproductions. required_remediation: - - Extend the committed MATLAB-derived fixtures beyond analytic lambda/gradient/Jacobian - outputs, nonlinear polynomial surfaces, and the deterministic recursive lambda - prefix to cover additional thinning and seeded simulation summaries. - - Add MATLAB/Simulink comparison fixtures for recursive CIF simulation trajectories - when the random-stream alignment question is resolved. + - Extend the committed MATLAB-derived fixtures beyond deterministic recursive + traces and thinning-from-lambda proposals to additional seeded simulation summaries. + - Add MATLAB/Simulink comparison fixtures for the remaining seeded recursive + trajectories when the random-stream alignment question is resolved. plotting_report_parity: Simulation/report plotting is limited; downstream notebooks generate figures with helper code rather than a full MATLAB-equivalent CIF report API. @@ -392,14 +411,21 @@ items: tensors instead of only Python-specific dictionaries. symbol_presence_verified: yes known_remaining_differences: + - The linear `PP_fixedIntervalSmoother` path is now fixture-backed against + MATLAB on the canonical single-state example, including the MATLAB-specific + `x0`/`Pi0` initialization semantics for the first lagged update. + - The linear `PPHybridFilterLinear` path now follows MATLAB's per-step model + mixing recursion and is fixture-backed on the canonical 2-model, 1-state + example, including mixed-state bank outputs and posterior model probabilities. - The nonlinear `PPDecodeFilter` path is now fixture-backed against MATLAB on a deterministic polynomial-CIF example, but it still shows small symbolic/numeric drift at the `1e-4` level and remains high-fidelity rather than exact. - Target-estimation augmentation, EM routines, and some advanced symbolic-CIF workflows remain thinner than MATLAB. required_remediation: - - Extend the committed MATLAB-derived numerical fixtures from `PPDecode_predict` - and the deterministic nonlinear `PPDecodeFilter` case to DecodingExample, + - Extend the committed MATLAB-derived numerical fixtures from `PPDecode_predict`, + `PP_fixedIntervalSmoother`, `PPHybridFilterLinear`, and the deterministic nonlinear + `PPDecodeFilter` case to DecodingExample, DecodingExampleWithHist, and HybridFilterExample summaries. - Port the remaining target-estimation, EM, and symbolic-CIF branches from the MATLAB toolbox. @@ -463,12 +489,16 @@ items: python_public_name: nstat.ConfidenceInterval python_impl_path: nstat/confidence_interval.py status: high_fidelity - constructor_parity: Basic time-and-bounds construction aligns with MATLAB intent. - property_parity: lower and upper accessors plus color metadata are exposed. - method_parity: Color assignment, plotting, and arithmetic composition with scalar - signals and other confidence intervals are implemented for the MATLAB-facing workflows - used by Covariate. - defaults_parity: Default color and time/bounds normalization are close to MATLAB. + constructor_parity: Fixture-backed time-and-bounds construction, metadata defaults, + and SignalObj-style serialization now follow MATLAB much more closely. + property_parity: lower and upper accessors plus color/value metadata and SignalObj-style + structure fields are exposed. + method_parity: Color assignment, SignalObj-style dataToStructure/fromStructure + behavior, plotting, and arithmetic composition with scalar signals and other + confidence intervals are implemented for the MATLAB-facing workflows used by + Covariate. + defaults_parity: Default color/value behavior and the MATLAB fromStructure reset + to blue 95% CI are now fixture-backed. indexing_parity: Bounds are stored in MATLAB-style n x 2 lower/upper form. error_warning_parity: Core validation is present, though some MATLAB display/plotting edge cases remain lighter. @@ -476,12 +506,13 @@ items: the expected workflow positions. symbol_presence_verified: yes known_remaining_differences: - - Full MATLAB serialization/display semantics remain lighter than the original toolbox. + - Full MATLAB display/plot styling semantics are still lighter than the original + toolbox. required_remediation: - - Add MATLAB-derived fixtures for serialized confidence-interval payloads and plot - styling. - plotting_report_parity: Core CI plotting works; full MATLAB display/serialization - styling remains lighter. + - Add MATLAB-derived fixtures for exact plot styling before promoting ConfidenceInterval + from fixture-backed high_fidelity to exact. + plotting_report_parity: Core CI plotting works, including MATLAB's string-color + quirk in line mode; full display/styling parity remains lighter. - matlab_name: CovColl kind: class matlab_path: CovColl.m diff --git a/parity/report.md b/parity/report.md index 4ab42bd4..dfb68f22 100644 --- a/parity/report.md +++ b/parity/report.md @@ -60,11 +60,12 @@ Generated from `parity/manifest.yml`, `parity/class_fidelity.yml`, `tools/notebo - 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. +- Clean-room boundary: the installed Python package, pure-Python tests, notebooks, examples, and default CI no longer invoke MATLAB; cross-language execution is confined to the MATLAB-side `tests/python_port_fidelity` harness. - Notebook fidelity: all tracked MATLAB-helpfile notebook ports are marked high fidelity or exact. - 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. - Runtime symbol verification: every audited MATLAB-facing Python symbol marked present in `parity/class_fidelity.yml` resolves on the live public surface. -- Simulink fidelity: native Python coverage exists for the required published workflows, and 10 inventoried MATLAB assets remain reference-only. +- Simulink fidelity: native Python coverage exists for the required published workflows, deterministic injected-input fixtures now back the required native paths, and 10 inventoried MATLAB assets remain reference-only. ## Remaining Mapping Deltas @@ -90,7 +91,7 @@ No audit/runtime symbol mismatches were detected. - `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. +- `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-side fixture export. - `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. diff --git a/parity/simulink_fidelity.yml b/parity/simulink_fidelity.yml index 52a5b7c0..490557e2 100644 --- a/parity/simulink_fidelity.yml +++ b/parity/simulink_fidelity.yml @@ -17,13 +17,15 @@ items: matlab_usage: required_for_behavioral_parity python_strategy: native_python current_python_status: high_fidelity - chosen_interoperability_strategy: Native Python simulation through `nstat.cif.CIF.simulateCIF`, with optional MATLAB Engine reference execution through `nstat.matlab_reference.run_point_process_reference` when MATLAB is available. + deterministic_fixture_backed: yes + stochastic_tolerance_backed: yes + chosen_interoperability_strategy: Native Python simulation through `nstat.cif.CIF.simulateCIF`, with MATLAB-side `tests/python_port_fidelity/TestPythonSimulinkParity.m` serving as the cross-language reference checker. fidelity_risks: - Exact stochastic spike-count realizations differ between MATLAB and NumPy random generators even when the same seed is requested. - The native Python path mirrors the Simulink transfer-function semantics for the published help workflows, but not every internal Simulink block configuration has a one-to-one Python analogue. validation_plan: - - Compare deterministic lambda traces against MATLAB Engine reference runs when MATLAB is available. - - Keep committed MATLAB gold fixtures for the leading lambda trace, plus seeded Python regression tests for FIR filtering, recursive history terms, and PPSimExample outputs in CI. + - Compare deterministic lambda traces against MATLAB-generated fixtures and the MATLAB-side Simulink parity harness. + - Keep committed MATLAB gold fixtures for deterministic injected-uniform traces covering recursive history, eta, lambda-delta, and spike indicators, plus seeded Python regression tests for PPSimExample outputs in CI. - model_name: PointProcessSimulationCont model_path: PointProcessSimulationCont.slx purpose: Continuous-time companion model kept with the MATLAB toolbox for simulation/reference work. @@ -96,7 +98,7 @@ items: 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. + chosen_interoperability_strategy: Keep as a MATLAB reference asset while the Python port validates thinning behavior through `CIF.simulateCIFByThinning` and MATLAB-side fixture export. fidelity_risks: - The legacy thinning model may encode branch behavior that is not separately exercised by current Python CI. validation_plan: @@ -129,14 +131,16 @@ items: matlab_usage: required_for_example_execution python_strategy: native_python current_python_status: high_fidelity - chosen_interoperability_strategy: Native Python execution through `nstat.simulators.simulate_two_neuron_network`, with optional MATLAB Engine reference execution through `nstat.matlab_reference.run_simulated_network_reference` when MATLAB is available. + deterministic_fixture_backed: yes + stochastic_tolerance_backed: yes + chosen_interoperability_strategy: Native Python execution through `nstat.simulators.simulate_two_neuron_network`, with MATLAB-side `tests/python_port_fidelity/TestPythonSimulinkParity.m` serving as the reference checker. fidelity_risks: - Exact spike trains still differ from Simulink because MATLAB and NumPy do not share the same binomial random stream. - The native port mirrors the published NetworkTutorial parameterization and one-sample-delay semantics, but not every internal Simulink block detail is separately exposed. validation_plan: - - Keep deterministic regression tests for the native simulator parameters, probability traces, binary state traces, and estimated network layout in CI. + - Keep deterministic injected-uniform fixtures for probability traces, binary state traces, history/ensemble terms, and eta traces in CI. - Keep committed MATLAB gold fixtures for `prob_head` and `state_head`, and treat seeded spike-count summaries as tolerance-based because MATLAB and NumPy random streams are not identical. - - Run optional MATLAB Engine smoke comparisons for the actual connectivity layout when MATLAB is available. + - Keep MATLAB-side parity checks for the actual connectivity layout and deterministic probability/state traces. - model_name: SimulatedNetwork2Cache model_path: helpfiles/SimulatedNetwork2.slxc purpose: Simulink compiled cache artifact for `SimulatedNetwork2`. diff --git a/pyproject.toml b/pyproject.toml index cc3cf9f8..9f042f18 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,4 +39,3 @@ include = ["nstat*"] [project.scripts] nstat-paper-examples = "nstat.paper_examples:main" nstat-install = "nstat.install:main" -nstat-release-check = "nstat.release_check:main" diff --git a/tests/parity/fixtures/matlab_gold/analysis_exactness.mat b/tests/parity/fixtures/matlab_gold/analysis_exactness.mat index effc52b0..71a9c25f 100644 Binary files a/tests/parity/fixtures/matlab_gold/analysis_exactness.mat and b/tests/parity/fixtures/matlab_gold/analysis_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/cif_exactness.mat b/tests/parity/fixtures/matlab_gold/cif_exactness.mat index e3fe01b8..b4473c51 100644 Binary files a/tests/parity/fixtures/matlab_gold/cif_exactness.mat and b/tests/parity/fixtures/matlab_gold/cif_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/confidence_interval_exactness.mat b/tests/parity/fixtures/matlab_gold/confidence_interval_exactness.mat new file mode 100644 index 00000000..0fb75820 Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/confidence_interval_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/config_exactness.mat b/tests/parity/fixtures/matlab_gold/config_exactness.mat new file mode 100644 index 00000000..a1a9404d Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/config_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/covariate_exactness.mat b/tests/parity/fixtures/matlab_gold/covariate_exactness.mat index 45b3a6b6..39b03e28 100644 Binary files a/tests/parity/fixtures/matlab_gold/covariate_exactness.mat and b/tests/parity/fixtures/matlab_gold/covariate_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/decoding_predict_exactness.mat b/tests/parity/fixtures/matlab_gold/decoding_predict_exactness.mat index 9b06adf9..f47ee130 100644 Binary files a/tests/parity/fixtures/matlab_gold/decoding_predict_exactness.mat and b/tests/parity/fixtures/matlab_gold/decoding_predict_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/decoding_smoother_exactness.mat b/tests/parity/fixtures/matlab_gold/decoding_smoother_exactness.mat new file mode 100644 index 00000000..667e8252 Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/decoding_smoother_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/fit_summary_exactness.mat b/tests/parity/fixtures/matlab_gold/fit_summary_exactness.mat new file mode 100644 index 00000000..dfa93e19 Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/fit_summary_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/hybrid_filter_exactness.mat b/tests/parity/fixtures/matlab_gold/hybrid_filter_exactness.mat new file mode 100644 index 00000000..254474ff Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/hybrid_filter_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/ksdiscrete_exactness.mat b/tests/parity/fixtures/matlab_gold/ksdiscrete_exactness.mat new file mode 100644 index 00000000..1425f67f Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/ksdiscrete_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/nonlinear_decode_exactness.mat b/tests/parity/fixtures/matlab_gold/nonlinear_decode_exactness.mat index 8ca489eb..7b3c4a3d 100644 Binary files a/tests/parity/fixtures/matlab_gold/nonlinear_decode_exactness.mat and b/tests/parity/fixtures/matlab_gold/nonlinear_decode_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat b/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat index 8c5ae069..968146c5 100644 Binary files a/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat and b/tests/parity/fixtures/matlab_gold/nspiketrain_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/nstcoll_exactness.mat b/tests/parity/fixtures/matlab_gold/nstcoll_exactness.mat index d3c89e49..974f7f7d 100644 Binary files a/tests/parity/fixtures/matlab_gold/nstcoll_exactness.mat and b/tests/parity/fixtures/matlab_gold/nstcoll_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/point_process_exactness.mat b/tests/parity/fixtures/matlab_gold/point_process_exactness.mat index 181a2466..81ae3662 100644 Binary files a/tests/parity/fixtures/matlab_gold/point_process_exactness.mat and b/tests/parity/fixtures/matlab_gold/point_process_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat b/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat index b3128bbb..1720d972 100644 Binary files a/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat and b/tests/parity/fixtures/matlab_gold/signalobj_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/simulated_network_exactness.mat b/tests/parity/fixtures/matlab_gold/simulated_network_exactness.mat index 3cc2ec3e..6d30202f 100644 Binary files a/tests/parity/fixtures/matlab_gold/simulated_network_exactness.mat and b/tests/parity/fixtures/matlab_gold/simulated_network_exactness.mat differ diff --git a/tests/parity/fixtures/matlab_gold/thinning_exactness.mat b/tests/parity/fixtures/matlab_gold/thinning_exactness.mat new file mode 100644 index 00000000..8826856a Binary files /dev/null and b/tests/parity/fixtures/matlab_gold/thinning_exactness.mat differ diff --git a/tests/test_analysis_pipeline.py b/tests/test_analysis_pipeline.py index 1da1fb7f..ca65709b 100644 --- a/tests/test_analysis_pipeline.py +++ b/tests/test_analysis_pipeline.py @@ -3,6 +3,8 @@ import numpy as np from nstat import Analysis, CIFModel, ConfigCollection, Covariate, CovariateCollection, FitSummary, Trial, TrialConfig +from nstat.nspikeTrain import nspikeTrain +from nstat.nstColl import nstColl def test_trial_analysis_pipeline() -> None: @@ -22,4 +24,25 @@ def test_trial_analysis_pipeline() -> None: summary = FitSummary(fits) assert summary.numNeurons == 3 - assert summary.AIC.shape == (1,) + assert summary.AIC.shape == (3, 1) + assert summary.meanAIC.shape == (1,) + + +def test_analysis_helpers_accept_multi_trial_spike_inputs_like_matlab() -> None: + time = np.arange(0.0, 1.1, 0.1) + lam = Covariate(time, np.full(time.shape, 2.0), "\\lambda(t)", "time", "s", "Hz", ["lambda"]) + + st1 = nspikeTrain([0.1, 0.3], "1", 10.0, 0.0, 0.5, "time", "s", "", "", -1) + st2 = nspikeTrain([0.2], "1", 10.0, 0.0, 0.5, "time", "s", "", "", -1) + coll = nstColl([st1, st2]) + collapsed = coll.toSpikeTrain() + + multi_ks = Analysis.computeKSStats([st1, st2], lam, DTCorrection=0) + collapsed_ks = Analysis.computeKSStats(collapsed, lam, DTCorrection=0) + for left, right in zip(multi_ks, collapsed_ks, strict=False): + np.testing.assert_allclose(np.asarray(left, dtype=float), np.asarray(right, dtype=float), rtol=1e-8, atol=1e-10) + + multi_residual = Analysis.computeFitResidual([st1, st2], lam, windowSize=0.1) + collapsed_residual = Analysis.computeFitResidual(collapsed, lam, windowSize=0.1) + np.testing.assert_allclose(multi_residual.time, collapsed_residual.time, rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(multi_residual.data, collapsed_residual.data, rtol=1e-8, atol=1e-10) diff --git a/tests/test_cleanroom_boundary.py b/tests/test_cleanroom_boundary.py new file mode 100644 index 00000000..d31487df --- /dev/null +++ b/tests/test_cleanroom_boundary.py @@ -0,0 +1,46 @@ +from __future__ import annotations + +import re +from pathlib import Path + +import nstat + + +REPO_ROOT = Path(__file__).resolve().parents[1] + +FORBIDDEN_RUNTIME_PATTERNS = [ + re.compile(r"\bimport\s+matlab\b"), + re.compile(r"\bimport\s+matlab\.engine\b"), + re.compile(r"\bmatlab\.engine\b"), + re.compile(r"\bstart_matlab\b"), + re.compile(r"\bsubprocess\.(run|Popen)\([^\\n]*matlab"), + re.compile(r"\bshutil\.which\(['\"]matlab['\"]\)"), +] + + +def _assert_clean(paths: list[Path]) -> None: + violations: list[str] = [] + for path in paths: + text = path.read_text(encoding="utf-8", errors="ignore") + for pattern in FORBIDDEN_RUNTIME_PATTERNS: + if pattern.search(text): + violations.append(f"{path.relative_to(REPO_ROOT)} :: {pattern.pattern}") + assert not violations, "Disallowed MATLAB runtime references found:\n" + "\n".join(violations) + + +def test_installable_package_has_no_matlab_runtime_dependency() -> None: + package_paths = sorted((REPO_ROOT / "nstat").glob("**/*.py")) + _assert_clean(package_paths) + + +def test_notebooks_examples_and_ci_do_not_shell_out_to_matlab() -> None: + targets = sorted((REPO_ROOT / "examples").glob("**/*.py")) + targets += sorted((REPO_ROOT / "tools").glob("**/*.py")) + targets += sorted((REPO_ROOT / "notebooks").glob("**/*.ipynb")) + targets += sorted((REPO_ROOT / ".github" / "workflows").glob("**/*.yml")) + _assert_clean(targets) + + +def test_package_root_does_not_export_matlab_reference_helpers() -> None: + for name in ("matlab_engine_available", "run_point_process_reference", "run_simulated_network_reference"): + assert not hasattr(nstat, name), f"nstat should not export MATLAB reference helper {name}" diff --git a/tests/test_fitresult_diagnostics.py b/tests/test_fitresult_diagnostics.py index a44220f9..ad2e2859 100644 --- a/tests/test_fitresult_diagnostics.py +++ b/tests/test_fitresult_diagnostics.py @@ -28,8 +28,10 @@ def test_fitresult_diagnostics_populate_ks_and_residual_fields() -> None: assert "ks_stat" in ks assert fit.KSStats.shape == (1, 1) - assert residual.name == "fit residual" + assert residual.name == "M(t_k)" assert np.asarray(inv, dtype=float).ndim == 1 + assert np.all(np.isfinite(np.asarray(inv, dtype=float))) + np.testing.assert_allclose(np.asarray(inv, dtype=float), np.asarray(fit.X, dtype=float)) assert fit.Residual is not None @@ -66,7 +68,8 @@ def test_fitresult_matlab_style_helpers_expose_plot_params_and_subsets() -> None assert len(labels) == coeffs.size == se.size assert param_vals.size == param_se.size == param_sig.size == 1 assert subset.numResults == 1 - assert fit.KSStats[0, 0] == 0.5 + assert np.isfinite(fit.KSStats[0, 0]) + assert np.isfinite(fit.KSPvalues[0]) assert fit.invGausStats["X"].shape == (1,) assert fit.Residual == {"value": 1} @@ -113,6 +116,8 @@ def test_fitsummary_matlab_style_helpers_cover_ic_and_coeff_views() -> None: restored = FitSummary.fromStructure(summary.toStructure()) assert len(fig1.axes) == 3 + assert summary.AIC.shape == (1, 1) + assert summary.meanAIC.shape == (1,) assert hasattr(ax1, "boxplot") assert hasattr(ax2, "boxplot") assert hasattr(ax3, "boxplot") diff --git a/tests/test_matlab_gold_fixtures.py b/tests/test_matlab_gold_fixtures.py index 5b7f7dff..8fa8ea2a 100644 --- a/tests/test_matlab_gold_fixtures.py +++ b/tests/test_matlab_gold_fixtures.py @@ -13,6 +13,7 @@ CovColl, Covariate, DecodingAlgorithms, + FitResult, FitResSummary, SignalObj, Trial, @@ -46,11 +47,25 @@ def _string(payload: dict[str, np.ndarray], key: str) -> str: if isinstance(value, str): return value arr = np.asarray(value) + if arr.size == 0: + return "" if arr.shape == (): return str(arr.item()) return str(arr.reshape(-1)[0]) +def _string_list(payload: dict[str, np.ndarray], key: str) -> list[str]: + value = payload[key] + if isinstance(value, list): + return [str(item) for item in value] + if isinstance(value, tuple): + return [str(item) for item in value] + arr = np.asarray(value, dtype=object) + if arr.shape == (): + return [str(arr.item())] + return [str(item) for item in arr.reshape(-1)] + + 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"]) @@ -135,15 +150,76 @@ def test_covariate_and_confidence_interval_match_matlab_gold_fixture() -> None: np.testing.assert_allclose(cov_single.ci[0].bounds, np.asarray(payload["explicit_ci"], dtype=float), rtol=1e-8, atol=1e-10) +def test_confidence_interval_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("confidence_interval_exactness.mat") + ci = ConfidenceInterval( + _vector(payload, "time"), + np.asarray(payload["bounds"], dtype=float), + "CI", + "time", + "s", + "a.u.", + ["lo", "hi"], + ["-.k"], + ) + ci.setColor(_string(payload, "color")) + ci.setValue(_scalar(payload, "value")) + structure = ci.toStructure() + roundtrip = ConfidenceInterval.fromStructure(structure) + + np.testing.assert_allclose(ci.dataToMatrix(), np.asarray(payload["bounds"], dtype=float), rtol=1e-12, atol=1e-12) + assert ci.name == _string(payload, "name") + assert ci.color == _string(payload, "color") + assert ci.plotProps == _string_list(payload, "plotProps") + np.testing.assert_allclose(ci.value, _scalar(payload, "value"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(np.asarray(structure["time"], dtype=float), _vector(payload, "structure_time"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(np.asarray(structure["signals"]["values"], dtype=float), np.asarray(payload["structure_values"], dtype=float), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(roundtrip.dataToMatrix(), np.asarray(payload["roundtrip_bounds"], dtype=float), rtol=1e-12, atol=1e-12) + assert roundtrip.color == _string(payload, "roundtrip_color") + np.testing.assert_allclose(roundtrip.value, _scalar(payload, "roundtrip_value"), rtol=1e-12, atol=1e-12) + assert roundtrip.name == _string(payload, "roundtrip_name") + assert roundtrip.plotProps == _string_list(payload, "roundtrip_plotProps") + + def test_nstcoll_matches_matlab_gold_fixture() -> None: payload = _load_fixture("nstcoll_exactness.mat") n1 = nspikeTrain(_vector(payload, "firstSpikeTimes"), "1", 10.0, 0.0, 0.5, "time", "s", "spikes", "spk", -1) n2 = nspikeTrain(_vector(payload, "secondSpikeTimes"), "2", 10.0, 0.0, 0.5, "time", "s", "spikes", "spk", -1) coll = nstColl([n1, n2]) + collapsed = coll.toSpikeTrain() np.testing.assert_equal(coll.numSpikeTrains, int(_scalar(payload, "numSpikeTrains"))) assert coll.getNST(1).name == _string(payload, "firstName") np.testing.assert_allclose(coll.dataToMatrix([1, 2], 0.1, 0.0, 0.5), np.asarray(payload["dataMatrix"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(collapsed.spikeTimes, _vector(payload, "collapsedSpikeTimes"), rtol=1e-8, atol=1e-10) + assert collapsed.name == _string(payload, "collapsedName") + np.testing.assert_allclose(float(collapsed.minTime), _scalar(payload, "collapsedMinTime"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(float(collapsed.maxTime), _scalar(payload, "collapsedMaxTime"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(float(collapsed.sampleRate), _scalar(payload, "collapsedSampleRate"), rtol=1e-12, atol=1e-12) + + +def test_trialconfig_and_configcoll_match_matlab_gold_fixture() -> None: + payload = _load_fixture("config_exactness.mat") + cfg = TrialConfig([["Position", "x"], ["Stimulus"]], 2.0, [0.0, 0.5, 1.0], [], [], 0.5, "stim_pos") + cfg2 = TrialConfig([["Stimulus"]], 2.0, [], [], [], [], "manual") + structure = cfg.toStructure() + roundtrip = TrialConfig.fromStructure(structure) + coll = ConfigColl([cfg, cfg2]) + subset = coll.getSubsetConfigs([1, 2]) + rebuilt = ConfigColl.fromStructure(coll.toStructure()) + + assert cfg.name == _string(payload, "cfg_name") + np.testing.assert_allclose(float(cfg.sampleRate), _scalar(payload, "cfg_sampleRate"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(float(cfg.covLag), _scalar(payload, "cfg_covLag"), rtol=1e-12, atol=1e-12) + assert roundtrip.name == _string(payload, "roundtrip_name") + assert roundtrip.covLag == _string(payload, "roundtrip_covLag") + np.testing.assert_allclose(float(roundtrip.ensCovMask), _scalar(payload, "roundtrip_ensCovMask"), rtol=1e-12, atol=1e-12) + assert coll.getConfigNames() == _string_list(payload, "config_names") + assert subset.getConfigNames() == _string_list(payload, "subset_names") + assert rebuilt.getConfigNames() == _string_list(payload, "rebuilt_names") + assert rebuilt.getConfig(1).name == _string(payload, "rebuilt_first_name") + assert rebuilt.getConfig(1).covLag == _string(payload, "rebuilt_first_covLag") + np.testing.assert_allclose(float(rebuilt.getConfig(1).ensCovMask), _scalar(payload, "rebuilt_first_ensCovMask"), rtol=1e-12, atol=1e-12) def test_cif_eval_surface_matches_matlab_gold_fixture() -> None: @@ -196,13 +272,153 @@ def test_analysis_fit_surface_matches_matlab_gold_fixture() -> None: 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) - np.testing.assert_allclose(float(summary.AIC[0]), _scalar(payload, "summaryAIC"), rtol=1e-8, atol=1e-10) - np.testing.assert_allclose(float(summary.BIC[0]), _scalar(payload, "summaryBIC"), rtol=1e-8, atol=1e-10) - assert np.isfinite(float(fit.logLL[0])) - assert np.isfinite(float(summary.logLL[0])) + np.testing.assert_allclose(float(fit.logLL[0]), _scalar(payload, "logLL"), rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(float(summary.AIC[0, 0]), _scalar(payload, "summaryAIC"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(summary.BIC[0, 0]), _scalar(payload, "summaryBIC"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(summary.logLL[0, 0]), _scalar(payload, "summarylogLL"), rtol=1e-6, atol=1e-8) + ks_stats = fit.computeKSStats(1) + np.testing.assert_allclose(float(ks_stats["ks_stat"]), _scalar(payload, "ks_stat"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(ks_stats["ks_pvalue"]), _scalar(payload, "ks_pvalue"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(ks_stats["within_conf_int"]), _scalar(payload, "ks_within_conf_int"), rtol=1e-8, atol=1e-10) + residual = fit.computeFitResidual(1) + np.testing.assert_allclose(residual.time, _vector(payload, "residual_time"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(residual.data[:, 0], _vector(payload, "residual_data"), rtol=1e-6, atol=1e-8) assert fit.fitType[0] == _string(payload, "distribution") +def test_analysis_discrete_ks_arrays_match_matlab_gold_fixture() -> None: + payload = _load_fixture("ksdiscrete_exactness.mat") + + spike_train = nspikeTrain( + _vector(payload, "spike_times"), + "1", + 0.1, + 0.0, + 1.0, + "time", + "s", + "", + "", + -1, + ) + lambda_signal = Covariate( + _vector(payload, "lambda_time"), + _vector(payload, "lambda_data"), + "\\lambda(t)", + "time", + "s", + "Hz", + ["\\lambda_{1}"], + ) + + Z, U, xAxis, KSSorted, ks_stat = Analysis.computeKSStats( + spike_train, + lambda_signal, + 1, + random_values=_vector(payload, "uniform_draws"), + ) + + np.testing.assert_allclose(np.asarray(Z, dtype=float).reshape(-1), _vector(payload, "Z"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(np.asarray(U, dtype=float).reshape(-1), _vector(payload, "U"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(np.asarray(xAxis, dtype=float).reshape(-1), _vector(payload, "xAxis"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(np.asarray(KSSorted, dtype=float).reshape(-1), _vector(payload, "KSSorted"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(ks_stat), _scalar(payload, "compute_ks_stat"), rtol=1e-8, atol=1e-10) + + stim = Covariate(_vector(payload, "time"), _vector(payload, "stim_data"), "Stimulus", "time", "s", "", ["stim"]) + trial = Trial(nstColl([spike_train]), CovColl([stim])) + cfg = TrialConfig([["Stimulus", "stim"]], _scalar(payload, "sample_rate"), [], [], name="stim") + fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl([cfg]), makePlot=0) + fit.setKSStats(Z, U, xAxis, KSSorted, np.asarray([ks_stat], dtype=float)) + np.testing.assert_allclose(float(fit.KSStats[0, 0]), _scalar(payload, "ks_stat"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(fit.KSPvalues[0]), _scalar(payload, "ks_pvalue"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(float(fit.withinConfInt[0]), _scalar(payload, "ks_within_conf_int"), rtol=1e-8, atol=1e-10) + + +def test_fit_summary_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("fit_summary_exactness.mat") + time = np.arange(0.0, 1.0 + 0.1, 0.1) + lambda_signal = Covariate( + time, + np.column_stack( + [ + np.linspace(2.0, 7.0, time.size, dtype=float), + np.linspace(3.0, 8.0, time.size, dtype=float), + ] + ), + "\\lambda(t)", + "time", + "s", + "Hz", + ["stim", "stim_hist"], + ) + st1 = nspikeTrain([0.1, 0.4, 0.7], "1", 0.1, 0.0, 1.0, "time", "s", "", "", -1) + st2 = nspikeTrain([0.2, 0.5, 0.8], "2", 0.1, 0.0, 1.0, "time", "s", "", "", -1) + stim_cfg = TrialConfig([["Stimulus", "stim"]], 10.0, [], [], name="stim") + stim_hist_cfg = TrialConfig([["Stimulus", "stim"]], 10.0, [0.0, 0.1, 0.2], [], name="stim_hist") + config_coll = ConfigColl([stim_cfg, stim_hist_cfg]) + fit1 = FitResult( + st1, + [["stim"], ["stim", "hist1", "hist2"]], + [0, 2], + [None, None], + [None, None], + lambda_signal, + [np.array([0.5]), np.array([0.3, -0.1, -0.05])], + np.array([1.0, 2.0]), + [ + {"se": np.array([0.05]), "p": np.array([0.01])}, + {"se": np.array([0.04, 0.03, 0.02]), "p": np.array([0.02, 0.04, 0.06])}, + ], + np.array([11.0, 7.0]), + np.array([12.0, 8.0]), + np.array([3.0, 5.0]), + config_coll, + [], + [], + "poisson", + ) + fit2 = FitResult( + st2, + [["stim"], ["stim", "hist1", "hist2"]], + [0, 2], + [None, None], + [None, None], + lambda_signal, + [np.array([0.4]), np.array([0.25, -0.08, -0.02])], + np.array([1.5, 2.5]), + [ + {"se": np.array([0.06]), "p": np.array([0.03])}, + {"se": np.array([0.05, 0.04, 0.03]), "p": np.array([0.01, 0.03, 0.07])}, + ], + np.array([13.0, 9.0]), + np.array([14.0, 10.0]), + np.array([2.0, 4.0]), + config_coll, + [], + [], + "poisson", + ) + fit1.KSStats[:, 0] = np.array([0.25, 0.50], dtype=float) + fit1.KSPvalues[:] = np.array([0.90, 0.40], dtype=float) + fit1.withinConfInt[:] = np.array([1.0, 1.0], dtype=float) + fit2.KSStats[:, 0] = np.array([0.35, 0.55], dtype=float) + fit2.KSPvalues[:] = np.array([0.80, 0.30], dtype=float) + fit2.withinConfInt[:] = np.array([1.0, 0.0], dtype=float) + summary = FitResSummary([fit1, fit2]) + + assert summary.fitNames == _string_list(payload, "fitNames") + np.testing.assert_allclose(np.asarray(summary.neuronNumbers, dtype=float), _vector(payload, "neuronNumbers"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(summary.AIC, np.asarray(payload["AIC"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(summary.BIC, np.asarray(payload["BIC"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(summary.logLL, np.asarray(payload["logLL"], dtype=float), rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(summary.KSStats, np.asarray(payload["KSStats"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(summary.KSPvalues, np.asarray(payload["KSPvalues"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(summary.withinConfInt, np.asarray(payload["withinConfInt"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(summary.getDiffAIC(1), np.asarray(payload["diffAIC"], dtype=float).reshape(summary.getDiffAIC(1).shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(summary.getDiffBIC(1), np.asarray(payload["diffBIC"], dtype=float).reshape(summary.getDiffBIC(1).shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(summary.getDifflogLL(1), np.asarray(payload["difflogLL"], dtype=float).reshape(summary.getDifflogLL(1).shape), rtol=1e-6, atol=1e-8) + + def test_point_process_lambda_trace_matches_matlab_gold_fixture() -> None: payload = _load_fixture("point_process_exactness.mat") @@ -225,6 +441,65 @@ def test_point_process_lambda_trace_matches_matlab_gold_fixture() -> None: np.testing.assert_allclose(lambda_cov.data[: _vector(payload, 'lambda_head').shape[0], 0], _vector(payload, "lambda_head"), rtol=1e-8, atol=1e-10) +def test_point_process_deterministic_trace_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("point_process_exactness.mat") + time = _vector(payload, "det_time") + stim_values = _vector(payload, "det_stimulus") + uniforms = _vector(payload, "det_uniforms").reshape(-1, 1) + stim = Covariate(time, stim_values, "Stimulus", "time", "s", "Voltage", ["sin"]) + ens = Covariate(time, np.zeros_like(time), "Ensemble", "time", "s", "Spikes", ["n1"]) + _, lambda_cov, details = CIF.simulateCIF( + -3.0, + np.array([-1.0, -2.0, -4.0], dtype=float), + np.array([1.0], dtype=float), + np.array([0.0], dtype=float), + stim, + ens, + numRealizations=1, + simType="binomial", + random_values=uniforms, + return_lambda=True, + return_details=True, + ) + + np.testing.assert_allclose(lambda_cov.data[:, 0], _vector(payload, "det_rate_hz"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(details["lambda_delta"][:, 0], _vector(payload, "det_lambda_delta"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(details["history_effect"][:, 0], _vector(payload, "det_history_effect"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(details["eta"][:, 0], _vector(payload, "det_eta"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(details["spike_indicator"][:, 0], _vector(payload, "det_spike_indicator"), rtol=1e-8, atol=1e-10) + + +def test_cif_thinning_from_lambda_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("thinning_exactness.mat") + lambda_cov = Covariate( + _vector(payload, "time"), + _vector(payload, "lambda_data"), + "\\lambda(t)", + "time", + "s", + "Hz", + ["\\lambda"], + ) + spike_coll, details = CIF.simulateCIFByThinningFromLambda( + lambda_cov, + numRealizations=1, + maxTimeRes=_scalar(payload, "maxTimeRes"), + random_values=_vector(payload, "arrival_uniforms"), + thinning_values=_vector(payload, "thinning_uniforms"), + return_details=True, + ) + + np.testing.assert_allclose(float(details["lambda_bound"]), _scalar(payload, "lambda_bound"), rtol=1e-12, atol=1e-12) + assert int(details["proposal_count"]) == int(_scalar(payload, "proposal_count")) + np.testing.assert_allclose(details["arrival_uniforms"], _vector(payload, "arrival_uniforms"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(details["interarrival_times"], _vector(payload, "interarrival_times"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(details["candidate_spike_times"], _vector(payload, "candidate_spike_times"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(details["lambda_ratio"], _vector(payload, "lambda_ratio"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(details["thinning_uniforms"], _vector(payload, "thinning_uniforms"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(details["accepted_spike_times"], _vector(payload, "rounded_spike_times"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(spike_coll.getNST(1).spikeTimes, _vector(payload, "rounded_spike_times"), rtol=1e-8, atol=1e-10) + + def test_decoding_predict_matches_matlab_gold_fixture() -> None: payload = _load_fixture("decoding_predict_exactness.mat") x_p, W_p = DecodingAlgorithms.PPDecode_predict( @@ -238,6 +513,50 @@ def test_decoding_predict_matches_matlab_gold_fixture() -> None: np.testing.assert_allclose(W_p, np.asarray(payload["W_p"], dtype=float), rtol=1e-8, atol=1e-10) +def test_pp_fixed_interval_smoother_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("decoding_smoother_exactness.mat") + x_pLag, W_pLag, x_uLag, W_uLag = DecodingAlgorithms.PP_fixedIntervalSmoother( + _scalar(payload, "A"), + _scalar(payload, "Q"), + np.asarray(payload["dN"], dtype=float), + int(_scalar(payload, "lags")), + _scalar(payload, "mu"), + _scalar(payload, "beta"), + _string(payload, "fitType"), + _scalar(payload, "delta"), + ) + + np.testing.assert_allclose(x_pLag, np.asarray(payload["x_pLag"], dtype=float).reshape(x_pLag.shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(W_pLag, np.asarray(payload["W_pLag"], dtype=float).reshape(W_pLag.shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(x_uLag, np.asarray(payload["x_uLag"], dtype=float).reshape(x_uLag.shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(W_uLag, np.asarray(payload["W_uLag"], dtype=float).reshape(W_uLag.shape), rtol=1e-8, atol=1e-10) + + +def test_pp_hybrid_filter_linear_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("hybrid_filter_exactness.mat") + S_est, X, W, MU_u, X_s, W_s, pNGivenS = DecodingAlgorithms.PPHybridFilterLinear( + [np.array([[float(_scalar(payload, "A1"))]], dtype=float), np.array([[float(_scalar(payload, "A2"))]], dtype=float)], + [np.array([[float(_scalar(payload, "Q1"))]], dtype=float), np.array([[float(_scalar(payload, "Q2"))]], dtype=float)], + np.asarray(payload["p_ij"], dtype=float), + _vector(payload, "Mu0"), + np.asarray(payload["dN"], dtype=float), + float(_scalar(payload, "mu")), + float(_scalar(payload, "beta")), + _string(payload, "fitType"), + _scalar(payload, "binwidth"), + ) + + np.testing.assert_allclose(S_est, _vector(payload, "S_est"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(X, np.asarray(payload["X"], dtype=float).reshape(X.shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(W, np.asarray(payload["W"], dtype=float).reshape(W.shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(MU_u, np.asarray(payload["MU_u"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(X_s[0], np.asarray(payload["X_s_1"], dtype=float).reshape(X_s[0].shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(X_s[1], np.asarray(payload["X_s_2"], dtype=float).reshape(X_s[1].shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(W_s[0], np.asarray(payload["W_s_1"], dtype=float).reshape(W_s[0].shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(W_s[1], np.asarray(payload["W_s_2"], dtype=float).reshape(W_s[1].shape), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(pNGivenS, np.asarray(payload["pNGivenS"], dtype=float), rtol=1e-8, atol=1e-10) + + def test_nonlinear_ppdecodefilter_matches_matlab_gold_fixture() -> None: payload = _load_fixture("nonlinear_decode_exactness.mat") lambda_cifs = [ @@ -279,3 +598,21 @@ def test_simulated_network_matches_matlab_gold_fixture() -> None: matlab_counts = _vector(payload, "spike_counts") assert native_counts.shape == matlab_counts.shape assert np.all(np.abs(native_counts - matlab_counts) <= 64.0) + + +def test_simulated_network_deterministic_trace_matches_matlab_gold_fixture() -> None: + payload = _load_fixture("simulated_network_exactness.mat") + sim = simulate_two_neuron_network( + duration_s=float(_vector(payload, "det_time")[-1]), + dt=float(_vector(payload, "det_time")[1] - _vector(payload, "det_time")[0]), + seed=None, + uniform_values=np.asarray(payload["det_uniforms"], dtype=float), + ) + + np.testing.assert_allclose(sim.time, _vector(payload, "det_time"), rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(sim.latent_drive, _vector(payload, "det_stimulus"), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(sim.lambda_delta, np.asarray(payload["det_probability"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(sim.spike_indicator, np.asarray(payload["det_state"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(sim.eta, np.asarray(payload["det_eta"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(sim.history_effect, np.asarray(payload["det_history_effect"], dtype=float), rtol=1e-8, atol=1e-10) + np.testing.assert_allclose(sim.ensemble_effect, np.asarray(payload["det_ensemble_effect"], dtype=float), rtol=1e-8, atol=1e-10) diff --git a/tests/test_matlab_reference.py b/tests/test_matlab_reference.py index 49416780..f58b1fa7 100644 --- a/tests/test_matlab_reference.py +++ b/tests/test_matlab_reference.py @@ -3,46 +3,20 @@ from pathlib import Path import numpy as np -import pytest +from scipy.io import loadmat from nstat import Analysis, CIF, ConfigColl, CovColl, Covariate, FitResSummary, 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, -) REPO_ROOT = Path(__file__).resolve().parents[1] -MATLAB_REPO_ROOT = REPO_ROOT.parent / "nSTAT" +FIXTURE_ROOT = REPO_ROOT / "tests" / "parity" / "fixtures" / "matlab_gold" +def _load_fixture(name: str): + return loadmat(FIXTURE_ROOT / name, squeeze_me=True, struct_as_record=False) -def test_matlab_engine_detection_is_boolean() -> None: - assert isinstance(matlab_engine_available(), bool) - - -@pytest.mark.skipif(not MATLAB_REPO_ROOT.exists(), reason="MATLAB reference repo not available") -def test_matlab_reference_executes_only_when_engine_is_available() -> None: - if not matlab_engine_available(): - pytest.skip("MATLAB Engine for Python is not installed") - - point_process = run_point_process_reference(matlab_repo=MATLAB_REPO_ROOT) - network = run_simulated_network_reference(matlab_repo=MATLAB_REPO_ROOT) - - assert point_process["spike_counts"].shape == (5,) - assert point_process["lambda_head"].shape == (10,) - assert network["spike_counts"].shape == (2,) - assert network["prob_head"].shape == (5, 2) - assert network["state_head"].shape == (5, 2) - np.testing.assert_allclose(network["actual_network"], np.array([[0.0, 1.0], [-4.0, 0.0]], dtype=float)) - - -@pytest.mark.skipif(not MATLAB_REPO_ROOT.exists(), reason="MATLAB reference repo not available") -def test_native_point_process_simulation_matches_matlab_lambda_head_when_engine_is_available() -> None: - if not matlab_engine_available(): - pytest.skip("MATLAB Engine for Python is not installed") - +def test_point_process_fixture_remains_consumable_without_matlab_runtime() -> None: + payload = _load_fixture("point_process_exactness.mat") + lambda_head = np.asarray(payload["lambda_head"], dtype=float).reshape(-1) time = np.arange(0.0, 50.0 + 0.001, 0.001, dtype=float) stim = Covariate(time, np.sin(2 * np.pi * 1.0 * time), "Stimulus", "time", "s", "Voltage", ["sin"]) ens = Covariate(time, np.zeros_like(time), "Ensemble", "time", "s", "Spikes", ["n1"]) @@ -58,36 +32,25 @@ def test_native_point_process_simulation_matches_matlab_lambda_head_when_engine_ seed=5, return_lambda=True, ) - matlab_ref = run_point_process_reference(matlab_repo=MATLAB_REPO_ROOT, seed=5) - - np.testing.assert_allclose(lambda_cov.data[:10, 0], matlab_ref["lambda_head"], rtol=1e-6, atol=1e-8) - - -@pytest.mark.skipif(not MATLAB_REPO_ROOT.exists(), reason="MATLAB reference repo not available") -def test_native_network_simulation_preserves_matlab_connectivity_layout_when_engine_is_available() -> None: - if not matlab_engine_available(): - pytest.skip("MATLAB Engine for Python is not installed") + np.testing.assert_allclose(lambda_cov.data[: lambda_head.shape[0], 0], lambda_head, rtol=1e-8, atol=1e-10) +def test_network_fixture_remains_consumable_without_matlab_runtime() -> None: + payload = _load_fixture("simulated_network_exactness.mat") native = simulate_two_neuron_network(seed=4) - matlab_ref = run_simulated_network_reference(matlab_repo=MATLAB_REPO_ROOT, seed=4) - 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) + np.testing.assert_allclose(native.actual_network, np.asarray(payload["actual_network"], dtype=float)) + np.testing.assert_allclose(native.lambda_delta[:5], np.asarray(payload["prob_head"], dtype=float), rtol=1e-8, atol=1e-10) dt = float(native.time[1] - native.time[0]) native_state_head = np.column_stack([ native.spikes.getNST(1).getSigRep(dt, float(native.time[0]), float(native.time[-1])).data[:5, 0], native.spikes.getNST(2).getSigRep(dt, float(native.time[0]), float(native.time[-1])).data[:5, 0], ]) - np.testing.assert_allclose(native_state_head, matlab_ref["state_head"], rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(native_state_head, np.asarray(payload["state_head"], dtype=float), rtol=1e-8, atol=1e-10) native_counts = np.array([len(native.spikes.getNST(1).spikeTimes), len(native.spikes.getNST(2).spikeTimes)], dtype=float) - assert np.all(np.abs(native_counts - matlab_ref["spike_counts"]) <= 64.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") + assert np.all(np.abs(native_counts - np.asarray(payload["spike_counts"], dtype=float).reshape(-1)) <= 64.0) +def test_analysis_fixture_remains_consumable_without_matlab_runtime() -> None: + payload = _load_fixture("analysis_exactness.mat") 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) @@ -95,12 +58,15 @@ def test_native_analysis_fit_matches_matlab_reference_when_engine_is_available() cfg = TrialConfig([["Stimulus", "stim"]], 10.0, [], [], name="stim") fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl([cfg])) summary = FitResSummary([fit]) - 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) - np.testing.assert_allclose(np.asarray(summary.AIC, dtype=float)[:1], matlab_ref["summary_aic"], rtol=1e-6, atol=1e-8) - np.testing.assert_allclose(np.asarray(summary.BIC, dtype=float)[:1], matlab_ref["summary_bic"], rtol=1e-6, atol=1e-8) + assert np.asarray(payload["coeffs"], dtype=float).reshape(-1).size >= 1 + assert np.asarray(payload["lambda_data"], dtype=float).reshape(-1).size >= 5 + assert np.isfinite(np.asarray(payload["AIC"], dtype=float).reshape(-1)[0]) + assert np.isfinite(np.asarray(payload["BIC"], dtype=float).reshape(-1)[0]) + assert np.isfinite(np.asarray(payload["summaryAIC"], dtype=float).reshape(-1)[0]) + assert np.isfinite(np.asarray(payload["summaryBIC"], dtype=float).reshape(-1)[0]) + assert np.isfinite(np.asarray(fit.AIC, dtype=float)[0]) + assert np.isfinite(np.asarray(fit.BIC, dtype=float)[0]) + assert np.isfinite(np.asarray(summary.AIC, dtype=float).reshape(-1)[0]) + assert np.isfinite(np.asarray(summary.BIC, dtype=float).reshape(-1)[0]) assert np.isfinite(np.asarray(fit.logLL, dtype=float)[0]) diff --git a/tests/test_parity_report.py b/tests/test_parity_report.py index 6400baf5..4a07a569 100644 --- a/tests/test_parity_report.py +++ b/tests/test_parity_report.py @@ -25,6 +25,8 @@ def test_parity_report_highlights_current_constraints() -> None: 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 "Clean-room boundary" in text + assert "cross-language execution is confined to the MATLAB-side `tests/python_port_fidelity` harness" in text assert "No partial notebook-fidelity items remain in `tools/notebooks/parity_notes.yml`." in text assert "No partial or missing items remain in the mapping inventory." in text assert "Remaining Class-Fidelity Deltas" in text diff --git a/tests/test_release_check.py b/tests/test_release_check.py index 7d91b856..c65cef5f 100644 --- a/tests/test_release_check.py +++ b/tests/test_release_check.py @@ -1,24 +1,21 @@ from __future__ import annotations +import sys from pathlib import Path -from nstat.release_check import build_release_gate_commands - REPO_ROOT = Path(__file__).resolve().parents[1] +RELEASE_TOOLS = REPO_ROOT / "tools" / "release" +if str(RELEASE_TOOLS) not in sys.path: + sys.path.insert(0, str(RELEASE_TOOLS)) - -def test_release_gate_includes_fixture_generation_and_matlab_suite() -> None: - commands = build_release_gate_commands(REPO_ROOT, matlab_repo_root=REPO_ROOT.parent / "nSTAT") - flattened = [" ".join(command) for command in commands] - - assert any("export_matlab_gold_fixtures.py" in item for item in flattened) - assert any("tests/python_port_fidelity" in item for item in flattened) - assert any("addpath(fullfile(pwd,'helpfiles'))" in item for item in flattened) +from release_gate_lib import build_release_gate_commands -def test_release_gate_can_skip_matlab() -> None: - commands = build_release_gate_commands(REPO_ROOT, skip_matlab=True) +def test_release_gate_is_pure_python() -> None: + commands = build_release_gate_commands(REPO_ROOT) flattened = [" ".join(command) for command in commands] - assert not any(item.startswith("matlab ") for item in flattened) + assert any("tests/test_cleanroom_boundary.py" in item for item in flattened) + assert not any("matlab " in item.lower() for item in flattened) + assert not any("export_matlab_gold_fixtures.py" in item for item in flattened) diff --git a/tests/test_signalobj_fidelity.py b/tests/test_signalobj_fidelity.py index c59e30f5..d1b7e695 100644 --- a/tests/test_signalobj_fidelity.py +++ b/tests/test_signalobj_fidelity.py @@ -152,3 +152,21 @@ def test_signalobj_math_and_summary_methods_match_matlab_surface() -> None: np.testing.assert_allclose(min_vals, [1.0, 1.0]) np.testing.assert_array_equal(min_idx, [0, 1]) np.testing.assert_allclose(min_time, [0.0, 1.0]) + + +def test_confidence_interval_line_plot_ignores_string_color_like_matlab() -> None: + ci = ConfidenceInterval([0.0, 1.0], [[0.8, 1.2], [1.8, 2.2]], "CI", "time", "s", "a.u.", ["lo", "hi"], ["-.k"]) + + fig1, ax1 = plt.subplots() + lines_default = ci.plot(color="r", drawPatches=0, ax=ax1) + default_colors = [line.get_color() for line in lines_default] + + fig2, ax2 = plt.subplots() + lines_numeric = ci.plot(color=(0.2, 0.4, 0.6), drawPatches=0, ax=ax2) + numeric_colors = [line.get_color() for line in lines_numeric] + + assert default_colors != ["r", "r"] + assert numeric_colors == [(0.2, 0.4, 0.6), (0.2, 0.4, 0.6)] + + plt.close(fig1) + plt.close(fig2) diff --git a/tests/test_simulink_fidelity_audit.py b/tests/test_simulink_fidelity_audit.py index a26d47f6..114e57af 100644 --- a/tests/test_simulink_fidelity_audit.py +++ b/tests/test_simulink_fidelity_audit.py @@ -51,6 +51,19 @@ def test_simulink_fidelity_audit_records_required_execution_fields() -> None: assert row["validation_plan"] +def test_required_native_simulink_paths_record_deterministic_and_stochastic_backing() -> None: + payload = _load_audit() + required = { + row["model_name"]: row + for row in payload["items"] + if row["model_name"] in {"PointProcessSimulation", "SimulatedNetwork2"} + } + assert bool(required["PointProcessSimulation"]["deterministic_fixture_backed"]) is True + assert bool(required["PointProcessSimulation"]["stochastic_tolerance_backed"]) is True + assert bool(required["SimulatedNetwork2"]["deterministic_fixture_backed"]) is True + assert bool(required["SimulatedNetwork2"]["stochastic_tolerance_backed"]) is True + + def test_simulink_fidelity_audit_covers_required_model_inventory() -> None: payload = _load_audit() model_paths = {row["model_path"] for row in payload["items"]} diff --git a/tests/test_trial_fidelity.py b/tests/test_trial_fidelity.py index 4bcf886d..beeaec45 100644 --- a/tests/test_trial_fidelity.py +++ b/tests/test_trial_fidelity.py @@ -83,16 +83,27 @@ def test_trialconfig_and_configcoll_apply_and_roundtrip() -> None: assert trial.getCovLabelsFromMask() == ["x", "stim"] roundtrip = TrialConfig.fromStructure(cfg.toStructure()) - assert roundtrip.name == "stim_pos" + assert roundtrip.name == "" + assert roundtrip.covLag == "stim_pos" + assert roundtrip.ensCovMask == 0.5 assert roundtrip.covariate_names == ["Position", "Stimulus"] - configs = ConfigColl([cfg, "manual", None]) - assert configs.numConfigs == 3 - assert configs.getConfigNames() == ["stim_pos", "manual", "Empty Config"] + cfg2 = TrialConfig( + covMask=[["Stimulus"]], + sampleRate=2.0, + history=[], + covLag=[], + name="manual", + ) + configs = ConfigColl([cfg, cfg2]) + assert configs.numConfigs == 2 + assert configs.getConfigNames() == ["stim_pos", "manual"] subset = configs.getSubsetConfigs([1, 2]) assert subset.numConfigs == 2 rebuilt = ConfigColl.fromStructure(configs.toStructure()) - assert rebuilt.getConfigNames() == ["stim_pos", "manual", "Empty Config"] + assert rebuilt.getConfigNames() == ["Fit 1", "Fit 2"] + assert rebuilt.getConfig(1).name == "" + assert rebuilt.getConfig(1).covLag == "stim_pos" def test_trial_partition_history_design_matrix_and_spike_vector() -> None: @@ -115,6 +126,7 @@ def test_trial_partition_history_design_matrix_and_spike_vector() -> None: assert design.shape[1] == 5 spikes = trial.getSpikeVector() assert spikes.shape[1] == 2 + np.testing.assert_allclose(trial.getSpikeVector(1).reshape(-1), spikes[:, 0]) def test_events_validation_and_history_collection_output() -> None: diff --git a/tests/test_workflow_fidelity.py b/tests/test_workflow_fidelity.py index 05d9fcb7..e4507c6a 100644 --- a/tests/test_workflow_fidelity.py +++ b/tests/test_workflow_fidelity.py @@ -80,7 +80,7 @@ def test_fitresult_roundtrip_and_summary_preserve_core_metadata() -> None: assert summary.numNeurons == 2 assert summary.numResults == 1 assert summary.fitNames == ["stim_only"] - assert summary.AIC.shape == (1,) + assert summary.AIC.shape == (2, 1) def test_cif_instantiation_evaluation_and_simulate_from_lambda() -> None: @@ -95,6 +95,17 @@ def test_cif_instantiation_evaluation_and_simulate_from_lambda() -> None: sim = CIF.simulateCIFByThinningFromLambda(cov, numRealizations=2) assert sim.numSpikeTrains == 2 + sim_single, details = CIF.simulateCIFByThinningFromLambda( + cov.getSubSignal(1), + numRealizations=1, + random_values=np.array([0.8, 0.6, 0.4, 0.2], dtype=float), + thinning_values=np.array([0.1, 0.9, 0.3, 0.7], dtype=float), + return_details=True, + ) + assert sim_single.numSpikeTrains == 1 + assert int(details["proposal_count"]) >= 0 + assert details["candidate_spike_times"].ndim == 1 + model = CIFModel(time, np.array([5.0, 6.0, 7.0]), name="lambda") sim2 = model.simulate(num_realizations=1, seed=1) assert sim2.numSpikeTrains == 1 diff --git a/tools/notebooks/build_foundational_help_notebooks.py b/tools/notebooks/build_foundational_help_notebooks.py index 77f3c316..63be7b66 100644 --- a/tools/notebooks/build_foundational_help_notebooks.py +++ b/tools/notebooks/build_foundational_help_notebooks.py @@ -361,13 +361,13 @@ def _figure(label: str, *, figsize=(8.5, 4.5)): summary = FitResSummary(results) fig = _figure("Summary.plotSummary", figsize=(10.0, 4.5)) summary.plotSummary(handle=fig) - print({"fit_names": summary.fitNames, "mean_AIC": np.asarray(summary.AIC, dtype=float).round(3).tolist()}) + print({"fit_names": summary.fitNames, "mean_AIC": np.asarray(summary.meanAIC, dtype=float).round(3).tolist()}) """, """ # SECTION 17: Summarize model selection fig = _figure("bar(summary.AIC)", figsize=(8.0, 4.5)) ax = fig.subplots(1, 1) - ax.bar(np.arange(len(summary.fitNames)), np.asarray(summary.AIC, dtype=float), color=["0.6", "tab:blue", "tab:green"]) + ax.bar(np.arange(len(summary.fitNames)), np.asarray(summary.meanAIC, dtype=float), color=["0.6", "tab:blue", "tab:green"]) ax.set_xticks(np.arange(len(summary.fitNames)), summary.fitNames, rotation=20) ax.set_ylabel("mean AIC") ax.set_title("Model comparison across realizations") diff --git a/tools/notebooks/build_helpfile_fidelity_notebooks.py b/tools/notebooks/build_helpfile_fidelity_notebooks.py index e5595d10..1a26eb15 100644 --- a/tools/notebooks/build_helpfile_fidelity_notebooks.py +++ b/tools/notebooks/build_helpfile_fidelity_notebooks.py @@ -538,10 +538,10 @@ def _plot_isi_hist(ax, train, lambda_hz, *, title): fig = _prepare_figure("Summary.plotSummary", figsize=(8.5, 4.5)) axs = fig.subplots(1, 2) xloc = np.arange(len(summary.fitNames)) - axs[0].bar(xloc, summary.AIC, color=["tab:blue", "tab:green"]) + axs[0].bar(xloc, summary.meanAIC, color=["tab:blue", "tab:green"]) axs[0].set_xticks(xloc, summary.fitNames) axs[0].set_title("Mean AIC across neurons") - axs[1].bar(xloc, summary.BIC, color=["tab:blue", "tab:green"]) + axs[1].bar(xloc, summary.meanBIC, color=["tab:blue", "tab:green"]) axs[1].set_xticks(xloc, summary.fitNames) axs[1].set_title("Mean BIC across neurons") diff --git a/tools/parity/export_matlab_gold_fixtures.py b/tools/parity/export_matlab_gold_fixtures.py deleted file mode 100644 index 8313279f..00000000 --- a/tools/parity/export_matlab_gold_fixtures.py +++ /dev/null @@ -1,51 +0,0 @@ -from __future__ import annotations - -import argparse -import shutil -import subprocess -import sys -from pathlib import Path - - -THIS_DIR = Path(__file__).resolve().parent -REPO_ROOT = THIS_DIR.parents[1] -HELPER_DIR = THIS_DIR / "matlab" - - -def _matlab_quote(value: str) -> str: - return value.replace("'", "''") - - -def default_matlab_repo(repo_root: Path | None = None) -> Path: - base = REPO_ROOT if repo_root is None else repo_root.resolve() - return base.parent / "nSTAT" - - -def build_matlab_batch(repo_root: Path, matlab_repo: Path) -> str: - return ( - f"addpath('{_matlab_quote(str(HELPER_DIR))}'); " - f"export_matlab_gold_fixtures('{_matlab_quote(str(repo_root))}','{_matlab_quote(str(matlab_repo))}'); " - "exit" - ) - - -def main(argv: list[str] | None = None) -> int: - parser = argparse.ArgumentParser(description="Generate MATLAB-derived gold fixtures for Python parity tests.") - parser.add_argument("--repo-root", type=Path, default=REPO_ROOT) - parser.add_argument("--matlab-repo", type=Path, default=None) - args = parser.parse_args(argv) - - repo_root = args.repo_root.resolve() - matlab_repo = default_matlab_repo(repo_root) if args.matlab_repo is None else args.matlab_repo.resolve() - if not matlab_repo.exists(): - raise FileNotFoundError(f"MATLAB reference repo not found at {matlab_repo}") - if shutil.which("matlab") is None: - raise FileNotFoundError("`matlab` is not on PATH") - - command = ["matlab", "-batch", build_matlab_batch(repo_root, matlab_repo)] - subprocess.run(command, cwd=repo_root, check=True) - return 0 - - -if __name__ == "__main__": - raise SystemExit(main(sys.argv[1:])) diff --git a/tools/parity/matlab/export_matlab_gold_fixtures.m b/tools/parity/matlab/export_matlab_gold_fixtures.m index a62e99ac..f96d1b95 100644 --- a/tools/parity/matlab/export_matlab_gold_fixtures.m +++ b/tools/parity/matlab/export_matlab_gold_fixtures.m @@ -19,17 +19,55 @@ function export_matlab_gold_fixtures(repoRoot, matlabRepoRoot) end export_signalobj_fixture(fixtureRoot); +export_confidence_interval_fixture(fixtureRoot); export_covariate_fixture(fixtureRoot); export_nspiketrain_fixture(fixtureRoot); export_nstcoll_fixture(fixtureRoot); +export_config_fixture(fixtureRoot); export_cif_fixture(fixtureRoot); export_analysis_fixture(fixtureRoot); +export_ksdiscrete_fixture(fixtureRoot); +export_fit_summary_fixture(fixtureRoot); export_point_process_fixture(fixtureRoot); +export_thinning_fixture(fixtureRoot); export_decoding_predict_fixture(fixtureRoot); +export_decoding_smoother_fixture(fixtureRoot); +export_hybrid_filter_fixture(fixtureRoot); export_nonlinear_decode_fixture(fixtureRoot); export_simulated_network_fixture(fixtureRoot); end +function export_confidence_interval_fixture(fixtureRoot) +t = (0:0.1:0.4)'; +bounds = [0.9 1.1; 1.9 2.1; 2.9 3.1; 3.9 4.1; 4.9 5.1]; +ci = ConfidenceInterval(t, bounds, 'CI', 'time', 's', 'a.u.', {'lo','hi'}, {'-.k'}); +ci.setColor('r'); +ci.setValue(0.9); +structure = ci.dataToStructure; +roundtrip = ConfidenceInterval.fromStructure(structure); + +payload = struct(); +payload.time = ci.time; +payload.bounds = ci.data; +payload.color = ci.color; +payload.value = ci.value; +payload.name = ci.name; +payload.xlabelval = ci.xlabelval; +payload.xunits = ci.xunits; +payload.yunits = ci.yunits; +payload.dataLabels = ci.dataLabels; +payload.plotProps = ci.plotProps; +payload.structure_time = structure.time; +payload.structure_values = structure.signals.values; +payload.roundtrip_bounds = roundtrip.data; +payload.roundtrip_color = roundtrip.color; +payload.roundtrip_value = roundtrip.value; +payload.roundtrip_name = roundtrip.name; +payload.roundtrip_plotProps = roundtrip.plotProps; + +save(fullfile(fixtureRoot, 'confidence_interval_exactness.mat'), '-struct', 'payload'); +end + function export_signalobj_fixture(fixtureRoot) t = (0:0.1:0.4)'; data = [1 0; 2 1; 4 0; 8 -1; 16 0]; @@ -128,6 +166,7 @@ function export_nstcoll_fixture(fixtureRoot) n2 = nspikeTrain([0.2], '2', 10, 0.0, 0.5, 'time', 's', 'spikes', 'spk', -1); coll = nstColl({n1, n2}); dataMat = coll.dataToMatrix([1 2], 0.1, 0.0, 0.5); +collapsed = coll.toSpikeTrain; payload = struct(); payload.numSpikeTrains = coll.numSpikeTrains; @@ -135,10 +174,43 @@ function export_nstcoll_fixture(fixtureRoot) payload.dataMatrix = dataMat; payload.firstSpikeTimes = coll.getNST(1).spikeTimes; payload.secondSpikeTimes = coll.getNST(2).spikeTimes; +payload.collapsedSpikeTimes = collapsed.spikeTimes; +payload.collapsedName = collapsed.name; +payload.collapsedMinTime = collapsed.minTime; +payload.collapsedMaxTime = collapsed.maxTime; +payload.collapsedSampleRate = collapsed.sampleRate; save(fullfile(fixtureRoot, 'nstcoll_exactness.mat'), '-struct', 'payload'); end +function export_config_fixture(fixtureRoot) +cfg = TrialConfig({{'Position','x'},{'Stimulus'}}, 2.0, [0 0.5 1.0], [], [], 0.5, 'stim_pos'); +cfg2 = TrialConfig({{'Stimulus'}}, 2.0, [], [], [], [], 'manual'); +structure = cfg.toStructure; +roundtrip = TrialConfig.fromStructure(structure); +coll = ConfigColl({cfg, cfg2}); +subset = coll.getSubsetConfigs([1 2]); +rebuilt = ConfigColl.fromStructure(coll.toStructure); + +payload = struct(); +payload.cfg_name = cfg.name; +payload.cfg_sampleRate = cfg.sampleRate; +payload.cfg_covLag = cfg.covLag; +payload.cfg_covMask = cfg.covMask; +payload.roundtrip_name = roundtrip.name; +payload.roundtrip_covLag = roundtrip.covLag; +payload.roundtrip_ensCovMask = roundtrip.ensCovMask; +payload.roundtrip_covMask = roundtrip.covMask; +payload.config_names = coll.getConfigNames(); +payload.subset_names = subset.getConfigNames(); +payload.rebuilt_names = rebuilt.getConfigNames(); +payload.rebuilt_first_name = rebuilt.getConfig(1).name; +payload.rebuilt_first_covLag = rebuilt.getConfig(1).covLag; +payload.rebuilt_first_ensCovMask = rebuilt.getConfig(1).ensCovMask; + +save(fullfile(fixtureRoot, 'config_exactness.mat'), '-struct', 'payload'); +end + function export_cif_fixture(fixtureRoot) cif = CIF([0.1 0.5], {'stim1', 'stim2'}, {'stim1', 'stim2'}, 'binomial'); stimVal = [0.6; -0.2]; @@ -174,6 +246,8 @@ function export_analysis_fixture(fixtureRoot) cfg.setName('stim'); fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl({cfg})); summary = FitResSummary({fit}); +Analysis.KSPlot(fit, 1, 0); +Analysis.plotFitResidual(fit, 0.01, 0); payload = struct(); payload.time = t; @@ -190,10 +264,108 @@ function export_analysis_fixture(fixtureRoot) payload.summaryAIC = summary.AIC(1); payload.summaryBIC = summary.BIC(1); payload.summarylogLL = summary.logLL(1); +payload.Z = fit.Z(:,1); +payload.U = fit.U(:,1); +payload.ks_xAxis = fit.KSStats.xAxis(:,1); +payload.ks_sorted = fit.KSStats.KSSorted(:,1); +payload.ks_stat = fit.KSStats.ks_stat(1); +payload.ks_pvalue = fit.KSStats.pValue(1); +payload.ks_within_conf_int = fit.KSStats.withinConfInt(1); +payload.residual_time = fit.Residual.time; +payload.residual_data = fit.Residual.data(:,1); save(fullfile(fixtureRoot, 'analysis_exactness.mat'), '-struct', 'payload'); end +function export_ksdiscrete_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})); + +pk = fit.lambda.data(:,1) .* (1 / fit.lambda.sampleRate); +pk = min(max(pk, 1e-10), 1); +spikeSignal = spikeTrain.getSigRep.data(:,1); +uniformDraws = [0.25; 0.75]; +[Z, ~, xAxis, ~, ~] = ksdiscrete_with_draws(pk, spikeSignal, 'spiketrain', uniformDraws); +U = 1 - exp(-Z); +KSSorted = sort(U, 'ascend'); +[differentDists, pValue, ksStat] = kstest2(xAxis, KSSorted); + +payload = struct(); +payload.time = t; +payload.stim_data = stimData; +payload.spike_times = spikeTrain.spikeTimes; +payload.lambda_time = fit.lambda.time; +payload.lambda_data = fit.lambda.data(:,1); +payload.sample_rate = fit.lambda.sampleRate; +payload.pk = pk; +payload.spike_signal = spikeSignal; +payload.uniform_draws = uniformDraws; +payload.Z = Z; +payload.U = U; +payload.xAxis = xAxis; +payload.KSSorted = KSSorted; +payload.compute_ks_stat = max(abs(KSSorted - xAxis)); +payload.ks_stat = ksStat; +payload.ks_pvalue = pValue; +payload.ks_within_conf_int = ~differentDists; + +save(fullfile(fixtureRoot, 'ksdiscrete_exactness.mat'), '-struct', 'payload'); +end + +function export_fit_summary_fixture(fixtureRoot) +t = (0:0.1:1.0)'; +lambdaData = [2.0 3.0; 2.5 3.5; 3.0 4.0; 3.5 4.5; 4.0 5.0; 4.5 5.5; 5.0 6.0; 5.5 6.5; 6.0 7.0; 6.5 7.5; 7.0 8.0]; +lambda = Covariate(t, lambdaData, '\lambda(t)', 'time', 's', 'Hz', {'stim','stim_hist'}); +st1 = nspikeTrain([0.1 0.4 0.7], '1', 0.1, 0.0, 1.0, 'time', 's', '', '', -1); +st2 = nspikeTrain([0.2 0.5 0.8], '2', 0.1, 0.0, 1.0, 'time', 's', '', '', -1); +stimCfg = TrialConfig({{'Stimulus', 'stim'}}, 10, [], []); +stimCfg.setName('stim'); +stimHistCfg = TrialConfig({{'Stimulus', 'stim'}}, 10, [0 0.1 0.2], []); +stimHistCfg.setName('stim_hist'); +configColl = ConfigColl({stimCfg, stimHistCfg}); +covLabels = {{'stim'}, {'stim','hist1','hist2'}}; +numHist = [0 2]; +histObjects = {[], []}; +ensHistObj = {[], []}; +b1 = {[0.5], [0.3 -0.1 -0.05]}; +b2 = {[0.4], [0.25 -0.08 -0.02]}; +stats1 = {struct('se', [0.05], 'p', [0.01]), struct('se', [0.04 0.03 0.02], 'p', [0.02 0.04 0.06])}; +stats2 = {struct('se', [0.06], 'p', [0.03]), struct('se', [0.05 0.04 0.03], 'p', [0.01 0.03 0.07])}; +fit1 = FitResult(st1, covLabels, numHist, histObjects, ensHistObj, lambda, b1, [1.0 2.0], stats1, [11.0 7.0], [12.0 8.0], [3.0 5.0], configColl, {}, {}, 'poisson'); +fit2 = FitResult(st2, covLabels, numHist, histObjects, ensHistObj, lambda, b2, [1.5 2.5], stats2, [13.0 9.0], [14.0 10.0], [2.0 4.0], configColl, {}, {}, 'poisson'); +fit1.KSStats.ks_stat = [0.25 0.50]; +fit1.KSStats.pValue = [0.90 0.40]; +fit1.KSStats.withinConfInt = [1 1]; +fit2.KSStats.ks_stat = [0.35 0.55]; +fit2.KSStats.pValue = [0.80 0.30]; +fit2.KSStats.withinConfInt = [1 0]; +summary = FitResSummary({fit1, fit2}); +dAIC = summary.getDiffAIC(1, 0); +dBIC = summary.getDiffBIC(1, 0); +dlogLL = summary.getDifflogLL(1, 0); + +payload = struct(); +payload.fitNames = summary.fitNames; +payload.neuronNumbers = summary.neuronNumbers; +payload.AIC = summary.AIC; +payload.BIC = summary.BIC; +payload.logLL = summary.logLL; +payload.KSStats = summary.KSStats; +payload.KSPvalues = summary.KSPvalues; +payload.withinConfInt = summary.withinConfInt; +payload.diffAIC = dAIC; +payload.diffBIC = dBIC; +payload.difflogLL = dlogLL; + +save(fullfile(fixtureRoot, 'fit_summary_exactness.mat'), '-struct', 'payload'); +end + function export_point_process_fixture(fixtureRoot) rng(5); Ts = .001; @@ -216,9 +388,52 @@ function export_point_process_fixture(fixtureRoot) payload.lambda_head = lambda.data(1:8, 1); payload.spike_counts = spikeCounts; +detTime = (0:Ts:0.01)'; +detStim = sin(2*pi*1*detTime); +detUniforms = [0.99; 0.01; 0.99; 0.2; 0.8; 0.05; 0.95; 0.15; 0.6; 0.4; 0.3]; +[hNum, ~] = tfdata(H, 'v'); +[sNum, ~] = tfdata(S, 'v'); +[eNum, ~] = tfdata(E, 'v'); +[detRate, detDelta, detHist, detEta, detSpikes] = local_discrete_point_process(mu, hNum, sNum, eNum, detStim, zeros(size(detStim)), detUniforms, Ts, 0); +payload.det_time = detTime; +payload.det_uniforms = detUniforms; +payload.det_stimulus = detStim; +payload.det_rate_hz = detRate; +payload.det_lambda_delta = detDelta; +payload.det_history_effect = detHist; +payload.det_eta = detEta; +payload.det_spike_indicator = detSpikes; + save(fullfile(fixtureRoot, 'point_process_exactness.mat'), '-struct', 'payload'); end +function export_thinning_fixture(fixtureRoot) +t = (0:0.1:1.0)'; +lambdaData = [2.0; 4.0; 3.5; 5.0; 1.5; 2.5; 4.5; 3.0; 2.0; 1.0; 0.5]; +lambdaCov = Covariate(t, lambdaData, '\lambda(t)', 'time', 's', 'Hz', {'\lambda'}); +arrivalUniforms = [0.95; 0.5; 0.2; 0.8; 0.3; 0.7; 0.4; 0.6]; +thinningUniforms = [0.1; 0.9; 0.2; 0.7; 0.3; 0.6; 0.4; 0.5]; +maxTimeRes = 0.05; +[proposalCount, lambdaBound, interarrival, candidateTimes, lambdaRatio, thinningUsed, acceptedTimes, roundedTimes] = ... + local_thin_from_lambda(lambdaCov, arrivalUniforms, thinningUniforms, maxTimeRes); + +payload = struct(); +payload.time = t; +payload.lambda_data = lambdaData; +payload.lambda_bound = lambdaBound; +payload.proposal_count = proposalCount; +payload.arrival_uniforms = arrivalUniforms(1:proposalCount); +payload.interarrival_times = interarrival; +payload.candidate_spike_times = candidateTimes; +payload.lambda_ratio = lambdaRatio; +payload.thinning_uniforms = thinningUsed; +payload.accepted_spike_times = acceptedTimes; +payload.rounded_spike_times = roundedTimes; +payload.maxTimeRes = maxTimeRes; + +save(fullfile(fixtureRoot, 'thinning_exactness.mat'), '-struct', 'payload'); +end + function export_decoding_predict_fixture(fixtureRoot) x_u = [0.1; -0.2]; W_u = [1.0 0.1; 0.1 2.0]; @@ -237,6 +452,72 @@ function export_decoding_predict_fixture(fixtureRoot) save(fullfile(fixtureRoot, 'decoding_predict_exactness.mat'), '-struct', 'payload'); end +function export_decoding_smoother_fixture(fixtureRoot) +A = 1.0; +Q = 0.05; +dN = [0 1 0 1]; +lags = 2; +mu = -1.0; +beta = 0.75; +fitType = 'binomial'; +delta = 0.1; +[x_pLag, W_pLag, x_uLag, W_uLag] = DecodingAlgorithms.PP_fixedIntervalSmoother(A, Q, dN, lags, mu, beta, fitType, delta); + +payload = struct(); +payload.A = A; +payload.Q = Q; +payload.dN = dN; +payload.lags = lags; +payload.mu = mu; +payload.beta = beta; +payload.fitType = fitType; +payload.delta = delta; +payload.x_pLag = x_pLag; +payload.W_pLag = W_pLag; +payload.x_uLag = x_uLag; +payload.W_uLag = W_uLag; + +save(fullfile(fixtureRoot, 'decoding_smoother_exactness.mat'), '-struct', 'payload'); +end + +function export_hybrid_filter_fixture(fixtureRoot) +A = {1.0, 0.9}; +Q = {0.02, 0.05}; +p_ij = [0.95 0.05; 0.10 0.90]; +Mu0 = [0.6; 0.4]; +dN = [0 1 1 0 1]; +mu = -1.0; +beta = 0.75; +fitType = 'binomial'; +binwidth = 0.1; +[S_est, X, W, MU_u, X_s, W_s, pNGivenS] = DecodingAlgorithms.PPHybridFilterLinear( ... + A, Q, p_ij, Mu0, dN, mu, beta, fitType, binwidth); + +payload = struct(); +payload.A1 = A{1}; +payload.A2 = A{2}; +payload.Q1 = Q{1}; +payload.Q2 = Q{2}; +payload.p_ij = p_ij; +payload.Mu0 = Mu0; +payload.dN = dN; +payload.mu = mu; +payload.beta = beta; +payload.fitType = fitType; +payload.binwidth = binwidth; +payload.S_est = S_est; +payload.X = X; +payload.W = W; +payload.MU_u = MU_u; +payload.X_s_1 = X_s{1}; +payload.X_s_2 = X_s{2}; +payload.W_s_1 = W_s{1}; +payload.W_s_2 = W_s{2}; +payload.pNGivenS = pNGivenS; + +save(fullfile(fixtureRoot, 'hybrid_filter_exactness.mat'), '-struct', 'payload'); +end + function export_nonlinear_decode_fixture(fixtureRoot) A = eye(2); Q = 0.01 * eye(2); @@ -315,9 +596,197 @@ function export_simulated_network_fixture(fixtureRoot) payload.state_head = yout(1:5,1:2); payload.spike_counts = [sum(yout(:,1) > .5), sum(yout(:,2) > .5)]; +detTime = (0:Ts:0.009)'; +detStim = sin(2*pi*1*detTime); +detUniforms = [ + 0.99 0.99; + 0.01 0.99; + 0.99 0.02; + 0.20 0.80; + 0.60 0.10; + 0.05 0.95; + 0.90 0.40; + 0.15 0.30; + 0.70 0.20; + 0.25 0.85 +]; +[stateDet, probDet, etaDet, histDet, ensDet] = local_two_neuron_network(mu{1}, mu{2}, h1Num, h2Num, s1Num, s2Num, e1Num, e2Num, detStim, detUniforms); +payload.det_time = detTime; +payload.det_uniforms = detUniforms; +payload.det_stimulus = detStim; +payload.det_probability = probDet; +payload.det_state = stateDet; +payload.det_eta = etaDet; +payload.det_history_effect = histDet; +payload.det_ensemble_effect = ensDet; + save(fullfile(fixtureRoot, 'simulated_network_exactness.mat'), '-struct', 'payload'); end +function [rateHz, lambdaDelta, histEffect, eta, spikes] = local_discrete_point_process(mu, hNum, sNum, eNum, stimInput, ensInput, uniforms, Ts, simTypeSelect) +stimInput = stimInput(:); +ensInput = ensInput(:); +uniforms = uniforms(:); +n = length(stimInput); +rateHz = zeros(n,1); +lambdaDelta = zeros(n,1); +histEffect = zeros(n,1); +eta = zeros(n,1); +spikes = zeros(n,1); +for idx = 1:n + stimEffect = 0; + for lag = 1:length(sNum) + if idx-lag+1 >= 1 + stimEffect = stimEffect + sNum(lag) * stimInput(idx-lag+1); + end + end + ensEffect = 0; + for lag = 1:length(eNum) + if idx-lag+1 >= 1 + ensEffect = ensEffect + eNum(lag) * ensInput(idx-lag+1); + end + end + histVal = 0; + for lag = 1:length(hNum) + if idx-lag >= 1 + histVal = histVal + hNum(lag) * spikes(idx-lag); + end + end + histEffect(idx) = histVal; + eta(idx) = mu + stimEffect + ensEffect + histVal; + if simTypeSelect == 0 + lambdaDelta(idx) = exp(eta(idx)) / (1 + exp(eta(idx))); + rateHz(idx) = lambdaDelta(idx) / Ts; + spikes(idx) = uniforms(idx) < lambdaDelta(idx); + else + rateHz(idx) = exp(eta(idx)); + lambdaDelta(idx) = 1 - exp(-rateHz(idx) * Ts); + spikes(idx) = uniforms(idx) < lambdaDelta(idx); + end +end +end + +function [proposalCount, lambdaBound, interarrival, candidateTimes, lambdaRatio, thinningUsed, acceptedTimes, roundedTimes] = local_thin_from_lambda(lambdaCov, arrivalUniforms, thinningUniforms, maxTimeRes) +lambdaBound = max(lambdaCov); +Tmax = lambdaCov.maxTime; +proposalCount = ceil(lambdaBound * (1.5 * Tmax)); +arrivalUniforms = arrivalUniforms(:); +thinningUniforms = thinningUniforms(:); +u = arrivalUniforms(1:proposalCount); +interarrival = -log(u) ./ lambdaBound; +candidateTimes = cumsum(interarrival); +candidateTimes = candidateTimes(candidateTimes <= Tmax); +if isempty(candidateTimes) + lambdaRatio = []; + thinningUsed = []; + acceptedTimes = []; +else + lambdaRatio = lambdaCov.getValueAt(candidateTimes) ./ lambdaBound; + thinningUsed = thinningUniforms(1:length(lambdaRatio)); + acceptedTimes = candidateTimes(lambdaRatio >= thinningUsed); +end +if isempty(maxTimeRes) + roundedTimes = acceptedTimes; +else + roundedTimes = unique(ceil(acceptedTimes ./ maxTimeRes) * maxTimeRes); +end +end + +function [stateMat, probMat, etaMat, histMat, ensMat] = local_two_neuron_network(mu1, mu2, h1Num, h2Num, s1Num, s2Num, e1Num, e2Num, stimInput, uniforms) +stimInput = stimInput(:); +uniforms = double(uniforms); +n = length(stimInput); +stateMat = zeros(n,2); +probMat = zeros(n,2); +etaMat = zeros(n,2); +histMat = zeros(n,2); +ensMat = zeros(n,2); +for idx = 1:n + hist1 = 0; + hist2 = 0; + for lag = 1:length(h1Num) + if idx-lag >= 1 + hist1 = hist1 + h1Num(lag) * stateMat(idx-lag,1); + end + end + for lag = 1:length(h2Num) + if idx-lag >= 1 + hist2 = hist2 + h2Num(lag) * stateMat(idx-lag,2); + end + end + stim1 = s1Num(1) * stimInput(idx); + stim2 = s2Num(1) * stimInput(idx); + ens1 = 0; + ens2 = 0; + if idx > 1 + ens1 = e1Num(1) * stateMat(idx-1,2); + ens2 = e2Num(1) * stateMat(idx-1,1); + end + histMat(idx,:) = [hist1 hist2]; + ensMat(idx,:) = [ens1 ens2]; + etaMat(idx,:) = [mu1 + hist1 + stim1 + ens1, mu2 + hist2 + stim2 + ens2]; + probMat(idx,:) = exp(etaMat(idx,:)) ./ (1 + exp(etaMat(idx,:))); + stateMat(idx,1) = uniforms(idx,1) < probMat(idx,1); + stateMat(idx,2) = uniforms(idx,2) < probMat(idx,2); +end +end + +function [rst,varargout] = ksdiscrete_with_draws(pk,st,spikeflag,draws) +[m1,m2]=size(pk); +if (m1 ~=1 && m2 ~=1); error('pk must be a vector'); end +if (m2>m1); pk=pk'; end +[m1,~]=size(pk); +if any(pk<0) || any(pk>1); error('all values for pk must be within [0,1]'); end + +if strcmp(spikeflag,'spiketrain') + [n1,n2]=size(st); + if (n1 ~=1 && n2 ~=1); error('spike train must be a vector'); end + if (n2>n1); st=st'; end + if m1 ~= n1; error('pk and spike train must be same length'); end + spikeindicies=find(st==1); + Nspikes=length(spikeindicies); +elseif strcmp(spikeflag,'spikeind') + [n1,n2]=size(st); + if (n1 ~=1 && n2 ~=1); error('spike indicies must be a vector'); end + if (n2>n1); st=st'; end + spikeindicies=unique(st); + Nspikes=length(spikeindicies); +else + error('invalid spikeflag'); +end + +if isempty(spikeindicies) + rst = pk; + return; +end +if spikeindicies(1)<1; error('There is at least one spike with index less than 0'); end +if spikeindicies(Nspikes)>length(pk); error('There is at least one spike with a index greater than the length of pk'); end + +qk=-log(1-pk); +rst=zeros(Nspikes-1,1); +rstold=zeros(Nspikes-1,1); +for r=1:Nspikes-1 + ind1=spikeindicies(r); + ind2=spikeindicies(r+1); + total=sum(qk(ind1+1:ind2-1)); + delta=-(1/qk(ind2))*log(1-draws(r)*(1-exp(-qk(ind2)))); + if(delta~=0) + total=total+qk(ind2)*delta; + end + rst(r)=total; + rstold(r)=sum(qk(ind1+1:ind2)); +end + +rstsort=sort(rst); +varargout{1}=rstsort; +inrst=1/(Nspikes-1); +xrst=(0.5*inrst:inrst:1-0.5*inrst)'; +varargout{2}=xrst; +cb=1.36*sqrt(inrst); +varargout{3}=cb; +varargout{4}=sort(rstold); +end + function cifObj = build_polynomial_binomial_cif(beta) beta = beta(:)'; syms x y real diff --git a/tools/release/release_gate_lib.py b/tools/release/release_gate_lib.py new file mode 100644 index 00000000..aeb33556 --- /dev/null +++ b/tools/release/release_gate_lib.py @@ -0,0 +1,56 @@ +from __future__ import annotations + +import argparse +import shutil +import subprocess +import sys +from pathlib import Path + + +def repo_root() -> Path: + return Path(__file__).resolve().parents[2] + + +def build_release_gate_commands(repo: Path | None = None) -> list[list[str]]: + base = repo_root() if repo is None else repo.resolve() + return [ + [ + sys.executable, + "-m", + "pytest", + "-q", + "tests/test_zernike.py", + "tests/test_matlab_gold_fixtures.py", + "tests/test_signalobj_fidelity.py", + "tests/test_nspiketrain_fidelity.py", + "tests/test_workflow_fidelity.py", + "tests/test_class_fidelity_audit.py", + "tests/test_matlab_symbol_surface.py", + "tests/test_simulink_fidelity_audit.py", + "tests/test_parity_report.py", + "tests/test_cleanroom_boundary.py", + ], + [sys.executable, str(base / "tools" / "notebooks" / "run_notebooks.py"), "--group", "parity_core", "--timeout", "900"], + [sys.executable, str(base / "tools" / "parity" / "build_report.py")], + ] + + +def run_release_gate(repo: Path | None = None) -> None: + base = repo_root() if repo is None else repo.resolve() + for command in build_release_gate_commands(base): + subprocess.run(command, cwd=base, check=True) + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description="Run the pure-Python fidelity release gate.") + parser.add_argument("--repo-root", type=Path, default=None) + args = parser.parse_args(argv) + + if shutil.which(sys.executable) is None: + raise FileNotFoundError(f"Python executable is not on PATH: {sys.executable}") + run_release_gate(args.repo_root) + return 0 + + +__all__ = ["build_release_gate_commands", "main", "repo_root", "run_release_gate"] + diff --git a/tools/release/run_fidelity_gate.py b/tools/release/run_fidelity_gate.py index 7121d0c7..fa1b80b8 100644 --- a/tools/release/run_fidelity_gate.py +++ b/tools/release/run_fidelity_gate.py @@ -9,7 +9,7 @@ if str(REPO_ROOT) not in sys.path: sys.path.insert(0, str(REPO_ROOT)) -from nstat.release_check import main +from tools.release.release_gate_lib import main if __name__ == "__main__":