sinusoidal_gamma_generator

sinusoidal_gamma_generator#

class brainpy.state.sinusoidal_gamma_generator(in_size=1, rate=Quantity(0., 'Hz'), amplitude=Quantity(0., 'Hz'), frequency=Quantity(0., 'Hz'), phase=0.0, order=1.0, individual_spike_trains=True, start=Quantity(0., 'ms'), stop=None, origin=Quantity(0., 'ms'), rng_seed=0, name=None)#

Sinusoidally modulated gamma spike generator compatible with NEST.

Description

sinusoidal_gamma_generator re-implements NEST’s stimulation device of the same name. It emits binary spikes from an inhomogeneous gamma renewal process whose instantaneous rate is sinusoidally modulated.

1. Instantaneous-rate model

The internal rate in spikes/ms is

\[\lambda(t) = r + a \sin(\omega t + \phi),\]

with parameter-to-symbol conversion:

  • \(r = \mathrm{rate}/1000\),

  • \(a = \mathrm{amplitude}/1000\),

  • \(\omega = 2\pi \cdot \mathrm{frequency}/1000\) (rad/ms),

  • \(\phi = \mathrm{phase}\cdot\pi/180\) (rad).

The validated constraint 0 <= amplitude <= rate guarantees \(\lambda(t) \ge 0\) for all \(t\).

2. Renewal integral, closed-form increment, and hazard

For gamma order \(k = \mathrm{order}\) and train-specific renewal origin \(t_0\), define the scaled integrated hazard as

\[\Lambda(t) = k \int_{t_0}^{t} \lambda(s)\,ds.\]

The implementation keeps t0_ms and Lambda_t0 as per-train state variables and advances \(\Lambda\) each step via the closed-form increment computed in _delta_lambda():

\[\Delta\Lambda = k r (t_b - t_a) - \frac{k a}{\omega}\Bigl[ \cos(\omega t_b + \phi) - \cos(\omega t_a + \phi) \Bigr].\]

When amplitude == 0 or frequency == 0 (i.e. \(\omega = 0\)), the cosine term is omitted and \(\Delta\Lambda = k r (t_b - t_a)\), which avoids division by zero and recovers the homogeneous Poisson limit (\(k = 1\)) or homogeneous gamma limit (\(k > 1\)).

The per-step hazard (already multiplied by dt) evaluated at time \(t\) is

\[h(t) = \Delta t \cdot \frac{k\,\lambda(t)\,\Lambda(t)^{k-1}\,e^{-\Lambda(t)}} {\Gamma(k,\,\Lambda(t))},\]

where \(\Gamma(k, \Lambda)\) is the upper incomplete gamma function evaluated via jax.lax.igammac and math.gamma. The ratio \(h(t)\) approximates \(\Pr(\text{spike in } [t, t+\Delta t))\) under the gamma renewal model.

3. Update ordering and activity-window semantics

Each call to update() mirrors the ordering in NEST models/sinusoidal_gamma_generator.cpp:

  1. Evaluate time at the right edge of the current step: t_eval = (step + 1) * dt.

  2. Compute \(\lambda(t_{\mathrm{eval}})\) and cache the value as the step-end instantaneous rate in spikes/s (accessible via get_recorded_rate()).

  3. If the generator is active and \(\lambda(t_{\mathrm{eval}}) > 0\), compute the per-train hazard and sample Bernoulli decisions.

  4. Reset t0_ms and Lambda_t0 to t_eval / 0 for any train that fired.

  5. Return binary spike outputs as int64 with shape self.varshape.

NEST spike-generator activity semantics use the half-open-right window

\[t_{\min} < n \le t_{\max},\]

where \(n\) is the current integer step index, t_min = round((origin + start) / dt), and t_max = round((origin + stop) / dt) after projection to grid steps.

4. Piecewise-integral semantics on parameter changes

When set() is called after initialization, the existing per-train renewal state is first advanced to the change time \(t_c\) using the previous process parameters, then future increments use the new parameters:

\[\Lambda(t) = \Lambda_{\mathrm{old}}(t_c) + k_{\mathrm{new}} \int_{t_c}^{t} \lambda_{\mathrm{new}}(s)\,ds.\]

