From d3f7e76e91d2d863ab64ea6420ec4a4bb87b260a Mon Sep 17 00:00:00 2001 From: chaoming Date: Fri, 19 Jun 2026 02:35:31 +0800 Subject: [PATCH] fix(dyn/synapses): DualExp equal-tau handling and STP per-neuron reset - DualExpon raised ZeroDivisionError / produced NaN when tau_rise == tau_decay; compute the normalization with the e/tau limit element-wise (High) - DualExponV2 silently output zeros/NaN for equal taus; raise a clear ValueError in the auto-normalizer pointing to bp.dyn.Alpha (High) - STP.reset_state crashed for a per-neuron array U (Variable.fill_ needs a scalar); broadcast U instead (High) Findings recorded in docs/issues-found-20260619-dyn-synapses.md --- brainpy/dyn/synapses/abstract_models.py | 58 ++++++- brainpy/dyn/synapses/abstract_models_test.py | 73 +++++++++ docs/issues-found-20260619-dyn-synapses.md | 152 +++++++++++++++++++ 3 files changed, 280 insertions(+), 3 deletions(-) create mode 100644 docs/issues-found-20260619-dyn-synapses.md diff --git a/brainpy/dyn/synapses/abstract_models.py b/brainpy/dyn/synapses/abstract_models.py index 96d666d18..8f17633f6 100644 --- a/brainpy/dyn/synapses/abstract_models.py +++ b/brainpy/dyn/synapses/abstract_models.py @@ -159,11 +159,57 @@ def return_info(self): def _format_dual_exp_A(self, A): A = parameter(A, sizes=self.varshape, allow_none=True, sharding=self.sharding) if A is None: + # The peak normalizer ``A = tau_decay / (tau_decay - tau_rise) * ...`` divides + # by zero when ``tau_rise == tau_decay``. For ``DualExponV2`` (the only consumer + # of this auto-normalizer that stores ``A`` directly and returns + # ``A * (g_decay - g_rise)``) the two gates collapse to identical trajectories, + # so ``g_decay - g_rise`` is identically zero and no finite ``A`` can recover a + # non-zero waveform. Fail loudly with an actionable message rather than emitting + # a ZeroDivisionError / silent NaN. ``DualExpon`` does not reach this branch for + # the equal-tau case: it computes its ``a`` coefficient directly (see below). + if bm.any(self.tau_rise == self.tau_decay): + raise ValueError( + 'The dual-exponential peak normalizer "A" is undefined when ' + '"tau_rise == tau_decay" (the rise and decay gates collapse to the ' + 'same trajectory). Use brainpy.dyn.Alpha for a single-time-constant ' + 'alpha synapse, or pass an explicit "A".' + ) A = (self.tau_decay / (self.tau_decay - self.tau_rise) * bm.float_power(self.tau_rise / self.tau_decay, self.tau_rise / (self.tau_rise - self.tau_decay))) return A +def _dual_exp_a(self, A): + r"""Compute the input-scaling coefficient ``a`` for :class:`DualExpon`. + + The conductance jump per spike is ``a``. With the auto peak normalizer + (``A is None``) the closed form simplifies to + + .. math:: + + a = \frac{1}{\tau_{rise}} + \left(\frac{\tau_{rise}}{\tau_{decay}}\right)^{\tau_{rise}/(\tau_{rise}-\tau_{decay})} + + which is finite even when :math:`\tau_{rise} = \tau_{decay}` (the + dual-exponential degenerates to the normalized alpha function). In that limit + L'Hôpital gives :math:`a = e/\tau`, evaluated element-wise so heterogeneous + time constants are supported. An explicitly supplied ``A`` is honoured as + ``a = (tau_decay - tau_rise) / (tau_rise * tau_decay) * A``. + """ + A = parameter(A, sizes=self.varshape, allow_none=True, sharding=self.sharding) + if A is not None: + return (self.tau_decay - self.tau_rise) / self.tau_rise / self.tau_decay * A + equal = self.tau_rise == self.tau_decay + ratio = self.tau_rise / self.tau_decay + # ``where(equal, ...)`` on the exponent avoids a 0/0 in the unused branch so the + # gradient/value stays finite; the equal-tau entries take the L'Hôpital value e/tau. + exponent = self.tau_rise / bm.where(equal, bm.ones_like(self.tau_decay), + self.tau_rise - self.tau_decay) + a_general = bm.float_power(ratio, exponent) / self.tau_rise + a_limit = bm.exp(bm.ones_like(self.tau_rise)) / self.tau_rise + return bm.where(equal, a_limit, a_general) + + class DualExpon(SynDyn): r"""Dual exponential synapse model. @@ -262,8 +308,11 @@ def __init__( # parameters self.tau_rise = self.init_param(tau_rise) self.tau_decay = self.init_param(tau_decay) - A = _format_dual_exp_A(self, A) - self.a = (self.tau_decay - self.tau_rise) / self.tau_rise / self.tau_decay * A + # Compute the conductance-jump coefficient ``a`` directly. This avoids the + # ``(tau_decay - tau_rise)`` cancellation that produced a ZeroDivisionError / + # NaN for equal time constants, and supports the alpha-function limit + # ``tau_rise == tau_decay`` (a = e/tau) element-wise. + self.a = _dual_exp_a(self, A) # integrator self.integral = odeint(JointEq(self.dg, self.dh), method=method) @@ -854,8 +903,11 @@ def __init__( def reset_state(self, batch_or_mode=None, **kwargs): self.x = self.init_variable(bm.ones, batch_or_mode) + # Initialise ``u`` to the release probability ``U`` by broadcasting rather than + # ``Variable.fill_`` (which only accepts a scalar). This supports a per-neuron + # array ``U`` and batched modes alike. self.u = self.init_variable(bm.ones, batch_or_mode) - self.u.fill_(self.U) + self.u.value = self.u.value * self.U @property def derivative(self): diff --git a/brainpy/dyn/synapses/abstract_models_test.py b/brainpy/dyn/synapses/abstract_models_test.py index b226ba538..9b336b407 100644 --- a/brainpy/dyn/synapses/abstract_models_test.py +++ b/brainpy/dyn/synapses/abstract_models_test.py @@ -15,6 +15,7 @@ import unittest import matplotlib.pyplot as plt +import numpy as np import brainpy as bp import brainpy.math as bm @@ -99,3 +100,75 @@ def update(self): bp.visualize.line_plot(indices * bm.get_dt(), gs, legend='g', show=show) plt.close('all') + + +class TestDualExpEqualTau(unittest.TestCase): + """Regression tests for P9-H1 / P9-H2: ``tau_rise == tau_decay``.""" + + def test_equal_tau_no_crash(self): + # P9-H1: scalar equal taus previously raised ZeroDivisionError at construction. + syn = bp.dyn.DualExpon(2, tau_rise=10., tau_decay=10.) + self.assertTrue(np.all(np.isfinite(np.asarray(syn.a)))) + # array taus with one equal pair previously yielded a NaN coefficient. + syn2 = bp.dyn.DualExpon(2, + tau_rise=bm.asarray([10., 5.]), + tau_decay=bm.asarray([10., 50.])) + self.assertTrue(np.all(np.isfinite(np.asarray(syn2.a)))) + + def test_equal_tau_matches_alpha(self): + # P9-H1: the equal-tau dual-exponential is the normalized alpha function; + # a single unit spike drives a response that peaks at ~1.0. + bm.set(dt=0.05) + + class Net(bp.DynSysGroup): + def __init__(self, tau): + super().__init__() + self.inp = bp.dyn.SpikeTimeGroup(1, bm.zeros(1, dtype=int), bm.asarray([1.])) + self.syn = bp.dyn.DualExpon(1, tau_rise=tau, tau_decay=tau) + + def update(self): + return self.syn(self.inp()) + + net = Net(10.) + indices = bm.as_numpy(bm.arange(4000)) + gs = bm.for_loop(net.step_run, indices) + peak = float(np.max(np.asarray(gs))) + self.assertTrue(np.isfinite(peak)) + self.assertAlmostEqual(peak, 1.0, places=2) + bm.set(dt=0.1) + + def test_v2_equal_tau_raises(self): + # P9-H2: DualExponV2 is structurally singular for equal taus; the auto + # normalizer must raise a clear error rather than crash / output zeros. + with self.assertRaises(ValueError): + bp.dyn.DualExponV2(2, tau_rise=10., tau_decay=10.) + + +class TestSTPReset(unittest.TestCase): + """Regression tests for P9-H3: per-neuron array ``U``.""" + + def test_array_U_reset(self): + # P9-H3: array U previously crashed in reset_state via Variable.fill_. + U = bm.asarray([0.1, 0.2, 0.3]) + syn = bp.dyn.STP(3, U=U) + np.testing.assert_allclose(np.asarray(syn.u.value), np.asarray(U)) + np.testing.assert_allclose(np.asarray(syn.x.value), np.ones(3)) + + def test_array_U_run(self): + # The synapse must also integrate without error for heterogeneous U. + from brainpy.context import share + bm.set(dt=0.1) + syn = bp.dyn.STP(3, U=bm.asarray([0.1, 0.2, 0.3])) + out = [] + for i in range(20): + share.save(t=i * 0.1, dt=0.1, i=i) + spk = bm.asarray([1., 0., 1.]) if i == 5 else bm.asarray([0., 0., 0.]) + out.append(np.asarray(syn.update(spk))) + out = np.asarray(out) + self.assertTrue(np.all(np.isfinite(out))) + + def test_scalar_U_batched(self): + # Scalar U with batching must keep working (no regression). + syn = bp.dyn.STP(3, U=0.15, mode=bm.BatchingMode(4)) + self.assertEqual(syn.u.shape, (4, 3)) + np.testing.assert_allclose(np.asarray(syn.u.value), np.full((4, 3), 0.15)) diff --git a/docs/issues-found-20260619-dyn-synapses.md b/docs/issues-found-20260619-dyn-synapses.md new file mode 100644 index 000000000..a29c850b7 --- /dev/null +++ b/docs/issues-found-20260619-dyn-synapses.md @@ -0,0 +1,152 @@ +# P9 — dyn/synapses + dyn/projections audit (2026-06-19) + +Branch: `fix/audit-20260619-dyn-synapses` +Scope: `brainpy/dyn/synapses/{abstract_models,bio_models,delay_couplings}.py`, +`brainpy/dyn/projections/{align_post,align_pre,base,conn,delta,inputs,plasticity,utils,vanilla}.py` +(+ co-located `*_test.py`). + +Note on prior audit (`dev/issues-found-20260618.md`): several previously-reported +synapse/projection bugs are **already fixed in this tree** and were re-verified as +not-present: C-06/H-39 (STP facilitation `u` uses decayed locals — see the +documented comment + `u = u + pre_spike*U*(1-u); x = x - pre_spike*u*x`), C-17 +(PoissonInput uses `scale=sqrt(b*p)` std, not variance), C-18 +(`HalfProjAlignPost.update` reuses `current` instead of calling `comm` twice), +C-19/H-41 (`STDP_Song2000` passes `W_min/W_max=None` through unchanged and uses +`bm.as_jax` on traces), H-40 (`projections/base.py` now re-exports the real +`Projection`/`SynConn`). These are noted for the record only. + +--- + +### P9-H1 — DualExpon raises ZeroDivisionError / yields NaN when `tau_rise == tau_decay` [High] +- File: brainpy/dyn/synapses/abstract_models.py:159-164,265-266 +- Category: numerics / edge / api-drift +- What: `_format_dual_exp_A` computes the peak normalizer + `A = tau_decay/(tau_decay - tau_rise) * (tau_rise/tau_decay)**(tau_rise/(tau_rise-tau_decay))`. + When `tau_rise == tau_decay` the leading factor divides by zero. For Python-float + taus this raises `ZeroDivisionError` at construction; for array taus it silently + produces `inf`, and `DualExpon`'s `a = (tau_decay-tau_rise)/(tau_rise*tau_decay)*A` + then becomes `0*inf = nan`. +- Why it's a bug: equal rise/decay time constants is a legitimate, common request + (the dual-exponential degenerates to the normalized alpha function). It is also the + parameterization used by the `dynold` `AlphaCUBA/AlphaCOBA` compat classes (C-20), + which route through `DualExpon`. A crash / silent NaN on a realistic default is wrong. +- Repro: `bp.dyn.DualExpon(2, tau_rise=10., tau_decay=10.)` -> `ZeroDivisionError`; + `bp.dyn.DualExpon(2, tau_rise=bm.asarray([10.,5.]), tau_decay=bm.asarray([10.,50.]))` + -> `a = [nan, 0.258...]`. +- Fix: compute `DualExpon.a` via the closed form `a = (1/tau_rise) * + (tau_rise/tau_decay)**(tau_rise/(tau_rise-tau_decay))` (algebraically equal to the + old expression but free of the `(tau_decay-tau_rise)` cancellation), and special-case + the equal-tau limit element-wise to its L'Hôpital value `a = e/tau` (verified: a + single-spike response then peaks at exactly 1.0). The `A is None` auto-normalizer path + only; an explicitly supplied `A` is honoured unchanged. +- Tests: `abstract_models_test.py::TestDualExpon::test_equal_tau_no_crash`, + `::test_equal_tau_matches_alpha` +- Status: fixed + +### P9-H2 — DualExponV2 silently outputs all-zeros (or NaN) when `tau_rise == tau_decay` [High] +- File: brainpy/dyn/synapses/abstract_models.py:159-164,413 +- Category: numerics / edge +- What: `DualExponV2` stores `self.a = A` and returns `a * (g_decay - g_rise)`. With + equal taus the two gates obey identical ODEs with identical inputs, so + `g_decay - g_rise ≡ 0` for all time; the synapse is structurally singular. With the + default `A=None` normalizer this additionally hits the same div-by-zero/`inf` as P9-H1, + giving `inf*0 = nan`. +- Why it's a bug: the model produces a silently-dead (identically zero) or NaN synapse + for an innocuous parameter choice, with no diagnostic. Unlike `DualExpon`, no finite + `A` can recover a non-zero waveform, so the only correct behaviour is a clear error. +- Repro: `bp.dyn.DualExponV2(2, tau_rise=10., tau_decay=10.)` -> `ZeroDivisionError`; + with an explicit `A=1.` the output is identically 0 over a full simulation. +- Fix: in `_format_dual_exp_A`, when `A is None` and `tau_rise == tau_decay`, raise a + clear `ValueError` telling the user the dual-exponential normalizer is undefined for + equal time constants and to use `brainpy.dyn.Alpha` (single-tau alpha synapse) instead. + (Only the V2/auto-normalizer path reaches this branch; `DualExpon` computes `a` + directly per P9-H1 and never calls the helper for the equal-tau case.) +- Tests: `abstract_models_test.py::TestDualExpon::test_v2_equal_tau_raises` +- Status: fixed + +### P9-H3 — STP.reset_state crashes for per-neuron (array) `U` [High] +- File: brainpy/dyn/synapses/abstract_models.py:855-858 +- Category: correctness / edge +- What: `reset_state` does `self.u = self.init_variable(bm.ones, ...)` then + `self.u.fill_(self.U)`. `fill_` requires a scalar fill value (`shape == ()`), so when + `U` is a per-neuron array the call raises + `MathError: The shape of the fill value must be ()`. +- Why it's a bug: heterogeneous release probability `U` across synapses is a standard, + realistic short-term-plasticity configuration; the model crashes on construction. +- Repro: `bp.dyn.STP(3, U=bm.asarray([0.1, 0.2, 0.3]))` -> `MathError`. +- Fix: initialise `u` by broadcasting `U` into the (possibly batched) state: + `self.u = self.init_variable(bm.ones, batch_or_mode); self.u.value = self.u.value * self.U`. + Works for scalar `U`, array `U`, and batched modes. +- Tests: `abstract_models_test.py::TestSTP::test_array_U_reset`, + `::test_array_U_run` +- Status: fixed + +### P9-M1 — STP/STD discrete jumps assume binary `pre_spike`; graded inputs scaled wrongly [Medium] +- File: brainpy/dyn/synapses/abstract_models.py:800-801,883-884 +- Category: numerics +- What: the "simplified" updates `x = x - pre_spike*U*self.x` (STD) and + `u = u + pre_spike*U*(1-u); x = x - pre_spike*u*x` (STP) reproduce the original + `bm.where(pre_spike, ...)` exactly only for `pre_spike in {0,1}`. For graded + (non-binary) presynaptic signals the depression/facilitation magnitude is scaled by + the graded value, which is not the documented Tsodyks–Markram jump. +- Why it's a bug: callers feeding graded "spikes" get a different (linearly scaled) + release, with no warning. (Matches prior-audit M-22.) +- Repro: static — compare `bm.where(pre>0, x-U*x, x)` vs `x-pre*U*x` for `pre=0.5`. +- Fix: recorded only. The simplified form is a defensible generalization, is faithful to + the historical `dynold` STD/STP convention for the binary case, and changing it risks + altering long-standing numeric behaviour. Documenting/asserting binary input is a + cross-cutting API decision left to maintainers. +- Tests: none +- Status: recorded-only + +### P9-M2 — STD applies depression jump to the pre-decay state, inconsistent with the STP fix [Medium] +- File: brainpy/dyn/synapses/abstract_models.py:795-801 +- Category: numerics +- What: `STD.update` integrates `x -> x_decayed` but then applies the release jump using + the **pre-decay** Variable: `self.x.value = x_decayed - pre_spike*U*self.x`. The + companion `STP` model was deliberately changed (prior-audit H-39, documented comment) + to apply jumps to the **decayed** local `x`. STD was left on the pre-decay form, so the + two short-term-plasticity models now discretize the spike jump inconsistently (off by + one decay step in the jump term). +- Why it's a bug: an asymmetry/correctness smell; the decayed-local form is the cleaner + discretization. +- Repro: static. +- Fix: recorded only. The pre-decay form matches the original commented code AND the + `dynold` reference (`short_term_plasticity.py` STD), and the per-step difference is + `O(dt/tau)` confined to the jump term. Changing it alters historical STD numerics for + every user; flagged for a maintainer decision rather than fixed unilaterally. +- Tests: none +- Status: recorded-only + +### P9-L1 — DelayCoupling docstrings reference a non-existent gain `g`; malformed `Parameters::` [Low] +- File: brainpy/dyn/synapses/delay_couplings.py:131-163,234-256 +- Category: style +- What: `DiffusiveCoupling`/`AdditiveCoupling` docstrings describe + `coupling = g * (...)` but no `g` gain parameter exists (the gain is folded into + `conn_mat`). They also use the literal-block `Parameters::` marker instead of a + NumPy-doc underlined `Parameters` section (CLAUDE.md mandates NumPy-doc). Matches + prior-audit L-09/L-14. +- Why it's a bug: misleading docs / Sphinx rendering. +- Fix: recorded only (Low). +- Tests: none +- Status: recorded-only + +### P9-L2 — GABAa docstring lists wrong default alpha/beta [Low] +- File: brainpy/dyn/synapses/bio_models.py:286-287 +- Category: style +- What: the `Args` block says `alpha: Default 0.062` and `beta: Default 3.57`, but the + constructor defaults are `alpha=0.53`, `beta=0.18`. +- Why it's a bug: documentation does not match code. +- Fix: recorded only (Low). +- Tests: none +- Status: recorded-only + +### P9-L3 — SynConn.update raises with non-f-string message; conn_mat shape check on raw input [Low] +- File: brainpy/dyn/projections/conn.py:119; delay_couplings.py:81 +- Category: style / edge +- What: minor robustness/style items (e.g. `AdditiveCoupling.update`'s final + `else: raise ValueError` has no message; `conn_mat.shape` is read before confirming + the object exposes `.shape`). +- Fix: recorded only (Low). +- Tests: none +- Status: recorded-only