Point Processes¶

The generic simulation framework: core abstractions and the Gillespie engine.

Core abstractions¶

process ¶

Generic continuous-time point process abstractions.

Defines the interface between a stochastic process (which knows its rates and transitions) and the Gillespie engine (which knows how to sample events).

Trajectory `dataclass` ¶

Trajectory(
    t: NDArray[float64],
    compartments: dict[str, NDArray[int_]],
)

Result from a single stochastic simulation run.

Stores time points and named compartment counts, allowing any process to record its own state variables without the engine needing to know what they represent.

Parameters:

Name	Type	Description	Default
`t`	`NDArray[float64]`	Monotonically increasing event times.	required
`compartments`	`dict[str, NDArray[int_]]`	Mapping from compartment name to count time series. Each array has the same length as `t`.	required

names `property` ¶

names: list[str]

Compartment names in insertion order.

getitem ¶

__getitem__(name: str) -> NDArray[int_]

Access a compartment by name, e.g. trajectory["infected"].

ConvergenceConfig `dataclass` ¶

ConvergenceConfig(
    t_end: float = 15.0,
    convergence_threshold: float = 0.02,
    min_realizations: int = 30,
    max_realizations: int = 10000,
    num_workers: int = 1,
)

Configuration for convergence-based ensemble stopping.

Controls when the Gillespie engine stops launching new realisations. Process-agnostic: knows nothing about SIR, Hawkes, etc.

ContinuousTimeProcess ¶

Bases: ABC

A continuous-time Markov chain on a network.

Implementations define the dynamics (rates and transitions) while the Gillespie engine handles the mechanics (exponential waiting times, event selection, convergence).

The process owns its mutable state internally. The engine interacts with it only through this interface.

Lifecycle (managed by the engine): 1. Engine calls rates() to get current event rates. 2. Engine samples dt ~ Exp(1 / total_rate) and selects an event index. 3. Engine calls execute(event_index, dt) to apply the transition. 4. Repeat until is_absorbing() returns True or t >= t_end. 5. Engine calls trajectory() to extract the recorded time series.

To define a new process, implement all abstract methods and provide a factory that constructs fresh instances (one per realisation).

time `abstractmethod` `property` ¶

time: float

Current simulation time.

rates `abstractmethod` ¶

rates() -> NDArray[float64]

Compute all event rates given the current state.

Returns:

Type	Description
`NDArray[float64]`	1-D array of non-negative rates. Length and interpretation
`NDArray[float64]`	are process-specific. The Gillespie engine selects event `i`
`NDArray[float64]`	with probability `rates[i] / rates.sum()`.

execute `abstractmethod` ¶

execute(event_index: int, dt: float) -> None

Apply the selected event and advance time by dt.

Parameters:

Name	Type	Description	Default
`event_index`	`int`	Index into the array returned by `rates()`.	required
`dt`	`float`	Time elapsed since the previous event.	required

is_absorbing `abstractmethod` ¶

is_absorbing() -> bool

Return True if no further events can occur.

The engine checks this before computing rates, so returning True here guarantees rates() will not be called on a dead state.

scalar_output `abstractmethod` ¶

scalar_output() -> float

Scalar summary for convergence monitoring.

The engine feeds this value into a Welford running-mean tracker after each realisation. Choose the quantity whose relative SEM you want to converge (e.g. final epidemic size, total event count).

should_accept ¶

should_accept() -> bool

Whether this realisation should be included in ensemble statistics.

Override to implement filtering (e.g. discard sub-critical epidemics). Discarded runs still count toward max_realizations but not toward min_realizations or convergence.

Default: always accept.

trajectory `abstractmethod` ¶

trajectory() -> Trajectory

Extract the trajectory recorded during simulation.

Called once after the simulation loop terminates. The returned Trajectory is frozen — safe to store, aggregate, or serialise.

ProcessFactory ¶

Bases: ABC

Creates fresh process instances for each realisation.

The Gillespie engine calls create() once per realisation rather than resetting state, ensuring clean isolation between runs.

Implementations hold the immutable configuration and network topology; create() builds a new mutable process from them.

convergence_config `abstractmethod` `property` ¶

convergence_config: ConvergenceConfig

Convergence parameters for the ensemble.

create `abstractmethod` ¶

create(rng: Generator) -> ContinuousTimeProcess

Build a fresh process instance with the given RNG.

Parameters:

Name	Type	Description	Default
`rng`	`Generator`	Random number generator for this realisation. Each realisation should receive an independent generator.	required

Gillespie engine¶

gillespie ¶

Generic Gillespie algorithm for continuous-time point processes.

Exact stochastic simulation engine that works with any process implementing the ContinuousTimeProcess interface. Provides convergence-based ensemble execution, trajectory aggregation, and parallel dispatch.

ProgressCallback `module-attribute` ¶

ProgressCallback = Callable[[int, int, float, bool], None]

EnsembleResult `dataclass` ¶

