gamma_sup_generator

gamma_sup_generator#

class brainpy.state.gamma_sup_generator(in_size=1, rate=Quantity(0., 'Hz'), gamma_shape=1, n_proc=1, start=Quantity(0., 'ms'), stop=None, origin=Quantity(0., 'ms'), rng_seed=0, name=None)#

Superposition of independent gamma processes (NEST-compatible).

Description

gamma_sup_generator re-implements NEST’s stimulation device of the same name. It emits, per output train and simulation step, the multiplicity of spikes produced by superimposing n_proc independent component renewal processes with gamma-distributed inter-spike intervals.

1. State-space model, derivation, and update equations

Let \(k = \mathrm{gamma\_shape}\) and \(r = \mathrm{rate}\) in Hz. Each component gamma process is represented as a cyclic chain of k exponential phases. Using an occupation vector per output train,

\[\mathbf{occ} = (occ_0, \dots, occ_{k-1}), \qquad \sum_{i=0}^{k-1} occ_i = n_{\mathrm{proc}},\]

a process in phase i transitions to phase i+1 (mod k) with per-step probability

\[p = r \cdot k \cdot \Delta t / 1000.\]

This is the discrete-time hazard form used by NEST. For each bin i,

\[n_i \sim \mathrm{Binomial}(occ_i, p),\]

except for NEST’s sparse/high-count approximation branch:

if occ_i >= 100 and p <= 0.01, or
if occ_i >= 500 and p * occ_i <= 0.1,

use Poisson(p * occ_i) and clip to occ_i.

After sampling, all n_i are moved simultaneously to preserve integer mass and avoid order-dependent updates. The emitted spike multiplicity for one train is

\[K_n = n_{k-1},\]

i.e., transitions leaving the last phase and re-entering phase 0. This allows per-step counts larger than 1, matching NEST SpikeEvent multiplicity semantics.

2. Timing semantics and activity window

Activity follows NEST StimulationDevice::is_active for spike generators:

\[t_{\min} < t \le t_{\max}, \qquad t_{\min} = origin + start,\quad t_{\max} = origin + stop.\]

Therefore start is exclusive and stop is inclusive. Internally, finite times are projected to integer steps with round(time_ms / dt_ms) and checked as t_min_step < curr_step <= t_max_step.

3. Assumptions, constraints, and computational implications

Parameters are scalarized to float64/int before simulation. Enforced constraints are rate >= 0, gamma_shape >= 1, n_proc >= 1, and stop >= start. If dt is available, finite origin, start, and stop must lie on the simulation grid (absolute tolerance 1e-12 in time/dt ratio).

Runtime complexity of update() is \(O(\prod \mathrm{varshape} \cdot \mathrm{gamma\_shape})\), with state memory \(O(\prod \mathrm{varshape} \cdot \mathrm{gamma\_shape})\) from the occupation matrix. RNG sampling uses numpy.random.Generator (seeded by rng_seed), so stochastic draws are CPU NumPy based rather than JAX-key based.

Parameters:

in_size (Size, optional) – Output size specification consumed by brainstate.nn.Dynamics. The exact shape returned by update() is self.varshape derived from in_size; each element corresponds to one independent output train. Default is 1.
rate (ArrayLike, optional) – Scalar component-process rate in spikes/s (Hz), shape () after conversion. Accepts a single-element numeric ArrayLike or a saiunit.Quantity convertible to u.Hz. Must satisfy rate >= 0. Default is 0.0 * u.Hz.
gamma_shape (ArrayLike, optional) – Scalar integer gamma shape \(k\), shape () after conversion. Parsed via nearest-integer check with absolute tolerance 1e-12. Must satisfy gamma_shape >= 1. Default is 1.
n_proc (ArrayLike, optional) – Scalar integer number of independent component processes per output train, shape () after conversion. Parsed by nearest-integer check with absolute tolerance 1e-12. Must satisfy n_proc >= 1. Default is 1.
start (ArrayLike, optional) – Scalar relative activation time in ms, shape () after conversion. Effective lower activity bound is origin + start and is exclusive. Must be grid-representable when dt is available. Default is 0.0 * u.ms.
stop (ArrayLike or None, optional) – Scalar relative deactivation time in ms, shape () after conversion. Effective upper activity bound is origin + stop and is inclusive. None maps to +inf. Must satisfy stop >= start and be grid-representable when finite and dt is available. Default is None.
origin (ArrayLike, optional) – Scalar time-origin offset in ms, shape () after conversion, added to start and stop to compute absolute active bounds. Must be grid-representable when finite and dt is available. Default is 0.0 * u.ms.
rng_seed (int, optional) – Seed used to initialize numpy.random.default_rng in init_state(). Default is 0.
name (str or None, optional) – Optional node name passed to brainstate.nn.Dynamics.

Parameter Mapping

Table 35 Parameter mapping to model symbols#
Parameter	Default	Math symbol	Semantics
`rate`	`0.0 * u.Hz`	\(r\)	Component-process rate in spikes/s.
`gamma_shape`	`1`	\(k\)	Number of cyclic exponential phases per component process.
`n_proc`	`1`	\(n_{\mathrm{proc}}\)	Number of superimposed component processes per output train.
`start`	`0.0 * u.ms`	\(t_{\mathrm{start,rel}}\)	Relative exclusive lower bound of activity.
`stop`	`None`	\(t_{\mathrm{stop,rel}}\)	Relative inclusive upper bound; `None` maps to `+\infty`.
`origin`	`0.0 * u.ms`	\(t_0\)	Global offset added to `start` and `stop`.
`in_size`	`1`		Defines `self.varshape` for independent output trains.
`rng_seed`	`0`		Seed of the NumPy generator used for transition draws.

Raises:

ValueError – If scalar conversion fails due to non-scalar inputs; if rate < 0; if gamma_shape < 1; if n_proc < 1; if stop < start; if integer-valued inputs are non-integral beyond tolerance; or if finite origin/start/stop are not multiples of simulation resolution when dt is available.
TypeError – If unit conversion to u.Hz or u.ms fails for supplied inputs.
KeyError – At runtime, if required simulation-context fields (for example dt used by brainstate.environ.get_dt()) are unavailable.

Notes

Initial occupation is the NEST equilibrium approximation used in pre_run_hook(): floor(n_proc / gamma_shape) in all bins, with the remainder added to the last bin.
As in NEST, each output train maintains independent internal occupation states.
In the direct binomial branch, transition_prob is numerically clamped to [0, 1] before calling NumPy RNG to avoid invalid probability inputs in edge cases.

Examples

>>> import brainpy
>>> import brainstate
>>> import saiunit as u
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     gen = brainpy.state.gamma_sup_generator(
...         in_size=(2, 3),
...         rate=20.0 * u.Hz,
...         gamma_shape=3,
...         n_proc=50,
...         start=5.0 * u.ms,
...         stop=40.0 * u.ms,
...         rng_seed=7,
...     )
...     with brainstate.environ.context(t=12.0 * u.ms):
...         counts = gen.update()
...     _ = counts.shape

>>> import brainpy
>>> import saiunit as u
>>> gen = brainpy.state.gamma_sup_generator(rate=15.0 * u.Hz, gamma_shape=2)
>>> gen.set(n_proc=20, stop=None, origin=1.0 * u.ms)
>>> params = gen.get()
>>> _ = params['gamma_shape'], params['n_proc']

gamma_sup_generator

Contents

gamma_sup_generator#