This preserves renewal history across parameter updates, matching NEST SetStatus behavior.

5. Assumptions, constraints, and computational implications

  • Public parameters are scalarized to float64 (or int for rng_seed); non-scalar inputs raise ValueError.

  • Enforced constraints: order >= 1, 0 <= amplitude <= rate, and stop >= start.

  • When dt is available at construction time, finite origin / start / stop must be representable on the simulation grid (absolute tolerance 1e-12 in time / dt ratio).

  • individual_spike_trains=True allocates one independent renewal state per element of self.varshape and draws independent Bernoulli samples; individual_spike_trains=False maintains one shared renewal state and broadcasts the single Bernoulli draw to all outputs.

  • Per-step runtime is \(O(n_{\mathrm{trains}})\) for hazard evaluation and sampling, with memory \(O(n_{\mathrm{trains}})\) for t0_ms and Lambda_t0.

  • At most one spike per train can be emitted per step because spike decisions are independent Bernoulli trials against the per-step hazard.

Parameters:

  • in_size (Size, optional) – Output size specification for brainstate.nn.Dynamics. self.varshape derived from in_size is the exact output shape of update(); each element corresponds to one emitted train. Default is 1.

  • rate (ArrayLike, optional) – Scalar mean firing rate in spikes/s (Hz), shape () after conversion. Accepted as a scalar ArrayLike or a saiunit.Quantity convertible to u.Hz. Must satisfy 0 <= amplitude <= rate. Default is 0.0 * u.Hz.

  • amplitude (ArrayLike, optional) – Scalar sinusoidal modulation amplitude in spikes/s (Hz), shape () after conversion. Must satisfy 0 <= amplitude <= rate after conversion; the constraint ensures \(\lambda(t) \ge 0\). Default is 0.0 * u.Hz.

  • frequency (ArrayLike, optional) – Scalar modulation frequency in Hz, shape () after conversion. Internally mapped to angular frequency \(\omega = 2\pi f / 1000\) (rad/ms). Default is 0.0 * u.Hz.

  • phase (ArrayLike, optional) – Scalar initial phase in degrees, shape () after conversion; internally converted to radians as \(\phi = \phi_{\deg} \pi / 180\). Default is 0.0.

  • order (ArrayLike, optional) – Scalar gamma renewal order \(k\), shape () after conversion. Must satisfy order >= 1; order 1 recovers an inhomogeneous Poisson process. Default is 1.0.

  • individual_spike_trains (bool, optional) – Spike-generation mode selector. If True, each output index in self.varshape maintains independent renewal state (t0_ms, Lambda_t0) and receives independent Bernoulli draws. If False, one shared renewal process is maintained and the same binary spike decision is broadcast to all outputs. Default is True.

  • start (ArrayLike, optional) – Scalar relative activation start time in ms, shape () after conversion. Effective lower activity bound is origin + start; the bound is exclusive in step space (t_min_step < curr_step). Default is 0.0 * u.ms.

  • stop (ArrayLike or None, optional) – Scalar relative deactivation stop time in ms, shape () after conversion. None maps to +inf. Effective upper activity bound is origin + stop; the bound is inclusive in step space (curr_step <= t_max_step). Must satisfy stop >= start after conversion. Default is None.

  • origin (ArrayLike, optional) – Scalar origin offset in ms, shape () after conversion. Added to start and stop to obtain the absolute activity bounds t_min and t_max. Default is 0.0 * u.ms.

  • rng_seed (int, optional) – Seed used to initialize jax.random.PRNGKey during init_state() and lazy initialization in update(). Default is 0.

  • name (str or None, optional) – Optional node name passed to brainstate.nn.Dynamics. Default is None.

Parameter Mapping

Table 36 Parameter mapping to model symbols#

Parameter

Default

Math symbol

Semantics

rate

0.0 * u.Hz

\(r\)

Baseline firing-rate term in spikes/ms after division by 1000.

amplitude

0.0 * u.Hz

\(a\)

Sinusoidal modulation amplitude in spikes/ms after division by 1000.

frequency

0.0 * u.Hz

\(f\)

