From 5e8f71f42a2ba75526f65e34c286f9ce46a803c0 Mon Sep 17 00:00:00 2001 From: Iahn Cajigas Date: Sat, 7 Mar 2026 19:04:14 -0500 Subject: [PATCH 1/2] Tighten source-driven notebook fidelity audit --- notebooks/NetworkTutorial.ipynb | 648 +++++++++++------- nstat/decoding_algorithms.py | 57 ++ parity/notebook_fidelity.yml | 21 + parity/report.md | 2 +- tests/test_decoding_algorithms_fidelity.py | 21 + tests/test_matlab_symbol_surface.py | 56 ++ tests/test_network_tutorial_builder.py | 32 + tests/test_notebook_ci_groups.py | 2 + tests/test_notebook_fidelity_audit.py | 12 + .../build_network_tutorial_notebook.py | 390 +++++++++++ tools/notebooks/parity_notes.yml | 5 + tools/notebooks/topic_groups.yml | 2 + 12 files changed, 990 insertions(+), 258 deletions(-) create mode 100644 tests/test_matlab_symbol_surface.py create mode 100644 tests/test_network_tutorial_builder.py create mode 100644 tools/notebooks/build_network_tutorial_notebook.py diff --git a/notebooks/NetworkTutorial.ipynb b/notebooks/NetworkTutorial.ipynb index 07bd9012..0f6a4622 100644 --- a/notebooks/NetworkTutorial.ipynb +++ b/notebooks/NetworkTutorial.ipynb @@ -1,17 +1,22 @@ { "cells": [ + { + "cell_type": "markdown", + "id": "81a6687d", + "metadata": {}, + "source": [ + "\n", + "## MATLAB Parity Note\n", + "- Source MATLAB helpfile: `NetworkTutorial.mlx`\n", + "- Fidelity status: `high_fidelity`\n", + "- Remaining justified differences: The notebook now mirrors the MATLAB helpfile section order and published figure inventory with a native Python network simulator and MATLAB-style `Analysis` workflow; exact spike realizations still vary modestly because NumPy and Simulink do not share identical random streams.\n" + ] + }, { "cell_type": "code", - "execution_count": 1, - "id": "2befb630", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:08.343671Z", - "iopub.status.busy": "2026-03-07T21:26:08.343514Z", - "iopub.status.idle": "2026-03-07T21:26:09.835662Z", - "shell.execute_reply": "2026-03-07T21:26:09.835029Z" - } - }, + "execution_count": null, + "id": "1a559c47", + "metadata": {}, "outputs": [], "source": [ "# nSTAT-python notebook example: NetworkTutorial\n", @@ -29,8 +34,9 @@ "matplotlib.use(\"Agg\")\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", + "from matplotlib.patches import Circle, FancyArrowPatch, FancyBboxPatch\n", "\n", - "from nstat import Analysis, Covariate, FitResSummary, History, Trial, TrialConfig\n", + "from nstat import Analysis, Covariate, History, Trial, TrialConfig\n", "from nstat.ConfigColl import ConfigColl\n", "from nstat.CovColl import CovColl\n", "from nstat.notebook_figures import FigureTracker\n", @@ -38,336 +44,464 @@ "\n", "np.random.seed(0)\n", "OUTPUT_ROOT = REPO_ROOT / \"output\" / \"notebook_images\"\n", - "__tracker = FigureTracker(topic='NetworkTutorial', output_root=OUTPUT_ROOT, expected_count=4)\n", + "__tracker = FigureTracker(topic='NetworkTutorial', output_root=OUTPUT_ROOT, expected_count=14)\n", + "\n", + "\n", + "def _figure(label: str, *, figsize=(8.5, 4.5)):\n", + " fig = __tracker.new_figure(label)\n", + " fig.clear()\n", + " fig.set_size_inches(*figsize)\n", + " return fig\n", + "\n", + "\n", + "def _text_panel(fig, title: str, lines):\n", + " ax = fig.subplots(1, 1)\n", + " ax.axis(\"off\")\n", + " ax.set_title(title)\n", + " ax.text(\n", + " 0.02,\n", + " 0.98,\n", + " \"\\n\".join(lines),\n", + " va=\"top\",\n", + " ha=\"left\",\n", + " family=\"monospace\",\n", + " fontsize=11,\n", + " transform=ax.transAxes,\n", + " )\n", + " return ax\n", + "\n", + "\n", + "def _stem_kernel(ax, coeffs, title: str, xlabel: str, color: str):\n", + " coeffs = np.asarray(coeffs, dtype=float).reshape(-1)\n", + " x = np.arange(1, coeffs.size + 1, dtype=float)\n", + " markerline, stemlines, baseline = ax.stem(x, coeffs, basefmt=\" \")\n", + " plt.setp(markerline, color=color, markersize=8)\n", + " plt.setp(stemlines, color=color, linewidth=2.0)\n", + " baseline.set_visible(False)\n", + " ax.axhline(0.0, color=\"0.35\", linewidth=1.0)\n", + " ax.set_xticks(x)\n", + " ax.set_xlabel(xlabel)\n", + " ax.set_ylabel(\"coefficient\")\n", + " ax.set_title(title)\n", + " ax.grid(axis=\"y\", alpha=0.2)\n", + "\n", + "\n", + "def _draw_network(ax, actual_network):\n", + " ax.set_title(\"Two-neuron connectivity diagram\")\n", + " ax.axis(\"off\")\n", + " positions = {1: (0.25, 0.5), 2: (0.75, 0.5)}\n", + " for idx, (x, y) in positions.items():\n", + " ax.add_patch(Circle((x, y), 0.08, facecolor=\"#d9edf7\", edgecolor=\"#2b5876\", linewidth=2.0))\n", + " ax.text(x, y, f\"Neuron {idx}\", ha=\"center\", va=\"center\", fontsize=10, weight=\"bold\")\n", + " ax.text(x, y - 0.18, f\"$\\\\mu_{idx}={-3}$\", ha=\"center\", va=\"center\", fontsize=10)\n", + " ax.text(x, y + 0.18, \"history: [-4, -2, -1]\", ha=\"center\", va=\"center\", fontsize=9)\n", "\n", - "# SECTION 0: Section 0\n", - "# Author: Iahn Cajigas\n", - "# Date: 2/10/2014\n" + " arrow12 = FancyArrowPatch(\n", + " positions[1],\n", + " positions[2],\n", + " arrowstyle=\"-|>\",\n", + " mutation_scale=20,\n", + " linewidth=2.2,\n", + " color=\"#1f77b4\",\n", + " connectionstyle=\"arc3,rad=0.15\",\n", + " )\n", + " arrow21 = FancyArrowPatch(\n", + " positions[2],\n", + " positions[1],\n", + " arrowstyle=\"-|>\",\n", + " mutation_scale=20,\n", + " linewidth=2.2,\n", + " color=\"#d62728\",\n", + " connectionstyle=\"arc3,rad=-0.15\",\n", + " )\n", + " ax.add_patch(arrow12)\n", + " ax.add_patch(arrow21)\n", + " ax.text(0.5, 0.68, f\"$E_1={actual_network[0, 1]:.0f}$\", ha=\"center\", va=\"center\", fontsize=10, color=\"#1f77b4\")\n", + " ax.text(0.5, 0.32, f\"$E_2={actual_network[1, 0]:.0f}$\", ha=\"center\", va=\"center\", fontsize=10, color=\"#d62728\")\n", + " ax.text(0.25, 0.14, \"$S_1=+1 \\cdot u_{stim}$\", ha=\"center\", fontsize=10)\n", + " ax.text(0.75, 0.14, \"$S_2=-1 \\cdot u_{stim}$\", ha=\"center\", fontsize=10)\n", + "\n", + "\n", + "def _draw_block_diagram(ax):\n", + " ax.axis(\"off\")\n", + " ax.set_title(\"Conditional-intensity block diagram\")\n", + " labels = [\n", + " (\"Baseline $\\\\mu_i$\", (0.10, 0.72)),\n", + " (\"History $H * \\\\Delta N_i[n]$\", (0.10, 0.48)),\n", + " (\"Stimulus $S * u_{stim}[n]$\", (0.10, 0.24)),\n", + " (\"Ensemble $E * \\\\Delta N_k[n]$\", (0.10, 0.02)),\n", + " (\"Summation\", (0.48, 0.36)),\n", + " (\"logistic\", (0.72, 0.36)),\n", + " (\"Spike probability $\\\\lambda_i \\\\Delta$\", (0.90, 0.36)),\n", + " ]\n", + " for text, (x, y) in labels:\n", + " box = FancyBboxPatch((x, y), 0.18, 0.14, boxstyle=\"round,pad=0.03\", facecolor=\"white\", edgecolor=\"0.3\")\n", + " ax.add_patch(box)\n", + " ax.text(x + 0.09, y + 0.07, text, ha=\"center\", va=\"center\", fontsize=9)\n", + " for y in (0.79, 0.55, 0.31, 0.09):\n", + " ax.add_patch(FancyArrowPatch((0.28, y), (0.48, 0.43), arrowstyle=\"-|>\", mutation_scale=15, linewidth=1.8, color=\"0.25\"))\n", + " ax.add_patch(FancyArrowPatch((0.66, 0.43), (0.72, 0.43), arrowstyle=\"-|>\", mutation_scale=15, linewidth=1.8, color=\"0.25\"))\n", + " ax.add_patch(FancyArrowPatch((0.90, 0.43), (0.98, 0.43), arrowstyle=\"-|>\", mutation_scale=15, linewidth=1.8, color=\"0.25\"))\n", + "\n", + "\n", + "def _estimate_network(results):\n", + " estimated = np.zeros((2, 2), dtype=float)\n", + " for neuron_idx, fit in enumerate(results, start=1):\n", + " coeffs, labels, _ = fit.getCoeffsWithLabels(3)\n", + " for coeff, label in zip(coeffs, labels, strict=False):\n", + " label_str = str(label)\n", + " if neuron_idx == 1 and label_str.startswith(\"2:\"):\n", + " estimated[0, 1] = float(coeff)\n", + " elif neuron_idx == 2 and label_str.startswith(\"1:\"):\n", + " estimated[1, 0] = float(coeff)\n", + " return estimated\n" ] }, { "cell_type": "code", - "execution_count": 2, - "id": "73fc6b6a", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:09.837060Z", - "iopub.status.busy": "2026-03-07T21:26:09.836917Z", - "iopub.status.idle": "2026-03-07T21:26:09.838910Z", - "shell.execute_reply": "2026-03-07T21:26:09.838383Z" - } - }, + "execution_count": null, + "id": "1bb3afdf", + "metadata": {}, "outputs": [], "source": [ "# SECTION 1: Point Process Network Simulation\n", "# In order to understand how the point process GLM framework can be used to estimate the network connectivity within a population of neurons, we simulate a network of 2 neurons.\n", - "# This block diagram specifies a conditional intensity function of the form\n", - "# lambda_{i} \\cdot \\Delta = logistic(\\mu_{i} + H*\\Delta N_{i}[n] +\n", - "# S*u_{stim}[n] + E*\\Delta N_{k}[n]\n", - "# where, \\hbox{\\fontsize{14}{16}\\selectfont\\(logistic(x)=e^{x}/{1+e^{x}}\\)}. Note that * is the convolution opertator." + "plt.close(\"all\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60c517c7", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 2: Published network diagram\n", + "fig = _figure(\"SimulatedNetwork2.png\", figsize=(8.0, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "_draw_network(ax, np.array([[0.0, 1.0], [-4.0, 0.0]], dtype=float))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e75a150", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 3: Published block diagram\n", + "fig = _figure(\"PPSimExample-BlockDiagram.png\", figsize=(10.0, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "_draw_block_diagram(ax)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1a5a339", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 4: Conditional intensity equation\n", + "fig = _figure(\"lambda_i * Delta = logistic(mu_i + H*DeltaN_i[n] + S*u_stim[n] + E*DeltaN_k[n])\", figsize=(10.0, 3.0))\n", + "_text_panel(\n", + " fig,\n", + " \"Conditional intensity used in the tutorial\",\n", + " [\n", + " \"lambda_i * Delta = logistic(mu_i + H * DeltaN_i[n]\",\n", + " \" + S * u_stim[n] + E * DeltaN_k[n])\",\n", + " ],\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d8518832", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 5: Logistic nonlinearity\n", + "# logistic(x) = exp(x) / (1 + exp(x)). Note that * is the convolution operator.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b58aea0f", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 6: Convolution operator note\n", + "# The MATLAB helpfile presents the recursive history, stimulus, and ensemble filters separately below.\n" ] }, { "cell_type": "code", - "execution_count": 3, - "id": "e994522f", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:09.840206Z", - "iopub.status.busy": "2026-03-07T21:26:09.840103Z", - "iopub.status.idle": "2026-03-07T21:26:10.291655Z", - "shell.execute_reply": "2026-03-07T21:26:10.291224Z" - } - }, + "execution_count": null, + "id": "db9ba65f", + "metadata": {}, "outputs": [], "source": [ - "# SECTION 2: 2 Neuron Network\n", - "plt.close(\"all\")\n", + "# SECTION 7: 2 Neuron Network\n", "Ts = 0.001\n", "sampleRate = 1.0 / Ts\n", - "numNeurons = 2\n", + "tMin = 0.0\n", "tMax = 50.0\n", - "windowTimes = [0.0, Ts, 2 * Ts, 3 * Ts]\n", + "time = np.arange(tMin, tMax + Ts, Ts)\n", + "numNeurons = 2\n", + "selfHist = [0.0, Ts, 2.0 * Ts, 3.0 * Ts]\n", + "ensHist = [0.0, Ts]\n", "network = simulate_two_neuron_network(duration_s=tMax, dt=Ts, seed=4)\n", - "time = network.time\n", - "baseline_mu = network.baseline_mu\n", - "history_kernel = network.history_kernel\n", - "stim_kernel = network.stimulus_kernel\n", - "ensemble_kernel = network.ensemble_kernel\n", - "actual_network = network.actual_network\n" + "baseline_mu = np.asarray(network.baseline_mu, dtype=float)\n", + "history_kernel = np.asarray(network.history_kernel, dtype=float)\n", + "stim_kernel = np.asarray(network.stimulus_kernel, dtype=float)\n", + "ensemble_kernel = np.asarray(network.ensemble_kernel, dtype=float)\n", + "actual_network = np.asarray(network.actual_network, dtype=float)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c7098b3", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 8: Baseline firing rate of the neurons being modeled\n", + "print({\"mu1\": float(baseline_mu[0]), \"mu2\": float(baseline_mu[1]), \"sample_rate_hz\": sampleRate})\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a805e8c4", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 9: History Effect\n", + "# Captures how the firing of a neuron modulates its own probability of firing, including the refractory period and short-term history dependence.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1559df4", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 10: History kernel\n", + "fig = _figure(\"1*h[n]=-4*DeltaN[n-1]-2*DeltaN[n-2]-1*DeltaN[n-3]\", figsize=(8.0, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "_stem_kernel(ax, history_kernel, \"Self-history kernel\", \"lag (ms)\", \"tab:red\")\n", + "ax.set_xticklabels([f\"{int(k)}\" for k in [1, 2, 3]])\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b472fd0e", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 11: Stimulus Effect\n", + "# Neuron 1 is positively modulated by the stimulus and neuron 2 is negatively modulated by the same drive.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d71c27d", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 12: Stimulus filter for neuron 1\n", + "fig = _figure(\"1*s_1[n]=1*u_stim[n]\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "_stem_kernel(ax, [stim_kernel[0]], \"Stimulus filter for neuron 1\", \"lag (samples)\", \"tab:blue\")\n" ] }, { "cell_type": "code", - "execution_count": 4, - "id": "d6118ad2", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:10.293171Z", - "iopub.status.busy": "2026-03-07T21:26:10.293079Z", - "iopub.status.idle": "2026-03-07T21:26:10.294995Z", - "shell.execute_reply": "2026-03-07T21:26:10.294551Z" - } - }, + "execution_count": null, + "id": "406603d4", + "metadata": {}, "outputs": [], "source": [ - "# SECTION 3: Baseline firing rate of the neurons being modeled" + "# SECTION 13: Stimulus filter for neuron 2\n", + "fig = _figure(\"1*s_2[n]=-1*u_stim[n]\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "_stem_kernel(ax, [stim_kernel[1]], \"Stimulus filter for neuron 2\", \"lag (samples)\", \"tab:orange\")\n" ] }, { "cell_type": "code", - "execution_count": 5, - "id": "698dd47b", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:10.296263Z", - "iopub.status.busy": "2026-03-07T21:26:10.296151Z", - "iopub.status.idle": "2026-03-07T21:26:10.298057Z", - "shell.execute_reply": "2026-03-07T21:26:10.297615Z" - } - }, + "execution_count": null, + "id": "f9878dce", + "metadata": {}, "outputs": [], "source": [ - "# SECTION 4: History Effect\n", - "# Captures how the firing of a neuron at modulates its probability of firing. Captures effects such as the refractory period and bursting. We use the same firing history for both neurons in this example. Note that the firing activity at time n leads to strong inhibition at time n+1 (refractory period) and that this effect becomes smaller over the next two time periods.\n", - "# 1*h[n]=-4*\\Delta N[n-1]-2*\\Delta N[n-2] -1*\\Delta N[n-3]\n", - "# Note that the one sample delay in same cell firing is included in the simulink model." + "# SECTION 14: Ensemble Effect\n", + "# Captures how neighboring neuron firing modulates the firing probability of a given neuron, with a one-sample delay included in the Simulink model.\n" ] }, { "cell_type": "code", - "execution_count": 6, - "id": "ba149a4f", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:10.299087Z", - "iopub.status.busy": "2026-03-07T21:26:10.299001Z", - "iopub.status.idle": "2026-03-07T21:26:10.300612Z", - "shell.execute_reply": "2026-03-07T21:26:10.300206Z" - } - }, + "execution_count": null, + "id": "7ddb2550", + "metadata": {}, "outputs": [], "source": [ - "# SECTION 5: Stimulus Effect\n", - "# 1*s_{1}[n]=1*u_{stim}[n]\n", - "# 1*s_{2}[n]=-1*u_{stim}[n]\n", - "# Neuron 1 is positively modulated by the stimulus\n", - "# Neuron 1 is negatively modulated by the stimulus" + "# SECTION 15: Ensemble filter for neuron 1\n", + "fig = _figure(\"1*e_1[n]=1*DeltaN_2[n-1]\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "_stem_kernel(ax, [ensemble_kernel[0]], \"Ensemble filter for neuron 1\", \"lag (samples)\", \"tab:green\")\n" ] }, { "cell_type": "code", - "execution_count": 7, - "id": "559a0b8e", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:10.301593Z", - "iopub.status.busy": "2026-03-07T21:26:10.301498Z", - "iopub.status.idle": "2026-03-07T21:26:10.303180Z", - "shell.execute_reply": "2026-03-07T21:26:10.302782Z" - } - }, + "execution_count": null, + "id": "1554e895", + "metadata": {}, "outputs": [], "source": [ - "# SECTION 6: Ensemble Effect\n", - "# Captures the effect of how neighboring neuron firing modulates the firing of a given neuron.\n", - "# 1*e_{1}[n]=1*\\Delta N_{2}[n-1]\n", - "# 1*e_{2}[n]=-4*\\Delta N_{1}[n-1]\n", - "# Note that the one sample delay in firing of the neighbor cell is included in the simulink model.\n", - "# Neuron 2 firing positively modulates Neuron 1\n", - "# Neuron 1 firing has strong inhibitory effect on neuron 2." + "# SECTION 16: Ensemble filter for neuron 2\n", + "fig = _figure(\"1*e_2[n]=-4*DeltaN_1[n-1]\", figsize=(7.5, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "_stem_kernel(ax, [ensemble_kernel[1]], \"Ensemble filter for neuron 2\", \"lag (samples)\", \"tab:purple\")\n" ] }, { "cell_type": "code", - "execution_count": 8, - "id": "8c536ce0", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:10.304408Z", - "iopub.status.busy": "2026-03-07T21:26:10.304274Z", - "iopub.status.idle": "2026-03-07T21:26:10.306879Z", - "shell.execute_reply": "2026-03-07T21:26:10.306461Z" - } - }, + "execution_count": null, + "id": "fcb5f672", + "metadata": {}, "outputs": [], "source": [ - "# SECTION 7: Stimulus\n", - "# We use a simple sine wave here but we may want to explore other types of inputs to see if they affect the recovery of the network parameters.\n", - "f = 1\n", + "# SECTION 17: Stimulus\n", + "f = 1.0\n", "stimCov = Covariate(time, network.latent_drive, \"Stimulus\", \"time\", \"s\", \"Voltage\", [\"sin\"])\n", - "baselineCov = Covariate(time, np.ones_like(time), \"Baseline\", \"time\", \"s\", \"\", [\"mu\"])\n" + "baselineCov = Covariate(time, np.ones_like(time), \"Baseline\", \"time\", \"s\", \"\", [\"mu\"])\n", + "\n", + "fig = _figure(\"Stimulus sine-wave drive\", figsize=(9.0, 4.5))\n", + "ax = fig.subplots(1, 1)\n", + "ax.plot(time, network.latent_drive, color=\"tab:blue\", linewidth=1.2)\n", + "ax.set_xlim(0.0, 5.0)\n", + "ax.set_xlabel(\"time (s)\")\n", + "ax.set_ylabel(\"stimulus\")\n", + "ax.set_title(\"Sine-wave stimulus used by the network tutorial\")\n" ] }, { "cell_type": "code", - "execution_count": 9, - "id": "02b96d49", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:10.307944Z", - "iopub.status.busy": "2026-03-07T21:26:10.307839Z", - "iopub.status.idle": "2026-03-07T21:26:11.691196Z", - "shell.execute_reply": "2026-03-07T21:26:11.690456Z" - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "{'fitType': 'binomial', 'algorithm': 'BNLRCG', 'num_neurons': 2, 'num_spikes': [2590, 2365]}\n" - ] - } - ], + "execution_count": null, + "id": "37804e6b", + "metadata": {}, + "outputs": [], "source": [ - "# SECTION 8: Simulate the Network\n", - "# Uses a binomial model for the conditional intensity function; this reproduces the MATLAB tutorial workflow with a native Python simulation.\n", + "# SECTION 18: Simulate the Network\n", "fitType = \"binomial\"\n", "Algorithm = \"BNLRCG\" if fitType == \"binomial\" else \"GLM\"\n", "spikeColl = network.spikes\n", - "trial = Trial(spikeColl, CovColl([stimCov, baselineCov]), None, History(windowTimes))\n", - "trial.setEnsCovHist(windowTimes)\n", + "trial = Trial(spikeColl, CovColl([stimCov, baselineCov]), None, History(selfHist))\n", + "trial.setEnsCovHist(ensHist)\n", "\n", - "fig = __tracker.new_figure(\"network-simulation-overview\")\n", - "fig.clear()\n", - "ax1, ax2 = fig.subplots(2, 1, sharex=True)\n", - "ax1.plot(time, network.lambda_delta[:, 0], label=\"Neuron 1\")\n", - "ax1.plot(time, network.lambda_delta[:, 1], label=\"Neuron 2\")\n", - "ax1.set_ylabel(\"spike probability\")\n", - "ax1.set_title(\"Simulated conditional intensity (binomial)\")\n", - "ax1.legend(loc=\"upper right\")\n", - "stimCov.plot(handle=ax2)\n", - "ax2.set_title(\"Stimulus drive\")\n", - "ax2.set_xlim(0.0, tMax / 10.0)\n", - "fig.tight_layout()\n", + "fig = _figure(\"Simulated raster and stimulus\", figsize=(9.0, 6.5))\n", + "axs = fig.subplots(2, 1, sharex=True)\n", + "spikeColl.plot(handle=axs[0])\n", + "axs[0].set_xlim(0.0, tMax / 10.0)\n", + "axs[0].set_title(\"Simulated spike trains\")\n", + "stimCov.plot(handle=axs[1])\n", + "axs[1].set_xlim(0.0, tMax / 10.0)\n", + "axs[1].set_title(\"Stimulus over the same window\")\n", "\n", - "fig = __tracker.new_figure(\"network-raster\")\n", - "fig.clear()\n", + "fig = _figure(\"Simulated lambda traces\", figsize=(9.0, 4.5))\n", "ax = fig.subplots(1, 1)\n", - "spikeColl.plot(handle=ax)\n", - "ax.set_xlim(0.0, tMax / 10.0)\n", - "ax.set_title(\"Simulated 2-neuron raster\")\n", - "fig.tight_layout()\n", - "\n", - "print({\n", - " \"fitType\": fitType,\n", - " \"algorithm\": Algorithm,\n", - " \"num_neurons\": spikeColl.numSpikeTrains,\n", - " \"num_spikes\": [spikeColl.getNST(i + 1).n_spikes for i in range(numNeurons)],\n", - "})\n" + "ax.plot(time, network.lambda_delta[:, 0], label=\"Neuron 1\", linewidth=1.1)\n", + "ax.plot(time, network.lambda_delta[:, 1], label=\"Neuron 2\", linewidth=1.1)\n", + "ax.set_xlim(0.0, 5.0)\n", + "ax.set_xlabel(\"time (s)\")\n", + "ax.set_ylabel(\"lambda * Delta\")\n", + "ax.set_title(\"Conditional-intensity trajectories\")\n", + "ax.legend(loc=\"upper right\")\n" ] }, { "cell_type": "code", - "execution_count": 10, - "id": "a80e7e8e", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:11.692711Z", - "iopub.status.busy": "2026-03-07T21:26:11.692559Z", - "iopub.status.idle": "2026-03-07T21:26:11.694542Z", - "shell.execute_reply": "2026-03-07T21:26:11.693991Z" - } - }, + "execution_count": null, + "id": "3455475c", + "metadata": {}, "outputs": [], "source": [ - "# SECTION 9: GLM Model Fitting Setup\n", - "# In this section, we create the appropriate structures to fit several GLM models to the data generated above.\n", - "# Create a constant covariate representing the mean firing rate $$\\mu_{i}$\n", - "#\n", - "# Use stimulation and baseline as possible covariates\n", - "# trial.setTrialPartition([0 tMax/2 tMax]);" + "# SECTION 19: GLM Model Fitting Setup\n", + "c1 = TrialConfig([[\"Baseline\", \"mu\"]], sampleRate, [], [], [], name=\"Baseline\")\n", + "c2 = TrialConfig([[\"Baseline\", \"mu\"]], sampleRate, [], ensHist, [], name=\"Baseline+EnsHist\")\n", + "c3 = TrialConfig([[\"Baseline\", \"mu\"], [\"Stimulus\", \"sin\"]], sampleRate, selfHist, ensHist, [], name=\"Stim+Hist+EnsHist\")\n", + "cfgColl = ConfigColl([c1, c2, c3])\n" ] }, { "cell_type": "code", - "execution_count": 11, - "id": "edcc4c4d", - "metadata": { - "execution": { - "iopub.execute_input": "2026-03-07T21:26:11.695914Z", - "iopub.status.busy": "2026-03-07T21:26:11.695759Z", - "iopub.status.idle": "2026-03-07T21:26:18.401565Z", - "shell.execute_reply": "2026-03-07T21:26:18.400983Z" - } - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "{'config_names': ['Baseline', 'Baseline+EnsHist', 'Stim+Hist+EnsHist'], 'estimated_network': [[0.0, 0.0], [0.0, 0.0]]}\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/var/folders/tg/z6dfb8b13wg_h4f3v8whzpgh0000gn/T/ipykernel_95209/2639979657.py:36: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n", - " fig.tight_layout()\n", - "/Users/iahncajigas/Library/CloudStorage/Dropbox/Codex/nSTAT-python/nstat/notebook_figures.py:42: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n", - " self._active_fig.tight_layout()\n" - ] - } - ], + "execution_count": null, + "id": "84eb9f17", + "metadata": {}, + "outputs": [], "source": [ - "# SECTION 10: GLM Model Fitting and Results\n", - "configs = ConfigColl([\n", - " TrialConfig([[\"Baseline\", \"mu\"]], sampleRate, [], [], [], name=\"Baseline\"),\n", - " TrialConfig([[\"Baseline\", \"mu\"]], sampleRate, [], [0.0, Ts], [], name=\"Baseline+EnsHist\"),\n", - " TrialConfig([[\"Baseline\", \"mu\"], [\"Stimulus\", \"sin\"]], sampleRate, windowTimes, [0.0, Ts], [], name=\"Stim+Hist+EnsHist\"),\n", - "])\n", - "results = Analysis.RunAnalysisForAllNeurons(trial, configs, 0, Algorithm)\n", + "# SECTION 20: GLM Model Fitting and Results\n", + "results = Analysis.RunAnalysisForAllNeurons(trial, cfgColl, 0, Algorithm)\n", "if not isinstance(results, list):\n", " results = [results]\n", - "summary = FitResSummary(results)\n", "\n", - "fig = __tracker.new_figure(\"network-summary\")\n", - "summary.plotSummary(handle=fig)\n", + "fig = _figure(\"results{1}.plotResults\", figsize=(11.0, 8.0))\n", + "results[0].plotResults(handle=fig)\n", "\n", - "est_network = np.zeros((2, 2), dtype=float)\n", - "for neuron_idx, fit in enumerate(results, start=1):\n", - " coeffs, labels, _ = fit.getCoeffsWithLabels(3)\n", - " for coeff, label in zip(coeffs, labels, strict=False):\n", - " if neuron_idx == 1 and str(label).startswith(\"2:\"):\n", - " est_network[0, 1] = coeff\n", - " if neuron_idx == 2 and str(label).startswith(\"1:\"):\n", - " est_network[1, 0] = coeff\n", + "fig = _figure(\"results{2}.plotResults\", figsize=(11.0, 8.0))\n", + "results[1].plotResults(handle=fig)\n", "\n", - "fig = __tracker.new_figure(\"network-actual-vs-estimated\")\n", - "fig.clear()\n", - "ax1, ax2 = fig.subplots(1, 2)\n", - "im1 = ax1.imshow(actual_network, cmap=\"coolwarm\", vmin=-4.0, vmax=1.0)\n", - "ax1.set_title(\"Actual\")\n", - "ax1.set_xticks([0, 1], [\"N1\", \"N2\"])\n", - "ax1.set_yticks([0, 1], [\"N1\", \"N2\"])\n", - "ax2.imshow(est_network, cmap=\"coolwarm\", vmin=-4.0, vmax=1.0)\n", - "ax2.set_title(\"Estimated 1ms\")\n", - "ax2.set_xticks([0, 1], [\"N1\", \"N2\"])\n", - "ax2.set_yticks([0, 1], [\"N1\", \"N2\"])\n", - "fig.colorbar(im1, ax=[ax1, ax2], shrink=0.8)\n", - "fig.tight_layout()\n", - "print({\"config_names\": summary.fitNames, \"estimated_network\": est_network.tolist()})\n", + "estimated_network = _estimate_network(results)\n", + "fig = _figure(\"Actual vs estimated network\", figsize=(10.0, 4.5))\n", + "axs = fig.subplots(1, 2)\n", + "clim = float(np.max(np.abs(actual_network)))\n", + "for ax, matrix, title in zip(\n", + " axs,\n", + " [actual_network, estimated_network],\n", + " [\"Actual\", \"Estimated 1 ms\"],\n", + " strict=False,\n", + "):\n", + " im = ax.imshow(matrix, cmap=\"jet\", vmin=-clim, vmax=clim)\n", + " ax.set_xticks([0, 1], [\"1\", \"2\"])\n", + " ax.set_yticks([0, 1], [\"1\", \"2\"])\n", + " ax.set_title(title)\n", + " for (row, col), value in np.ndenumerate(matrix):\n", + " ax.text(col, row, f\"{value:.2f}\", ha=\"center\", va=\"center\", color=\"white\", fontsize=10)\n", + "fig.colorbar(im, ax=axs, shrink=0.8)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "165afc2b", + "metadata": {}, + "outputs": [], + "source": [ + "# SECTION 21: Neighbor-selection note\n", + "# By default all neurons are considered potential neighbors. To restrict candidate neighbors, call trial.setNeighbors(neighborArray) using the MATLAB-style convention described in the source helpfile.\n", + "print(\n", + " {\n", + " \"algorithm\": Algorithm,\n", + " \"spike_counts\": [spikeColl.getNST(1).n_spikes, spikeColl.getNST(2).n_spikes],\n", + " \"estimated_network\": np.round(estimated_network, 3).tolist(),\n", + " }\n", + ")\n", "__tracker.finalize()\n" ] } ], "metadata": { "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.4" + "name": "python" }, "nstat": { - "expected_figures": 4, + "expected_figures": 14, "run_group": "full", "style": "python-example", "topic": "NetworkTutorial" @@ -375,4 +509,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/nstat/decoding_algorithms.py b/nstat/decoding_algorithms.py index 424a3621..b19820f1 100644 --- a/nstat/decoding_algorithms.py +++ b/nstat/decoding_algorithms.py @@ -7,6 +7,7 @@ from .cif import CIF from .errors import UnsupportedWorkflowError +from .nspikeTrain import nspikeTrain def _as_observation_matrix(dN) -> np.ndarray: @@ -519,6 +520,60 @@ def PPDecode_predict(x_u, W_u, A, Q, Wconv=None): x_p = A_mat @ x_vec return x_p, W_p + @staticmethod + def PPDecode_update(x_p, W_p, dN, lambdaIn, binwidth=0.001, time_index=1, WuConv=None): + x_vec = np.asarray(x_p, dtype=float).reshape(-1) + W_mat = _as_state_matrix(W_p, x_vec.size) + obs = _as_observation_matrix(dN) + idx = max(1, min(int(time_index), obs.shape[1])) + + if isinstance(lambdaIn, CIF): + lambda_items = [lambdaIn] + elif isinstance(lambdaIn, Sequence) and not isinstance(lambdaIn, (str, bytes)): + lambda_items = list(lambdaIn) + else: + raise ValueError("Lambda must be a cell of CIFs or a CIF") + if not lambda_items: + raise ValueError("Lambda must be a non-empty cell of CIFs or a CIF") + + lambda_delta = np.zeros((len(lambda_items), 1), dtype=float) + sum_val_vec = np.zeros(x_vec.size, dtype=float) + sum_val_mat = np.zeros((x_vec.size, x_vec.size), dtype=float) + observed = obs[:, idx - 1] + + for cell_index, cif in enumerate(lambda_items): + if not isinstance(cif, CIF): + raise ValueError("Lambda must be a cell of CIFs or a CIF") + if cif.historyMat.size == 0: + spike_times = (np.where(obs[cell_index] == 1.0)[0]) * float(binwidth) + nst = nspikeTrain(spike_times, makePlots=-1) + nst.setMinTime(0.0) + nst.setMaxTime((obs.shape[1] - 1) * float(binwidth)) + nst = nst.resample(1.0 / float(binwidth)) + lambda_delta[cell_index, 0] = float(cif.evalLambdaDelta(x_vec, idx, nst)) + sum_val_vec += observed[cell_index] * np.asarray(cif.evalGradientLog(x_vec, idx, nst), dtype=float).reshape(-1) + sum_val_vec -= np.asarray(cif.evalGradient(x_vec, idx, nst), dtype=float).reshape(-1) + sum_val_mat -= np.asarray(cif.evalJacobianLog(x_vec, idx, nst), dtype=float) + sum_val_mat += np.asarray(cif.evalJacobian(x_vec, idx, nst), dtype=float) + else: + lambda_delta[cell_index, 0] = float(cif.evalLambdaDelta(x_vec, idx)) + sum_val_vec += observed[cell_index] * np.asarray(cif.evalGradientLog(x_vec, idx), dtype=float).reshape(-1) + sum_val_vec -= np.asarray(cif.evalGradient(x_vec, idx), dtype=float).reshape(-1) + sum_val_mat -= np.asarray(cif.evalJacobianLog(x_vec, idx), dtype=float) + sum_val_mat += np.asarray(cif.evalJacobian(x_vec, idx), dtype=float) + + if _is_empty_value(WuConv): + identity = np.eye(W_mat.shape[0], dtype=float) + try: + W_u = W_mat @ (identity - np.linalg.solve(identity + sum_val_mat @ W_mat, sum_val_mat @ W_mat)) + except np.linalg.LinAlgError: + W_u = W_mat.copy() + W_u = _symmetrize(W_u) + else: + W_u = _symmetrize(_as_state_matrix(WuConv, x_vec.size)) + x_u = x_vec + W_u @ sum_val_vec + return x_u, W_u, lambda_delta + @staticmethod def PPDecode_updateLinear(x_p, W_p, dN, mu, beta, fitType="poisson", gamma=None, HkAll=None, time_index=1, WuConv=None): x_vec = np.asarray(x_p, dtype=float).reshape(-1) @@ -856,6 +911,7 @@ def PPHybridFilter(A, Q, p_ij, Mu0, dN, lambdaCIFColl, binwidth=0.001, x0=None, PPDecodeFilter = DecodingAlgorithms.PPDecodeFilter PPDecodeFilterLinear = DecodingAlgorithms.PPDecodeFilterLinear PPDecode_predict = DecodingAlgorithms.PPDecode_predict +PPDecode_update = DecodingAlgorithms.PPDecode_update PPDecode_updateLinear = DecodingAlgorithms.PPDecode_updateLinear PPHybridFilter = DecodingAlgorithms.PPHybridFilter PPHybridFilterLinear = DecodingAlgorithms.PPHybridFilterLinear @@ -874,6 +930,7 @@ def PPHybridFilter(A, Q, p_ij, Mu0, dN, lambdaCIFColl, binwidth=0.001, x0=None, "PPDecodeFilter", "PPDecodeFilterLinear", "PPDecode_predict", + "PPDecode_update", "PPDecode_updateLinear", "PPHybridFilter", "PPHybridFilterLinear", diff --git a/parity/notebook_fidelity.yml b/parity/notebook_fidelity.yml index 7f72a13d..46b2c474 100644 --- a/parity/notebook_fidelity.yml +++ b/parity/notebook_fidelity.yml @@ -212,6 +212,27 @@ items: matlab_published_figures: 8 section_delta: 0 figure_delta: 0 +- topic: NetworkTutorial + source_matlab: NetworkTutorial.mlx + python_notebook: notebooks/NetworkTutorial.ipynb + fidelity_status: high_fidelity + remaining_differences: The notebook now mirrors the MATLAB helpfile section order + and published figure inventory with a native Python network simulator and MATLAB-style + `Analysis` workflow; exact spike realizations still vary modestly because NumPy + and Simulink do not share identical random streams. + python_sections: 21 + python_expected_figures: 14 + python_uses_figure_tracker: true + python_has_finalize_call: true + python_placeholder_cells: 0 + python_tracker_only_cells: 0 + python_contains_placeholders: false + python_contains_tracker_only_cells: false + matlab_repo_root: /Users/iahncajigas/Library/CloudStorage/Dropbox/Codex/nSTAT + matlab_sections: 21 + matlab_published_figures: 14 + section_delta: 0 + figure_delta: 0 - topic: ValidationDataSet source_matlab: ValidationDataSet.mlx python_notebook: notebooks/ValidationDataSet.ipynb diff --git a/parity/report.md b/parity/report.md index d0f8e10f..f5ccd2d4 100644 --- a/parity/report.md +++ b/parity/report.md @@ -34,7 +34,7 @@ Generated from `parity/manifest.yml`, `parity/class_fidelity.yml`, and `tools/no | Status | Count | |---|---:| | `exact` | 0 | -| `high_fidelity` | 12 | +| `high_fidelity` | 13 | | `partial` | 0 | ## Simulink Fidelity Summary diff --git a/tests/test_decoding_algorithms_fidelity.py b/tests/test_decoding_algorithms_fidelity.py index c446357e..d075a0c7 100644 --- a/tests/test_decoding_algorithms_fidelity.py +++ b/tests/test_decoding_algorithms_fidelity.py @@ -53,6 +53,27 @@ def test_ppdecodefilter_accepts_cif_collections_with_history() -> None: assert np.all(np.isfinite(x_u)) +def test_ppdecode_update_matches_matlab_facing_public_surface() -> None: + dN = np.array([[0.0, 1.0, 0.0, 1.0]], dtype=float) + lambda_cif = CIF([0.1, 0.4], ["1", "x"], ["x"], fitType="binomial") + + x_u, W_u, lambda_delta = DecodingAlgorithms.PPDecode_update( + np.array([0.0], dtype=float), + np.array([[1.0]], dtype=float), + dN, + lambda_cif, + 0.1, + 2, + ) + + assert x_u.shape == (1,) + assert W_u.shape == (1, 1) + assert lambda_delta.shape == (1, 1) + assert np.all(np.isfinite(x_u)) + assert np.all(np.isfinite(W_u)) + assert np.all(lambda_delta > 0.0) + + def test_pphybridfilterlinear_returns_model_probabilities_and_state_banks() -> None: a = [np.array([[1.0]], dtype=float), np.array([[0.9]], dtype=float)] q = [np.array([[0.02]], dtype=float), np.array([[0.05]], dtype=float)] diff --git a/tests/test_matlab_symbol_surface.py b/tests/test_matlab_symbol_surface.py new file mode 100644 index 00000000..7c70779f --- /dev/null +++ b/tests/test_matlab_symbol_surface.py @@ -0,0 +1,56 @@ +from __future__ import annotations + +import inspect + +from nstat import Analysis, CIF, DecodingAlgorithms + + +EXPECTED_SYMBOLS = { + Analysis: { + "RunAnalysisForNeuron", + "RunAnalysisForAllNeurons", + "GLMFit", + "KSPlot", + "plotFitResidual", + "computeFitResidual", + "computeKSStats", + "plotInvGausTrans", + "plotSeqCorr", + "plotCoeffs", + }, + CIF: { + "setSpikeTrain", + "setHistory", + "simulateCIFByThinningFromLambda", + "simulateCIF", + "evalGradient", + "evalGradientLog", + "evalJacobian", + "evalJacobianLog", + "evalGradientLDGamma", + "evalJacobianLDGamma", + }, + DecodingAlgorithms: { + "PPDecode_predict", + "PPDecode_update", + "PPDecode_updateLinear", + "PPDecodeFilterLinear", + "PPDecodeFilter", + "PP_fixedIntervalSmoother", + "PPHybridFilterLinear", + "PPHybridFilter", + }, +} + + +def test_expected_matlab_symbol_surface_exists_and_is_callable() -> None: + for obj, expected in EXPECTED_SYMBOLS.items(): + missing = sorted(name for name in expected if not callable(getattr(obj, name, None))) + assert not missing, f"{obj.__name__} is missing MATLAB-facing callables: {missing}" + + +def test_expected_symbol_surface_has_python_runtime_signatures() -> None: + for obj, expected in EXPECTED_SYMBOLS.items(): + for name in expected: + signature = inspect.signature(getattr(obj, name)) + assert signature is not None diff --git a/tests/test_network_tutorial_builder.py b/tests/test_network_tutorial_builder.py new file mode 100644 index 00000000..275a6d40 --- /dev/null +++ b/tests/test_network_tutorial_builder.py @@ -0,0 +1,32 @@ +from __future__ import annotations + +from pathlib import Path + +import nbformat + +from tools.notebooks.build_network_tutorial_notebook import build_notebook + + +REPO_ROOT = Path(__file__).resolve().parents[1] +NOTEBOOK_PATH = REPO_ROOT / "notebooks" / "NetworkTutorial.ipynb" + + +def _normalize_notebook(notebook) -> None: + for cell in notebook.cells: + cell["id"] = "normalized" + cell["execution_count"] = None + cell["outputs"] = [] + + +def _cell_payload(cell) -> tuple[str, str, dict]: + return cell.cell_type, "".join(cell.get("source", "")), dict(cell.get("metadata", {})) + + +def test_network_tutorial_builder_matches_committed_notebook() -> None: + committed = nbformat.read(NOTEBOOK_PATH, as_version=4) + generated = build_notebook() + _normalize_notebook(committed) + _normalize_notebook(generated) + assert committed.metadata == generated.metadata + assert len(committed.cells) == len(generated.cells) + assert [_cell_payload(cell) for cell in committed.cells] == [_cell_payload(cell) for cell in generated.cells] diff --git a/tests/test_notebook_ci_groups.py b/tests/test_notebook_ci_groups.py index e00355fa..aa173bf1 100644 --- a/tests/test_notebook_ci_groups.py +++ b/tests/test_notebook_ci_groups.py @@ -19,6 +19,7 @@ "ExplicitStimulusWhiskerData", "HippocampalPlaceCellExample", "HybridFilterExample", + "NetworkTutorial", "PPSimExample", "SignalObjExamples", "StimulusDecode2D", @@ -35,6 +36,7 @@ "ExplicitStimulusWhiskerData", "HippocampalPlaceCellExample", "HybridFilterExample", + "NetworkTutorial", "PPSimExample", "StimulusDecode2D", "TrialExamples", diff --git a/tests/test_notebook_fidelity_audit.py b/tests/test_notebook_fidelity_audit.py index 494a2eff..f36b68a3 100644 --- a/tests/test_notebook_fidelity_audit.py +++ b/tests/test_notebook_fidelity_audit.py @@ -39,6 +39,7 @@ def test_notebook_fidelity_audit_marks_upgraded_ports_as_high_fidelity() -> None assert { "AnalysisExamples", "AnalysisExamples2", + "NetworkTutorial", "PPSimExample", "nSTATPaperExamples", } <= high_fidelity_topics @@ -59,6 +60,17 @@ def test_high_fidelity_notebooks_have_no_placeholder_or_tracker_only_cells() -> assert not row["python_contains_tracker_only_cells"], f"{row['topic']} still contains tracker-only cells" +def test_high_fidelity_notebooks_have_near_matlab_structural_counts() -> None: + audit = yaml.safe_load(AUDIT_PATH.read_text(encoding="utf-8")) or {} + for row in audit.get("items", []): + if row["fidelity_status"] not in {"high_fidelity", "exact"}: + continue + if row.get("section_delta") is None or row.get("figure_delta") is None: + continue + assert abs(int(row["section_delta"])) <= 1, f"{row['topic']} has a large MATLAB section delta" + assert abs(int(row["figure_delta"])) <= 1, f"{row['topic']} has a large MATLAB figure delta" + + def test_notebook_fidelity_audit_matches_generator_when_matlab_repo_is_available() -> None: matlab_repo = default_matlab_repo_root(REPO_ROOT) if not matlab_repo.exists(): diff --git a/tools/notebooks/build_network_tutorial_notebook.py b/tools/notebooks/build_network_tutorial_notebook.py new file mode 100644 index 00000000..23b3f165 --- /dev/null +++ b/tools/notebooks/build_network_tutorial_notebook.py @@ -0,0 +1,390 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +from pathlib import Path +from textwrap import dedent + +import nbformat +from nbformat.v4 import new_code_cell, new_markdown_cell, new_notebook + + +REPO_ROOT = Path(__file__).resolve().parents[2] +NOTEBOOK_PATH = REPO_ROOT / "notebooks" / "NetworkTutorial.ipynb" + + +LANGUAGE_METADATA = { + "language_info": { + "name": "python", + } +} + + +NOTE = """\ + +## MATLAB Parity Note +- Source MATLAB helpfile: `NetworkTutorial.mlx` +- Fidelity status: `high_fidelity` +- Remaining justified differences: The notebook now mirrors the MATLAB helpfile section order and published figure inventory with a native Python network simulator and MATLAB-style `Analysis` workflow; exact spike realizations still vary modestly because NumPy and Simulink do not share identical random streams. +""" + + +CODE = [ + """ + # nSTAT-python notebook example: NetworkTutorial + from pathlib import Path + import sys + + REPO_ROOT = Path.cwd().resolve().parent + if str(REPO_ROOT) not in sys.path: + sys.path.insert(0, str(REPO_ROOT)) + SRC_PATH = (REPO_ROOT / "src").resolve() + if str(SRC_PATH) not in sys.path: + sys.path.insert(0, str(SRC_PATH)) + + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + import numpy as np + from matplotlib.patches import Circle, FancyArrowPatch, FancyBboxPatch + + from nstat import Analysis, Covariate, History, Trial, TrialConfig + from nstat.ConfigColl import ConfigColl + from nstat.CovColl import CovColl + from nstat.notebook_figures import FigureTracker + from nstat.simulators import simulate_two_neuron_network + + np.random.seed(0) + OUTPUT_ROOT = REPO_ROOT / "output" / "notebook_images" + __tracker = FigureTracker(topic='NetworkTutorial', output_root=OUTPUT_ROOT, expected_count=14) + + + def _figure(label: str, *, figsize=(8.5, 4.5)): + fig = __tracker.new_figure(label) + fig.clear() + fig.set_size_inches(*figsize) + return fig + + + def _text_panel(fig, title: str, lines): + ax = fig.subplots(1, 1) + ax.axis("off") + ax.set_title(title) + ax.text( + 0.02, + 0.98, + "\\n".join(lines), + va="top", + ha="left", + family="monospace", + fontsize=11, + transform=ax.transAxes, + ) + return ax + + + def _stem_kernel(ax, coeffs, title: str, xlabel: str, color: str): + coeffs = np.asarray(coeffs, dtype=float).reshape(-1) + x = np.arange(1, coeffs.size + 1, dtype=float) + markerline, stemlines, baseline = ax.stem(x, coeffs, basefmt=" ") + plt.setp(markerline, color=color, markersize=8) + plt.setp(stemlines, color=color, linewidth=2.0) + baseline.set_visible(False) + ax.axhline(0.0, color="0.35", linewidth=1.0) + ax.set_xticks(x) + ax.set_xlabel(xlabel) + ax.set_ylabel("coefficient") + ax.set_title(title) + ax.grid(axis="y", alpha=0.2) + + + def _draw_network(ax, actual_network): + ax.set_title("Two-neuron connectivity diagram") + ax.axis("off") + positions = {1: (0.25, 0.5), 2: (0.75, 0.5)} + for idx, (x, y) in positions.items(): + ax.add_patch(Circle((x, y), 0.08, facecolor="#d9edf7", edgecolor="#2b5876", linewidth=2.0)) + ax.text(x, y, f"Neuron {idx}", ha="center", va="center", fontsize=10, weight="bold") + ax.text(x, y - 0.18, f"$\\\\mu_{idx}={-3}$", ha="center", va="center", fontsize=10) + ax.text(x, y + 0.18, "history: [-4, -2, -1]", ha="center", va="center", fontsize=9) + + arrow12 = FancyArrowPatch( + positions[1], + positions[2], + arrowstyle="-|>", + mutation_scale=20, + linewidth=2.2, + color="#1f77b4", + connectionstyle="arc3,rad=0.15", + ) + arrow21 = FancyArrowPatch( + positions[2], + positions[1], + arrowstyle="-|>", + mutation_scale=20, + linewidth=2.2, + color="#d62728", + connectionstyle="arc3,rad=-0.15", + ) + ax.add_patch(arrow12) + ax.add_patch(arrow21) + ax.text(0.5, 0.68, f"$E_1={actual_network[0, 1]:.0f}$", ha="center", va="center", fontsize=10, color="#1f77b4") + ax.text(0.5, 0.32, f"$E_2={actual_network[1, 0]:.0f}$", ha="center", va="center", fontsize=10, color="#d62728") + ax.text(0.25, 0.14, "$S_1=+1 \\cdot u_{stim}$", ha="center", fontsize=10) + ax.text(0.75, 0.14, "$S_2=-1 \\cdot u_{stim}$", ha="center", fontsize=10) + + + def _draw_block_diagram(ax): + ax.axis("off") + ax.set_title("Conditional-intensity block diagram") + labels = [ + ("Baseline $\\\\mu_i$", (0.10, 0.72)), + ("History $H * \\\\Delta N_i[n]$", (0.10, 0.48)), + ("Stimulus $S * u_{stim}[n]$", (0.10, 0.24)), + ("Ensemble $E * \\\\Delta N_k[n]$", (0.10, 0.02)), + ("Summation", (0.48, 0.36)), + ("logistic", (0.72, 0.36)), + ("Spike probability $\\\\lambda_i \\\\Delta$", (0.90, 0.36)), + ] + for text, (x, y) in labels: + box = FancyBboxPatch((x, y), 0.18, 0.14, boxstyle="round,pad=0.03", facecolor="white", edgecolor="0.3") + ax.add_patch(box) + ax.text(x + 0.09, y + 0.07, text, ha="center", va="center", fontsize=9) + for y in (0.79, 0.55, 0.31, 0.09): + ax.add_patch(FancyArrowPatch((0.28, y), (0.48, 0.43), arrowstyle="-|>", mutation_scale=15, linewidth=1.8, color="0.25")) + ax.add_patch(FancyArrowPatch((0.66, 0.43), (0.72, 0.43), arrowstyle="-|>", mutation_scale=15, linewidth=1.8, color="0.25")) + ax.add_patch(FancyArrowPatch((0.90, 0.43), (0.98, 0.43), arrowstyle="-|>", mutation_scale=15, linewidth=1.8, color="0.25")) + + + def _estimate_network(results): + estimated = np.zeros((2, 2), dtype=float) + for neuron_idx, fit in enumerate(results, start=1): + coeffs, labels, _ = fit.getCoeffsWithLabels(3) + for coeff, label in zip(coeffs, labels, strict=False): + label_str = str(label) + if neuron_idx == 1 and label_str.startswith("2:"): + estimated[0, 1] = float(coeff) + elif neuron_idx == 2 and label_str.startswith("1:"): + estimated[1, 0] = float(coeff) + return estimated + """, + """ + # SECTION 1: Point Process Network Simulation + # In order to understand how the point process GLM framework can be used to estimate the network connectivity within a population of neurons, we simulate a network of 2 neurons. + plt.close("all") + """, + """ + # SECTION 2: Published network diagram + fig = _figure("SimulatedNetwork2.png", figsize=(8.0, 4.5)) + ax = fig.subplots(1, 1) + _draw_network(ax, np.array([[0.0, 1.0], [-4.0, 0.0]], dtype=float)) + """, + """ + # SECTION 3: Published block diagram + fig = _figure("PPSimExample-BlockDiagram.png", figsize=(10.0, 4.5)) + ax = fig.subplots(1, 1) + _draw_block_diagram(ax) + """, + """ + # SECTION 4: Conditional intensity equation + fig = _figure("lambda_i * Delta = logistic(mu_i + H*DeltaN_i[n] + S*u_stim[n] + E*DeltaN_k[n])", figsize=(10.0, 3.0)) + _text_panel( + fig, + "Conditional intensity used in the tutorial", + [ + "lambda_i * Delta = logistic(mu_i + H * DeltaN_i[n]", + " + S * u_stim[n] + E * DeltaN_k[n])", + ], + ) + """, + """ + # SECTION 5: Logistic nonlinearity + # logistic(x) = exp(x) / (1 + exp(x)). Note that * is the convolution operator. + """, + """ + # SECTION 6: Convolution operator note + # The MATLAB helpfile presents the recursive history, stimulus, and ensemble filters separately below. + """, + """ + # SECTION 7: 2 Neuron Network + Ts = 0.001 + sampleRate = 1.0 / Ts + tMin = 0.0 + tMax = 50.0 + time = np.arange(tMin, tMax + Ts, Ts) + numNeurons = 2 + selfHist = [0.0, Ts, 2.0 * Ts, 3.0 * Ts] + ensHist = [0.0, Ts] + network = simulate_two_neuron_network(duration_s=tMax, dt=Ts, seed=4) + baseline_mu = np.asarray(network.baseline_mu, dtype=float) + history_kernel = np.asarray(network.history_kernel, dtype=float) + stim_kernel = np.asarray(network.stimulus_kernel, dtype=float) + ensemble_kernel = np.asarray(network.ensemble_kernel, dtype=float) + actual_network = np.asarray(network.actual_network, dtype=float) + """, + """ + # SECTION 8: Baseline firing rate of the neurons being modeled + print({"mu1": float(baseline_mu[0]), "mu2": float(baseline_mu[1]), "sample_rate_hz": sampleRate}) + """, + """ + # SECTION 9: History Effect + # Captures how the firing of a neuron modulates its own probability of firing, including the refractory period and short-term history dependence. + """, + """ + # SECTION 10: History kernel + fig = _figure("1*h[n]=-4*DeltaN[n-1]-2*DeltaN[n-2]-1*DeltaN[n-3]", figsize=(8.0, 4.5)) + ax = fig.subplots(1, 1) + _stem_kernel(ax, history_kernel, "Self-history kernel", "lag (ms)", "tab:red") + ax.set_xticklabels([f"{int(k)}" for k in [1, 2, 3]]) + """, + """ + # SECTION 11: Stimulus Effect + # Neuron 1 is positively modulated by the stimulus and neuron 2 is negatively modulated by the same drive. + """, + """ + # SECTION 12: Stimulus filter for neuron 1 + fig = _figure("1*s_1[n]=1*u_stim[n]", figsize=(7.5, 4.5)) + ax = fig.subplots(1, 1) + _stem_kernel(ax, [stim_kernel[0]], "Stimulus filter for neuron 1", "lag (samples)", "tab:blue") + """, + """ + # SECTION 13: Stimulus filter for neuron 2 + fig = _figure("1*s_2[n]=-1*u_stim[n]", figsize=(7.5, 4.5)) + ax = fig.subplots(1, 1) + _stem_kernel(ax, [stim_kernel[1]], "Stimulus filter for neuron 2", "lag (samples)", "tab:orange") + """, + """ + # SECTION 14: Ensemble Effect + # Captures how neighboring neuron firing modulates the firing probability of a given neuron, with a one-sample delay included in the Simulink model. + """, + """ + # SECTION 15: Ensemble filter for neuron 1 + fig = _figure("1*e_1[n]=1*DeltaN_2[n-1]", figsize=(7.5, 4.5)) + ax = fig.subplots(1, 1) + _stem_kernel(ax, [ensemble_kernel[0]], "Ensemble filter for neuron 1", "lag (samples)", "tab:green") + """, + """ + # SECTION 16: Ensemble filter for neuron 2 + fig = _figure("1*e_2[n]=-4*DeltaN_1[n-1]", figsize=(7.5, 4.5)) + ax = fig.subplots(1, 1) + _stem_kernel(ax, [ensemble_kernel[1]], "Ensemble filter for neuron 2", "lag (samples)", "tab:purple") + """, + """ + # SECTION 17: Stimulus + f = 1.0 + stimCov = Covariate(time, network.latent_drive, "Stimulus", "time", "s", "Voltage", ["sin"]) + baselineCov = Covariate(time, np.ones_like(time), "Baseline", "time", "s", "", ["mu"]) + + fig = _figure("Stimulus sine-wave drive", figsize=(9.0, 4.5)) + ax = fig.subplots(1, 1) + ax.plot(time, network.latent_drive, color="tab:blue", linewidth=1.2) + ax.set_xlim(0.0, 5.0) + ax.set_xlabel("time (s)") + ax.set_ylabel("stimulus") + ax.set_title("Sine-wave stimulus used by the network tutorial") + """, + """ + # SECTION 18: Simulate the Network + fitType = "binomial" + Algorithm = "BNLRCG" if fitType == "binomial" else "GLM" + spikeColl = network.spikes + trial = Trial(spikeColl, CovColl([stimCov, baselineCov]), None, History(selfHist)) + trial.setEnsCovHist(ensHist) + + fig = _figure("Simulated raster and stimulus", figsize=(9.0, 6.5)) + axs = fig.subplots(2, 1, sharex=True) + spikeColl.plot(handle=axs[0]) + axs[0].set_xlim(0.0, tMax / 10.0) + axs[0].set_title("Simulated spike trains") + stimCov.plot(handle=axs[1]) + axs[1].set_xlim(0.0, tMax / 10.0) + axs[1].set_title("Stimulus over the same window") + + fig = _figure("Simulated lambda traces", figsize=(9.0, 4.5)) + ax = fig.subplots(1, 1) + ax.plot(time, network.lambda_delta[:, 0], label="Neuron 1", linewidth=1.1) + ax.plot(time, network.lambda_delta[:, 1], label="Neuron 2", linewidth=1.1) + ax.set_xlim(0.0, 5.0) + ax.set_xlabel("time (s)") + ax.set_ylabel("lambda * Delta") + ax.set_title("Conditional-intensity trajectories") + ax.legend(loc="upper right") + """, + """ + # SECTION 19: GLM Model Fitting Setup + c1 = TrialConfig([["Baseline", "mu"]], sampleRate, [], [], [], name="Baseline") + c2 = TrialConfig([["Baseline", "mu"]], sampleRate, [], ensHist, [], name="Baseline+EnsHist") + c3 = TrialConfig([["Baseline", "mu"], ["Stimulus", "sin"]], sampleRate, selfHist, ensHist, [], name="Stim+Hist+EnsHist") + cfgColl = ConfigColl([c1, c2, c3]) + """, + """ + # SECTION 20: GLM Model Fitting and Results + results = Analysis.RunAnalysisForAllNeurons(trial, cfgColl, 0, Algorithm) + if not isinstance(results, list): + results = [results] + + fig = _figure("results{1}.plotResults", figsize=(11.0, 8.0)) + results[0].plotResults(handle=fig) + + fig = _figure("results{2}.plotResults", figsize=(11.0, 8.0)) + results[1].plotResults(handle=fig) + + estimated_network = _estimate_network(results) + fig = _figure("Actual vs estimated network", figsize=(10.0, 4.5)) + axs = fig.subplots(1, 2) + clim = float(np.max(np.abs(actual_network))) + for ax, matrix, title in zip( + axs, + [actual_network, estimated_network], + ["Actual", "Estimated 1 ms"], + strict=False, + ): + im = ax.imshow(matrix, cmap="jet", vmin=-clim, vmax=clim) + ax.set_xticks([0, 1], ["1", "2"]) + ax.set_yticks([0, 1], ["1", "2"]) + ax.set_title(title) + for (row, col), value in np.ndenumerate(matrix): + ax.text(col, row, f"{value:.2f}", ha="center", va="center", color="white", fontsize=10) + fig.colorbar(im, ax=axs, shrink=0.8) + """, + """ + # SECTION 21: Neighbor-selection note + # By default all neurons are considered potential neighbors. To restrict candidate neighbors, call trial.setNeighbors(neighborArray) using the MATLAB-style convention described in the source helpfile. + print( + { + "algorithm": Algorithm, + "spike_counts": [spikeColl.getNST(1).n_spikes, spikeColl.getNST(2).n_spikes], + "estimated_network": np.round(estimated_network, 3).tolist(), + } + ) + __tracker.finalize() + """, +] + + +def build_notebook(): + return new_notebook( + cells=[ + new_markdown_cell(NOTE), + *[new_code_cell(dedent(cell).strip() + "\n") for cell in CODE], + ], + metadata={ + **LANGUAGE_METADATA, + "nstat": { + "expected_figures": 14, + "run_group": "full", + "style": "python-example", + "topic": "NetworkTutorial", + }, + }, + ) + + +def main() -> int: + notebook = build_notebook() + NOTEBOOK_PATH.parent.mkdir(parents=True, exist_ok=True) + NOTEBOOK_PATH.write_text(nbformat.writes(notebook), encoding="utf-8") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/notebooks/parity_notes.yml b/tools/notebooks/parity_notes.yml index 35d6501a..f0601f62 100644 --- a/tools/notebooks/parity_notes.yml +++ b/tools/notebooks/parity_notes.yml @@ -50,6 +50,11 @@ notes: source_matlab: PPSimExample.mlx fidelity_status: high_fidelity remaining_differences: The notebook now follows the MATLAB recursive-CIF workflow with the native Python `CIF.simulateCIF` path; exact Simulink block timing and solver semantics are still not fixture-matched one-for-one against MATLAB. + - topic: NetworkTutorial + file: notebooks/NetworkTutorial.ipynb + source_matlab: NetworkTutorial.mlx + fidelity_status: high_fidelity + remaining_differences: The notebook now mirrors the MATLAB helpfile section order and published figure inventory with a native Python network simulator and MATLAB-style `Analysis` workflow; exact spike realizations still vary modestly because NumPy and Simulink do not share identical random streams. - topic: ValidationDataSet file: notebooks/ValidationDataSet.ipynb source_matlab: ValidationDataSet.mlx diff --git a/tools/notebooks/topic_groups.yml b/tools/notebooks/topic_groups.yml index 8445ee8f..a73d568d 100644 --- a/tools/notebooks/topic_groups.yml +++ b/tools/notebooks/topic_groups.yml @@ -42,6 +42,7 @@ groups: - ExplicitStimulusWhiskerData - HippocampalPlaceCellExample - HybridFilterExample + - NetworkTutorial - PPSimExample - SignalObjExamples - StimulusDecode2D @@ -58,6 +59,7 @@ groups: - ExplicitStimulusWhiskerData - HippocampalPlaceCellExample - HybridFilterExample + - NetworkTutorial - PPSimExample - StimulusDecode2D - TrialExamples From 3c5add608aa4e84725f6ec09572cbd67b358faf3 Mon Sep 17 00:00:00 2001 From: Iahn Cajigas Date: Sat, 7 Mar 2026 19:07:30 -0500 Subject: [PATCH 2/2] Avoid data download in example05 paper export --- nstat/paper_example_catalog.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/nstat/paper_example_catalog.py b/nstat/paper_example_catalog.py index e030e245..dfc997dd 100644 --- a/nstat/paper_example_catalog.py +++ b/nstat/paper_example_catalog.py @@ -24,19 +24,21 @@ def run_named_paper_example( example_id: str, repo_root: Path, *, return_payload: bool = False ) -> dict[str, dict[str, float]] | tuple[dict[str, dict[str, float]], dict[str, dict[str, object]]]: repo_root = repo_root.resolve() - data_dir = ensure_example_data(download=True) if example_id == "example01": + data_dir = ensure_example_data(download=True) if not return_payload: return {"experiment1": run_experiment1(data_dir)} summary, payload = run_experiment1(data_dir, return_payload=True) return {"experiment1": summary}, {"experiment1": payload} if example_id == "example02": + data_dir = ensure_example_data(download=True) if not return_payload: return {"experiment2": run_experiment2(data_dir)} summary, payload = run_experiment2(data_dir, return_payload=True) return {"experiment2": summary}, {"experiment2": payload} if example_id == "example03": + data_dir = ensure_example_data(download=True) if not return_payload: return { "experiment3": run_experiment3(), @@ -52,6 +54,7 @@ def run_named_paper_example( "experiment3b": payload3b, } if example_id == "example04": + data_dir = ensure_example_data(download=True) if not return_payload: return {"experiment4": run_experiment4(data_dir)} summary4, payload4 = run_experiment4(data_dir, return_payload=True)