molfoundry.Stochsim

class molfoundry.Stochsim(*, graphDatabase, expandStrategy, initialState, labelSettings=None, reactionRate=None, inputRate=None, outputRate=None, drawTime=None, retainHistory=True, onStep=None, seed=42)

Bases: TSObject

Rule-based Gillespie stochastic simulation with on-demand network expansion.

Implements Herrera Machado et al., “Rule-Based Gillespie Simulation of Chemical Systems” (2025): the reaction network is not enumerated up front but grown lazily as new species appear, while an exact SSA evolves the species counts.

Rates may be constants (evaluated once) or callbacks. Following MØD’s DrawMassAction, a callback returns either a rate or a (rate, cache) pair and – unless it caches – is re-evaluated every step, so rates may depend on the current time or state(). The underlying simulator applies the Law of Mass Action on top of the returned rate constant. The waiting-time distribution is customizable via drawTime (default: exponential), mirroring MØD’s DrawTimeExponential.

Create a rule-based Gillespie simulation.

retainHistory (default True) keeps the whole trajectory in memory, so trajectory() / series() and print() work directly. Set it to False for long / large runs to bound memory to O(species): the engine then holds only the start and current step. The full trajectory is instead streamed to a history file under ./out (adjacent to ./summary) so print() can still rebuild the interactive playback; trajectory() / series(), which need the in-memory log, are unavailable in that mode.

onStep is an optional callback invoked once for the start step and then for every step as it is produced, receiving a StepInfo (index, time, reaction). It fires synchronously as the step is recorded, so it may read the live state via this simulation (e.g. state()) — the streaming hook that pairs with retainHistory=False.

Parameters:

• graphDatabase (List[Graph])

• expandStrategy (Any)

• initialState (Dict[Graph, int])

• labelSettings (LabelSettings | None)

• reactionRate (float | int | Tuple[float, bool] | Callable[[...], Any] | None)

• inputRate (float | int | Tuple[float, bool] | Callable[[...], Any] | None)

• outputRate (float | int | Tuple[float, bool] | Callable[[...], Any] | None)

• drawTime (Callable[[float, Callable[[], float]], float] | None)

• retainHistory (bool)

• onStep (Callable[[Stochsim.StepInfo], Any] | None)

• seed (int)

setReactionRate(rate)

Set the rate for each reaction (default: 1 per hyperedge).

rate is a constant (float, or a (rate, cache) pair) or a callback (DG.HyperEdge) -> float | (float, bool). A callback is re-evaluated every step unless it returns cache=True, so reaction rates may depend on time or state(). The engine applies the Law of Mass Action on top of the returned rate constant.

Parameters:

rate (float | int | Tuple[float, bool] | Callable[[...], Any])

Return type:

None

setInputRate(rate)

Set the input-flow (∅ -> species) rate (default: 0, closed system).

rate is a constant or a callback (DG.Vertex) -> float | (float, bool) with the same re-evaluation semantics as setReactionRate().

Parameters:

rate (float | int | Tuple[float, bool] | Callable[[...], Any])

Return type:

None

setOutputRate(rate)

Set the output-flow (species -> ∅) rate (default: 0, closed system).

rate is a constant or a callback (DG.Vertex) -> float | (float, bool) with the same re-evaluation semantics as setReactionRate().

Parameters:

rate (float | int | Tuple[float, bool] | Callable[[...], Any])

Return type:

None

setDrawTime(drawTime)

Set the waiting-time strategy, or None to restore the default Gillespie exponential draw ln(1 / uniform()) / totalPropensity.

drawTime is called as drawTime(totalPropensity, uniform), where uniform() yields independent draws in [0, 1) from the simulator’s own seeded RNG, and must return the time increment until the next reaction. Must be set before the simulation starts.

Parameters:

drawTime (Callable[[float, Callable[[], float]], float] | None)

Return type:

None

setOnStep(onStep)

Set a per-step callback, or None to clear it.

onStep receives a StepInfo once for the start step and then for every step produced. Set it before simulate() so the start step is seen; the callback runs synchronously as the step is recorded and may read the live state through this simulation.

Parameters:

onStep (Callable[[StepInfo], Any] | None)

Return type:

None

simulate(time=None, iterations=None, advanceToEndTime=False)

Advance the simulation up to a time bound and/or iteration bound.

Both bounds are relative to the current state, so repeated calls continue. With neither bound the simulation runs until it deadlocks.

With a time bound, advanceToEndTime parks the clock at exactly that time when no reaction fires before it (the drawn event overshot, or the system is momentarily dead), appending a no-reaction marker step at the bound. This lets a later simulate segment re-evaluate a time-dependent rate at that instant – reviving an otherwise dead system – or a scheduled intervention (e.g. state()-dependent rate, or a follow-up action) land at exactly that time. Off by default, so the clock otherwise stops at the last event.

Parameters:

• time (float | None)

• iterations (int | None)

• advanceToEndTime (bool)

Return type:

None

property dg: DG

property iteration: int

property time: float

state(graph)

Current molecule count of a species graph.

Parameters:

graph (Graph)

Return type:

int

trajectory()

Full event-by-event trajectory of the simulation.

Return type:

Trajectory

series(graph)

Per-step count series of graph aligned with trajectory().times.

Parameters:

graph (Graph)

Return type:

List[int]

print(printer=None, *, name=None, dark='dynamic')

Write an interactive HTML playback of this simulation to the summary folder.

Produces one standalone, seekable page (<name>.html next to summary.html) that lays out the final derivation graph as oval-framed molecule depictions, shows each molecule’s current count below it, and reveals nodes/hyperedges at the step they were first expanded. Play / pause / stop and a speed control step through the whole simulation.

printer is an optional GraphPrinter – exactly as for DG.print() – threaded to every molecule depiction so its drawMode / collapseHydrogens settings control how the molecules are drawn. name sets the filename (default stochsim). Like the rest of the summary the file is flushed when the process exits; returns the target filename.

With retainHistory=False the trajectory is rebuilt from the streamed ./out history file rather than the (absent) in-memory log.

Parameters:

• printer (GraphPrinter | None)

• name (str | None)

• dark (Any)

Return type:

str

class StepInfo(index, time, reaction)

Bases: object

One step handed to an onStep callback.

Parameters:

• index (int)

• time (float)

• reaction (str | None)

index

Absolute step index in the trajectory (0 is the start step).

time

Simulation time at this step.

reaction

Name of the reaction that fired, or None for the start step and no-reaction time markers.

class Trajectory(times, vertexIds, counts)

Bases: object

Parameters:

• times (List[float])

• vertexIds (List[int])

• counts (List[List[int]])

times

Time stamp of each recorded step.

vertexIds

DG vertex index of each column of counts.

counts

counts[step][column] molecule counts.

series(vertexId)

Parameters:

vertexId (int)

Return type:

List[int]