Frequency in Hz mapped to \(\omega = 2\pi f/1000\) (rad/ms).

phase

0.0

\(\phi\)

Phase in degrees mapped to radians as \(\phi_{\deg}\pi/180\).

order

1.0

\(k\)

Gamma renewal order; 1 = inhomogeneous Poisson.

start

0.0 * u.ms

\(t_{\mathrm{start,rel}}\)

Relative exclusive lower bound of activity window.

stop

None

\(t_{\mathrm{stop,rel}}\)

Relative inclusive upper bound; None maps to \(+\infty\).

origin

0.0 * u.ms

\(t_{\mathrm{origin}}\)

Global time offset added to start and stop.

in_size

1

—

Defines self.varshape and the total output train count.

individual_spike_trains

True

—

Independent per-output renewal states vs. shared broadcast process.

rng_seed

0

—

Seed for JAX random key initialization and splitting.
Raises:

  • ValueError – If scalar-conversion fails due to non-scalar inputs; if 0 <= amplitude <= rate is violated; if order < 1; if stop < start; or if finite origin / start / stop are not multiples of the simulation resolution when dt is available.

  • TypeError – If provided values cannot be converted to numeric values or to the required units (e.g. a non-convertible u.Hz or u.ms quantity).

  • KeyError – At runtime in update(), if required simulation-context entries (notably dt) are unavailable from brainstate.environ.

Notes

  • Hazard values are computed in float64; tiny negative \(\Lambda\) values arising from floating-point roundoff are clamped to zero before hazard evaluation.

  • The value returned by get_recorded_rate() is the step-end instantaneous rate in spikes/s, matching NEST’s rate recordable.

  • Renewal state is revalidated against the timing grid whenever dt changes between update() calls.

Examples

Simulate a 2×3 array of independent sinusoidally modulated gamma trains:

>>> import brainpy
>>> import brainstate
>>> import saiunit as u
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     gen = brainpy.state.sinusoidal_gamma_generator(
...         in_size=(2, 3),
...         rate=50.0 * u.Hz,
...         amplitude=20.0 * u.Hz,
...         frequency=8.0 * u.Hz,
...         phase=30.0,
...         order=3.0,
...         start=5.0 * u.ms,
...         stop=80.0 * u.ms,
...         rng_seed=9,
...     )
...     with brainstate.environ.context(t=12.0 * u.ms):
...         spikes = gen.update()
...     _ = spikes.shape

Use individual_spike_trains=False and update parameters at runtime:

>>> import brainpy
>>> import brainstate
>>> import saiunit as u
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     gen = brainpy.state.sinusoidal_gamma_generator(
...         individual_spike_trains=False
...     )
...     gen.set(rate=40.0 * u.Hz, amplitude=10.0 * u.Hz, order=2.0)
...     params = gen.get()
...     _ = params['rate'], params['order']

See also

sinusoidal_poisson_generator

Sinusoidally modulated Poisson generator.

gamma_sup_generator

Superposition of independent gamma-renewal processes.

References

get()[source]#

Return current public parameters as plain Python scalars.

Returns:

out – Mapping of parameter names to their current values. Keys and value types are:

  • 'rate' : float — mean firing rate in spikes/s.

  • 'amplitude' : float — sinusoidal modulation amplitude in spikes/s.

  • 'frequency' : float — modulation frequency in Hz.

  • 'phase' : float — modulation phase in degrees.

  • 'order' : float — gamma renewal order \(k\).

  • 'individual_spike_trains' : bool — spike-generation mode flag.

  • 'start' : float — relative activation start in ms.

  • 'stop' : float — relative deactivation stop in ms (float('inf') when no stop was set).

  • 'origin' : float — time-origin offset in ms.

Return type:

dict

get_recorded_rate()[source]#

Return the latest step-end instantaneous rate in spikes/s.

The value is updated by update() at each simulation step to \(\lambda(t_{\mathrm{eval}}) \times 1000\) spikes/s, where \(t_{\mathrm{eval}} = (\mathrm{step} + 1) \times dt\) is the right edge of the current step. This matches NEST’s rate recordable quantity.

Returns:

out – Most recently cached instantaneous rate in spikes/s (float64 scalar). Returns 0.0 if init_state() has not been called yet.

Return type:

float

init_state(batch_size=None, **kwargs)[source]#

Initialize random key and per-train renewal state.

Allocates three brainstate.ShortTermState variables:

  • rng_key: a JAX PRNGKey seeded from self.rng_seed.

  • t0_ms: per-train renewal origin, initialized to the current simulation time (float64 array of length self._num_trains).

  • Lambda_t0: per-train accumulated scaled hazard at t0_ms, initialized to zero (float64 array of length self._num_trains).

  • _recorded_rate_hz: cached step-end instantaneous rate in spikes/s, initialized to 0.0.

The timing cache (_t_min_step, _t_max_step) is also refreshed if dt is available in the current brainstate.environ context.

Parameters:

  • batch_size (int or None, optional) – Unused placeholder kept for brainstate.nn.Dynamics signature compatibility. Default is None.

  • **kwargs – Unused extra keyword arguments; silently ignored.

Raises:

ValueError – If finite origin / start / stop do not lie on the simulation grid when dt is available.

set(*, rate=<object object>, amplitude=<object object>, frequency=<object object>, phase=<object object>, order=<object object>, individual_spike_trains=<object object>, start=<object object>, stop=<object object>, origin=<object object>)[source]#

Update public parameters and refresh the process and timing caches.

Any parameter left at its sentinel value _UNSET retains its current value. When called after init_state(), the internal renewal state is first advanced to the current simulation time using the previous process parameters before switching to the new ones, preserving the piecewise-integral semantics described in the class docstring. If individual_spike_trains changes in a way that alters the required number of trains, t0_ms and Lambda_t0 are grown (new trains start fresh) or truncated accordingly.

Parameters:

  • rate (ArrayLike or None, optional) – New scalar mean rate in spikes/s (Hz), shape () after conversion. _UNSET retains the current value.

  • amplitude (ArrayLike or None, optional) – New scalar modulation amplitude in spikes/s (Hz), shape () after conversion. _UNSET retains the current value.

  • frequency (ArrayLike or None, optional) – New scalar modulation frequency in Hz, shape () after conversion. _UNSET retains the current value.

  • phase (ArrayLike or None, optional) – New scalar modulation phase in degrees, shape () after conversion. _UNSET retains the current value.

  • order (ArrayLike or None, optional) – New scalar gamma order \(k\), shape () after conversion. _UNSET retains the current value.

  • individual_spike_trains (bool or None, optional) – New spike-generation mode. _UNSET retains the current value.

  • start (ArrayLike or None, optional) – New scalar relative activation start time in ms, shape () after conversion. _UNSET retains the current value.

  • stop (ArrayLike, None, or sentinel, optional) – New scalar relative stop time in ms, shape () after conversion; explicit None maps to +inf. _UNSET retains the current value.

  • origin (ArrayLike or None, optional) – New scalar origin offset in ms, shape () after conversion. _UNSET retains the current value.

Raises:

  • ValueError – If scalar conversion fails due to non-scalar inputs; if the constraints order >= 1, 0 <= amplitude <= rate, or stop >= start are violated after resolving new values; or if finite timing parameters do not lie on the simulation grid when dt is available.

  • TypeError – If numeric or unit conversion fails for any supplied input.

update()[source]#

Advance one simulation step and emit binary spike events.

Reads the current time t and resolution dt from brainstate.environ, lazily calls init_state() if state has not been allocated, and refreshes the timing cache if dt has changed since the last call. The step-end time t_eval = (step + 1) * dt is used for rate evaluation and \(\Lambda\) accumulation. Trains outside the active window or with \(\lambda(t_{\mathrm{eval}}) \le 0\) receive zero spikes without consuming random draws.

Returns:

out – int64 JAX array with shape self.varshape. Each element is 0 or 1, giving the binary spike decision for the corresponding output train in the current step. When individual_spike_trains=False, all elements share the same value.

Return type:

jax.Array

Raises:

  • KeyError – If dt is not available in the current brainstate.environ context.

  • ValueError – If timing-parameter grid validation fails after the simulation resolution changes between calls.