Skip to content

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

Closed
RaulSimpetru wants to merge 17 commits intomainfrom
reviewer_suggestions
Closed

RNG accessor refactor, FR_std NaN guard, optional elephant, figure distribution plots#5
RaulSimpetru wants to merge 17 commits intomainfrom
reviewer_suggestions

Conversation

@RaulSimpetru
Copy link
Copy Markdown
Member

@RaulSimpetru RaulSimpetru commented Apr 18, 2026

Summary

Maintenance pass on the MyoGen package: RNG architecture, NaN-guarding a firing-rate statistic, moving elephant to an optional extra, aligning NEURON version references, and reshaping several example figures.

Key changes

  • RNG refactormyogen.get_random_generator() / myogen.get_random_seed() / myogen.derive_subseed(*labels) replace module-level RANDOM_GENERATOR / SEED imports. Legacy names still resolve via __getattr__ with a DeprecationWarning. All 6 internal modules and 8 examples/ scripts migrated. derive_subseed replaces the old non-injective SEED + (class_id+1)*(global_id+1) pattern at 3 cells.py sites and the motor_unit_sim.py KMeans site.
  • Deterministic tests — new tests/test_determinism.py (10 cases) and tests/test_firing_rate_statistics.py (4 cases) covering byte-level reproducibility, seed propagation, deprecated-attribute warnings, sub-seed collision avoidance, and single-unit NaN guards. Added pytest>=8.0 to the dev group.
  • FR_std NaN guardcalculate_firing_rate_statistics now returns FR_std = 0.0 when fewer than two units are active, rather than NaN from np.std(..., ddof=1). A parallel guard in the per-neuron CV_ISI path handles the n=1 ISI edge case.
  • elephant optional — moved from [project.dependencies] to [project.optional-dependencies].elephant. The five internal sites that import it were already guarded with try/except ImportError; their "Install with: pip install myogen[elephant]" message is now accurate.
  • NEURON 8.2.7 — README, docs, and setup.py all point at 8.2.7 with the correct Windows installer filename (py-39-310-311-312-313).
  • 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 converted to dumbbell / violin + box + jittered scatter (all data points shown when n<10, box/violin when n≥10). New plot_cv_vs_fr_per_muscle() in 02_finetune/05_plot_isi_cv_multi_muscle_comparison.py emits a three-panel VL / VM / FDI breakdown of the pooled ISI/CV comparison; the pooled plot itself is unchanged.
  • Docs buildsphinx<9 pinned in the docs dep-group to work around a sphinx-hoverxref regression on Sphinx 9; the [tool.pytest.ini_options] table was moved so it no longer re-parents the docs group.

The previous design bound RANDOM_GENERATOR and SEED at module import time
via `from myogen import RANDOM_GENERATOR, SEED`. Python's from-import copies
the current binding into the importing module's namespace — so when
set_random_seed() later rebinds myogen.RANDOM_GENERATOR, existing importers
keep their stale reference. The SEED constant was never updated at all,
so seeds derived from it (cells.py Cython Mersenne seeds, motor_unit_sim.py
KMeans random_state) were frozen to 180319 for the life of the process.

Fix: introduce get_random_generator() and get_random_seed() accessors in
myogen/__init__.py. A function lookup always reads the current module
global, so subsequent set_random_seed() calls propagate correctly. The
RANDOM_GENERATOR and SEED module attributes still resolve via a module-
level __getattr__ that emits DeprecationWarning and returns the current
state, preserving backwards compatibility for external callers.

Updated sites:
 - myogen/simulator/neuron/network.py (5 RNG usages)
 - myogen/simulator/neuron/cells.py (3 SEED-derived Cython seeds + 1 RNG draw)
 - myogen/simulator/core/muscle/muscle.py (3 RNG usages)
 - myogen/simulator/core/emg/surface/surface_emg.py (3 RNG usages)
 - myogen/simulator/core/emg/intramuscular/motor_unit_sim.py
   (1 KMeans random_state, 2 rng=RANDOM_GENERATOR caches, 1 draw)
 - myogen/simulator/core/emg/intramuscular/intramuscular_emg.py (1 RNG draw)

Cell-specific seeds now use `get_random_seed() + (class_id+1)*(global_id+1)`
— same mathematical structure as before, but the base seed tracks
set_random_seed().

Tests: added tests/test_determinism.py with 8 regression tests covering
byte-level reproducibility, seed propagation to downstream seed-derivation
sites, and DeprecationWarning emission for the legacy attribute names.
Added pytest>=8.0 to the dev dependency group.
Follow-up to the previous commit. Codex review flagged:

1) MEDIUM: cells.py seed derivation `get_random_seed() + (class_id+1)*(global_id+1)`
   is not injective — e.g. (0, 5) and (1, 2) both add 6. Added a
   collision-free `derive_subseed(*labels)` helper in myogen/__init__.py
   built on `numpy.random.SeedSequence`. All three cells.py sites and the
   motor_unit_sim.py KMeans `random_state` now use this helper.