EnsembleResult(
    t: NDArray[float64],
    means: dict[str, NDArray[float64]],
    stds: dict[str, NDArray[float64]],
    scalar_outputs: NDArray[float64],
    convergence: ConvergenceInfo,
)

Aggregated statistics from a stochastic simulation ensemble.

Compartment time series are stored in dicts keyed by compartment name, making this type process-agnostic.

scalar_output_mean `property` ¶

scalar_output_mean: float

Mean of the scalar convergence quantity.

scalar_output_std `property` ¶

scalar_output_std: float

Standard deviation of the scalar convergence quantity.

ConvergenceInfo `dataclass` ¶

ConvergenceInfo(
    converged: bool,
    n_realizations: int,
    relative_sem: float,
    threshold: float,
    n_discarded: int = 0,
)

Metadata about convergence of a stochastic ensemble.

ConvergenceMonitor ¶

ConvergenceMonitor(threshold: float)

Track running statistics for convergence detection using Welford's algorithm.

Welford's algorithm computes running mean and variance in a single pass with numerical stability, avoiding catastrophic cancellation.

n `property` ¶

n: int

Number of observations.

mean `property` ¶

mean: float

Current running mean.

variance `property` ¶

variance: float

Current sample variance (n-1 denominator).

std `property` ¶

std: float

Current sample standard deviation.

sem `property` ¶

sem: float

Standard error of the mean.

relative_sem `property` ¶

relative_sem: float

Relative standard error (SEM / mean).

update ¶

update(value: float) -> None

Add a new observation using Welford's online algorithm.

is_converged ¶

is_converged(min_n: int) -> bool

Check if convergence criterion is met.

run_once ¶

run_once(
    factory: ProcessFactory, t_end: float, rng: Generator
) -> tuple[Trajectory, float, bool]

Run a single realisation of the Gillespie algorithm.

Parameters:

Name	Type	Description	Default
`factory`	`ProcessFactory`	Creates a fresh process instance for this realisation.	required
`t_end`	`float`	Maximum simulation time.	required
`rng`	`Generator`	Random number generator for event sampling.	required

Returns:

Type	Description
`tuple[Trajectory, float, bool]`	(trajectory, scalar_output, should_accept) tuple.

simulate ¶

simulate(
    factory: ProcessFactory,
    rng: Generator | None = None,
    n_points: int = 200,
    progress_callback: ProgressCallback | None = None,
) -> EnsembleResult

Run stochastic simulation ensemble with convergence-based stopping.

Uses the Gillespie algorithm for exact stochastic simulation. Runs until the relative standard error of the scalar output drops below the configured threshold.

Parameters:

Name	Type	Description	Default
`factory`	`ProcessFactory`	Creates fresh process instances for each realisation.	required
`rng`	`Generator \| None`	Random number generator (used for seed generation in parallel mode).	`None`
`n_points`	`int`	Number of time points for output grid.	`200`
`progress_callback`	`ProgressCallback \| None`	Optional callback for progress updates. Called with (n_completed, max_n, relative_sem, converged).	`None`

Returns:

Type	Description
`EnsembleResult`	EnsembleResult with ensemble statistics and convergence info.

Point Processes¶

Core abstractions¶

process ¶

Trajectory dataclass ¶

names property ¶

__getitem__ ¶

ConvergenceConfig dataclass ¶

ContinuousTimeProcess ¶

time abstractmethod property ¶

rates abstractmethod ¶

execute abstractmethod ¶

is_absorbing abstractmethod ¶

scalar_output abstractmethod ¶

should_accept ¶

trajectory abstractmethod ¶

ProcessFactory ¶

convergence_config abstractmethod property ¶

create abstractmethod ¶

Gillespie engine¶

gillespie ¶

ProgressCallback module-attribute ¶

EnsembleResult dataclass ¶

scalar_output_mean property ¶

scalar_output_std property ¶

ConvergenceInfo dataclass ¶

ConvergenceMonitor ¶

n property ¶

mean property ¶

variance property ¶

std property ¶

sem property ¶

relative_sem property ¶

update ¶

is_converged ¶

run_once ¶

simulate ¶

Trajectory `dataclass` ¶

names `property` ¶

getitem ¶

ConvergenceConfig `dataclass` ¶

time `abstractmethod` `property` ¶

rates `abstractmethod` ¶

execute `abstractmethod` ¶

is_absorbing `abstractmethod` ¶

scalar_output `abstractmethod` ¶

trajectory `abstractmethod` ¶

convergence_config `abstractmethod` `property` ¶

create `abstractmethod` ¶

ProgressCallback `module-attribute` ¶

EnsembleResult `dataclass` ¶

scalar_output_mean `property` ¶

scalar_output_std `property` ¶

ConvergenceInfo `dataclass` ¶

n `property` ¶

mean `property` ¶

variance `property` ¶

std `property` ¶

sem `property` ¶

relative_sem `property` ¶