diff --git a/notebooks/AnalysisExamples.ipynb b/notebooks/AnalysisExamples.ipynb index 64b86ee9..68f00805 100644 --- a/notebooks/AnalysisExamples.ipynb +++ b/notebooks/AnalysisExamples.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "98fc6fc8", + "id": "4ba2084f", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `AnalysisExamples.mlx`\n", "- Fidelity status: `exact`\n", - "- Remaining justified differences: Complete MATLAB standard-GLM workflow with the canonical glm_data.mat dataset and real KS/model-visualization figures. Only inherent GLM solver numerics and matplotlib styling differ." + "- Remaining justified differences: The notebook now follows the MATLAB standard-GLM workflow with the canonical `glm_data.mat` dataset and real KS/model-visualization figures; coefficient values and styling still vary modestly because the Python GLM backend and plotting defaults differ from MATLAB.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "7807842d", + "id": "1bf27b7c", "metadata": {}, "outputs": [], "source": [ @@ -44,12 +44,14 @@ "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", "__tracker = FigureTracker(topic=\"AnalysisExamples\", output_root=OUTPUT_ROOT, expected_count=4)\n", "\n", + "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", " fig.clear()\n", " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _poisson_standard_errors(design_matrix, result):\n", " x = np.asarray(design_matrix, dtype=float)\n", " if x.ndim == 1:\n", @@ -60,6 +62,7 @@ " cov = np.linalg.pinv(x_aug.T @ (lam[:, None] * x_aug))\n", " return np.sqrt(np.clip(np.diag(cov), 0.0, None))\n", "\n", + "\n", "T = np.asarray(GLM_DATA[\"T\"], dtype=float).reshape(-1)\n", "xN = np.asarray(GLM_DATA[\"xN\"], dtype=float).reshape(-1)\n", "yN = np.asarray(GLM_DATA[\"yN\"], dtype=float).reshape(-1)\n", @@ -68,25 +71,25 @@ "x_at_spiketimes = np.asarray(GLM_DATA[\"x_at_spiketimes\"], dtype=float).reshape(-1)\n", "y_at_spiketimes = np.asarray(GLM_DATA[\"y_at_spiketimes\"], dtype=float).reshape(-1)\n", "sample_rate = 1.0 / float(np.median(np.diff(T)))\n", - "nst = nspikeTrain(spiketimes, name=\"1\", minTime=float(T[0]), maxTime=float(T[-1]), makePlots=-1)" + "nst = nspikeTrain(spiketimes, name=\"1\", minTime=float(T[0]), maxTime=float(T[-1]), makePlots=-1)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "37ac20c9", + "id": "39e80ef9", "metadata": {}, "outputs": [], "source": [ "# SECTION 1: Analysis Examples\n", "plt.close(\"all\")\n", - "print({\"n_samples\": int(T.shape[0]), \"n_spikes\": int(spiketimes.shape[0]), \"sample_rate_hz\": round(sample_rate, 3)})" + "print({\"n_samples\": int(T.shape[0]), \"n_spikes\": int(spiketimes.shape[0]), \"sample_rate_hz\": round(sample_rate, 3)})\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "dbdc74f9", + "id": "b2e4c2dc", "metadata": {}, "outputs": [], "source": [ @@ -104,13 +107,13 @@ "x_quadratic = np.column_stack([xN, yN, xN**2, yN**2, xN * yN])\n", "linear_fit = fit_poisson_glm(x_linear, spikes_binned)\n", "quadratic_fit = fit_poisson_glm(x_quadratic, spikes_binned)\n", - "centered_fit = fit_poisson_glm(x_quadratic_centered, spikes_binned)" + "centered_fit = fit_poisson_glm(x_quadratic_centered, spikes_binned)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "5c38cff1", + "id": "15ae5117", "metadata": {}, "outputs": [], "source": [ @@ -122,13 +125,13 @@ "ax.set_aspect(\"equal\", adjustable=\"box\")\n", "ax.set_xlabel(\"x position (m)\")\n", "ax.set_ylabel(\"y position (m)\")\n", - "ax.set_title(\"Rat trajectory with spike locations\")" + "ax.set_title(\"Rat trajectory with spike locations\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "5af52914", + "id": "a8839159", "metadata": {}, "outputs": [], "source": [ @@ -141,13 +144,13 @@ "ax.errorbar(xpos, centered_beta, yerr=centered_se, fmt=\".\", color=\"tab:blue\", capsize=3)\n", "ax.set_xticks(xpos, [\"baseline\", \"x\", \"y\", \"x^2\", \"y^2\", \"x*y\"])\n", "ax.set_ylabel(\"coefficient value\")\n", - "ax.set_title(\"Quadratic GLM coefficients\")" + "ax.set_title(\"Quadratic GLM coefficients\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "5cac7309", + "id": "6bda2bd9", "metadata": {}, "outputs": [], "source": [ @@ -166,13 +169,13 @@ "ax.set_xlabel(\"x position (m)\")\n", "ax.set_ylabel(\"y position (m)\")\n", "ax.set_zlabel(\"lambda\")\n", - "ax.set_title(\"Quadratic GLM spatial intensity\")" + "ax.set_title(\"Quadratic GLM spatial intensity\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "36dd8e70", + "id": "0d526433", "metadata": {}, "outputs": [], "source": [ @@ -186,13 +189,13 @@ " \"linear_mean_rate_hz\": round(float(np.mean(lambda_linear_hz)), 4),\n", " \"quadratic_mean_rate_hz\": round(float(np.mean(lambda_quadratic_hz)), 4),\n", " }\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "2d8b81fd", + "id": "59472d54", "metadata": {}, "outputs": [], "source": [ @@ -217,7 +220,7 @@ "ax.set_ylabel(\"Empirical CDF of Rescaled ISIs\")\n", "ax.set_title(\"KS Plot with 95% Confidence Intervals\")\n", "ax.legend(loc=\"lower right\", frameon=False)\n", - "__tracker.finalize()" + "__tracker.finalize()\n" ] } ], @@ -234,4 +237,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/AnalysisExamples2.ipynb b/notebooks/AnalysisExamples2.ipynb index 426a9bd8..6a21b2a0 100644 --- a/notebooks/AnalysisExamples2.ipynb +++ b/notebooks/AnalysisExamples2.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "66d56086", + "id": "04a1108d", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `AnalysisExamples2.mlx`\n", "- Fidelity status: `exact`\n", - "- Remaining justified differences: Complete MATLAB toolbox workflow on the canonical glm_data.mat dataset with executable Trial, ConfigColl, and Analysis calls. Only inherent GLM solver numerics and plot styling differ." + "- Remaining justified differences: The notebook now follows the MATLAB toolbox workflow on the canonical `glm_data.mat` dataset with executable `Trial`, `ConfigColl`, and `Analysis` calls; exact coefficients and plot styling still vary modestly because the Python GLM backend differs from MATLAB.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "62e21501", + "id": "f4ecd812", "metadata": {}, "outputs": [], "source": [ @@ -42,7 +42,8 @@ "\n", "GLM_DATA = load_glm_data_for_notebook()\n", "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", - "__tracker = FigureTracker(topic=\"AnalysisExamples2\", output_root=OUTPUT_ROOT, expected_count=4)\n", + "__tracker = FigureTracker(topic=\"AnalysisExamples2\", output_root=OUTPUT_ROOT, expected_count=5)\n", + "\n", "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", @@ -50,6 +51,7 @@ " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "T = np.asarray(GLM_DATA[\"T\"], dtype=float).reshape(-1)\n", "xN = np.asarray(GLM_DATA[\"xN\"], dtype=float).reshape(-1)\n", "yN = np.asarray(GLM_DATA[\"yN\"], dtype=float).reshape(-1)\n", @@ -65,36 +67,36 @@ "velocity = Covariate(T, np.column_stack([vxN, vyN]), \"Velocity\", \"time\", \"s\", \"m/s\", [\"v_x\", \"v_y\"])\n", "radial = Covariate(T, np.column_stack([xN, yN, xN**2, yN**2, xN * yN]), \"Radial\", \"time\", \"s\", \"m\", [\"x\", \"y\", \"x^2\", \"y^2\", \"x*y\"])\n", "values_at_spiketimes = position.getValueAt(spiketimes)\n", - "values_at_spiketimes_upsampled = position.resample(1.0 / np.min(np.diff(spiketimes))).getValueAt(spiketimes)" + "values_at_spiketimes_upsampled = position.resample(1.0 / np.min(np.diff(spiketimes))).getValueAt(spiketimes)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "1836e297", + "id": "4dd70d35", "metadata": {}, "outputs": [], "source": [ "# SECTION 1: Analysis Examples 2\n", "plt.close(\"all\")\n", - "print({\"n_samples\": int(T.shape[0]), \"n_spikes\": int(spiketimes.shape[0]), \"analysis_sample_rate_hz\": sample_rate})" + "print({\"n_samples\": int(T.shape[0]), \"n_spikes\": int(spiketimes.shape[0]), \"analysis_sample_rate_hz\": sample_rate})\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "bf657cd0", + "id": "da89b88a", "metadata": {}, "outputs": [], "source": [ "# SECTION 2: load the rat trajectory and spiking data\n", - "print({\"position_shape\": list(position.data.shape), \"velocity_shape\": list(velocity.data.shape), \"radial_shape\": list(radial.data.shape)})" + "print({\"position_shape\": list(position.data.shape), \"velocity_shape\": list(velocity.data.shape), \"radial_shape\": list(radial.data.shape)})\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "fe47aacc", + "id": "4346a9a7", "metadata": {}, "outputs": [], "source": [ @@ -104,13 +106,13 @@ " \"direct_spike_position_head\": np.asarray(values_at_spiketimes[:3], dtype=float).round(4).tolist(),\n", " \"upsampled_spike_position_head\": np.asarray(values_at_spiketimes_upsampled[:3], dtype=float).round(4).tolist(),\n", " }\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "2f28e3be", + "id": "d5785121", "metadata": {}, "outputs": [], "source": [ @@ -122,13 +124,13 @@ "ax.set_aspect(\"equal\", adjustable=\"box\")\n", "ax.set_xlabel(\"x position (m)\")\n", "ax.set_ylabel(\"y position (m)\")\n", - "ax.set_title(\"Trajectory and interpolated spike locations\")" + "ax.set_title(\"Trajectory and interpolated spike locations\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "d40c40c9", + "id": "6cf72911", "metadata": {}, "outputs": [], "source": [ @@ -141,13 +143,13 @@ " TrialConfig([[\"Baseline\", \"mu\"], [\"Radial\", \"x\", \"y\", \"x^2\", \"y^2\", \"x*y\"]], sampleRate=sample_rate, history=[], name=\"Quadratic\"),\n", " TrialConfig([[\"Baseline\", \"mu\"], [\"Radial\", \"x\", \"y\", \"x^2\", \"y^2\", \"x*y\"]], sampleRate=sample_rate, history=[0.0, 1.0 / sample_rate], name=\"Quadratic+Hist\"),\n", "]\n", - "tcc = ConfigColl(tc)" + "tcc = ConfigColl(tc)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "f9160bdb", + "id": "e3b67df9", "metadata": {}, "outputs": [], "source": [ @@ -155,13 +157,13 @@ "fitResults = Analysis.RunAnalysisForAllNeurons(trial, tcc, 0)\n", "fig = _prepare_figure(\"fitResults.plotResults\", figsize=(11.0, 8.0))\n", "fitResults.plotResults(handle=fig)\n", - "print({\"config_names\": fitResults.configNames, \"aic\": np.asarray(fitResults.AIC, dtype=float).round(3).tolist()})" + "print({\"config_names\": fitResults.configNames, \"aic\": np.asarray(fitResults.AIC, dtype=float).round(3).tolist()})\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "a2fafcc8", + "id": "79879555", "metadata": {}, "outputs": [], "source": [ @@ -178,21 +180,34 @@ "ax.set_xlabel(\"x position (m)\")\n", "ax.set_ylabel(\"y position (m)\")\n", "ax.set_zlabel(\"lambda\")\n", - "ax.set_title(\"Toolbox-model spatial intensity comparison\")" + "ax.set_title(\"Toolbox-model spatial intensity comparison\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "64cdac30", + "id": "eef1351f", "metadata": {}, "outputs": [], - "source": "# SECTION 8: Toolbox vs. Standard GLM comparison\n# Compare the nSTAT fit with a standalone glmfit using the same Quadratic covariates\n# MATLAB: [b,dev,stats] = glmfit([xN yN xN.^2 yN.^2 xN.*yN], spikes_binned, 'poisson');\n# b - fitResults.b{2} % should be close to zero\nX_quad = np.column_stack([xN, yN, xN**2, yN**2, xN * yN])\nglm_result = fit_poisson_glm(X_quad, spikes_binned)\nb = np.concatenate([[glm_result.intercept], glm_result.coefficients])\nb_diff = b - fitResults.getCoeffs(2)\nprint(\"b - fitResults.b{2} =\", b_diff)" + "source": [ + "# SECTION 8: Toolbox vs. Standard GLM comparison\n", + "standard_fit = fit_poisson_glm(np.column_stack([np.ones_like(xN), xN, yN, xN**2, yN**2, xN * yN]), spikes_binned, include_intercept=False)\n", + "coeff_diff = np.asarray(standard_fit.coefficients - fitResults.getCoeffs(2)[0], dtype=float)\n", + "fig = _prepare_figure(\"b-fitResults.b{2}\", figsize=(7.0, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "labels = [\"mu\", \"x\", \"y\", \"x^2\", \"y^2\", \"x*y\"]\n", + "ax.bar(np.arange(coeff_diff.size), coeff_diff, color=\"tab:blue\")\n", + "ax.axhline(0.0, color=\"0.3\", linestyle=\"--\", linewidth=1.0)\n", + "ax.set_xticks(np.arange(coeff_diff.size), labels, rotation=20)\n", + "ax.set_ylabel(\"standard minus toolbox\")\n", + "ax.set_title(\"Coefficient agreement between workflows\")\n", + "print({\"quadratic_coeff_diff_max_abs\": round(float(np.max(np.abs(coeff_diff))), 6)})\n" + ] }, { "cell_type": "code", "execution_count": null, - "id": "8782d383", + "id": "04cd0c92", "metadata": {}, "outputs": [], "source": [ @@ -208,7 +223,7 @@ "ax.set_ylabel(\"AIC\")\n", "ax.set_title(\"History-lag model comparison\")\n", "print({\"history_config_names\": histConfigs.getConfigNames(), \"summary_fit_names\": histSummary.fitNames})\n", - "__tracker.finalize()" + "__tracker.finalize()\n" ] } ], @@ -218,7 +233,7 @@ }, "nstat": { "expected_figures": 5, - "run_group": "full", + "run_group": "smoke", "style": "python-example", "topic": "AnalysisExamples2" } diff --git a/notebooks/DecodingExample.ipynb b/notebooks/DecodingExample.ipynb index bbb4eab6..8b2e695c 100644 --- a/notebooks/DecodingExample.ipynb +++ b/notebooks/DecodingExample.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "0813e1e0", + "id": "6ff93288", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `DecodingExample.mlx`\n", - "- Fidelity status: `exact`\n", - "- Remaining justified differences: Workflow, model fitting, and decoded-stimulus figures follow the MATLAB helpfile. Only stochastic simulation draws and Python plotting defaults cause trace-level variation." + "- Fidelity status: `high_fidelity`\n", + "- Remaining justified differences: Workflow, model fitting, and decoded-stimulus figures now follow the MATLAB helpfile closely; exact traces still depend on stochastic simulation draws and Python plotting defaults.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "44d88ec9", + "id": "bb9c9790", "metadata": {}, "outputs": [], "source": [ @@ -42,12 +42,14 @@ "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", "__tracker = FigureTracker(topic=\"DecodingExample\", output_root=OUTPUT_ROOT, expected_count=5)\n", "\n", + "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", " fig.clear()\n", " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _plot_raster(ax, spike_coll):\n", " for row in range(1, spike_coll.numSpikeTrains + 1):\n", " train = spike_coll.getNST(row)\n", @@ -57,6 +59,7 @@ " ax.set_ylabel(\"Neuron\")\n", " ax.set_ylim(0.5, spike_coll.numSpikeTrains + 0.5)\n", "\n", + "\n", "def _plot_decoded_ci(ax, time, decoded, cov, stim, title):\n", " center = np.asarray(decoded, dtype=float).reshape(-1)\n", " variance = np.asarray(cov, dtype=float).reshape(-1)\n", @@ -72,14 +75,15 @@ " ax.set_xlabel(\"time (s)\")\n", " ax.legend(loc=\"upper right\", frameon=False, fontsize=8)\n", "\n", + "\n", "# SECTION 0: STIMULUS DECODING\n", - "# In this example we decode a univariate stimulus from simulated point-process observations by following the MATLAB DecodingExample workflow." + "# In this example we decode a univariate stimulus from simulated point-process observations by following the MATLAB DecodingExample workflow.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "2f1ec431", + "id": "c9032274", "metadata": {}, "outputs": [], "source": [ @@ -106,13 +110,13 @@ "axs[1].plot(time, lambda_cov.data[:, 0], color=\"b\", linewidth=2.0)\n", "axs[1].set_title(\"Conditional intensity λ(t)\")\n", "axs[1].set_xlabel(\"time (s)\")\n", - "axs[1].set_ylabel(\"Hz\")" + "axs[1].set_ylabel(\"Hz\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "a7048402", + "id": "282af22b", "metadata": {}, "outputs": [], "source": [ @@ -142,7 +146,7 @@ ")\n", "results = Analysis.RunAnalysisForAllNeurons(trial, cfgColl, 0)\n", "\n", - "paramEst = np.column_stack([fit.getCoeffs(2)[:2] for fit in results])\n", + "paramEst = np.column_stack([fit.getCoeffs(2)[0][:2] for fit in results])\n", "meanParams = np.mean(paramEst, axis=1)\n", "aic_matrix = np.vstack([fit.AIC for fit in results])\n", "logll_matrix = np.vstack([fit.logLL for fit in results])\n", @@ -170,13 +174,13 @@ "axs[0].set_title(\"Mean AIC across neurons\")\n", "axs[1].bar(xloc, np.mean(logll_matrix, axis=0), color=[\"0.6\", \"0.3\"])\n", "axs[1].set_xticks(xloc, config_names, rotation=15)\n", - "axs[1].set_title(\"Mean log-likelihood across neurons\")" + "axs[1].set_title(\"Mean log-likelihood across neurons\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "df079e59", + "id": "49257a4c", "metadata": {}, "outputs": [], "source": [ @@ -194,7 +198,7 @@ "fig = _prepare_figure(\"figure\", figsize=(8.0, 4.5))\n", "ax = fig.subplots(1, 1)\n", "_plot_decoded_ci(ax, time, x_u, W_u, stim.data[:, 0], f\"Decoded stimulus using {numRealizations} cells\")\n", - "__tracker.finalize()" + "__tracker.finalize()\n" ] } ], @@ -211,4 +215,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/DecodingExampleWithHist.ipynb b/notebooks/DecodingExampleWithHist.ipynb index eb2ecffc..32e9549d 100644 --- a/notebooks/DecodingExampleWithHist.ipynb +++ b/notebooks/DecodingExampleWithHist.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "b9af2115", + "id": "a2b60473", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `DecodingExampleWithHist.mlx`\n", - "- Fidelity status: `exact`\n", - "- Remaining justified differences: Mirrors the MATLAB history-aware decoding workflow. Only inherent stochastic trajectories and figure styling differ under Python execution." + "- Fidelity status: `high_fidelity`\n", + "- Remaining justified differences: The notebook now mirrors the MATLAB history-aware decoding workflow closely; exact stochastic trajectories and figure styling still vary slightly under Python execution.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "a847f096", + "id": "0a4bc8bf", "metadata": {}, "outputs": [], "source": [ @@ -42,12 +42,14 @@ "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", "__tracker = FigureTracker(topic=\"DecodingExampleWithHist\", output_root=OUTPUT_ROOT, expected_count=2)\n", "\n", + "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", " fig.clear()\n", " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _plot_raster(ax, spike_coll):\n", " for row in range(1, spike_coll.numSpikeTrains + 1):\n", " train = spike_coll.getNST(row)\n", @@ -57,6 +59,7 @@ " ax.set_ylabel(\"Neuron\")\n", " ax.set_ylim(0.5, spike_coll.numSpikeTrains + 0.5)\n", "\n", + "\n", "def _plot_decoded_ci(ax, time, decoded, cov, stim, title):\n", " center = np.asarray(decoded, dtype=float).reshape(-1)\n", " spread = np.asarray(cov, dtype=float).reshape(-1)\n", @@ -71,6 +74,7 @@ " ax.set_xlabel(\"time (s)\")\n", " ax.legend(loc=\"upper right\", frameon=False, fontsize=8)\n", "\n", + "\n", "def _simulate_history_spike_train(time, stim_data, baseline, hist_coeffs, window_times):\n", " spikes = []\n", " for idx in range(1, len(time)):\n", @@ -89,14 +93,15 @@ " spikes.append(t)\n", " return np.asarray(spikes, dtype=float)\n", "\n", + "\n", "# SECTION 0: 1-D Stimulus Decode with History Effect\n", - "# We simulate neurons with refractory-history effects and compare point-process decoding with and without the correct history terms." + "# We simulate neurons with refractory-history effects and compare point-process decoding with and without the correct history terms.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "b9bfa418", + "id": "7a31db6f", "metadata": {}, "outputs": [], "source": [ @@ -118,7 +123,7 @@ "trains = []\n", "for idx in range(numRealizations):\n", " spikes = _simulate_history_spike_train(time, stimData, b0, histCoeffs, windowTimes)\n", - " trains.append(nspikeTrain(spikes, str(idx + 1), delta, 0.0, Tmax, makePlots=-1))\n", + " trains.append(nspikeTrain(spikes, str(idx + 1), 1.0 / delta, 0.0, Tmax, makePlots=-1))\n", "sC = nstColl(trains)\n", "\n", "fig = _prepare_figure(\"figure\", figsize=(8.0, 5.5))\n", @@ -152,7 +157,7 @@ "axs = fig.subplots(2, 1, sharex=True)\n", "_plot_decoded_ci(axs[0], time, x_u, W_u, stim.data[:, 0], f\"Decoded stimulus with history using {numRealizations} cells\")\n", "_plot_decoded_ci(axs[1], time, x_u_no_hist, W_u_no_hist, stim.data[:, 0], f\"Decoded stimulus without history using {numRealizations} cells\")\n", - "__tracker.finalize()" + "__tracker.finalize()\n" ] } ], @@ -169,4 +174,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/ExplicitStimulusWhiskerData.ipynb b/notebooks/ExplicitStimulusWhiskerData.ipynb index 06b995a6..4e4996d0 100644 --- a/notebooks/ExplicitStimulusWhiskerData.ipynb +++ b/notebooks/ExplicitStimulusWhiskerData.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "199165d0", + "id": "30c2134a", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `ExplicitStimulusWhiskerData.mlx`\n", - "- Fidelity status: `exact`\n", - "- Remaining justified differences: Reproduces the dataset-backed lag search, stimulus-effect, and history-effect workflow with real figures. Only inherent GLM solver numerics and plotting defaults differ." + "- Fidelity status: `high_fidelity`\n", + "- Remaining justified differences: The notebook now reproduces the dataset-backed lag search, stimulus-effect, and history-effect workflow with real figures; exact KS traces and coefficient values still vary modestly from MATLAB because the Python GLM backend and plotting defaults are different.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "e102e8bb", + "id": "5f00236d", "metadata": {}, "outputs": [], "source": [ @@ -42,7 +42,8 @@ "np.random.seed(0)\n", "DATA_DIR = notebook_example_data_dir(allow_synthetic=True)\n", "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", - "__tracker = FigureTracker(topic='ExplicitStimulusWhiskerData', output_root=OUTPUT_ROOT, expected_count=10)\n", + "__tracker = FigureTracker(topic='ExplicitStimulusWhiskerData', output_root=OUTPUT_ROOT, expected_count=9)\n", + "\n", "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", @@ -50,6 +51,7 @@ " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _plot_spike_indicator(ax, time_s, spike_indicator):\n", " spike_times = np.asarray(time_s, dtype=float)[np.asarray(spike_indicator, dtype=float) > 0.5]\n", " if spike_times.size:\n", @@ -57,11 +59,12 @@ " ax.set_ylim(0.0, 1.0)\n", " ax.set_ylabel(\"spikes\")\n", "\n", + "\n", "def _plot_ks(ax, ideal, empirical, ci, *, label, color):\n", " ideal_arr = np.asarray(ideal, dtype=float)\n", " empirical_arr = np.asarray(empirical, dtype=float)\n", " ci_arr = np.asarray(ci, dtype=float)\n", - " ax.plot(ideal_arr, ideal_arr, color=\"0.2\", linewidth=1.0, linestyle=\"--\", label=\"45\u00b0 line\")\n", + " ax.plot(ideal_arr, ideal_arr, color=\"0.2\", linewidth=1.0, linestyle=\"--\", label=\"45° line\")\n", " ax.plot(ideal_arr, empirical_arr, color=color, linewidth=1.5, label=label)\n", " ax.fill_between(\n", " ideal_arr,\n", @@ -74,13 +77,13 @@ " ax.set_xlabel(\"Theoretical quantiles\")\n", " ax.set_ylabel(\"Empirical quantiles\")\n", " ax.set_xlim(0.0, 1.0)\n", - " ax.set_ylim(0.0, 1.0)" + " ax.set_ylim(0.0, 1.0)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "45734023", + "id": "ed25adb7", "metadata": {}, "outputs": [], "source": [ @@ -97,13 +100,13 @@ " \"peak_lag_ms\": round(float(summary[\"peak_lag_seconds\"]) * 1000.0, 1),\n", " \"best_history_window_bins\": best_history_window,\n", " }\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "0e41ab95", + "id": "c780cff1", "metadata": {}, "outputs": [], "source": [ @@ -125,13 +128,13 @@ "axs[1].plot(payload[\"time_s\"], payload[\"velocity\"], color=\"tab:orange\", linewidth=1.2)\n", "axs[1].set_title(\"Stimulus derivative\")\n", "axs[1].set_ylabel(\"d(stimulus)/dt\")\n", - "axs[1].set_xlabel(\"time (s)\")" + "axs[1].set_xlabel(\"time (s)\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "6c2191fc", + "id": "6768c8ae", "metadata": {}, "outputs": [], "source": [ @@ -140,13 +143,13 @@ "ax = fig.subplots(1, 1)\n", "_plot_ks(ax, payload[\"ks_ideal\"], payload[\"ks_const_empirical\"], payload[\"ks_ci\"], label=\"Baseline model\", color=\"tab:blue\")\n", "ax.set_title(\"Baseline model KS plot\")\n", - "ax.legend(loc=\"lower right\", frameon=False, fontsize=8)" + "ax.legend(loc=\"lower right\", frameon=False, fontsize=8)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "da4b0be4", + "id": "e2145f7b", "metadata": {}, "outputs": [], "source": [ @@ -161,13 +164,13 @@ "ax.scatter([lags_ms[peak_idx]], [xcorr_vals[peak_idx]], color=\"tab:red\", zorder=3)\n", "ax.set_title(\"Cross-covariance used to identify the stimulus lag\")\n", "ax.set_xlabel(\"lag (ms)\")\n", - "ax.set_ylabel(\"cross-covariance\")" + "ax.set_ylabel(\"cross-covariance\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "334952df", + "id": "ae00cde0", "metadata": {}, "outputs": [], "source": [ @@ -189,13 +192,13 @@ "_plot_ks(ax, payload[\"ks_ideal\"], payload[\"ks_const_empirical\"], payload[\"ks_ci\"], label=\"Baseline\", color=\"tab:blue\")\n", "ax.plot(np.asarray(payload[\"ks_ideal\"], dtype=float), np.asarray(payload[\"ks_stim_empirical\"], dtype=float), color=\"tab:orange\", linewidth=1.5, label=\"Baseline+Stimulus\")\n", "ax.set_title(\"Baseline vs stimulus-augmented model\")\n", - "ax.legend(loc=\"lower right\", frameon=False, fontsize=8)" + "ax.legend(loc=\"lower right\", frameon=False, fontsize=8)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "eb6dc162", + "id": "8a22e891", "metadata": {}, "outputs": [], "source": [ @@ -209,10 +212,10 @@ "axs[0].set_title(\"History-window scan\")\n", "axs[1].plot(history_windows, payload[\"delta_aic\"], marker=\"o\", color=\"tab:green\", linewidth=1.2)\n", "axs[1].scatter([history_windows[best_history_idx]], [payload[\"delta_aic\"][best_history_idx]], color=\"tab:red\", zorder=3)\n", - "axs[1].set_ylabel(\"\u0394AIC\")\n", + "axs[1].set_ylabel(\"ΔAIC\")\n", "axs[2].plot(history_windows, payload[\"delta_bic\"], marker=\"o\", color=\"tab:brown\", linewidth=1.2)\n", "axs[2].scatter([history_windows[best_history_idx]], [payload[\"delta_bic\"][best_history_idx]], color=\"tab:red\", zorder=3)\n", - "axs[2].set_ylabel(\"\u0394BIC\")\n", + "axs[2].set_ylabel(\"ΔBIC\")\n", "axs[2].set_xlabel(\"history window count\")\n", "\n", "fig = _prepare_figure(\"plot(x,dBIC,'.')\", figsize=(8.0, 4.5))\n", @@ -221,13 +224,13 @@ "ax.axvline(history_windows[best_history_idx], color=\"tab:red\", linestyle=\"--\", linewidth=1.0)\n", "ax.set_title(\"BIC improvement across history-window choices\")\n", "ax.set_xlabel(\"history window count\")\n", - "ax.set_ylabel(\"\u0394BIC relative to first history model\")" + "ax.set_ylabel(\"ΔBIC relative to first history model\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "09de73d5", + "id": "5b152fd9", "metadata": {}, "outputs": [], "source": [ @@ -270,7 +273,7 @@ "ax.plot(np.asarray(payload[\"ks_ideal\"], dtype=float), np.asarray(payload[\"ks_hist_empirical\"], dtype=float), color=\"tab:green\", linewidth=1.5, label=\"Baseline+Stimulus+History\")\n", "ax.set_title(\"Final KS comparison across the three models\")\n", "ax.legend(loc=\"lower right\", frameon=False, fontsize=8)\n", - "__tracker.finalize()" + "__tracker.finalize()\n" ] } ], @@ -280,7 +283,7 @@ }, "nstat": { "expected_figures": 9, - "run_group": "full", + "run_group": "smoke", "style": "python-example", "topic": "ExplicitStimulusWhiskerData" } diff --git a/notebooks/HippocampalPlaceCellExample.ipynb b/notebooks/HippocampalPlaceCellExample.ipynb index d9ef8901..b3ad9d87 100644 --- a/notebooks/HippocampalPlaceCellExample.ipynb +++ b/notebooks/HippocampalPlaceCellExample.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "74d6cdfa", + "id": "f563f4e9", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `HippocampalPlaceCellExample.mlx`\n", - "- Fidelity status: `exact`\n", - "- Remaining justified differences: Reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with the same normalized 10-term Zernike basis used by MATLAB. Only inherent GLM solver numerics and surface styling differ." + "- Fidelity status: `high_fidelity`\n", + "- Remaining justified differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with the same normalized 10-term Zernike basis used by MATLAB; exact AIC/BIC values and surface styling still vary modestly because the Python GLM solver and plotting backend are not byte-identical to MATLAB.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "15ef53db", + "id": "e0302c60", "metadata": {}, "outputs": [], "source": [ @@ -42,7 +42,8 @@ "np.random.seed(0)\n", "DATA_DIR = notebook_example_data_dir(allow_synthetic=True)\n", "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", - "__tracker = FigureTracker(topic='HippocampalPlaceCellExample', output_root=OUTPUT_ROOT, expected_count=10)\n", + "__tracker = FigureTracker(topic='HippocampalPlaceCellExample', output_root=OUTPUT_ROOT, expected_count=11)\n", + "\n", "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", @@ -50,6 +51,7 @@ " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _interp_spike_positions(time_s, x_pos, y_pos, spike_times):\n", " spike_times = np.asarray(spike_times, dtype=float)\n", " return (\n", @@ -57,6 +59,7 @@ " np.interp(spike_times, np.asarray(time_s, dtype=float), np.asarray(y_pos, dtype=float)),\n", " )\n", "\n", + "\n", "def _plot_field_grid(fig, animal_key, field_key, title):\n", " animal = payload[animal_key]\n", " grid_x = np.asarray(animal[\"grid_x\"], dtype=float)\n", @@ -76,13 +79,13 @@ " ax.set_xticks([])\n", " ax.set_yticks([])\n", " fig.suptitle(title)\n", - " fig.colorbar(image, ax=axs.ravel().tolist(), shrink=0.78)" + " fig.colorbar(image, ax=axs.ravel().tolist(), shrink=0.78)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "30094bfc", + "id": "270e5d9d", "metadata": {}, "outputs": [], "source": [ @@ -96,13 +99,13 @@ " \"mean_delta_aic\": round(float(summary[\"mean_delta_aic_gaussian_minus_zernike\"]), 3),\n", " \"mean_delta_bic\": round(float(summary[\"mean_delta_bic_gaussian_minus_zernike\"]), 3),\n", " }\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "33671840", + "id": "8079dbfd", "metadata": {}, "outputs": [], "source": [ @@ -116,64 +119,142 @@ "ax.set_title(f\"Animal 1, Cell {int(mesh['cell_index']) + 1}\")\n", "ax.set_xlabel(\"x\")\n", "ax.set_ylabel(\"y\")\n", - "ax.set_aspect(\"equal\", adjustable=\"box\")" + "ax.set_aspect(\"equal\", adjustable=\"box\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "081e6179", + "id": "cd6e5704", "metadata": {}, "outputs": [], "source": [ "# SECTION 2: Analyze All Cells\n", - "fig = _prepare_figure(\"Summary.plotSummary\", figsize=(12.0, 4.5))\n", - "axs = fig.subplots(1, 2)\n", + "fig = _prepare_figure(\"Summary.plotSummary\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", "animal1 = payload[\"animal1\"]\n", "labels = [f\"Cell {int(idx) + 1}\" for idx in np.asarray(animal1[\"selected_indices\"], dtype=int)]\n", - "axs[0].bar(np.arange(len(labels)), animal1[\"delta_aic\"], color=\"tab:purple\")\n", - "axs[0].axhline(0.0, color=\"0.2\", linewidth=1.0)\n", - "axs[0].set_xticks(np.arange(len(labels)), labels, rotation=20)\n", - "axs[0].set_ylabel(\"Gaussian - Zernike AIC\")\n", - "axs[0].set_title(\"Animal 1 AIC\")\n", - "axs[1].bar(np.arange(len(labels)), animal1[\"delta_bic\"], color=\"tab:green\")\n", - "axs[1].axhline(0.0, color=\"0.2\", linewidth=1.0)\n", - "axs[1].set_xticks(np.arange(len(labels)), labels, rotation=20)\n", - "axs[1].set_ylabel(\"Gaussian - Zernike BIC\")\n", - "axs[1].set_title(\"Animal 1 BIC\")\n" + "ax.bar(np.arange(len(labels)), animal1[\"delta_aic\"], color=\"tab:purple\")\n", + "ax.axhline(0.0, color=\"0.2\", linewidth=1.0)\n", + "ax.set_xticks(np.arange(len(labels)), labels, rotation=20)\n", + "ax.set_ylabel(\"Gaussian - Zernike AIC\")\n", + "ax.set_title(\"Animal 1 model comparison\")\n", + "\n", + "fig = _prepare_figure(\"Summary.plotSummary\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "ax.bar(np.arange(len(labels)), animal1[\"delta_bic\"], color=\"tab:green\")\n", + "ax.axhline(0.0, color=\"0.2\", linewidth=1.0)\n", + "ax.set_xticks(np.arange(len(labels)), labels, rotation=20)\n", + "ax.set_ylabel(\"Gaussian - Zernike BIC\")\n", + "ax.set_title(\"Animal 1 model comparison\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "aa269c6b", + "id": "2487ae0f", "metadata": {}, "outputs": [], "source": [ "# SECTION 3: View Summary Statistics\n", - "fig = _prepare_figure(\"Summary.plotSummary\", figsize=(12.0, 4.5))\n", - "axs = fig.subplots(1, 2)\n", + "fig = _prepare_figure(\"Summary.plotSummary\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", "animal2 = payload[\"animal2\"]\n", "labels = [f\"Cell {int(idx) + 1}\" for idx in np.asarray(animal2[\"selected_indices\"], dtype=int)]\n", - "axs[0].bar(np.arange(len(labels)), animal2[\"delta_aic\"], color=\"tab:purple\")\n", - "axs[0].axhline(0.0, color=\"0.2\", linewidth=1.0)\n", - "axs[0].set_xticks(np.arange(len(labels)), labels, rotation=20)\n", - "axs[0].set_ylabel(\"Gaussian - Zernike AIC\")\n", - "axs[0].set_title(\"Animal 2 AIC\")\n", - "axs[1].bar(np.arange(len(labels)), animal2[\"delta_bic\"], color=\"tab:green\")\n", - "axs[1].axhline(0.0, color=\"0.2\", linewidth=1.0)\n", - "axs[1].set_xticks(np.arange(len(labels)), labels, rotation=20)\n", - "axs[1].set_ylabel(\"Gaussian - Zernike BIC\")\n", - "axs[1].set_title(\"Animal 2 BIC\")\n" + "ax.bar(np.arange(len(labels)), animal2[\"delta_aic\"], color=\"tab:purple\")\n", + "ax.axhline(0.0, color=\"0.2\", linewidth=1.0)\n", + "ax.set_xticks(np.arange(len(labels)), labels, rotation=20)\n", + "ax.set_ylabel(\"Gaussian - Zernike AIC\")\n", + "ax.set_title(\"Animal 2 model comparison\")\n", + "\n", + "fig = _prepare_figure(\"Summary.plotSummary\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "ax.bar(np.arange(len(labels)), animal2[\"delta_bic\"], color=\"tab:green\")\n", + "ax.axhline(0.0, color=\"0.2\", linewidth=1.0)\n", + "ax.set_xticks(np.arange(len(labels)), labels, rotation=20)\n", + "ax.set_ylabel(\"Gaussian - Zernike BIC\")\n", + "ax.set_title(\"Animal 2 model comparison\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "26aafec5", + "id": "f92a805b", "metadata": {}, "outputs": [], - "source": "# SECTION 4: Visualize the results\nfig = _prepare_figure(\"h4=figure(4)\", figsize=(8.5, 8.0))\n_plot_field_grid(fig, \"animal1\", \"gaussian_fields\", \"Gaussian place fields - Animal 1\")\n\nfig = _prepare_figure(\"h5=figure(5)\", figsize=(8.5, 8.0))\n_plot_field_grid(fig, \"animal1\", \"zernike_fields\", \"Zernike place fields - Animal 1\")\n\nfig = _prepare_figure(\"h6=figure(6)\", figsize=(8.5, 8.0))\n_plot_field_grid(fig, \"animal2\", \"gaussian_fields\", \"Gaussian place fields - Animal 2\")\n\nfig = _prepare_figure(\"h7=figure(7)\", figsize=(8.5, 8.0))\n_plot_field_grid(fig, \"animal2\", \"zernike_fields\", \"Zernike place fields - Animal 2\")\n\nfig = _prepare_figure(\"figure(8)\", figsize=(7.0, 5.5))\nax = fig.subplots(1, 1)\nax.imshow(\n mesh[\"gaussian_field\"],\n origin=\"lower\",\n extent=[float(np.min(mesh[\"grid_x\"])), float(np.max(mesh[\"grid_x\"])), float(np.min(mesh[\"grid_y\"])), float(np.max(mesh[\"grid_y\"]))],\n aspect=\"equal\",\n cmap=\"viridis\",\n)\nax.plot(mesh[\"x_pos\"], mesh[\"y_pos\"], color=\"white\", linewidth=0.5, alpha=0.35)\nax.scatter(spike_x, spike_y, s=8, color=\"tab:red\", alpha=0.7)\nax.set_title(f\"Gaussian receptive field - Cell {int(mesh['cell_index']) + 1}\")\nax.set_xlabel(\"x\")\nax.set_ylabel(\"y\")\n\nfig = _prepare_figure(\"figure(9)\", figsize=(7.0, 5.5))\nax = fig.subplots(1, 1)\nax.imshow(\n mesh[\"zernike_field\"],\n origin=\"lower\",\n extent=[float(np.min(mesh[\"grid_x\"])), float(np.max(mesh[\"grid_x\"])), float(np.min(mesh[\"grid_y\"])), float(np.max(mesh[\"grid_y\"]))],\n aspect=\"equal\",\n cmap=\"viridis\",\n)\nax.plot(mesh[\"x_pos\"], mesh[\"y_pos\"], color=\"white\", linewidth=0.5, alpha=0.35)\nax.scatter(spike_x, spike_y, s=8, color=\"tab:red\", alpha=0.7)\nax.set_title(f\"Zernike receptive field - Cell {int(mesh['cell_index']) + 1}\")\nax.set_xlabel(\"x\")\nax.set_ylabel(\"y\")\n\n# figure(9) overlay — matches MATLAB hold-on composite (published as _10.png)\nfig = _prepare_figure(\"figure(9) overlay\", figsize=(10.0, 5.0))\naxs = fig.subplots(1, 2)\next = [float(np.min(mesh[\"grid_x\"])), float(np.max(mesh[\"grid_x\"])), float(np.min(mesh[\"grid_y\"])), float(np.max(mesh[\"grid_y\"]))]\nfor ax, field, label in zip(axs, [mesh[\"gaussian_field\"], mesh[\"zernike_field\"]], [\"Gaussian\", \"Zernike\"]):\n ax.imshow(field, origin=\"lower\", extent=ext, aspect=\"equal\", cmap=\"viridis\")\n ax.plot(mesh[\"x_pos\"], mesh[\"y_pos\"], color=\"white\", linewidth=0.5, alpha=0.35)\n ax.scatter(spike_x, spike_y, s=8, color=\"tab:red\", alpha=0.7)\n ax.set_title(f\"{label} - Cell {int(mesh['cell_index']) + 1}\")\n ax.set_xlabel(\"x\")\n ax.set_ylabel(\"y\")\n\n__tracker.finalize()" + "source": [ + "# SECTION 4: Visualize the results\n", + "fig = _prepare_figure(\"h4=figure(4)\", figsize=(8.5, 8.0))\n", + "_plot_field_grid(fig, \"animal1\", \"gaussian_fields\", \"Gaussian place fields - Animal 1\")\n", + "\n", + "fig = _prepare_figure(\"h5=figure(5)\", figsize=(8.5, 8.0))\n", + "_plot_field_grid(fig, \"animal1\", \"zernike_fields\", \"Zernike place fields - Animal 1\")\n", + "\n", + "fig = _prepare_figure(\"h6=figure(6)\", figsize=(8.5, 8.0))\n", + "_plot_field_grid(fig, \"animal2\", \"gaussian_fields\", \"Gaussian place fields - Animal 2\")\n", + "\n", + "fig = _prepare_figure(\"h7=figure(7)\", figsize=(8.5, 8.0))\n", + "_plot_field_grid(fig, \"animal2\", \"zernike_fields\", \"Zernike place fields - Animal 2\")\n", + "\n", + "fig = _prepare_figure(\"figure(8)\", figsize=(7.0, 5.5))\n", + "ax = fig.subplots(1, 1)\n", + "ax.imshow(\n", + " mesh[\"gaussian_field\"],\n", + " origin=\"lower\",\n", + " extent=[float(np.min(mesh[\"grid_x\"])), float(np.max(mesh[\"grid_x\"])), float(np.min(mesh[\"grid_y\"])), float(np.max(mesh[\"grid_y\"]))],\n", + " aspect=\"equal\",\n", + " cmap=\"viridis\",\n", + ")\n", + "ax.plot(mesh[\"x_pos\"], mesh[\"y_pos\"], color=\"white\", linewidth=0.5, alpha=0.35)\n", + "ax.scatter(spike_x, spike_y, s=8, color=\"tab:red\", alpha=0.7)\n", + "ax.set_title(f\"Gaussian receptive field - Cell {int(mesh['cell_index']) + 1}\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "\n", + "fig = _prepare_figure(\"figure(9)\", figsize=(7.0, 5.5))\n", + "ax = fig.subplots(1, 1)\n", + "ax.imshow(\n", + " mesh[\"zernike_field\"],\n", + " origin=\"lower\",\n", + " extent=[float(np.min(mesh[\"grid_x\"])), float(np.max(mesh[\"grid_x\"])), float(np.min(mesh[\"grid_y\"])), float(np.max(mesh[\"grid_y\"]))],\n", + " aspect=\"equal\",\n", + " cmap=\"viridis\",\n", + ")\n", + "ax.plot(mesh[\"x_pos\"], mesh[\"y_pos\"], color=\"white\", linewidth=0.5, alpha=0.35)\n", + "ax.scatter(spike_x, spike_y, s=8, color=\"tab:red\", alpha=0.7)\n", + "ax.set_title(f\"Zernike receptive field - Cell {int(mesh['cell_index']) + 1}\")\n", + "ax.set_xlabel(\"x\")\n", + "ax.set_ylabel(\"y\")\n", + "\n", + "fig = _prepare_figure(\"figure(10)\", figsize=(9.0, 4.5))\n", + "axs = fig.subplots(1, 2)\n", + "axs[0].hist(np.concatenate([payload[\"animal1\"][\"delta_aic\"], payload[\"animal2\"][\"delta_aic\"]]), bins=8, color=\"tab:purple\", alpha=0.8)\n", + "axs[0].axvline(0.0, color=\"0.2\", linewidth=1.0)\n", + "axs[0].set_title(\"Distribution of ΔAIC\")\n", + "axs[1].hist(np.concatenate([payload[\"animal1\"][\"delta_bic\"], payload[\"animal2\"][\"delta_bic\"]]), bins=8, color=\"tab:green\", alpha=0.8)\n", + "axs[1].axvline(0.0, color=\"0.2\", linewidth=1.0)\n", + "axs[1].set_title(\"Distribution of ΔBIC\")\n", + "\n", + "fig = _prepare_figure(\"figure(11)\", figsize=(6.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "ax.axis(\"off\")\n", + "ax.text(\n", + " 0.0,\n", + " 0.95,\n", + " \"\\n\".join(\n", + " [\n", + " f\"Cells analyzed: {int(summary['num_cells_fit'])}\",\n", + " f\"Mean Gaussian-Zernike ΔAIC: {summary['mean_delta_aic_gaussian_minus_zernike']:.2f}\",\n", + " f\"Mean Gaussian-Zernike ΔBIC: {summary['mean_delta_bic_gaussian_minus_zernike']:.2f}\",\n", + " \"Negative values favor the Zernike model.\",\n", + " ]\n", + " ),\n", + " va=\"top\",\n", + " family=\"monospace\",\n", + " fontsize=10,\n", + ")\n", + "__tracker.finalize()\n" + ] } ], "metadata": { @@ -182,7 +263,7 @@ }, "nstat": { "expected_figures": 11, - "run_group": "full", + "run_group": "smoke", "style": "python-example", "topic": "HippocampalPlaceCellExample" } diff --git a/notebooks/HybridFilterExample.ipynb b/notebooks/HybridFilterExample.ipynb index 575b57cb..2d84d178 100644 --- a/notebooks/HybridFilterExample.ipynb +++ b/notebooks/HybridFilterExample.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "d60b95a7", + "id": "031f2d59", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `HybridFilterExample.mlx`\n", - "- Fidelity status: `exact`\n", - "- Remaining justified differences: Reproduces the hybrid-filter simulation, single-run decoding, and averaged summary figures with real outputs. Only inherent stochastic trajectories and Python hybrid-filter implementation details differ.\n" + "- Fidelity status: `high_fidelity`\n", + "- Remaining justified differences: The notebook now reproduces the hybrid-filter simulation, single-run decoding, and averaged summary figures with real outputs; the Python port still uses the current hybrid-filter implementation instead of every MATLAB-specific reporting branch.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "d22dd216", + "id": "90518773", "metadata": {}, "outputs": [], "source": [ @@ -42,12 +42,14 @@ "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", "__tracker = FigureTracker(topic='HybridFilterExample', output_root=OUTPUT_ROOT, expected_count=3)\n", "\n", + "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", " fig.clear()\n", " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _plot_raster(ax, time_s, spikes, *, max_cells=18):\n", " n_cells = min(int(spikes.shape[1]), max_cells)\n", " for row in range(n_cells):\n", @@ -55,13 +57,13 @@ " if spike_times.size:\n", " ax.vlines(spike_times, row + 0.6, row + 1.4, color=\"k\", linewidth=0.35)\n", " ax.set_ylim(0.5, n_cells + 0.5)\n", - " ax.set_ylabel(\"cell\")" + " ax.set_ylabel(\"cell\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "cacac4c3", + "id": "1bccc86c", "metadata": {}, "outputs": [], "source": [ @@ -79,35 +81,35 @@ " \"num_cells\": int(summary[\"num_cells\"]),\n", " \"state_accuracy\": round(float(summary[\"state_accuracy\"]), 3),\n", " }\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "b031dd85", + "id": "5f1f9aee", "metadata": {}, "outputs": [], "source": [ "# SECTION 1: Problem Statement\n", - "# We infer both a discrete movement state and a continuous reach trajectory from point-process observations." + "# We infer both a discrete movement state and a continuous reach trajectory from point-process observations.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "fd11f1d4", + "id": "93db6f6c", "metadata": {}, "outputs": [], "source": [ "# SECTION 2: Hybrid state-space setup\n", - "# The Python port keeps the same two-state problem structure as MATLAB: a low-motion state and a movement state." + "# The Python port keeps the same two-state problem structure as MATLAB: a low-motion state and a movement state.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "fd690f04", + "id": "215fd0f4", "metadata": {}, "outputs": [], "source": [ @@ -153,24 +155,24 @@ " va=\"top\",\n", " family=\"monospace\",\n", " fontsize=9,\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "f9e4eb9d", + "id": "c7531b74", "metadata": {}, "outputs": [], "source": [ "# SECTION 4: Simulate Neural Firing\n", - "# The simulated spike population depends on the latent state and the movement dynamics." + "# The simulated spike population depends on the latent state and the movement dynamics.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "57220fb5", + "id": "90456226", "metadata": {}, "outputs": [], "source": [ @@ -229,7 +231,7 @@ " color=[\"tab:blue\", \"tab:orange\"],\n", ")\n", "axs[1, 1].set_title(\"Single-run decoding RMSE\")\n", - "__tracker.finalize()" + "__tracker.finalize()\n" ] } ], @@ -239,11 +241,11 @@ }, "nstat": { "expected_figures": 3, - "run_group": "full", + "run_group": "smoke", "style": "python-example", "topic": "HybridFilterExample" } }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/PointProcessSimulation.slxc b/notebooks/PointProcessSimulation.slxc new file mode 100644 index 00000000..83ec2739 Binary files /dev/null and b/notebooks/PointProcessSimulation.slxc differ diff --git a/notebooks/StimulusDecode2D.ipynb b/notebooks/StimulusDecode2D.ipynb index e9b01a4c..2f971d4b 100644 --- a/notebooks/StimulusDecode2D.ipynb +++ b/notebooks/StimulusDecode2D.ipynb @@ -2,14 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "4599bf89", + "id": "c20fb45a", "metadata": {}, - "source": "\n## MATLAB Parity Note\n- Source MATLAB helpfile: `StimulusDecode2D.mlx`\n- Fidelity status: `exact`\n- Remaining justified differences: Follows the MATLAB nonlinear-CIF decoding workflow with `DecodingAlgorithms.PPDecodeFilter` and all 4 published figures. Only inherent Python symbolic/numeric stack and random streams differ." + "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.\n" + ] }, { "cell_type": "code", "execution_count": null, - "id": "6f7a27a5", + "id": "1e9ccf67", "metadata": {}, "outputs": [], "source": [ @@ -34,7 +40,8 @@ "\n", "np.random.seed(0)\n", "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", - "__tracker = FigureTracker(topic='StimulusDecode2D', output_root=OUTPUT_ROOT, expected_count=4)\n", + "__tracker = FigureTracker(topic='StimulusDecode2D', output_root=OUTPUT_ROOT, expected_count=6)\n", + "\n", "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", @@ -42,11 +49,13 @@ " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _subplot_grid(count):\n", " rows = max(int(np.floor(np.sqrt(count))), 1)\n", " cols = int(np.ceil(count / rows))\n", " return rows, cols\n", "\n", + "\n", "def _simulate_decode(seed=0, *, num_realizations=80, delta=0.001, tmax=1.0):\n", " rng = np.random.default_rng(seed)\n", " time = np.arange(0.0, tmax + delta, delta)\n", @@ -140,6 +149,7 @@ " \"num_cells\": num_realizations,\n", " }\n", "\n", + "\n", "def _plot_raster(ax, time_s, spikes, *, max_cells=20):\n", " n_cells = min(int(spikes.shape[1]), max_cells)\n", " for row in range(n_cells):\n", @@ -147,17 +157,17 @@ " if spike_times.size:\n", " ax.vlines(spike_times, row + 0.6, row + 1.4, color=\"k\", linewidth=0.35)\n", " ax.set_ylim(0.5, n_cells + 0.5)\n", - " ax.set_ylabel(\"cell\")" + " ax.set_ylabel(\"cell\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "8809d377", + "id": "feb17f45", "metadata": {}, "outputs": [], "source": [ - "# SECTION 1: 2-D Stimulus Decode\n", + "# SECTION 0: 2-D Stimulus Decode\n", "# This notebook follows the MATLAB 2-D decoding workflow with simulated spatial receptive fields.\n", "plt.close(\"all\")\n", "payload = _simulate_decode()\n", @@ -168,17 +178,17 @@ " \"decode_rmse\": round(float(payload[\"decode_rmse\"]), 4),\n", " \"fallback_error\": payload[\"decode_error\"] or \"\",\n", " }\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "d87abf0b", + "id": "ad9b9210", "metadata": {}, "outputs": [], "source": [ - "# SECTION 2: Generate the random receptive fields to simulate different neurons\n", + "# SECTION 1: Generate the random receptive fields to simulate different neurons\n", "fig = _prepare_figure(\"figure; plot(px,py)\", figsize=(6.0, 6.0))\n", "ax = fig.subplots(1, 1)\n", "ax.plot(payload[\"px\"], payload[\"py\"], color=\"tab:blue\", linewidth=1.5)\n", @@ -214,13 +224,31 @@ " ax.set_xticks([])\n", " ax.set_yticks([])\n", "if image is not None:\n", - " fig.colorbar(image, ax=axs.ravel().tolist(), shrink=0.78)" + " fig.colorbar(image, ax=axs.ravel().tolist(), shrink=0.78)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33ba6aa5", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 2: Visualize the simulated neural activity\n", + "fig = _prepare_figure(\"spikeColl.plot\", figsize=(9.0, 5.0))\n", + "axs = fig.subplots(2, 1, sharex=True)\n", + "_plot_raster(axs[0], payload[\"time_s\"], payload[\"spikes\"])\n", + "axs[0].set_title(\"Population raster\")\n", + "axs[1].plot(payload[\"time_s\"], np.mean(payload[\"spikes\"], axis=1), color=\"tab:green\", linewidth=1.2)\n", + "axs[1].set_title(\"Population firing fraction\")\n", + "axs[1].set_xlabel(\"time (s)\")\n", + "axs[1].set_ylabel(\"mean spike/bin\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "5eddb87b", + "id": "f07dbd8c", "metadata": {}, "outputs": [], "source": [ @@ -235,6 +263,29 @@ "ax.legend(loc=\"best\", frameon=False, fontsize=8)\n", "ax.set_aspect(\"equal\", adjustable=\"box\")\n", "\n", + "fig = _prepare_figure(\"plot(decoded trajectories)\", figsize=(10.0, 5.5))\n", + "axs = fig.subplots(2, 1, sharex=True)\n", + "axs[0].plot(payload[\"time_s\"], payload[\"px\"], color=\"k\", linewidth=1.6, label=\"True x\")\n", + "axs[0].plot(payload[\"time_s\"][: payload[\"decoded_x\"].size], payload[\"decoded_x\"], color=\"tab:blue\", linewidth=1.2, label=\"Decoded x\")\n", + "axs[0].plot(payload[\"time_s\"][: payload[\"predicted_x\"].size], payload[\"predicted_x\"], color=\"tab:orange\", linewidth=1.0, label=\"Predicted x\")\n", + "axs[0].legend(loc=\"best\", frameon=False, fontsize=8)\n", + "axs[0].set_ylabel(\"x\")\n", + "axs[1].plot(payload[\"time_s\"], payload[\"py\"], color=\"k\", linewidth=1.6, label=\"True y\")\n", + "axs[1].plot(payload[\"time_s\"][: payload[\"decoded_y\"].size], payload[\"decoded_y\"], color=\"tab:blue\", linewidth=1.2, label=\"Decoded y\")\n", + "axs[1].plot(payload[\"time_s\"][: payload[\"predicted_y\"].size], payload[\"predicted_y\"], color=\"tab:orange\", linewidth=1.0, label=\"Predicted y\")\n", + "axs[1].set_ylabel(\"y\")\n", + "axs[1].set_xlabel(\"time (s)\")\n", + "\n", + "fig = _prepare_figure(\"decode_rmse\", figsize=(7.0, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "error_decode = np.sqrt((payload[\"decoded_x\"] - payload[\"px\"][: payload[\"decoded_x\"].size]) ** 2 + (payload[\"decoded_y\"] - payload[\"py\"][: payload[\"decoded_y\"].size]) ** 2)\n", + "error_predict = np.sqrt((payload[\"predicted_x\"] - payload[\"px\"][: payload[\"predicted_x\"].size]) ** 2 + (payload[\"predicted_y\"] - payload[\"py\"][: payload[\"predicted_y\"].size]) ** 2)\n", + "ax.plot(payload[\"time_s\"][: error_decode.size], error_decode, color=\"tab:blue\", linewidth=1.2, label=\"Filtered decode\")\n", + "ax.plot(payload[\"time_s\"][: error_predict.size], error_predict, color=\"tab:orange\", linewidth=1.0, label=\"Predicted decode\")\n", + "ax.set_title(f\"Pointwise decoding error (RMSE={payload['decode_rmse']:.3f})\")\n", + "ax.set_xlabel(\"time (s)\")\n", + "ax.set_ylabel(\"Euclidean error\")\n", + "ax.legend(loc=\"best\", frameon=False, fontsize=8)\n", "__tracker.finalize()\n" ] } @@ -245,7 +296,7 @@ }, "nstat": { "expected_figures": 6, - "run_group": "full", + "run_group": "smoke", "style": "python-example", "topic": "StimulusDecode2D" } diff --git a/notebooks/ValidationDataSet.ipynb b/notebooks/ValidationDataSet.ipynb index 766d9311..cb5a92f5 100644 --- a/notebooks/ValidationDataSet.ipynb +++ b/notebooks/ValidationDataSet.ipynb @@ -2,20 +2,20 @@ "cells": [ { "cell_type": "markdown", - "id": "ff1426a4", + "id": "b5c82c5b", "metadata": {}, "source": [ "\n", "## MATLAB Parity Note\n", "- Source MATLAB helpfile: `ValidationDataSet.mlx`\n", - "- Fidelity status: `exact`\n", - "- Remaining justified differences: Reproduces the constant-rate and piecewise-rate validation workflows with real Trial/Analysis objects and figure outputs. CI uses a documented shorter deterministic fast path for stability." + "- Fidelity status: `high_fidelity`\n", + "- Remaining justified differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` objects and figure outputs; local execution uses the MATLAB-scale simulation sizes, while CI switches to a documented shorter deterministic fast path for stability.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "5a19211e", + "id": "e1e838d2", "metadata": {}, "outputs": [], "source": [ @@ -44,12 +44,14 @@ "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", "__tracker = FigureTracker(topic='ValidationDataSet', output_root=OUTPUT_ROOT, expected_count=10)\n", "\n", + "\n", "def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)):\n", " fig = __tracker.new_figure(matlab_line)\n", " fig.clear()\n", " fig.set_size_inches(*figsize)\n", " return fig\n", "\n", + "\n", "def _lambda_columns(fit_result):\n", " time = np.asarray(fit_result.lambda_signal.time, dtype=float)\n", " data = np.asarray(fit_result.lambda_signal.data, dtype=float)\n", @@ -57,8 +59,10 @@ " data = data[:, None]\n", " return time, data\n", "\n", + "\n", "CI_FAST_PATH = os.environ.get(\"CI\", \"\").strip().lower() in {\"1\", \"true\", \"yes\"}\n", "\n", + "\n", "def _simulate_constant_case(seed=0, *, p=0.01, n_samples=None, delta=0.001):\n", " if n_samples is None:\n", " n_samples = 20001 if CI_FAST_PATH else 100001\n", @@ -71,7 +75,7 @@ " for idx in range(2):\n", " spike_mask = rng.random(n_samples) < p\n", " spike_times = time[spike_mask]\n", - " train = nspikeTrain(spike_times, str(idx + 1), delta, 0.0, total_time, makePlots=-1)\n", + " train = nspikeTrain(spike_times, str(idx + 1), 1.0 / delta, 0.0, total_time, makePlots=-1)\n", " trains.append(train)\n", " spike_coll = nstColl(trains)\n", " cov = Covariate(time, np.ones((time.shape[0], 1), dtype=float), \"Baseline\", \"time\", \"s\", \"\", [\"mu\"])\n", @@ -87,6 +91,7 @@ " \"trains\": trains,\n", " }\n", "\n", + "\n", "def _simulate_piecewise_case(seed=1, *, p1=0.001, p2=0.01, n1=None, n2=None, delta=0.001):\n", " if n1 is None:\n", " n1 = 20000 if CI_FAST_PATH else 100000\n", @@ -104,7 +109,7 @@ " spikes1 = t1[:-1][rng.random(n1) < p1]\n", " spikes2 = t2[rng.random(n2) < p2]\n", " spike_times = np.concatenate([spikes1, spikes2])\n", - " train = nspikeTrain(spike_times, str(idx + 1), delta, 0.0, total_time, makePlots=-1)\n", + " train = nspikeTrain(spike_times, str(idx + 1), 1.0 / delta, 0.0, total_time, makePlots=-1)\n", " trains.append(train)\n", " time = np.concatenate([t1[:-1], t2])\n", " cov_data = np.column_stack(\n", @@ -134,6 +139,7 @@ " \"trains\": trains,\n", " }\n", "\n", + "\n", "def _plot_isi_hist(ax, train, lambda_hz, *, title):\n", " isi = np.asarray(train.getISIs(), dtype=float)\n", " if isi.size:\n", @@ -142,13 +148,13 @@ " ax.plot(x, lambda_hz * np.exp(-lambda_hz * x), color=\"tab:red\", linewidth=1.5)\n", " ax.set_title(title)\n", " ax.set_xlabel(\"ISI (s)\")\n", - " ax.set_ylabel(\"density\")" + " ax.set_ylabel(\"density\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "de07a751", + "id": "a0dd7789", "metadata": {}, "outputs": [], "source": [ @@ -164,36 +170,36 @@ " \"piecewise_lambda1_hz\": round(float(piecewise_case[\"lambda1_hz\"]), 4),\n", " \"piecewise_lambda2_hz\": round(float(piecewise_case[\"lambda2_hz\"]), 4),\n", " }\n", - ")" + ")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "4c326ba4", + "id": "d247fcbf", "metadata": {}, "outputs": [], "source": [ "# SECTION 1: Case #1: Constant Rate Poisson Process\n", - "# First we verify that the analysis recovers a constant Poisson rate from simulated spike trains." + "# First we verify that the analysis recovers a constant Poisson rate from simulated spike trains.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "0348ff8b", + "id": "14ae875f", "metadata": {}, "outputs": [], "source": [ "# SECTION 2: Generate constant-rate neural firing activity\n", "constant_time = np.asarray(constant_case[\"time_s\"], dtype=float)\n", - "constant_trains = list(constant_case[\"trains\"])" + "constant_trains = list(constant_case[\"trains\"])\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "127d7e94", + "id": "50be0fae", "metadata": {}, "outputs": [], "source": [ @@ -201,19 +207,19 @@ "fig = _prepare_figure(\"nst{1}.plotISIHistogram\", figsize=(10.0, 4.0))\n", "axs = fig.subplots(1, 2)\n", "_plot_isi_hist(axs[0], constant_trains[0], constant_case[\"lambda_hz\"], title=\"Neuron 1 ISI histogram\")\n", - "_plot_isi_hist(axs[1], constant_trains[1], constant_case[\"lambda_hz\"], title=\"Neuron 2 ISI histogram\")" + "_plot_isi_hist(axs[1], constant_trains[1], constant_case[\"lambda_hz\"], title=\"Neuron 2 ISI histogram\")\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "ce2f5297", + "id": "0719a6e2", "metadata": {}, "outputs": [], "source": [ "# SECTION 4: Setup the constant-rate analysis\n", "constant_results = Analysis.RunAnalysisForAllNeurons(constant_case[\"trial\"], constant_case[\"cfg\"], 0)\n", - "constant_intercepts = np.asarray([fit.getCoeffs(1)[0] for fit in constant_results], dtype=float)\n", + "constant_intercepts = np.asarray([fit.getCoeffs(1)[0][0] for fit in constant_results], dtype=float)\n", "\n", "fig = _prepare_figure(\"plot(mu,'ro', 'MarkerSize',10)\", figsize=(7.5, 4.5))\n", "ax = fig.subplots(1, 1)\n", @@ -223,13 +229,13 @@ "ax.set_xticks(xloc, [f\"Neuron {idx}\" for idx in xloc])\n", "ax.set_ylabel(\"μ coefficient\")\n", "ax.set_title(\"Estimated constant-rate coefficient\")\n", - "ax.legend(loc=\"best\", frameon=False)" + "ax.legend(loc=\"best\", frameon=False)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "e4a199c4", + "id": "3a611243", "metadata": {}, "outputs": [], "source": [ @@ -245,26 +251,26 @@ " ax.set_xlabel(\"time (s)\")\n", " ax.grid(alpha=0.25)\n", "axs[0].set_ylabel(\"rate (Hz)\")\n", - "axs[1].legend(loc=\"best\", frameon=False, fontsize=8)" + "axs[1].legend(loc=\"best\", frameon=False, fontsize=8)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "26380744", + "id": "0930472e", "metadata": {}, "outputs": [], "source": [ "# SECTION 6: Case #2: Piece-wise Constant Rate Poisson Process\n", "# Next we compare a single-rate model against a two-epoch rate model.\n", "piecewise_time = np.asarray(piecewise_case[\"time_s\"], dtype=float)\n", - "piecewise_trains = list(piecewise_case[\"trains\"])" + "piecewise_trains = list(piecewise_case[\"trains\"])\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "cf2c6402", + "id": "c5d601a1", "metadata": {}, "outputs": [], "source": [ @@ -288,24 +294,24 @@ "ax.set_title(\"Ground-truth rates for the two-epoch simulation\")\n", "ax.set_xlabel(\"time (s)\")\n", "ax.set_ylabel(\"rate (Hz)\")\n", - "ax.legend(loc=\"best\", frameon=False, fontsize=8)" + "ax.legend(loc=\"best\", frameon=False, fontsize=8)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "25dfb2e5", + "id": "6c4195eb", "metadata": {}, "outputs": [], "source": [ "# SECTION 8: Setup the piecewise-rate analysis\n", - "piecewise_results = Analysis.RunAnalysisForAllNeurons(piecewise_case[\"trial\"], piecewise_case[\"cfg\"], 0)" + "piecewise_results = Analysis.RunAnalysisForAllNeurons(piecewise_case[\"trial\"], piecewise_case[\"cfg\"], 0)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "113d3272", + "id": "e5968267", "metadata": {}, "outputs": [], "source": [ @@ -334,13 +340,13 @@ " ax.set_xlabel(\"time (s)\")\n", " ax.grid(alpha=0.25)\n", "axs[0].set_ylabel(\"rate (Hz)\")\n", - "axs[1].legend(loc=\"best\", frameon=False, fontsize=8)" + "axs[1].legend(loc=\"best\", frameon=False, fontsize=8)\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "1d9d4cce", + "id": "989a4bd1", "metadata": {}, "outputs": [], "source": [ @@ -367,7 +373,7 @@ "ax.set_ylabel(\"log-likelihood\")\n", "ax.set_title(\"Per-neuron log-likelihood comparison\")\n", "ax.legend(loc=\"best\", frameon=False, fontsize=8)\n", - "__tracker.finalize()" + "__tracker.finalize()\n" ] } ], @@ -377,11 +383,11 @@ }, "nstat": { "expected_figures": 10, - "run_group": "full", + "run_group": "smoke", "style": "python-example", "topic": "ValidationDataSet" } }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/notebooks/slprj/sim/varcache/PointProcessSimulation/checksumOfCache.mat b/notebooks/slprj/sim/varcache/PointProcessSimulation/checksumOfCache.mat new file mode 100644 index 00000000..3785973a Binary files /dev/null and b/notebooks/slprj/sim/varcache/PointProcessSimulation/checksumOfCache.mat differ diff --git a/notebooks/slprj/sim/varcache/PointProcessSimulation/tmwinternal/simulink_cache.xml b/notebooks/slprj/sim/varcache/PointProcessSimulation/tmwinternal/simulink_cache.xml new file mode 100644 index 00000000..8534d00d --- /dev/null +++ b/notebooks/slprj/sim/varcache/PointProcessSimulation/tmwinternal/simulink_cache.xml @@ -0,0 +1,6 @@ + + + + D8jFo51WLz3upfIvn7uPG3TiMPEGU99lrUkoFpzEw9xmzA+++Ess23EmQAq77nvBZV26wlcL+MuRQoZL7Bz1BQ== + + \ No newline at end of file diff --git a/notebooks/slprj/sim/varcache/PointProcessSimulation/varInfo.mat b/notebooks/slprj/sim/varcache/PointProcessSimulation/varInfo.mat new file mode 100644 index 00000000..74e22f18 Binary files /dev/null and b/notebooks/slprj/sim/varcache/PointProcessSimulation/varInfo.mat differ diff --git a/nstat/confidence_interval.py b/nstat/confidence_interval.py index 3d43aef7..aa59ea91 100644 --- a/nstat/confidence_interval.py +++ b/nstat/confidence_interval.py @@ -184,6 +184,77 @@ def __rsub__(self, other): def __neg__(self): return ConfidenceInterval(self.time, np.column_stack([-self.upper, -self.lower]), self.color) + # ------------------------------------------------------------------ + # SignalObj-compatible interface methods + # In MATLAB, ConfidenceInterval < SignalObj, inheriting ~70 methods. + # We provide the most commonly used ones here for API parity. + # ------------------------------------------------------------------ + + def getSigInTimeWindow(self, wMin: float, wMax: float, holdVals: int = 0): + """Return a new ConfidenceInterval restricted to [wMin, wMax].""" + mask = (self.time >= wMin) & (self.time <= wMax) + return ConfidenceInterval(self.time[mask], self.bounds[mask], self.color, value=self.value) + + def windowedSignal(self, windowTimes): + """Return a ConfidenceInterval restricted to the given time window.""" + wt = np.asarray(windowTimes, dtype=float).reshape(-1) + return self.getSigInTimeWindow(float(wt[0]), float(wt[-1])) + + def shift(self, deltaT: float, updateLabels: bool = False): + """Return a time-shifted copy.""" + return ConfidenceInterval(self.time + deltaT, self.bounds.copy(), self.color, value=self.value) + + def resample(self, sample_rate: float): + """Resample the CI bounds to a new sample rate via linear interpolation.""" + new_time = np.arange(self.minTime, self.maxTime, 1.0 / sample_rate) + new_lower = np.interp(new_time, self.time, self.lower) + new_upper = np.interp(new_time, self.time, self.upper) + return ConfidenceInterval(new_time, np.column_stack([new_lower, new_upper]), self.color, value=self.value) + + @property + def derivative(self): + """Numerical derivative of both bounds.""" + dt = np.diff(self.time) + dt[dt == 0] = 1.0 + d_lower = np.concatenate([np.diff(self.lower) / dt, [0.0]]) + d_upper = np.concatenate([np.diff(self.upper) / dt, [0.0]]) + return ConfidenceInterval(self.time, np.column_stack([d_lower, d_upper]), self.color, value=self.value) + + def merge(self, other, holdVals: int = 0): + """Merge with another ConfidenceInterval along the time axis.""" + new_time = np.concatenate([self.time, other.time]) + new_bounds = np.vstack([self.bounds, other.bounds]) + order = np.argsort(new_time) + return ConfidenceInterval(new_time[order], new_bounds[order], self.color, value=self.value) + + def copySignal(self): + """Return a deep copy.""" + return ConfidenceInterval(self.time.copy(), self.bounds.copy(), self.color, value=self.value) + + def getSubSignal(self, identifier): + """Return a single column (lower=0 or upper=1) as a 1-column CI.""" + idx = int(identifier) + col = self.bounds[:, idx : idx + 1] + return ConfidenceInterval(self.time, np.column_stack([col, col]), self.color, value=self.value) + + def setMinTime(self, minTime=None, holdVals: int = 0): + """Restrict to times >= minTime.""" + if minTime is None: + return + mask = self.time >= minTime + self.time = self.time[mask] + self.bounds = self.bounds[mask] + self.minTime = float(self.time[0]) if self.time.size else float("nan") + + def setMaxTime(self, maxTime=None, holdVals: int = 0): + """Restrict to times <= maxTime.""" + if maxTime is None: + return + mask = self.time <= maxTime + self.time = self.time[mask] + self.bounds = self.bounds[mask] + self.maxTime = float(self.time[-1]) if self.time.size else float("nan") + def toStructure(self) -> dict: """Alias for :meth:`dataToStructure`.""" return self.dataToStructure() diff --git a/nstat/decoding_algorithms.py b/nstat/decoding_algorithms.py index a8245eaa..8c7b783c 100644 --- a/nstat/decoding_algorithms.py +++ b/nstat/decoding_algorithms.py @@ -390,28 +390,26 @@ def linear_decode(spike_counts: np.ndarray, stimulus: np.ndarray) -> dict[str, n @staticmethod def kalman_filter( - observations: np.ndarray, - transition: np.ndarray, - observation_matrix: np.ndarray, - q_cov: np.ndarray, - r_cov: np.ndarray, - x0: np.ndarray, - p0: np.ndarray, - ) -> dict[str, np.ndarray]: - """Discrete-time Kalman filter — public Python API. - - Runs a Kalman filter on time-major observations and returns a - dict with updated state estimates and covariances. + A=None, C=None, Pv=None, Pw=None, Px0=None, x0=None, y=None, GnConv=None, + *, + observations=None, transition=None, observation_matrix=None, + q_cov=None, r_cov=None, p0=None, + ): + """Discrete-time Kalman filter. - Parameters - ---------- - observations : (N, Dy) — observation time-series, one row per step. - transition : (Dx, Dx) — state-transition matrix A. - observation_matrix : (Dy, Dx) — observation matrix C (H). - q_cov : (Dx, Dx) — process-noise covariance. - r_cov : (Dy, Dy) — observation-noise covariance. - x0 : (Dx,) — initial state estimate. - p0 : (Dx, Dx) — initial error covariance. + Accepts **both** the MATLAB-compatible signature:: + + kalman_filter(A, C, Pv, Pw, Px0, x0, y, GnConv) + + and the Pythonic keyword signature:: + + kalman_filter(observations=y, transition=A, + observation_matrix=C, q_cov=Pv, + r_cov=Pw, x0=x0, p0=Px0) + + When MATLAB-style positional args are provided, delegates to + ``_kalman_filter_matlab``. Otherwise runs a simple Pythonic + Kalman filter returning ``dict(state=..., cov=...)``. Returns ------- @@ -419,15 +417,56 @@ def kalman_filter( ``state`` : (N, Dx) — updated (posterior) state estimates. ``cov`` : (N, Dx, Dx) — updated covariances. """ - y = np.asarray(observations, dtype=float) - a = np.asarray(transition, dtype=float) - h = np.asarray(observation_matrix, dtype=float) - q = np.asarray(q_cov, dtype=float) - r = np.asarray(r_cov, dtype=float) - x_prev = np.asarray(x0, dtype=float).reshape(-1) - p_prev = np.asarray(p0, dtype=float) - - n_t = y.shape[0] + # Detect MATLAB-style positional call: kalman_filter(A, C, Pv, Pw, Px0, x0, y, GnConv) + # MATLAB call has 7-8 positional args where y (arg 7) is the observation matrix. + # Pythonic call has 7 args where the first is observations (2D time-series). + # Distinguish by checking: in MATLAB style, A is square (Dx,Dx) and y is (Dy,N); + # in Pythonic style, the first arg ('A' slot) is observations (N,Dy). + _all_positional = A is not None and C is not None and Pv is not None and Pw is not None + if _all_positional and y is not None and GnConv is None: + # 7 positional args: could be MATLAB (A,C,Pv,Pw,Px0,x0,y) or Pythonic (obs,A,C,Q,R,x0,p0) + # MATLAB: A is square state matrix, y is observation time-series + # Pythonic: first arg (A slot) IS the observation time-series + a_arr = np.asarray(A, dtype=float) + y_arr = np.asarray(y, dtype=float) + # If A is 2D square and y has more rows than A, it's MATLAB style + if a_arr.ndim == 2 and a_arr.shape[0] == a_arr.shape[1] and y_arr.ndim >= 2 and y_arr.shape[1] > a_arr.shape[0]: + return DecodingAlgorithms._kalman_filter_matlab(A, C, Pv, Pw, Px0, x0, y, GnConv) + elif _all_positional and y is not None and GnConv is not None: + # 8 positional args: must be MATLAB (A,C,Pv,Pw,Px0,x0,y,GnConv) + return DecodingAlgorithms._kalman_filter_matlab(A, C, Pv, Pw, Px0, x0, y, GnConv) + + # Pythonic positional call: kalman_filter(observations, transition, obs_matrix, q, r, x0, p0) + # When 7 positional args are given and it's NOT MATLAB-style, treat as: + # A=observations, C=transition, Pv=obs_matrix, Pw=q_cov, Px0=r_cov, x0=x0, y=p0 + if observations is not None: + # Keyword call: use keyword values + obs = observations + a = transition if transition is not None else A + h = observation_matrix if observation_matrix is not None else C + q = q_cov if q_cov is not None else Pv + r = r_cov if r_cov is not None else Pw + x0_vec = x0 + p0_mat = p0 if p0 is not None else Px0 + else: + # Positional Pythonic call: A=obs, C=A, Pv=C, Pw=Q, Px0=R, x0=x0, y=p0 + obs = A + a = C + h = Pv + q = Pw + r = Px0 + x0_vec = x0 + p0_mat = y + + obs = np.asarray(obs, dtype=float) + a = np.asarray(a, dtype=float) + h = np.asarray(h, dtype=float) + q = np.asarray(q, dtype=float) + r = np.asarray(r, dtype=float) + x_prev = np.asarray(x0_vec, dtype=float).reshape(-1) + p_prev = np.asarray(p0_mat, dtype=float) + + n_t = obs.shape[0] n_x = x_prev.shape[0] xs = np.zeros((n_t, n_x), dtype=float) ps = np.zeros((n_t, n_x, n_x), dtype=float) @@ -436,7 +475,7 @@ def kalman_filter( x_pred = a @ x_prev p_pred = a @ p_prev @ a.T + q - innovation = y[t] - h @ x_pred + innovation = obs[t] - h @ x_pred s_cov = h @ p_pred @ h.T + r k_gain = p_pred @ h.T @ np.linalg.pinv(s_cov) diff --git a/nstat/fit.py b/nstat/fit.py index 6cb51e2e..310d1cca 100644 --- a/nstat/fit.py +++ b/nstat/fit.py @@ -9,6 +9,9 @@ matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np + +# matplotlib >= 3.9 renamed boxplot 'labels' to 'tick_labels' +_MPL_TICK_LABELS_KEY = "tick_labels" if tuple(int(x) for x in matplotlib.__version__.split(".")[:2]) >= (3, 9) else "labels" from scipy.stats import norm, pearsonr from .core import Covariate, nspikeTrain @@ -728,26 +731,19 @@ def _rawCoeffs(self, fit_num: int = 1) -> np.ndarray: """Return the raw coefficient vector for *fit_num* (1-based).""" return self.b[fit_num - 1].copy() - def getCoeffs(self, fit_num: int = 1) -> np.ndarray: - """Return the coefficient vector for *fit_num* (1-based). + def getCoeffs(self, fit_num: int = 1) -> tuple[np.ndarray, list[str], np.ndarray]: + """Return ``(coeffMat, labels, SEMat)`` for *fit_num* (1-based). - In Matlab ``[coeffMat, labels, SEMat] = getCoeffs(fitObj, fitNum)`` - returns multiple outputs. Use :meth:`getCoeffsWithLabels` to obtain - the full ``(coeffMat, labels, SEMat)`` tuple. + Matches MATLAB: ``[coeffMat, labels, SEMat] = getCoeffs(fitObj, fitNum)``. """ - return self._rawCoeffs(fit_num) + return self.getCoeffsWithLabels(fit_num) - def getHistCoeffs(self, fit_num: int = 1) -> np.ndarray: - """Return the history-coefficient vector for *fit_num* (1-based). + def getHistCoeffs(self, fit_num: int = 1) -> tuple[np.ndarray, list[str], np.ndarray]: + """Return ``(histMat, labels, SEMat)`` for *fit_num* (1-based). - In Matlab ``[histMat, labels, SEMat] = getHistCoeffs(fitObj, fitNum)`` - returns multiple outputs. Use :meth:`getHistCoeffsWithLabels` to - obtain the full ``(histMat, labels, SEMat)`` tuple. + Matches MATLAB: ``[histMat, labels, SEMat] = getHistCoeffs(fitObj, fitNum)``. """ - num_hist = int(self.numHist[fit_num - 1]) if fit_num - 1 < len(self.numHist) else 0 - if num_hist <= 0: - return np.array([], dtype=float) - return self._rawCoeffs(fit_num)[-num_hist:] + return self.getHistCoeffsWithLabels(fit_num) def getHistCoeffsWithLabels(self, fit_num: int = 1) -> tuple[np.ndarray, list[str], np.ndarray]: """Return ``(histMat, labels, SEMat)`` — Matlab multi-output form. @@ -1840,7 +1836,7 @@ def plotIC(self, handle=None): def plotAIC(self, handle=None): """Box-plot of AIC across neurons (Matlab ``plotAIC``).""" ax = handle if handle is not None else plt.subplots(1, 1, figsize=(5.0, 3.5))[1] - ax.boxplot(self.AIC, tick_labels=self.fitNames) + ax.boxplot(self.AIC, **{_MPL_TICK_LABELS_KEY: self.fitNames}) ax.set_ylabel("AIC") ax.set_title("AIC Across Neurons") return ax @@ -1848,7 +1844,7 @@ def plotAIC(self, handle=None): def plotBIC(self, handle=None): """Box-plot of BIC across neurons (Matlab ``plotBIC``).""" ax = handle if handle is not None else plt.subplots(1, 1, figsize=(5.0, 3.5))[1] - ax.boxplot(self.BIC, tick_labels=self.fitNames) + ax.boxplot(self.BIC, **{_MPL_TICK_LABELS_KEY: self.fitNames}) ax.set_ylabel("BIC") ax.set_title("BIC Across Neurons") return ax @@ -1856,7 +1852,7 @@ def plotBIC(self, handle=None): def plotlogLL(self, handle=None): """Box-plot of log-likelihood across neurons (Matlab ``plotlogLL``).""" ax = handle if handle is not None else plt.subplots(1, 1, figsize=(5.0, 3.5))[1] - ax.boxplot(self.logLL, tick_labels=self.fitNames) + ax.boxplot(self.logLL, **{_MPL_TICK_LABELS_KEY: self.fitNames}) ax.set_ylabel("log likelihood") ax.set_title("log likelihood Across Neurons") return ax @@ -1910,7 +1906,7 @@ def boxPlot(self, X, diffIndex: int = 1, h=None, dataLabels=None, **kwargs): 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, tick_labels=labels) + ax.boxplot(values, **{_MPL_TICK_LABELS_KEY: labels}) return ax # ------------------------------------------------------------------ diff --git a/nstat/trial.py b/nstat/trial.py index b990942d..4eac39e8 100644 --- a/nstat/trial.py +++ b/nstat/trial.py @@ -1367,22 +1367,39 @@ def psth( time = (window_times[1:] + window_times[:-1]) * 0.5 return Covariate(time, psth_data, "PSTH", "time", "s", "Hz", ["psth"]) - def psthGLM(self, binwidth: float, windowTimes=None, fitType: str = "poisson", - *, alphaVal: float = 0.05, Mc: int = 1000): + def psthGLM(self, basisWidth: float = None, history=None, fitType: str = "poisson", + selectorArray=None, minTime=None, maxTime=None, sampleRate=None, + *, binwidth: float = None, windowTimes=None, + alphaVal: float = 0.05, Mc: int = 1000): """GLM-based PSTH estimation (Matlab ``nstColl.psthGLM``). + Parameters match MATLAB: + ``psthGLM(basisWidth, history, fitType, selectorArray, minTime, maxTime, sampleRate)`` + + The Python keyword-only ``binwidth`` and ``windowTimes`` are accepted + as aliases for ``basisWidth`` and ``history`` respectively. + Returns ``(psth_covariate, histSignal, psthFitResult)`` matching the - Matlab signature. The PSTH and history covariates carry Monte Carlo - confidence intervals matching the Matlab implementation (1000 draws - from the normal approximation to the coefficient posterior, transformed - through the link function, with empirical quantile CIs). + Matlab signature. """ + # Resolve aliases: binwidth ↔ basisWidth, windowTimes ↔ history + if basisWidth is None and binwidth is not None: + basisWidth = binwidth + elif basisWidth is None and binwidth is None: + basisWidth = 0.100 # MATLAB default + if windowTimes is not None and history is None: + history = windowTimes + windowTimes = history # use unified name internally from .analysis import Analysis from .confidence_interval import ConfidenceInterval from .glm import fit_poisson_glm + # Use MATLAB-compatible param names (selectorArray/minTime/maxTime/sampleRate unused in current impl) + _sr = float(sampleRate) if sampleRate is not None else float(self.sampleRate) + _minT = float(minTime) if minTime is not None else float(self.minTime) + _maxT = float(maxTime) if maxTime is not None else float(self.maxTime) basis = self.generateUnitImpulseBasis( - float(binwidth), float(self.minTime), float(self.maxTime), float(self.sampleRate) + float(basisWidth), _minT, _maxT, _sr ) trial = Trial( SpikeTrainCollection([train.nstCopy() for train in self.nstrain]), @@ -1833,7 +1850,7 @@ def ssglm( cfgColl = ConfigCollection([cfg]) psthResult = Analysis.RunAnalysisForAllNeurons(trial, cfgColl, 0, "GLM", [], 1) fit = psthResult[0] if isinstance(psthResult, list) else psthResult - gamma0 = np.asarray(fit.getHistCoeffs(1), dtype=float).reshape(-1) + gamma0 = np.asarray(fit.getHistCoeffs(1)[0], dtype=float).reshape(-1) gamma0 = np.where(np.isnan(gamma0), -5.0, gamma0) except Exception: numHist = len(windowTimes) - 1 @@ -1927,7 +1944,7 @@ def ssglmFB( cfgColl = ConfigCollection([cfg]) psthResult = Analysis.RunAnalysisForAllNeurons(trial, cfgColl, 0, "GLM", [], 1) fit = psthResult[0] if isinstance(psthResult, list) else psthResult - gamma0 = np.asarray(fit.getHistCoeffs(1), dtype=float).reshape(-1) + gamma0 = np.asarray(fit.getHistCoeffs(1)[0], dtype=float).reshape(-1) gamma0 = np.where(np.isnan(gamma0), -5.0, gamma0) except Exception: numHist = len(windowTimes) - 1 diff --git a/parity/notebook_fidelity.yml b/parity/notebook_fidelity.yml index 049995bb..2aa9a945 100644 --- a/parity/notebook_fidelity.yml +++ b/parity/notebook_fidelity.yml @@ -1,5 +1,5 @@ version: 1 -generated_on: '2026-03-12' +generated_on: '2026-03-22' source_repositories: matlab: https://github.com/cajigaslab/nSTAT python: https://github.com/cajigaslab/nSTAT-python @@ -8,7 +8,7 @@ status_legend: - high_fidelity - partial - missing -matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT +matlab_repo_root: /private/tmp/nSTAT items: - topic: nSTATPaperExamples source_matlab: nSTATPaperExamples.mlx @@ -29,10 +29,10 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT - matlab_sections: 40 + matlab_repo_root: /private/tmp/nSTAT + matlab_sections: 37 matlab_published_figures: 30 - section_delta: 0 + section_delta: 3 figure_delta: 0 - topic: TrialExamples source_matlab: TrialExamples.mlx @@ -52,7 +52,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 9 matlab_published_figures: 6 section_delta: 0 @@ -76,7 +76,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 7 matlab_published_figures: 4 section_delta: 0 @@ -100,7 +100,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 9 matlab_published_figures: 4 section_delta: 0 @@ -124,7 +124,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 4 matlab_published_figures: 7 section_delta: 0 @@ -147,7 +147,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 2 matlab_published_figures: 2 section_delta: 0 @@ -171,7 +171,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 7 matlab_published_figures: 10 section_delta: 0 @@ -195,7 +195,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 5 matlab_published_figures: 10 section_delta: 0 @@ -219,7 +219,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 6 matlab_published_figures: 3 section_delta: 0 @@ -243,7 +243,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 17 matlab_published_figures: 9 section_delta: 0 @@ -267,7 +267,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 21 matlab_published_figures: 13 section_delta: 0 @@ -291,7 +291,7 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT + matlab_repo_root: /private/tmp/nSTAT matlab_sections: 11 matlab_published_figures: 10 section_delta: 0 @@ -315,8 +315,8 @@ items: python_tracker_only_cells: 0 python_contains_placeholders: false python_contains_tracker_only_cells: false - matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Claude/nSTAT - matlab_sections: 3 + matlab_repo_root: /private/tmp/nSTAT + matlab_sections: 4 matlab_published_figures: 4 - section_delta: 0 + section_delta: -1 figure_delta: 0 diff --git a/parity/report.md b/parity/report.md index 2afa1de9..d59f02f2 100644 --- a/parity/report.md +++ b/parity/report.md @@ -41,8 +41,8 @@ Generated from `parity/manifest.yml`, `parity/class_fidelity.yml`, `tools/notebo | Status | Count | |---|---:| -| `exact` | 13 | -| `high_fidelity` | 0 | +| `exact` | 6 | +| `high_fidelity` | 7 | | `partial` | 0 | ## Simulink Fidelity Summary diff --git a/tests/test_matlab_gold_fixtures.py b/tests/test_matlab_gold_fixtures.py index 36baa2f7..eb58bc04 100644 --- a/tests/test_matlab_gold_fixtures.py +++ b/tests/test_matlab_gold_fixtures.py @@ -588,7 +588,7 @@ def test_analysis_fit_surface_matches_matlab_gold_fixture() -> None: fit = Analysis.RunAnalysisForNeuron(trial, 1, ConfigColl([cfg])) summary = FitResSummary([fit]) - np.testing.assert_allclose(fit.getCoeffs(1), _vector(payload, "coeffs"), rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(fit.getCoeffs(1)[0], _vector(payload, "coeffs"), rtol=1e-6, atol=1e-8) np.testing.assert_allclose(fit.lambdaSignal.time, _vector(payload, "lambda_time"), rtol=1e-12, atol=1e-12) np.testing.assert_allclose(fit.lambdaSignal.data[:, 0], _vector(payload, "lambda_data"), rtol=1e-8, atol=1e-10) np.testing.assert_allclose(float(fit.AIC[0]), _scalar(payload, "AIC"), rtol=1e-8, atol=1e-10) @@ -620,8 +620,8 @@ def test_analysis_multineuron_surface_matches_matlab_gold_fixture() -> None: assert isinstance(fits, list) assert len(fits) == int(_scalar(payload, "num_fits")) - np.testing.assert_allclose(fits[0].getCoeffs(1), _vector(payload, "fit1_coeffs"), rtol=1e-6, atol=1e-8) - np.testing.assert_allclose(fits[1].getCoeffs(1), _vector(payload, "fit2_coeffs"), rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(fits[0].getCoeffs(1)[0], _vector(payload, "fit1_coeffs"), rtol=1e-6, atol=1e-8) + np.testing.assert_allclose(fits[1].getCoeffs(1)[0], _vector(payload, "fit2_coeffs"), rtol=1e-6, atol=1e-8) np.testing.assert_allclose(float(fits[0].AIC[0]), _scalar(payload, "fit1_AIC"), rtol=1e-8, atol=1e-10) np.testing.assert_allclose(float(fits[1].AIC[0]), _scalar(payload, "fit2_AIC"), rtol=1e-8, atol=1e-10) np.testing.assert_allclose(float(fits[0].BIC[0]), _scalar(payload, "fit1_BIC"), rtol=1e-8, atol=1e-10) diff --git a/tests/test_workflow_fidelity.py b/tests/test_workflow_fidelity.py index 67d42f93..277d53c5 100644 --- a/tests/test_workflow_fidelity.py +++ b/tests/test_workflow_fidelity.py @@ -62,8 +62,10 @@ def test_analysis_returns_matlab_style_fitresult_surface() -> None: assert fit.neuronNumber == 1.0 assert len(fit.covLabels) == 2 assert "stim" in fit.uniqueCovLabels - assert fit.getCoeffs(1).shape[0] >= 2 - assert fit.getHistCoeffs(1).shape[0] == 2 + coeffs, labels, se = fit.getCoeffs(1) + assert coeffs.shape[0] >= 2 + hist_coeffs, hist_labels, hist_se = fit.getHistCoeffs(1) + assert hist_coeffs.shape[0] == 2 def test_fitresult_roundtrip_and_summary_preserve_core_metadata() -> None: diff --git a/tools/notebooks/build_analysis_help_notebooks.py b/tools/notebooks/build_analysis_help_notebooks.py index 659dc666..a3ea0d2d 100644 --- a/tools/notebooks/build_analysis_help_notebooks.py +++ b/tools/notebooks/build_analysis_help_notebooks.py @@ -49,7 +49,7 @@ def _write_notebook( ## MATLAB Parity Note - Source MATLAB helpfile: `AnalysisExamples.mlx` -- Fidelity status: `high_fidelity` +- Fidelity status: `exact` - Remaining justified differences: The notebook now follows the MATLAB standard-GLM workflow with the canonical `glm_data.mat` dataset and real KS/model-visualization figures; coefficient values and styling still vary modestly because the Python GLM backend and plotting defaults differ from MATLAB. """ @@ -217,7 +217,7 @@ def _poisson_standard_errors(design_matrix, result): ## MATLAB Parity Note - Source MATLAB helpfile: `AnalysisExamples2.mlx` -- Fidelity status: `high_fidelity` +- Fidelity status: `exact` - Remaining justified differences: The notebook now follows the MATLAB toolbox workflow on the canonical `glm_data.mat` dataset with executable `Trial`, `ConfigColl`, and `Analysis` calls; exact coefficients and plot styling still vary modestly because the Python GLM backend differs from MATLAB. """ @@ -341,7 +341,7 @@ def _prepare_figure(matlab_line: str, *, figsize=(8.0, 4.5)): """ # SECTION 8: Toolbox vs. Standard GLM comparison standard_fit = fit_poisson_glm(np.column_stack([np.ones_like(xN), xN, yN, xN**2, yN**2, xN * yN]), spikes_binned, include_intercept=False) - coeff_diff = np.asarray(standard_fit.coefficients - fitResults.getCoeffs(2), dtype=float) + coeff_diff = np.asarray(standard_fit.coefficients - fitResults.getCoeffs(2)[0], dtype=float) fig = _prepare_figure("b-fitResults.b{2}", figsize=(7.0, 4.5)) ax = fig.subplots(1, 1) labels = ["mu", "x", "y", "x^2", "y^2", "x*y"] diff --git a/tools/notebooks/build_decoding_fidelity_notebooks.py b/tools/notebooks/build_decoding_fidelity_notebooks.py index 5efff3aa..ed0e9a7b 100644 --- a/tools/notebooks/build_decoding_fidelity_notebooks.py +++ b/tools/notebooks/build_decoding_fidelity_notebooks.py @@ -159,7 +159,7 @@ def _plot_decoded_ci(ax, time, decoded, cov, stim, title): ) results = Analysis.RunAnalysisForAllNeurons(trial, cfgColl, 0) - paramEst = np.column_stack([fit.getCoeffs(2)[:2] for fit in results]) + paramEst = np.column_stack([fit.getCoeffs(2)[0][:2] for fit in results]) meanParams = np.mean(paramEst, axis=1) aic_matrix = np.vstack([fit.AIC for fit in results]) logll_matrix = np.vstack([fit.logLL for fit in results]) diff --git a/tools/notebooks/build_helpfile_fidelity_notebooks.py b/tools/notebooks/build_helpfile_fidelity_notebooks.py index 38d94345..66fd99de 100644 --- a/tools/notebooks/build_helpfile_fidelity_notebooks.py +++ b/tools/notebooks/build_helpfile_fidelity_notebooks.py @@ -444,7 +444,7 @@ def _plot_isi_hist(ax, train, lambda_hz, *, title): """ # SECTION 4: Setup the constant-rate analysis constant_results = Analysis.RunAnalysisForAllNeurons(constant_case["trial"], constant_case["cfg"], 0) - constant_intercepts = np.asarray([fit.getCoeffs(1)[0] for fit in constant_results], dtype=float) + constant_intercepts = np.asarray([fit.getCoeffs(1)[0][0] for fit in constant_results], dtype=float) fig = _prepare_figure("plot(mu,'ro', 'MarkerSize',10)", figsize=(7.5, 4.5)) ax = fig.subplots(1, 1) diff --git a/tools/notebooks/parity_notes.yml b/tools/notebooks/parity_notes.yml index ade32e96..9d8f65d0 100644 --- a/tools/notebooks/parity_notes.yml +++ b/tools/notebooks/parity_notes.yml @@ -16,44 +16,49 @@ notes: file: notebooks/AnalysisExamples.ipynb source_matlab: AnalysisExamples.mlx fidelity_status: exact - remaining_differences: Complete MATLAB standard-GLM workflow with the canonical glm_data.mat dataset and real KS/model-visualization - figures. Only inherent GLM solver numerics and matplotlib styling differ. + remaining_differences: The notebook now follows the MATLAB standard-GLM workflow with the canonical `glm_data.mat` dataset + and real KS/model-visualization figures; coefficient values and styling still vary modestly because the Python GLM backend + and plotting defaults differ from MATLAB. - topic: AnalysisExamples2 file: notebooks/AnalysisExamples2.ipynb source_matlab: AnalysisExamples2.mlx fidelity_status: exact - remaining_differences: Complete MATLAB toolbox workflow on the canonical glm_data.mat dataset with executable Trial, ConfigColl, - and Analysis calls. Only inherent GLM solver numerics and plot styling differ. + remaining_differences: The notebook now follows the MATLAB toolbox workflow on the canonical `glm_data.mat` dataset with + executable `Trial`, `ConfigColl`, and `Analysis` calls; exact coefficients and plot styling still vary modestly because + the Python GLM backend differs from MATLAB. - topic: DecodingExample file: notebooks/DecodingExample.ipynb source_matlab: DecodingExample.mlx - fidelity_status: exact - remaining_differences: Workflow, model fitting, and decoded-stimulus figures follow the MATLAB helpfile. Only stochastic - simulation draws and Python plotting defaults cause trace-level variation. + fidelity_status: high_fidelity + remaining_differences: Workflow, model fitting, and decoded-stimulus figures now follow the MATLAB helpfile closely; exact + traces still depend on stochastic simulation draws and Python plotting defaults. - topic: DecodingExampleWithHist file: notebooks/DecodingExampleWithHist.ipynb source_matlab: DecodingExampleWithHist.mlx - fidelity_status: exact - remaining_differences: Mirrors the MATLAB history-aware decoding workflow. Only inherent stochastic trajectories and figure - styling differ under Python execution. + fidelity_status: high_fidelity + remaining_differences: The notebook now mirrors the MATLAB history-aware decoding workflow closely; exact stochastic trajectories + and figure styling still vary slightly under Python execution. - topic: ExplicitStimulusWhiskerData file: notebooks/ExplicitStimulusWhiskerData.ipynb source_matlab: ExplicitStimulusWhiskerData.mlx - fidelity_status: exact - remaining_differences: Reproduces the dataset-backed lag search, stimulus-effect, and history-effect workflow with real - figures. Only inherent GLM solver numerics and plotting defaults differ. + fidelity_status: high_fidelity + remaining_differences: The notebook now reproduces the dataset-backed lag search, stimulus-effect, and history-effect workflow + with real figures; exact KS traces and coefficient values still vary modestly from MATLAB because the Python GLM backend + and plotting defaults are different. - topic: HippocampalPlaceCellExample file: notebooks/HippocampalPlaceCellExample.ipynb source_matlab: HippocampalPlaceCellExample.mlx - fidelity_status: exact - remaining_differences: Reproduces the dataset-backed place-cell model-comparison and field-visualization workflow with the - same normalized 10-term Zernike basis used by MATLAB. Only inherent GLM solver numerics and surface styling differ. + fidelity_status: high_fidelity + remaining_differences: The notebook now reproduces the dataset-backed place-cell model-comparison and field-visualization + workflow with the same normalized 10-term Zernike basis used by MATLAB; exact AIC/BIC values and surface styling still + vary modestly because the Python GLM solver and plotting backend are not byte-identical to MATLAB. - topic: HybridFilterExample file: notebooks/HybridFilterExample.ipynb source_matlab: HybridFilterExample.mlx - fidelity_status: exact - remaining_differences: Reproduces the hybrid-filter simulation, single-run decoding, and averaged summary figures with real - outputs. Only inherent stochastic trajectories and Python hybrid-filter implementation details differ. + fidelity_status: high_fidelity + remaining_differences: The notebook now reproduces the hybrid-filter simulation, single-run decoding, and averaged summary + figures with real outputs; the Python port still uses the current hybrid-filter implementation instead of every MATLAB-specific + reporting branch. - topic: PPSimExample file: notebooks/PPSimExample.ipynb source_matlab: PPSimExample.mlx @@ -69,12 +74,14 @@ notes: - topic: ValidationDataSet file: notebooks/ValidationDataSet.ipynb source_matlab: ValidationDataSet.mlx - fidelity_status: exact - remaining_differences: Reproduces the constant-rate and piecewise-rate validation workflows with real Trial/Analysis objects - and figure outputs. CI uses a documented shorter deterministic fast path for stability. + fidelity_status: high_fidelity + remaining_differences: The notebook now reproduces the constant-rate and piecewise-rate validation workflows with real `Trial`/`Analysis` + objects and figure outputs; local execution uses the MATLAB-scale simulation sizes, while CI switches to a documented + shorter deterministic fast path for stability. - topic: StimulusDecode2D file: notebooks/StimulusDecode2D.ipynb source_matlab: StimulusDecode2D.mlx - fidelity_status: exact - remaining_differences: Follows the MATLAB nonlinear-CIF decoding workflow with `DecodingAlgorithms.PPDecodeFilter` and all - 4 published figures. Only inherent Python symbolic/numeric stack and random streams differ. + fidelity_status: high_fidelity + remaining_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.