2) MEDIUM: examples/02_finetune/03_optimize_dd_for_target_force.py and
   examples/03_papers/watanabe/02_optimize_oscillating_dc.py imported
   RANDOM_GENERATOR *and* later called set_random_seed(42) — the exact
   stale-reference bug leaking into examples. Migrated both to
   get_random_generator().

3) LOW: migrated the remaining six example scripts and their narrative
   docstring blocks from RANDOM_GENERATOR / SEED to get_random_generator(),
   so examples no longer teach the deprecated API or emit
   DeprecationWarnings at import.

4) LOW: added two regression tests:
   - test_derive_subseed_is_collision_free_for_distinct_labels — verifies
     the 8x8 sub-seed grid is injective, catching any future regression
     to a non-injective mixing function.
   - test_derive_subseed_tracks_set_random_seed — verifies derived
     sub-seeds change with the global seed.
Codex re-review noted that derive_subseed compresses inputs to a uint32,
so collisions remain possible in principle (birthday probability
N^2 / 2^33, ~1e-7 at 1000 cells). Reworded the docstring to promise
"deterministic, seed-tracking" rather than "collision-free", documented
the quantitative collision bound, and renamed the regression test to
test_derive_subseed_avoids_old_collision_pattern with an explicit check
on the (0,5) vs (1,2) pair that the old formula mangled. Also extended
coverage to a 16x16 label grid.
ddof=1 is undefined for n<2 (division by n-1 = 0) and returns NaN,
which could silently corrupt downstream ensemble-statistics analysis
when a pool has exactly one active unit after firing-rate filtering.

Guards the sample-std branch in myogen/utils/helper.py:175 to return
0.0 when fewer than two units are active. The empty-population path
(n=0) already returned 0.0 and is unchanged.

Adds tests/test_firing_rate_statistics.py with three cases:
 - single active unit -> FR_std is 0.0 (regression for the NaN case)
 - empty population -> FR_std stays 0.0
 - two distinct rates -> FR_std is finite and positive (sanity)
Codex noted a second NaN landmine in the same function: when a caller
sets return_per_neuron=True with min_spikes_for_cv=2, a neuron with
exactly two spikes yields one ISI, and np.std(..., ddof=1) over that
single element is NaN. Guarded with len(isis_array) > 1; falls back
to CV_ISI = 0.0 otherwise. Added a regression test
test_per_neuron_cv_is_finite_with_min_spikes_for_cv_2.
1) elephant is now a proper optional extra.

   - Removed elephant>=1.1.1 from [project.dependencies].
   - Added elephant = ["elephant>=1.1.1"] to [project.optional-dependencies].
   - The five code sites that import elephant
     (myogen/utils/helper.py, surface_emg.py, intramuscular_emg.py,
     force_model.py, force_model_vectorized.py) already wrap the import
     in a try/except ImportError and print a helpful
     "Install with: pip install myogen[elephant]" message when missing;
     those guards now match reality.

2) NEURON version unified to 8.2.7 in the Windows install docs.

   - README.md lines 50, 73: bumped 8.2.6 -> 8.2.7 (both installer
     link and version label).
   - docs/source/README.md line 42 and docs/source/index.md lines
     50, 73: same bump.
   - Corrected the Python-version stub in the Windows installer
     filename to match the 8.2.7 release: py-39-310-311-312-313
     (8.2.7 dropped 3.8 and added 3.13 compared to 8.2.6's
     py-38-39-310-311-312). Verified via
     `gh api repos/neuronsimulator/nrn/releases/tags/8.2.7`.
   - pyproject.toml and the CI workflow already pinned 8.2.7 for
     Linux/macOS; no change needed there.
Codex spotted remaining 8.2.6 references in setup.py (lines 204, 208)
in the Windows-install failure message. Bumped version string and
Windows installer URL to 8.2.7 (py-39-310-311-312-313 stub).
 - "extandable" -> "extensible" in README.md, docs/source/index.md,
   docs/source/README.md.
 - muscle.py docstring (two places) replaces the informal
   "as determined by no one" placeholders with concrete rationale:
   the 30 mm default length and the 0.25 max territory-to-muscle
   ratio are nominal starting points for the FDI; users are told to
   revisit them for other muscles.

Broader repo-wide grep for common English typos (recieve, teh,
seperate, occured, lenght, etc.) found no additional hits.
the underlying data rather than summary-only bars .

 - examples/02_finetune/01_optimize_dd_for_target_firing_rate.py
   Target-vs-Achieved comparison is two scalar endpoints per category
   (not a distribution), so converted to a labelled dumbbell plot that
   explicitly marks each data point.

 - examples/02_finetune/02_compute_force_from_optimized_dd.py
 - examples/03_papers/watanabe/01_compute_baseline_force.py
   Per-motor-unit firing rate bar chart replaced with a violin +
   box-and-whisker shell (only when n >= 10, matching the journal's
   threshold) plus a jittered scatter of every active unit.
   Centrality (Mean) and dispersion (SD) annotated in the legend.

