ppd_sup_generator

ppd_sup_generator#

class brainpy.state.ppd_sup_generator(in_size=1, rate=Quantity(0., 'Hz'), dead_time=Quantity(0., 'ms'), n_proc=1, frequency=Quantity(0., 'Hz'), relative_amplitude=0.0, start=Quantity(0., 'ms'), stop=None, origin=Quantity(0., 'ms'), rng_seed=0, name=None)#

Superposition of Poisson processes with dead time (NEST-compatible).

Description

ppd_sup_generator re-implements NEST’s stimulation device with the same name. For each output train, it emits the per-step multiplicity of a superposition of n_proc independent Poisson-like component processes with absolute dead time.

1. State model, derivation, and update equations

Let \(r=\mathrm{rate}\) (Hz), \(\tau_d=\mathrm{dead\_time}\) (ms), and \(\Delta t\) be the simulation resolution in ms. For each output train, the internal state is an age-discretized occupancy model:

occ_active: number of currently active component processes.
occ_refractory[a] for \(a=0,\dots,\lfloor\tau_d/\Delta t\rfloor-1\): number of processes in refractory age bin a.
activate: rotating pointer indicating the bin whose occupants become active at the current step.

NEST’s per-step hazard for one active process is

\[h_{\mathrm{step}} = \frac{\Delta t}{1000/r-\tau_d}.\]

This discretization is valid under NEST’s model constraint \(1000/r>\tau_d\) (or rate == 0). If sinusoidal modulation is enabled, the instantaneous hazard becomes

\[h_t = h_{\mathrm{step}} \left(1 + A \sin\left(2\pi f t / 1000\right)\right),\]

where \(A=\mathrm{relative\_amplitude}\in[0,1]\) and \(f=\mathrm{frequency}\) in Hz.

For each train and step, emitted multiplicity n_spikes is sampled from the active pool using NEST’s branch logic:

Binomial branch: Binomial(occ_active, h_t).
Poisson approximation branch when (occ_active >= 100 and h_t <= 0.01) or (occ_active >= 500 and h_t * occ_active <= 0.1): sample Poisson(h_t * occ_active) and clip to occ_active.

State transition for nonzero refractory bins is

\[occ\_active' = occ\_active + occ\_refractory[p] - n\_spikes, \quad occ\_refractory[p]' = n\_spikes,\]

with pointer update \(p'=(p+1)\bmod B\), \(B=\lfloor\tau_d/\Delta t\rfloor\). If B == 0 (zero dead time), the active pool is not decremented by refractory bookkeeping and each component can contribute at most one spike per step through the binomial/Poisson draw.

2. Timing semantics and activity window

Activity follows NEST StimulationDevice semantics for 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 steps with round(time_ms / dt_ms) and checked as t_min_step < curr_step <= t_max_step.

3. Assumptions, constraints, and computational implications

All physical parameters are scalarized to host-side float64 or int before simulation. Enforced constraints are:

dead_time >= 0.
n_proc >= 1.
relative_amplitude in [0, 1].
stop >= start.
1000 / rate > dead_time (or rate == 0).

If dt is available, finite origin, start, and stop must be exact grid multiples (absolute tolerance 1e-12 in time/dt ratio). Runtime of update() is \(O(\prod \mathrm{varshape})\) per step; memory is \(O(\prod \mathrm{varshape} \cdot \lfloor\tau_d/\Delta t\rfloor)\). Random draws are produced by numpy.random.Generator (seeded by rng_seed), so stochasticity is NumPy-host based rather than JAX-key based.

Parameters:

in_size (Size, optional) – Output size specification consumed by brainstate.nn.Dynamics. self.varshape derived from this value is the exact shape returned by update(), and 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 1000 / rate > dead_time when rate > 0. Default is 0.0 * u.Hz.
dead_time (ArrayLike, optional) – Scalar absolute refractory time in ms, shape () after conversion. Accepts a single-element numeric ArrayLike or a saiunit.Quantity convertible to u.ms. Must satisfy dead_time >= 0. Default is 0.0 * u.ms.
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.
frequency (ArrayLike, optional) – Scalar sinusoidal modulation frequency in Hz, shape () after conversion. frequency == 0 disables sinusoidal variation even when relative_amplitude > 0. Default is 0.0 * u.Hz.
relative_amplitude (ArrayLike, optional) – Scalar dimensionless modulation amplitude \(A\), shape () after conversion. Must satisfy 0 <= relative_amplitude <= 1. Default is 0.0.
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 38 Parameter mapping to model symbols#
Parameter	Default	Math symbol	Semantics
`rate`	`0.0 * u.Hz`	\(r\)	Component-process rate in spikes/s.
`dead_time`	`0.0 * u.ms`	\(\tau_d\)	Absolute refractory duration in ms.
`n_proc`	`1`	\(n_{\mathrm{proc}}\)	Number of component processes per output train.
`frequency`	`0.0 * u.Hz`	\(f\)	Modulation frequency in Hz.
`relative_amplitude`	`0.0`	\(A\)	Relative sinusoidal modulation amplitude.
`start`	`0.0 * u.ms`	\(t_{\mathrm{start,rel}}\)	Relative exclusive lower activity bound.
`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 for NumPy RNG used by stochastic transition draws.

Raises:

ValueError – If scalar conversion fails due to non-scalar inputs; if dead_time is negative; if n_proc < 1; if relative_amplitude is outside [0, 1]; if stop < start; if 1000 / rate <= dead_time for nonzero rate; 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 conversion to u.Hz/u.ms or numeric casting fails for provided parameter types.
KeyError – At runtime, if required simulation-context fields (for example dt used by brainstate.environ.get_dt()) are unavailable.

Notes

Initial occupancy matches NEST pre_run_hook(): floor(rate / 1000 * n_proc * dt) in each refractory bin and the remainder in occ_active.
NEST does not initialize to sinusoidal equilibrium, so modulation can show startup transients.
Stimulation-backend parameter order in NEST is [dead_time, rate, n_proc, frequency, relative_amplitude].

Examples

>>> import brainpy
>>> import brainstate
>>> import saiunit as u
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     gen = brainpy.state.ppd_sup_generator(
...         in_size=(2, 2),
...         rate=20.0 * u.Hz,
...         dead_time=2.0 * u.ms,
...         n_proc=80,
...         frequency=8.0 * u.Hz,
...         relative_amplitude=0.25,
...         start=5.0 * u.ms,
...         stop=50.0 * u.ms,
...         rng_seed=3,
...     )
...     with brainstate.environ.context(t=12.0 * u.ms):
...         counts = gen.update()
...     _ = counts.shape

>>> import brainpy
>>> import saiunit as u
>>> gen = brainpy.state.ppd_sup_generator(rate=15.0 * u.Hz, n_proc=30)
>>> gen.set(dead_time=1.5 * u.ms, stop=None, origin=2.0 * u.ms)
>>> params = gen.get()
>>> _ = params['dead_time'], params['stop']

ppd_sup_generator

Contents

ppd_sup_generator#