RNG accessor refactor, FR_std NaN guard, optional elephant, figure distribution plots

.github/workflows/tests.yml

@@ -0,0 +1,39 @@
+name: tests
+
+on:
+  pull_request:
+    branches: [main, improved-semg]
+  push:
+    branches: [main, improved-semg]
+  workflow_dispatch:
+
+jobs:
+  pytest:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python 3.12
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.12.8"
+
+      - name: System packages for NEURON / MPI
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y libopenmpi-dev openmpi-bin libgtk2.0-dev
+
+      - name: Install uv
+        run: |
+          python -m pip install --upgrade pip
+          python -m pip install uv
+
+      - name: Install project with dev group
+        run: |
+          uv sync --group dev
+
+      - name: Compile NEURON mechanisms
+        run: uv run poe setup_myogen
+
+      - name: Run pytest
+        run: uv run pytest tests/ -v

.gitignore

@@ -175,10 +175,14 @@ test_*.png
 
 myogen/simulator/neuron/_cython/*.c
 myogen/simulator/nmodl_files/x86_64/
+myogen/simulator/nmodl_files/arm64/
 myogen/simulator/x86_64/
+myogen/simulator/arm64/
 myogen/simulator/nmodl_files/*.c
 myogen/simulator/nmodl_files/*.o
 myogen/simulator/nmodl_files/*.dll
+myogen/simulator/nmodl_files/*.dylib
+myogen/simulator/nmodl_files/*.so
 
 data/
 # Note: examples/data/ is NOT ignored - it contains pre-generated static data for docs

CHANGELOG.md

@@ -7,6 +7,34 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+- **RNG accessors**: new public API for reproducible random draws
+  - `myogen.get_random_generator()` — always returns the current global RNG (tracks the latest `set_random_seed` call)
+  - `myogen.get_random_seed()` — returns the seed currently in effect
+  - `myogen.derive_subseed(*labels)` — deterministic sub-seed helper for seeding non-NumPy generators (Cython spike generators, sklearn `random_state`, etc.) so they also track `set_random_seed`
+- **Optional `elephant` extra**: install via `pip install myogen[elephant]`. The extra bundles both `elephant>=1.1.1` and `viziphant>=0.4.0` (viziphant was moved out of core dependencies because its only transitive import chain in core came from `elephant`). Core install is now genuinely elephant-free — verified with `importlib.util.find_spec('elephant') is None` on a fresh `uv sync`. The five modules that import elephant (`utils/helper.py`, `surface_emg.py`, `intramuscular_emg.py`, `force_model.py`, `force_model_vectorized.py`) were already guarded with `try/except ImportError` and now emit an accurate install hint when the extra is missing. `neo>=0.14.0` is now an explicit core dependency — it was previously only reached transitively through `elephant`/`viziphant`, so removing those would otherwise break `import myogen`
+- **Test suite** under `tests/`
+  - `test_determinism.py` — 10 regressions covering seed propagation, sub-seed collision avoidance, and deprecation warnings for the legacy RNG names
+  - `test_firing_rate_statistics.py` — 4 regressions covering `FR_std` and `CV_ISI` NaN edge cases
+  - Added `pytest>=8.0` to the `dev` dependency-group and a `[tool.pytest.ini_options]` section
+- **Per-muscle ISI/CV figure**: new `plot_cv_vs_fr_per_muscle()` in `examples/02_finetune/05_plot_isi_cv_multi_muscle_comparison.py` emits a 3-panel VL / VM / FDI breakdown alongside the existing pooled plot
+
+### Changed
+- **RNG architecture** (internal refactor): the six internal modules and eight example scripts that previously did `from myogen import RANDOM_GENERATOR, SEED` at module scope now call the accessor at each site, eliminating the stale-reference bug where `set_random_seed` didn't reach downstream modules or seed-derived generators
+- **Per-cell seed derivation**: the three sites in `myogen/simulator/neuron/cells.py` that built a per-cell seed as `SEED + (class_id+1)*(global_id+1)` — a formula that collided on swapped factors, e.g. `(0, 5)` and `(1, 2)` both resolved to `+6` — now use `derive_subseed(class_id, global_id)` (collision-free with respect to label order, built on `numpy.random.SeedSequence`). The `motor_unit_sim.py` `KMeans` `random_state` uses the same helper. **Consequence**: for a given global seed, RNG output differs from earlier releases; byte-identical reproduction of pre-`Unreleased` outputs is not possible without matching code
+- **Example figures**: bar graphs in `02_finetune/01_optimize_dd_for_target_firing_rate.py`, `02_finetune/02_compute_force_from_optimized_dd.py` and `03_papers/watanabe/01_compute_baseline_force.py` replaced with dumbbell / violin + box + jittered-scatter plots that show the underlying distribution (all points when n<10, violin + box when n≥10)
+- **NEURON version alignment**: `README.md`, `docs/source/README.md`, `docs/source/index.md`, and `setup.py` Windows-install messaging now reference NEURON **8.2.7** (matching the Linux/macOS pip pin in `pyproject.toml` and the CI workflow) with the correct installer filename (`py-39-310-311-312-313`)
+- **Docs build**: pinned `sphinx<9` in the `docs` dependency-group to work around a `sphinx-hoverxref` regression on Sphinx 9, and moved `[tool.pytest.ini_options]` in `pyproject.toml` so it no longer re-parents the `docs` dependency-group
+
+### Deprecated
+- **`myogen.RANDOM_GENERATOR`** and **`myogen.SEED`** module attributes. They remain accessible via a module-level `__getattr__` that emits a `DeprecationWarning` and returns the current RNG / current seed. External code should migrate to `get_random_generator()` / `get_random_seed()`; module-level `from myogen import RANDOM_GENERATOR` captures a stale reference and does not reflect later `set_random_seed` calls
+
+### Fixed
+- **`FR_std = NaN` with a single active unit**: `myogen/utils/helper.py::calculate_firing_rate_statistics` now returns `FR_std = 0.0` when fewer than two units pass the firing-rate filter, instead of propagating the `np.std(..., ddof=1)` NaN into downstream ensemble statistics
+- **`CV_ISI = NaN` for n=1 ISI**: same function's per-neuron branch now returns `0.0` when `min_spikes_for_cv` is set to 2 and a neuron has exactly two spikes
+- **`set_random_seed` did not propagate**: module-level `from myogen import RANDOM_GENERATOR, SEED` imports in the internal modules captured a stale reference at import time; reseeding the package rebinding did not affect them. The accessor refactor above is the fix
+- **Typos / placeholder docstrings**: `extandable` → `extensible` in README and docs; `"as determined by no one"` docstring comments in `myogen/simulator/core/muscle/muscle.py` replaced with concrete physiological rationale
+
 ## [0.8.5] - 2026-01-15
 
 ### Fixed

README.md

@@ -4,7 +4,7 @@
     <img src="https://raw.githubusercontent.com/NsquaredLab/MyoGen/main/docs/source/_static/myogen_logo.png" height="100" alt="MyoGen Logo">
   </h1>
 
-  <h2>The modular and extandable simulation toolkit for neurophysiology</h2>
+  <h2>The modular and extensible simulation toolkit for neurophysiology</h2>
 
   [![Documentation](https://img.shields.io/badge/docs-latest-blue.svg)](https://nsquaredlab.github.io/MyoGen/)
   [![Python 3.12+](https://img.shields.io/badge/python-3.12+-blue.svg)](https://www.python.org/downloads/)
@@ -47,7 +47,7 @@ MyoGen is designed for algorithm validation, hypothesis-driven research, and edu
 
 | Platform | Before Installing MyoGen |
 |----------|--------------------------|
-| **Windows** | [NEURON 8.2.6](https://github.com/neuronsimulator/nrn/releases/download/8.2.6/nrn-8.2.6.w64-mingw-py-38-39-310-311-312-setup.exe) - Download, run installer, select "Add to PATH" |
+| **Windows** | [NEURON 8.2.7](https://github.com/neuronsimulator/nrn/releases/download/8.2.7/nrn-8.2.7.w64-mingw-py-39-310-311-312-313-setup.exe) - Download, run installer, select "Add to PATH" |
 | **Linux** | `sudo apt install libopenmpi-dev` (Ubuntu/Debian) or `sudo dnf install openmpi-devel` (Fedora) |
 | **macOS** | `brew install open-mpi` |
 
@@ -70,7 +70,7 @@ MyoGen is designed for algorithm validation, hypothesis-driven research, and edu
 >
 > ### 2. NEURON Simulator
 >
-> 1. **Download**: [NEURON 8.2.6 Installer](https://github.com/neuronsimulator/nrn/releases/download/8.2.6/nrn-8.2.6.w64-mingw-py-38-39-310-311-312-setup.exe)
+> 1. **Download**: [NEURON 8.2.7 Installer](https://github.com/neuronsimulator/nrn/releases/download/8.2.7/nrn-8.2.7.w64-mingw-py-39-310-311-312-313-setup.exe)
 > 2. **Run the installer** and select **"Add to PATH"** when prompted
 > 3. **Restart your terminal** (close and reopen)
 > 4. Then continue with the installation below

docs/source/README.md

@@ -13,7 +13,7 @@
   [Paper](#citation)
 </div>
 
-# MyoGen - The modular and extandable simulation toolkit for neurophysiology
+# MyoGen - The modular and extensible simulation toolkit for neurophysiology
 
 MyoGen is a **modular and extensible neuromuscular simulation framework** for generating physiologically grounded motor-unit activity, muscle force, and surface EMG signals.  
 
@@ -39,7 +39,7 @@ MyoGen is designed for algorithm validation, hypothesis-driven research, and edu
 **Prerequisites**: Python ≥3.12, Linux/Windows/macOS
 
 > [!IMPORTANT]
-> **Windows users**: Install [NEURON 8.2.6](https://github.com/neuronsimulator/nrn/releases/download/8.2.6/nrn-8.2.6.w64-mingw-py-38-39-310-311-312-setup.exe) before running `uv sync`
+> **Windows users**: Install [NEURON 8.2.7](https://github.com/neuronsimulator/nrn/releases/download/8.2.7/nrn-8.2.7.w64-mingw-py-39-310-311-312-313-setup.exe) before running `uv sync`
 
 ```bash
 # Clone and install

docs/source/index.md

@@ -4,7 +4,7 @@
     <img src="https://raw.githubusercontent.com/NsquaredLab/MyoGen/main/docs/source/_static/myogen_logo.png" height="100" alt="MyoGen Logo">
   </h1>
 
-  <h2>The modular and extandable simulation toolkit for neurophysiology</h2>
+  <h2>The modular and extensible simulation toolkit for neurophysiology</h2>
 
   [![Documentation](https://img.shields.io/badge/docs-latest-blue.svg)](https://nsquaredlab.github.io/MyoGen/)
   [![Python 3.12+](https://img.shields.io/badge/python-3.12+-blue.svg)](https://www.python.org/downloads/)
@@ -47,7 +47,7 @@ MyoGen is designed for algorithm validation, hypothesis-driven research, and edu
 
 | Platform | Before Installing MyoGen |
 |----------|--------------------------|
-| **Windows** | [NEURON 8.2.6](https://github.com/neuronsimulator/nrn/releases/download/8.2.6/nrn-8.2.6.w64-mingw-py-38-39-310-311-312-setup.exe) - Download, run installer, select "Add to PATH" |
+| **Windows** | [NEURON 8.2.7](https://github.com/neuronsimulator/nrn/releases/download/8.2.7/nrn-8.2.7.w64-mingw-py-39-310-311-312-313-setup.exe) - Download, run installer, select "Add to PATH" |
 | **Linux** | `sudo apt install libopenmpi-dev` (Ubuntu/Debian) or `sudo dnf install openmpi-devel` (Fedora) |
 | **macOS** | `brew install open-mpi` |
 
@@ -70,7 +70,7 @@ MyoGen is designed for algorithm validation, hypothesis-driven research, and edu
 >
 > ### 2. NEURON Simulator
 >
-> 1. **Download**: [NEURON 8.2.6 Installer](https://github.com/neuronsimulator/nrn/releases/download/8.2.6/nrn-8.2.6.w64-mingw-py-38-39-310-311-312-setup.exe)
+> 1. **Download**: [NEURON 8.2.7 Installer](https://github.com/neuronsimulator/nrn/releases/download/8.2.7/nrn-8.2.7.w64-mingw-py-39-310-311-312-313-setup.exe)
 > 2. **Run the installer** and select **"Add to PATH"** when prompted
 > 3. **Restart your terminal** (close and reopen)
 > 4. Then continue with the installation below

examples/01_basic/02_simulate_spike_trains_current_injection.py

@@ -20,17 +20,17 @@
 # ----------------
 #
 # .. important::
-#    In **MyoGen** all **random number generation** is handled by the :data:`~myogen.RANDOM_GENERATOR` object.
+#    In **MyoGen** all **random number generation** is handled by the RNG returned from
+#    :func:`~myogen.get_random_generator`, a thin wrapper around :mod:`numpy.random`.
 #
-#    This object is a wrapper around the :mod:`numpy.random` module and is used to generate random numbers.
-#
-#    It is intended to be used with the following API:
+#    Always fetch the generator at the call site so the current seed is honored:
 #
 #    .. code-block:: python
 #
-#       from myogen import simulator, RANDOM_GENERATOR
+#       from myogen import simulator, get_random_generator
+#       get_random_generator().normal(0, 1)
 #
-#    To change the default seed, use :func:`~myogen.set_random_seed`:
+#    To change the seed, use :func:`~myogen.set_random_seed`:
 #
 #    .. code-block:: python
 #
@@ -50,7 +50,7 @@
 from neuron import h
 from viziphant.rasterplot import rasterplot_rates
 
-from myogen import RANDOM_GENERATOR
+from myogen import get_random_generator
 from myogen.simulator.neuron.populations import AlphaMN__Pool
 from myogen.utils.currents import create_trapezoid_current
 from myogen.utils.neuron.inject_currents_into_populations import (
@@ -107,9 +107,9 @@
 timestep = 0.05 * pq.ms
 simulation_time = 4000 * pq.ms
 
-rise_time_ms = list(RANDOM_GENERATOR.uniform(100, 500, size=n_pools)) * pq.ms
-plateau_time_ms = list(RANDOM_GENERATOR.uniform(1000, 2000, size=n_pools)) * pq.ms
-fall_time_ms = list(RANDOM_GENERATOR.uniform(1000, 2000, size=n_pools)) * pq.ms
+rise_time_ms = list(get_random_generator().uniform(100, 500, size=n_pools)) * pq.ms
+plateau_time_ms = list(get_random_generator().uniform(1000, 2000, size=n_pools)) * pq.ms
+fall_time_ms = list(get_random_generator().uniform(1000, 2000, size=n_pools)) * pq.ms
 
 input_current__AnalogSignal = create_trapezoid_current(
     n_pools,

examples/01_basic/03_simulate_spike_trains_descending_drive.py

@@ -28,17 +28,17 @@
 # ----------------
 #
 # .. important::
-#    In **MyoGen** all **random number generation** is handled by the ``RANDOM_GENERATOR`` object.
+#    In **MyoGen** all **random number generation** is handled by the RNG returned from
+#    ``get_random_generator()``, a thin wrapper around ``numpy.random``.
 #
-#    This object is a wrapper around the ``numpy.random`` module and is used to generate random numbers.
-#
-#    It is intended to be used with the following API:
+#    Always fetch the generator at the call site so the current seed is honored:
 #
 #    .. code-block:: python
 #
-#       from myogen import simulator, RANDOM_GENERATOR
+#       from myogen import simulator, get_random_generator
+#       get_random_generator().normal(0, 1)
 #
-#    To change the default seed, use ``set_random_seed``:
+#    To change the seed, use ``set_random_seed``:
 #
 #    .. code-block:: python
 #
@@ -59,7 +59,7 @@
 from neuron import h
 from tqdm import tqdm
 
-from myogen import RANDOM_GENERATOR
+from myogen import get_random_generator
 from myogen.simulator.neuron import Network
 from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
 from myogen.utils.nmodl import load_nmodl_mechanisms
@@ -161,7 +161,7 @@
 
 # Add small noise for realism
 trapezoid_drive = (
-    trapezoid_drive + np.clip(RANDOM_GENERATOR.normal(0, 1.0, size=time_points), 0, None) * pps
+    trapezoid_drive + np.clip(get_random_generator().normal(0, 1.0, size=time_points), 0, None) * pps
 )
 
 # Create AnalogSignal

examples/02_finetune/01_optimize_dd_for_target_firing_rate.py

@@ -56,7 +56,7 @@
 from neo import Segment, SpikeTrain
 from neuron import h
 
-from myogen import RANDOM_GENERATOR
+from myogen import get_random_generator
 from myogen.simulator import RecruitmentThresholds
 from myogen.simulator.neuron import Network
 from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
@@ -224,7 +224,7 @@ def objective(trial):
         # Generate constant drive signal with small noise
         time_points = int(SIMULATION_TIME_MS / TIMESTEP_MS)
         drive_signal = np.ones(time_points) * dd_drive__Hz + np.clip(
-            RANDOM_GENERATOR.normal(0, 1.0, size=time_points), 0, None
+            get_random_generator().normal(0, 1.0, size=time_points), 0, None
         )
 
         # Run NEURON simulation
@@ -504,7 +504,9 @@ def trial_to_dict(t):
 
 fig, ax = plt.subplots(1, 1, figsize=(10, 6))
 
-# Bar plot comparing target vs achieved
+# Dumbbell plot comparing target vs achieved (two scalar endpoints per
+# category — not a distribution, so Editor 8's "show all data points"
+# rule is honoured by marking each value as an explicit point).
 categories = ["Mean FR (Hz)", "Std FR (Hz)"]
 targets = [TARGET_FR_MEAN__HZ, TARGET_FR_STD__HZ]
 achieved = [
@@ -513,22 +515,16 @@ def trial_to_dict(t):
 ]
 
 x = np.arange(len(categories))
-width = 0.35
-
-bars1 = ax.bar(x - width / 2, targets, width, label="Target")
-bars2 = ax.bar(x + width / 2, achieved, width, label="Achieved")
-
-# Add value labels on bars
-for bars in [bars1, bars2]:
-    for bar in bars:
-        height = bar.get_height()
-        ax.text(
-            bar.get_x() + bar.get_width() / 2.0,
-            height,
-            f"{height:.2f}",
-            ha="center",
-            va="bottom",
-        )
+
+for xi, t, a in zip(x, targets, achieved):
+    ax.plot([xi, xi], [t, a], color="gray", linestyle="--", alpha=0.6, zorder=1)
+
+ax.scatter(x, targets, s=120, marker="o", label="Target", zorder=3)
+ax.scatter(x, achieved, s=120, marker="D", label="Achieved", zorder=3)
+
+for xi, t, a in zip(x, targets, achieved):
+    ax.text(xi + 0.05, t, f"{t:.2f}", va="center", ha="left")
+    ax.text(xi + 0.05, a, f"{a:.2f}", va="center", ha="left")
 
 # Calculate percent errors
 mean_error = abs(achieved[0] - targets[0]) / targets[0] * 100

examples/02_finetune/02_compute_force_from_optimized_dd.py

@@ -58,7 +58,7 @@
 from neo import Block, Segment, SpikeTrain
 from neuron import h
 
-from myogen import RANDOM_GENERATOR
+from myogen import get_random_generator
 from myogen.simulator import RecruitmentThresholds
 from myogen.simulator.core.force.force_model import ForceModel
 from myogen.simulator.neuron import Network
@@ -215,7 +215,7 @@
 
 time_points = int(SIMULATION_TIME_MS / TIMESTEP_MS)
 drive_signal = np.ones(time_points) * dd_drive__Hz + np.clip(
-    RANDOM_GENERATOR.normal(0, 1.0, size=time_points), 0, None
+    get_random_generator().normal(0, 1.0, size=time_points), 0, None
 )
 
 ##############################################################################
@@ -405,15 +405,33 @@
         mu_ids.append(i)
 
 if firing_rates:
-    axes[2].bar(mu_ids, firing_rates, alpha=0.7)
-    axes[2].axhline(fr_mean, linestyle="--", label=f"Mean = {fr_mean:.1f} Hz")
-    axes[2].set_xlabel("Motor Unit ID")
+    # Distribution plot with all data points (Editor 8): violin shell for
+    # the density, box-and-whisker for the quartiles, jittered points for
+    # every active motor unit.
+    fr_array = np.asarray(firing_rates)
+    fr_sd = float(np.std(fr_array, ddof=1)) if fr_array.size > 1 else 0.0
+    if fr_array.size >= 10:
+        axes[2].violinplot(
+            fr_array, positions=[0], widths=0.7, showmeans=False, showmedians=False
+        )
+        axes[2].boxplot(
+            fr_array, positions=[0], widths=0.3, showfliers=False
+        )
+    jitter = np.random.default_rng(0).uniform(-0.15, 0.15, size=fr_array.size)
+    axes[2].scatter(jitter, fr_array, alpha=0.6, s=18, zorder=3)
+    axes[2].axhline(
+        fr_mean, linestyle="--", label=f"Mean = {fr_mean:.1f} Hz (SD = {fr_sd:.1f} Hz)"
+    )
+    axes[2].set_xticks([0])
+    axes[2].set_xticklabels([f"n = {fr_array.size}"])
+    axes[2].set_xlabel("Motor Units")
     axes[2].set_ylabel("Firing Rate (Hz)")
     axes[2].set_title("Firing Rate Distribution")
     axes[2].legend(framealpha=1.0, edgecolor="none")
 
-# Format all axes
-for ax in axes:
+# Format time-series axes (subplot 2 is the MU firing-rate distribution,
+# whose x-axis is categorical, not time)
+for ax in axes[:2]:
     ax.set_xlim(0, SIMULATION_TIME_MS / 1000)
 
 plt.tight_layout()

examples/02_finetune/03_optimize_dd_for_target_force.py

@@ -50,7 +50,7 @@
 from neo import Block, Segment, SpikeTrain
 from neuron import h
 
-from myogen import RANDOM_GENERATOR, set_random_seed
+from myogen import get_random_generator, set_random_seed
 from myogen.simulator import RecruitmentThresholds
 from myogen.simulator.core.force.force_model import ForceModel
 from myogen.simulator.neuron import Network
@@ -212,7 +212,7 @@ def run_simulation_and_compute_force(dd_drive__Hz, gamma_shape, recruitment_thre
     # Create constant drive signal
     time_points = int(SIMULATION_TIME_MS / TIMESTEP_MS)
     drive_signal = np.ones(time_points) * dd_drive__Hz + np.clip(
-        RANDOM_GENERATOR.normal(0, 1.0, size=time_points), 0, None
+        get_random_generator().normal(0, 1.0, size=time_points), 0, None
     )
 
     # Initialize simulation