No functional change to the underlying data or simulation pipeline.
Codex caught a visual bug in the two reshaped figures: the existing
`for ax in axes: ax.set_xlim(0, SIMULATION_TIME_MS / 1000)` loop now
squashed the new firing-rate distribution subplot onto the left edge
of the axis and clipped the negative jitter. Restricted the loop to
the two time-series subplots (axes[:2]); the categorical distribution
subplot sets its own ticks.

Applies to:
 - examples/02_finetune/02_compute_force_from_optimized_dd.py
 - examples/03_papers/watanabe/01_compute_baseline_force.py
The main Figure 4B pools experimental ISI/CV statistics across VL, VM
and FDI into one axis. A per-muscle breakdown was requested
so the simulation overlay can be compared against each individual
muscle's envelope. The pooled main figure stays untouched; this
commit adds a companion three-panel supplementary figure.

Added to examples/02_finetune/05_plot_isi_cv_multi_muscle_comparison.py:

 - plot_cv_vs_fr_per_muscle(all_muscle_data, exp_data, muscles=...)
   builds a 1 x N grid (default VL | VM | FDI) where each panel shows
   that muscle's experimental convex hull + scatter alongside the full
   simulation overlay (same colors, same force-level markers as the
   pooled figure) so the simulation-vs-muscle match is evaluated
   panel-local.

 - Main-script block that calls the new function after the pooled plot
   and saves the companion figure to
   `isi_cv_per_muscle_supplement.{OUTPUT_FORMAT}` in the existing
   results directory.
Codex's LOW suggestion: without a legend the per-panel force markers
were only decoded implicitly via the main Fig 4B. Added a shared
figure-level legend covering muscle types and force-level markers
(same mapping used by the pooled plot), positioned below the panels.
Chunk 1 inadvertently inserted [tool.pytest.ini_options] between the
dev and docs dependency-groups, which reparented the docs entries
under pytest config and removed them from uv's dep-group resolver.
Moved the pytest section below the dep-groups.
sphinx 9.1.0 broke sphinx-hoverxref 1.4.2 (cannot unpack non-iterable
_Opt object in deprecated_configs_warning). sphinx-hoverxref has no
released fix yet; pinning sphinx<9 restores the docs build. Matches
what CI must have been resolving to before sphinx 9 released.

Verified: full `make html` build completes (with the documented
`|| make html` retry for NEURON nrn_finitialize state drift) and
renders every example in 01_basic/, 02_finetune/, and
03_papers/watanabe/, including all files touched by this branch.
@RaulSimpetru RaulSimpetru force-pushed the reviewer_suggestions branch 2 times, most recently from 3a2d087 to e6cd862 Compare April 19, 2026 04:53
@RaulSimpetru RaulSimpetru changed the title Reviewer-driven code changes for MyoGen Nature Communications rebuttal RNG accessor refactor, FR_std NaN guard, optional elephant, figure distribution plots Apr 19, 2026
… core

Codex review flagged that the previous "elephant is an optional extra"
change was only half-true: the core install still pulled viziphant,
which in turn pulled elephant transitively. Users running
`pip install myogen` therefore still received elephant, contradicting
the documented intent.

Fix:
 - Moved `viziphant>=0.4.0` out of core `[project.dependencies]` into
   the `elephant` optional-dependency bundle, where it lives alongside
   `elephant>=1.1.1`. The single `from viziphant.rasterplot import ...`
   use site is `examples/01_basic/02_simulate_spike_trains_current_injection.py`,
   which also imports elephant directly, so it already requires the
   elephant extra.
 - Added `neo>=0.14.0` as an explicit core dependency. It was reached
   transitively via `viziphant -> neo`, so dropping viziphant from core
   would otherwise break `import myogen` (12 direct import sites in
   core, including `myogen/utils/continuous_saver.py` which is on the
   top-level import path).
 - Added `elephant` and `viziphant` to the docs dep-group so
   sphinx-gallery can execute the examples that depend on them.

Verification: `uv sync` followed by
`python -c "import importlib.util; print(importlib.util.find_spec('elephant'))"`
now prints `None`; `import myogen` succeeds; 14/14 tests pass.

Also:
 - Documented `derive_subseed`'s non-negative-label precondition in
   the docstring (SeedSequence rejects negatives).
 - Added a `tests` GitHub Actions workflow that runs pytest on every
   PR to main (the existing workflows only build wheels, smoke-test
   imports, or publish docs on release).

CHANGELOG updated to reflect the true scope.
@RaulSimpetru RaulSimpetru force-pushed the reviewer_suggestions branch from 4db734c to 4c9ff67 Compare April 19, 2026 05:19
Codex noted the tests workflow was scoped only to main, which means
merging reviewer_suggestions into improved-semg wouldn't trigger CI
on that integration branch. Added improved-semg to both the pull_request
and push triggers.
@RaulSimpetru
Copy link
Copy Markdown
Member Author

Superseded — the reviewer_suggestions branch has been merged into improved-semg (see 3d05085). The combined work will land on main via the improved-semg PR.

@RaulSimpetru RaulSimpetru deleted the reviewer_suggestions branch April 19, 2026 06:26
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant