API reference ¶

Raises:

ImportError – If the format’s optional dependency is not installed (antimony for .ant — pip install 'bngsim[antimony]').
FileNotFoundError – If the file does not exist.
ModelError – If the suffix is not a loadable format, the file cannot be parsed, compartment_sizes is given for a .net/.bngl model, or BNG2.pl is unavailable for a .bngl one.

classmethod from_antimony(path, *, compartment_sizes=None)[source]¶

Load a model from an Antimony .ant file.

Antimony is a human-readable model description language. Internally converts to SBML via libantimony, then loads via libsbml for correct SBML semantics. A text with more than one module is loaded at its main module, as libantimony names it (issue #1508); see the user guide’s “Antimony and SBML model loading”.

Requires: pip install antimony python-libsbml

Parameters:

path (str | Path) – Path to the .ant file.
compartment_sizes (dict[str, float] | None) – Compartment id → volume, applied to the converted SBML before interpretation (issue #164). See from_sbml().

Returns:

The loaded model.

Return type:

Raises:

ImportError – If antimony or libsbml is not installed.
FileNotFoundError – If the file does not exist.
ModelError – If the file cannot be parsed.

classmethod from_antimony_string(text, *, compartment_sizes=None)[source]¶

Load a model from an Antimony string.

A text with more than one module is loaded at its main module, as from_antimony() loads a file (issue #1508).

Parameters:

text (str) – Antimony model text.
compartment_sizes (dict[str, float] | None) – Compartment id → volume, applied to the converted SBML before interpretation (issue #164). See from_sbml().

Returns:

The loaded model.

Return type:

classmethod from_sbml(path, *, defer_jacobian=None, compartment_sizes=None)[source]¶

Load a model from an SBML .xml file.

Parameters:

path (str | Path) – Path to the SBML file.
defer_jacobian (bool | None) – GH #145 escape hatch. The analytical Functional Jacobian (GH #76) is derived lazily at the first ODE-solve setup by default (None); pass defer_jacobian=False to derive it eagerly at load instead (the pre-#145 behavior, for A/B and safety). BNGSIM_EAGER_JACOBIAN=1 forces eager for every load path.
compartment_sizes (dict[str, float] | None) –
Compartment id → volume, applied to the parsed document before bngsim interprets it (issue #164). This is how a compartment volume is changed: set_param() refuses a compartment-size write, because the size is folded at load into constants the write cannot reach — per-species volume factors, amount-declared initial conditions, mass-action rate constants, SSA propensity volumes, and the emitted RHS/sensitivity sources. Overriding here moves the size before any of those are derived, so the result is the model you would get by editing the size= attribute in the file. Scan or fit a volume by looping over loads:
```
>>> for v in [1.0, 2.0, 4.0]:
...     m = Model.from_sbml("pbpk.xml", compartment_sizes={"Liver": v})
```
Ids are SBML ids as the document carries them (before bngsim’s identifier mangling, after any comp flattening). An initialAssignment on an overridden compartment is dropped, since it would otherwise take precedence. A compartment whose size an assignment rule computes is refused — the rule, not the attribute, is its volume.

Returns:

The loaded model.

Return type:

Raises:

ImportError – If python-libsbml is not installed.
FileNotFoundError – If the file does not exist.
ModelError – If the file cannot be parsed, or compartment_sizes names an unknown compartment, a non-positive size, or an assignment-rule compartment.

classmethod from_sbml_string(text, *, defer_jacobian=None, compartment_sizes=None)[source]¶

Load a model from an SBML XML string.

Parameters:

text (str) – SBML XML text.
defer_jacobian (bool | None) – GH #145 escape hatch (see from_sbml()). Default lazy; pass defer_jacobian=False (or set BNGSIM_EAGER_JACOBIAN=1) to derive the analytical Functional Jacobian eagerly at load.
compartment_sizes (dict[str, float] | None) – Compartment id → volume, applied before interpretation — the supported way to change a volume, since set_param() refuses a compartment-size write (issue #164). See from_sbml().

Returns:

The loaded model.

Return type:

classmethod from_net(path, *, defer_jacobian=None)[source]¶

Load a model from a BNG .net file.

Parameters:

path (str | Path) – Path to the .net file.
defer_jacobian (bool | None) – GH #145 escape hatch. The analytical Functional Jacobian (GH #76) is derived lazily at the first ODE-solve setup by default (None); pass defer_jacobian=False (or set BNGSIM_EAGER_JACOBIAN=1) to derive it eagerly at load instead (pre-#145 behavior, for A/B).

Returns:

The loaded model.

Return type:

Raises:

ModelError – If the file cannot be parsed, or holds none of the parameters, species, functions, reactions and groups blocks, as a 0-byte file, one of comments only, or one whose lines end in a bare carriage return does (issue #1511); or if a line the loader reads is not UTF-8 (issue #1653).
FileNotFoundError – If the file does not exist.

classmethod from_bngl(path, *, protocol=False, bng2_pl=None, net_out=None, timeout=600, cache=True, defer_jacobian=None)[source]¶

Load a rule-based BNGL (.bngl) model, via BNG2.pl (GH #162).

Overloads:

path (str | Path), *, protocol (Literal[False]) = ..., bng2_pl (str | Path | None) = ..., net_out (str | Path | None) = ..., timeout (int) = ..., cache (bool) = ..., defer_jacobian (bool | None) = ... → Model
path (str | Path), *, protocol (Literal[True]), bng2_pl (str | Path | None) = ..., net_out (str | Path | None) = ..., timeout (int) = ..., cache (bool) = ..., defer_jacobian (bool | None) = ... → tuple[Model, ProtocolSpec]

Parameters:

path (str | Path)
protocol (bool)
bng2_pl (str | Path | None)
net_out (str | Path | None)
timeout (int)
cache (bool)
defer_jacobian (bool | None)

Return type:

Model | tuple[Model, ProtocolSpec]

bngsim simulates reaction networks, and BNGL describes rules, so this is not a parser: the model block is handed to BNG2.pl generate_network, and the .net it emits is loaded through from_net(). That makes BNG2.pl a requirement — pip install 'bngsim[bngl]' supplies it (PyBioNetGen bundles it), or point $BNGPATH/$BNG2_PL at a BioNetGen you already have. It is a Perl script, so perl must be on PATH too (macOS and most Linux distributions ship one; stock Windows does not).

The file’s experiment is not executed. A .bngl typically ends in simulate({...}) or parameter_scan({...}), and running the file as written would run the author’s entire experiment just to obtain a network. Its actions before the first generate_network are applied, except file writes, and a relative path one reads (readFile) resolves against the model’s directory. generate_network itself runs with the source’s own max_iter/max_agg/max_stoich, since those are what make an unbounded rule set finite. Everything after it is dropped. A file with no generate_network before its first run stops at that run, the state the run starts from. Pass protocol=True to get the actions back as a ProtocolSpec rather than lose them, and bngsim.run_protocol() to run them as BNG2.pl does (bngsim.run_bngl() loads and runs in one call).

A ConversionWarning names the dropped actions that set a parameter or a concentration the first simulate, parameter_scan or bifurcate starts from, where the loaded model holds another value (issue #1366). With protocol=True there is no warning: the spec returns the actions themselves.

Parameters:

path (str | Path) – Path to the .bngl file.
protocol (bool) – Return (model, ProtocolSpec) instead of just the model, the spec being the parsed actions block (bngsim.convert.parse_bngl_protocol() in non-strict mode, so a fidelity-affecting action is recorded in spec.lossy and warned about rather than raising).
bng2_pl (str | Path | None) – An explicit BioNetGen folder or BNG2.pl path, outranking $BNG2_PL, $BNGPATH, PATH and an installed PyBioNetGen.
net_out (str | Path | None) – Write the generated network here and load it from there, instead of from the cache. The way to keep the .net as an artifact.
timeout (int) – Seconds to allow generate_network (default 600). A rule set whose network is unbounded hits this.
cache (bool) – Reuse a previously generated network for unchanged BNGL (default). The cache is content-addressed under ~/.cache/bngsim/networks (override with $BNGSIM_BNGL_CACHE_DIR), keyed on the flattened model text and the BioNetGen that produced it (BNG2.pl and the Perl2 library it loads), so it cannot go stale. A model that reads another file before generate_network (readFile, readModel, readNetwork, readSBML, readNFspecies) is not cached: its network depends on that file, which the key cannot see (issue #1494). False regenerates into a per-process directory.
defer_jacobian (bool | None) – GH #145 escape hatch, forwarded to from_net().

Returns:

The loaded model, and the actions block when protocol=True.

Return type:

Model or tuple[Model, ProtocolSpec]

Raises:

FileNotFoundError – If the file does not exist.
ModelError – If BNG2.pl (or perl) cannot be found — with every mechanism that was consulted named — or if generate_network fails or times out, with the tail of BNG2.pl’s own output attached.

Notes

Compartmental BNGL loads: BNG2.pl bakes each compartment’s volume into the generated rate constants, exactly as for a hand-generated .net. For the same reason there is no compartment_sizes= here — the volume has to change in the BNGL source, before generation.

Examples

>>> model = bngsim.Model.from_bngl("egfr.bngl")
>>> model, spec = bngsim.Model.from_bngl(
...     "egfr.bngl", protocol=True)
>>> [e.t_end for e in spec.experiments]
[100.0]

time_discontinuity_conditions()[source]¶

The branch conditions this model switches on simulation time with.

What Simulator.run() resolves to crossing times and stops the integration step on, so an if(time() >= 100, ...) window cannot be stepped over (issue #305 for the mechanism, issue #440 for why a BNGL model needs it).

Two sources, in that order. The SBML loader walks the document at load time and registers each condition it finds as a CVODE root, and hands the registered set over verbatim — reading it back here rather than re-deriving one is what keeps the stops a subset of the roots. A .net/BNGL model is built entirely in C++, so its loader has no seam to register a root at and there is nothing to read back; the conditions are recovered from the built model’s own function bodies instead, and only when something asks.

The recovery is cached, keyed by what it reads from parameter values (issue #726): which species pass the unit-rate clock probes at the model’s current values, and the size of the live compartment each such amount-valued species is read through. A BNGL counter fed at k_t is a clock at k_t = 1 and not otherwise, so its thresholds are time conditions at one value and not at another; an SBML substance-only counter fed at its compartment’s size is one only while that size is 1. Cached once, a model first run at k_t = 0 and then at 1 got no stop at a clock-gated window, and one first run at 1 and then at 0.5 got neither a stop nor a root. Every recovery is kept by its key and shared with clones, so a batch whose rows sit at another verdict than the parent’s recovers that verdict once.

A recovered entry is not always a condition. Since issue #869 it can also be a bare step call of a clock, such as floor((time()-t0)/P), with no comparison in it, and since issue #904 a call that bends (max(0, 1 - abs(time()-t0)/w)). Such a call jumps or turns the rate law itself, and each jump and turn is a stop.

The recovery runs for a model that registered conditions too, and what it finds beside them is added. The loader registers comparisons, so a floor pulse or a max pulse standing beside a piecewise in an SBML model had no entry at all, and was stepped over where the same pulse alone was not.

Return type:: tuple[str, ...]

state_switch_root_conditions()[source]¶

The rate-law conditions over model state a plain run roots on (issue #897).

A rate law switched by a comparison that reads the state, such as if((time() >= 4*S) && (time() < 4*S + w), k, 0), has crossings at times the state decides, so no stop can be placed at them before the run. Simulator.run() registers each one as a CVODE root instead, which stops the step at the crossing rather than letting it span a narrow window and miss it. A step or bend call over state, such as floor(T/100) or max(0, 1 - abs(T-50)/w), is rooted where it jumps or bends, since it spells the same pulse with no comparison (issue #1068).

The crossings time_discontinuity_conditions() already covers are left out, and so is every condition in a function that no rate law reads. Past those, the answer depends on the model’s structure alone, so it is cached for the covered set it was derived beside, and a clone inherits it. The covered set moves with a clock’s rate (issue #726): a threshold on a counter is covered while the counter is a clock and has to be rooted once it is not.

Return type:: tuple[str, ...]

table_function_breaks()[source]¶

Where the model’s table functions break (issue #709): for each table with any, its index kind, the observable it is indexed by, the breakpoints at which a run restarts, and those at which it only ends a step.

Simulator.run() stops the integration on each one a table over time, or over a unit-rate counter, reaches in the run. Which breakpoints count is bngsim._switch_sensitivity.table_function_breaks()’s rule: a table of many samples of a smooth signal has none.

Read from the tables alone, so it is kept until a table is added or given new data (add_table_function() drops it), and a clone inherits it. A model with no table pays one attribute read.

Return type:: list

prepare_analytical_jacobian()[source]¶

Derive and attach the analytical Functional Jacobian (GH #76), at most once.

Idempotent (GH #145): the SymPy derivation runs only on the first call; later calls are no-ops guarded by the model’s once-only sentinel. Returns whether the model now carries a complete analytical Jacobian (False if it fell back to finite differences, or was already FD / all-Elementary with the closed-form C++ Jacobian).

This is the lazy-derivation entry point. The Jacobian is consumed only by ODE solves (CVODE’s dense Jacobian, the steady-state Newton solver, and codegen’s analytical-Jacobian emitter), so it is deferred off the model- load path (from_sbml / from_net no longer derive) and triggered at ODE-solve setup. Call it directly to warm a parent template before clone() fan-out: a warmed parent passes the derived terms to clones (which re-compile the derivative ExprTk strings with no SymPy), so parallel fitting derives once, not once per worker.

Return type:: bool

property analytical_jacobian_status: str¶

Whether the model carries a complete analytical Jacobian, in words.

"complete" when every reaction’s derivative is closed-form — a mass-action model from the start, a model with functional rate laws once prepare_analytical_jacobian() attached them; "pending" while a model with functional rate laws has not derived them yet (an ODE solve, or that call, derives them); and "declined: <reason>" when it runs on the finite-difference Jacobian, with the reason the derivation gave — the rate law and what stopped it, the derivation budget, the C++ attach gate’s verdict (the entry it refused, with the probe state and the analytical and finite-difference values; issue #534), or the environment switch. The same reason goes to the bngsim logger at INFO when it happens; this is the readable form a harness can record instead of only the boolean bngsim._bngsim_core.NetworkModel.analytical_jacobian_complete (issue #506).

property analytical_jacobian_decline: dict | None¶

Why the model runs on the finite-difference Jacobian, as a stable category a census can group by (issue #538).

None while analytical_jacobian_status is "complete" or "pending". Otherwise a dict with kind, one of bngsim._jacobian.DECLINE_KINDS, and reason, the status’s text after "declined: ". For the C++ attach gate’s kinds (fd_mismatch, nonfinite_entry, …) it also carries the fields that kind’s verdict sets — of rxn_idx, per_observable, target_idx, row, col, probe, analytical, finite_difference, n_nonfinite and detail — and the species row and col name as row_species / col_species. A field the kind does not set is absent, not a placeholder.

The kind is what to group by: the reason embeds the rate law, names and values, and its wording changes. A clone carries it, as it carries the status.

property last_libsbml_parse_sec: float¶: Wall seconds the SBML loader spent in the libSBML parse phase (readSBML* + document-level error check). 0.0 for a non-SBML model (e.g. Model.from_net). See Simulator.last_libsbml_parse_sec.

property last_interpret_sec: float¶: Wall seconds spent interpreting the parsed libSBML document into the internal _core model (excludes libSBML parse, Jacobian derivation, and codegen). 0.0 for a non-SBML model. See Simulator.last_interpret_sec.

property last_jacobian_sec: float¶: Wall seconds spent symbolically deriving this model’s analytical Functional Jacobian (GH #76). 0.0 until the derivation runs (it is lazy since GH #145, and never runs on the SSA/PSA/NFsim paths). See Simulator.last_jacobian_sec.

clone()[source]¶

Deep copy the model for parallel workers.

Each clone is fully independent — it has its own parameter values, species concentrations, and expression evaluator state. (What a model and its clones do hold together is a memo of analyses of their functions, which are the same for all of them: read-only once made, and kept under a lock.)

This is also the answer to “give me my own copy” in general: all of a Model’s state lives behind one handle into the compiled extension, so neither copy.copy nor copy.deepcopy can produce an independent model, and both refuse rather than hand back an alias (issue #643).

Returns:: An independent deep copy.
Return type:: Model

validate_for_ssa()[source]¶

Return SSA-compatibility issues detected by the SBML loader.

Returns:: One entry per detected construct; empty for SSA-clean models and for models loaded outside the SBML path (Model.from_net, builder).
Return type:: list

See also

bngsim.validate_for_ssa: module-level function with the same body.

set_param(name, value, *, force_override=False)[source]¶

Set a parameter value by name.

Parameters:

name (str) – Parameter name (e.g. “kf”, “Km”).
value (float) – New value.
force_override (bool) – Pin a derived parameter to value regardless of whether its expression currently produces that value, and permanently for this model object — neither reset(), clone(), nor a later write lifts it. The default (False) overrides only while the value differs from the expression’s, which is what makes an ordinary write round-trip (issue #188, see Notes). Use this only when the caller’s contract is “this parameter is an independent input”; the one such caller in bngsim is the function bngsim.jax.differentiable_solve() with flat=True. Has no effect on a parameter that has no defining expression.

Raises:

ParameterError – If the parameter name is not found, or if name is one of frozen_params and value differs from the value it holds.
ValueError – If name is an SBML compartment size this model cannot resolve to a live volume, or a synthesized slot, and value differs from the value it holds. See Notes.

Return type:

Notes

A compartment’s value is folded at load into constants — per-species volume factors, amount-declared initial conditions, mass-action rate constants, SSA propensity volumes, the emitted RHS, and every <initialAssignment> that reads it. Issue #170 put the volume back into each fold it can reach, so most sizes are now ordinary writable parameters; the residue is refused by name rather than half-honoured, because honoring only some folds leaves the model internally inconsistent rather than moving the volume. They are listed in unwritable_compartment_size_params; load at the size you want instead (from_sbml() accepts compartment_sizes={...}). Writing the value it already holds is allowed, so round-tripping a full parameter vector still works.

A synthesized slot — param_is_internal — is refused on the same unchanged-write rule. One of those is a function’s name: the engine keeps the function’s evaluated value in a parameter slot and rewrites it from the function’s own expression before every derivative evaluation, so a write here is discarded at the next one. It used to be accepted — get_param() echoed the new value back and the trajectory did not move (issue #227). Write the parameters the function’s expression reads instead.

Writing a parameter also re-derives every expression-valued parameter that reads it (BNG setParameter semantics), and — issue #79 — every species initial condition that names one of them. A() Stot in a .net species block, or an SBML initialAssignment that is a bare <ci>, declares that species’ initial condition to be the parameter, so a dose scan over a total amount moves it:

>>> model.set_param("Stot", 1e6)
>>> model.get_state()[0]
1000000.0

The declared initial condition always follows, so reset() rebuilds from current parameter values rather than from a snapshot taken at load. The live concentration follows only while the species is still sitting on that initial condition: a species a Simulator.run() has advanced, or that set_concentration() has assigned, keeps its value and picks the new initial condition up at the next reset().

Two cases where the write deliberately does not reach an initial condition: after save_concentrations() with no label (the baseline is now a captured state, which the declared initial condition no longer describes — dose such a protocol with set_concentration()), and for an SBML initialAssignment too complex to be a single parameter reference (2*init_X + offset), which is evaluated once at load.

Writing a derived parameter — one of the symbols in param_is_expression, d in d = d__FREE or a loader-synthesized _rateLaw{N} — overrides its expression rather than moving it: d stops tracking d__FREE, which is BNG’s setParameter semantics. The override is keyed on the value and so is reversible (issue #188): writing back the value the expression produces re-attaches it, and a write of the value the parameter already holds is not an override at all, so set_params(dict(zip(param_names, vec))) round-trips a full parameter vector unchanged. Read the current state off param_is_expression, which goes False for the duration of an override while param_expressions keeps the defining expression.

An override is a structural change and is meant to be: a derived parameter pinned to a literal no longer carries the chain rule from the primaries underneath it, so those primaries lose that reaction’s term from their sensitivity columns, and the generated RHS changes to match. That is the correct derivative of the model you asked for — but it is why an accidental override used to be so hard to see, and why the round-trip above is the one to rely on.

property frozen_params: list[str]¶

Parameters that were folded to a number when the model was built (issues #313, #695, #696, #711).

An SBML document may write a compartment’s size, a stoichiometry or a conversion factor over its parameters, or set an initial value, of a parameter or of a species, by an initialAssignment that bngsim cannot keep symbolic. Each is evaluated once, at load, and the model holds the number. Every parameter that number was read from is listed: the ones the expression names, and the ones that gave each symbol it names its value at load (through another assignment, a rule, a reaction’s kinetic law, a species’ initial value or its compartment’s size). So is what sized a compartment when a species’ declared value was converted with it: a species declared in the unit it is not held in, in a compartment that a rate rule or an event resizes and an initialAssignment starts, or that an assignment rule sizes (issue #742; every name such a rule reads).

A write to one does not move what it was folded into, so set_param() refuses to change one, and a forward-sensitivity column for one is refused, where both used to be silent: the write took the value and left the number, and the column lacked everything that goes through it (an exact 0 where nothing else reads the parameter). Edit the document and load it again to move one.

A compartmental BNGL model has them too (issue #711). BNG2.pl writes each reaction’s volume factor into its rate constant as a number, 0.1*kb for 1/Ve at Ve = 10, and leaves the expression in a comment, unit_conversion=1/Ve. Every parameter such a comment names is listed, with what a derived parameter or a function among them reads: the network holds the number, and it is the BNGL source that has to change.

Empty for most models.

Returns:: Parameter names, in model parameter order.
Return type:: list[str]

property compartment_size_params: list[str]¶

Names of the parameters that are SBML compartment sizes.

set_param() writes these (issue #170): the storage convention a volume decides — the amount↔concentration conversion, an amount-declared initial condition, the mass-action scalar’s Π V^n, the SSA propensity volume — is re-derived from the parameter rather than left at the size the model happened to load at. What they are still not is differentiable: forward sensitivity refuses a d/dV column (issue #170 stage 3), which is why this list is worth having.

A handful cannot be written even so; see unwritable_compartment_size_params.

Empty for .net models, and for any compartment the SBML loader promoted to a species (rate-rule or event-resized): that one is live state, written with set_concentration(), not a parameter.

Returns:: Parameter names, in model parameter order.
Return type:: list[str]

Examples

>>> model = bngsim.Model.from_sbml("pbpk.xml")
>>> gradient_free = [p for p in model.param_names
...                  if p not in set(model.compartment_size_params)]

property unwritable_compartment_size_params: list[str]¶

Compartment sizes set_param() still refuses to change (#170).

The residue of issue #164’s blanket refusal: a compartment whose size an assignment rule recomputes every step (a write would not survive the next evaluation); one whose storage divide a single mass-action scalar shares across two compartments that merely happen to have the same load-time size (that scalar stops being exact the moment they differ); and one an initialAssignment folds into a quantity no parameter holds — a species amount, a reaction rate — since SBML evaluates every such expression once, at load, against the sizes. All three are decided at load and named in the error.

A subset of compartment_size_params, and usually empty — reload with Model.from_sbml(..., compartment_sizes={...}) to move one.

Returns:: Parameter names, in model parameter order.
Return type:: list[str]

get_param(name)[source]¶

Get a parameter value by name.

Parameters:: name (str) – Parameter name.
Returns:: Current value.
Return type:: float
Raises:: ParameterError – If the parameter name is not found.

set_params(params)[source]¶

Set multiple parameters from a dict.

Parameters:

params (dict[str, float]) – Parameter name → value mapping.

Raises:

ParameterError – If any parameter name is not found, or any value cannot be converted to float. Atomic: either all succeed or none do.
ValueError – If any entry is an SBML compartment size being changed (issue #164 — see set_param()). Checked in the validation phase with the names and the values, so the atomicity above holds: a dict with one compartment write in it applies none of its entries.

Return type:

Examples

>>> model.set_params({"kf": 0.5, "kr": 0.1})

Notes

Each write goes through the same path as set_param(), so a parameter that names a species initial condition re-resolves it here too (issue #79).

reset()[source]¶

Reset all species to their initial concentrations.

Parameter values are not reset — only species concentrations.

Equivalent to restore_concentrations() with no label: it returns to the seed initial conditions, or — after an unlabeled save_concentrations() — to that saved snapshot. Named snapshots (save_concentrations(label=...)) are unaffected.

The seed initial conditions are the ones the current parameter values imply, not the ones the network was generated with: a species whose initial condition names a parameter (A() Stot) returns to that parameter’s value as of now (issue #79).

Return type:: None

save_concentrations(label=None)[source]¶

Snapshot the current species concentrations for later restore.

Implements BNG saveConcentrations() / saveConcentrations("name").

Parameters:: label (str | None) – Name for the snapshot. When omitted (or None), this preserves the historical single-slot behavior: the current concentrations become the new baseline initial state, so a subsequent reset() (or restore_concentrations() with no label) returns here rather than to the original .net seed. When a label is given, the snapshot is stored under that name in a separate multi-slot store and does not disturb the default slot; restore it later with restore_concentrations(label). Multiple named states coexist, so a multi-phase protocol (e.g. saveConcentrations("t=0") … saveConcentrations("start_competition")) round-trips faithfully.
Return type:: None

Notes

A named snapshot captures only the species concentrations (the bulk state vector, ordered like species_names); parameters and the current time are not part of it, matching BNG resetConcentrations.

A snapshot taken of a pre-equilibrated state also captures that state’s forward-sensitivity matrix dx/dθ when one is pending, so restore_concentrations() puts the θ-derivative back with the concentrations (issue #81). Without it, a restored equilibrated state would be re-seeded as a fresh start (∂x(0)/∂θ = 0) — wrong, not approximate, since the restored initial condition is a function of θ. The unlabeled form hands the derivative to the new IC baseline instead, so reset() restores it too.

Saving does not change whether a run advanced the state. A state a run advanced without tracking dx/dθ is still x(θ) with no derivative to carry, so a sensitivity run from it raises, from the live state, from the baseline after reset(), and from a restored snapshot alike. Run the equilibration with the sensitivity_params and pass carry_sensitivities=True to measure from it. It used to start from dx(0)/dθ = 0, with no error (issue #1528).

restore_concentrations(label=None)[source]¶

Restore species concentrations from a saved snapshot.

Implements BNG resetConcentrations() / resetConcentrations("name").

A snapshot saved with a forward-sensitivity seed (a pre-equilibrated state — see save_concentrations()) restores that dx/dθ along with the concentrations, so a following run(carry_sensitivities=True) measures from the right seed (issue #81). A snapshot of a state a run advanced restores its carried-over mark too, seed or not, so a sensitivity run from it raises without one (issue #1528). Restoring a snapshot never clears a mark the model already carries, as set_state() does not. Like the state a run left, a restored snapshot of an advanced state is not re-dosed by a later set_param() of an initial-condition parameter; dose it with set_concentration().

Parameters:: label (str | None) – Name of the snapshot to restore. When omitted (or None), restores the default slot — identical to reset() (the seed initial conditions, or the last unlabeled save_concentrations()). When a label is given, restores the named snapshot saved by save_concentrations(label).
Raises:: ModelError – If label is given but no snapshot was saved under that name.
Return type:: None

has_saved_concentrations(label=None)[source]¶

Whether a named concentration snapshot is available to restore.

Parameters:: label (str | None) – When given, reports whether a snapshot saved under that exact name exists. When omitted (or None), reports whether any named snapshot exists. The default (unlabeled) slot is always restorable via reset() and is not reflected here.
Return type:: bool

property saved_concentration_labels: list[str]¶

Sorted names of the currently saved named concentration snapshots.

Does not include the default (unlabeled) slot, which is restored via reset() / restore_concentrations() with no label.

set_concentration(name, value)[source]¶

Set a single species concentration by name.

Parameters:

name (str) – Species name (e.g. "A(b)"). A fixed species’ $ may be given or left out: the model names $A(b) in a .net A(b), and BNG2.pl’s setConcentration ignores the $ too.
value (float) – New concentration value.

Raises:

ModelError – If the species name is not found.

Return type:

Notes

Implements BNG setConcentration("name", value) action.

The value is the start of the next Simulator.run(), of every replicate of Simulator.run_replicates(), and of Simulator.parameter_scan(), whose default reset target is the model’s state at the call. reset() does not keep it, and neither do Simulator.run_batch() (issue #1362), Simulator.steady_state_batch() without states= (issue #1521) or bngsim.jax.differentiable_solve, which start every row from the initial conditions and warn that they do. Call save_concentrations() after it to make it the state reset() returns to.

For forward sensitivity, this is a literal initial condition: ∂x_k(0)/∂θ = 0. It therefore supersedes whatever the species’ declared initial condition depended on — if the .net set it from a parameter (R() R0), that parameter no longer reaches this species’ IC and its seed row is dropped (issue #113). When the value you assign does depend on a differentiated parameter, say so with declare_ic_sensitivity(). On a state a run advanced, the other species are still x(θ): a sensitivity run from it raises, and saving it does not change that (issue #1528).

get_concentration(name)[source]¶

Get a single species concentration by name.

Parameters:: name (str) – Species name, a fixed species’ $ given or left out, as for set_concentration().
Returns:: Current concentration. For an assignment-rule species after a run, the value the run reports for its last row, in the model’s storage convention (amount over the compartment’s static size), before a variable-volume Result rescales it to the live size (issue #1510). Simulator.steady_state() writes no species, so it leaves a rule species as the last run did.
Return type:: float
Raises:: ModelError – If the species name is not found.

get_state()[source]¶

Bulk-copy the full live species-concentration vector (GH #102).

Returns a fresh float64 array of length n_species, ordered like species_names. This is the low-overhead per-step state-exchange primitive for driving bngsim as a reaction kernel from an external orchestrator (e.g. a hybrid SSA/ODE splitting loop): one Python call marshals the entire state, so per-step exchange cost stays negligible next to the ODE solve even at ~100K species.

See also

set_state: the inverse bulk assignment.
species_names: the ordering of the returned vector.

Return type:: ndarray

set_state(state)[source]¶

Bulk-assign the full live species-concentration vector (GH #102).

Parameters:: state (ndarray) – 1-D array of length n_species, ordered like species_names. Copied into the model’s live concentrations; observables and other derived state are recomputed on the next RHS or observable evaluation.
Raises:: ValueError – If state is not 1-D or its length differs from n_species.
Return type:: None

declare_ic_sensitivity(sens)[source]¶

Declare ∂x_k(0)/∂θ for initial conditions a per-point hook assigns.

Call this from a Simulator.parameter_scan() on_point hook when the value it assigns is computed from a parameter being differentiated — a ligand dose converted with a fitted volume, say. The scan then uses the declared row verbatim instead of measuring it through the hook, which is both exact and cheaper (a declared species is not probed).

Parameters:: sens (Mapping[str, Mapping[str, float]]) – {species_name: {param_name: d_ic_d_param}}. A declared species’ row is taken as fully specified: parameters not named get ∂x_k(0)/∂θ = 0. Declaring {} for a species therefore pins its whole row to zero (an explicitly θ-independent literal). Repeated calls merge per species; a species declared twice keeps the last row.
Raises:: ModelError – If a species or parameter name is unknown.
Return type:: None

Notes

Also honoured by a plain Simulator.run() — an initial condition assigned by hand is no longer described by the .net expression the parameter-graph seeding (issue #43) differentiates, so for a hand-assigned θ-dependent IC this declaration is the only way the engine can know ∂x_k(0)/∂θ.

A declaration lasts until the initial condition it describes changes: assigning the same species again (set_concentration()), or a wholesale reset() / set_state(), drops it. The scan primitive additionally clears declarations before each on_point call, so a hook declares afresh per point and points stay independent.

Nothing needs declaring for the ordinary dose: a literal set_concentration inside an on_point hook has ∂x_k(0)/∂θ = 0, which the scan measures on its own.

Examples

>>> def on_point(model, dose_nM):    # dose in molecules, fitted volume
...     v = dose_nM * 1e-9 * NA * model.get_param("Vecf")
...     model.set_concentration("L(r)", v)
...     model.declare_ic_sensitivity({"L(r)": {"Vecf": v / model.get_param("Vecf")}})

effective_ic_sensitivity(params=None)[source]¶

∂x(0)/∂θ — the initial-condition seeding a run from this state would use.

The reader paired with declare_ic_sensitivity(), and the single defining site for the seeding the forward-sensitivity solver is initialized with. Answers from model structure alone — the parameter graph, the live initial conditions and any declarations — with no integration and no simulation, so a frontend can build its gradient routing once at setup (issue #155).

Why a consumer needs it: Result.output_sensitivities() with axis="parameter" is a total derivative,

d_param[θ] = (right-hand-side path) + Σ_k (∂x_k(0)/∂θ)·d_ic[x_k]

so the seeding is already inside the parameter axis. A caller that routes a fitted parameter to several columns and sums them must add an ic-axis term only for the part of ∂x(0)/∂θ reported here as absent — anything present is carried already, and adding it again counts it twice.

This is the effective matrix, not the model file’s: a species retired because an assignment moved it off the initial condition its expression describes (issue #113) is gone, and a species the caller declared (issue #111) carries the declared row in place of the parameter-graph one. It is therefore state-dependent by design — call it on the configured model, before the run whose gradients it describes.

Parameters:

params (Iterable[str] | None) – Restrict to these parameter names. Default: every parameter. Pass the same list as sensitivity_params= to see exactly the rows that run’s parameter columns will carry.

Returns:

{species_name: {param_name: ∂x_k(0)/∂θ}}, keyed by the names species_names / param_names report — the same ids sensitivity_params= accepts. A compound <initialAssignment> lowered by issue #147 reports the original symbols it was written over, never the synthetic _ic_<species> carrier.

A present entry whose value is 0.0 means seeded, and the coefficient is zero at this state — a chain-rule factor that vanishes here but need not at another point. An absent entry means no seeding path at all, and is the one that says “this parameter’s initial-condition term is yours to supply”. The two are different answers; do not collapse them.

Return type:

dict[str, dict[str, float]]

Examples

>>> model.effective_ic_sensitivity(["R0", "kf"])
{'R(r)': {'R0': 1.0}}

property n_species: int¶: Number of species in the model.

property n_reactions: int¶: Number of reactions in the model.

property n_observables: int¶: Number of observable groups in the model.

property n_parameters: int¶: Number of parameters in the model.

property n_functions: int¶: Number of functions in the model.

property function_names: list[str]¶

Names of the model’s functions, in declaration order.

A function is not a parameter, even though each one also names a parameter slot — where the engine stores the value it last evaluated to, rewritten before every derivative evaluation. param_is_internal flags those slots and primary_param_names omits them (issue #227).

property n_events: int¶: Number of events (SBML/Antimony at (...) triggers) in the model.

property conservation_laws: dict¶

Conservation laws detected from the stoichiometry matrix.

Found for every input format (.net, Antimony, SBML, programmatic ModelBuilder). The returned dict has keys:

n_laws: number of independent conservation laws
n_species: number of species the laws are expressed over
dependent / independent: 0-based species index lists
constants: conservation constants evaluated from the initial conditions
coefficients: n_laws x n_species coefficient matrix

Each law is over the stored state, coefficients[k] @ y constant along a run, so coefficients @ rhs(y) is zero at any y. What a reaction between compartments of different size conserves is an amount, so a species’ coefficient there carries its compartment’s size: A in a compartment of size 1 exchanged with B in one of size 2 gives A + 2*B (issue #758). Such a model’s laws follow a size written with set_param(). A species whose share of such a reaction is divided by a size that moves, a concentration in a compartment that a rate rule or an event resizes, is in no law.

Consumed internally by the reduced-space Newton steady-state solver, which needs the independent subspace to sidestep the rank-deficient Jacobian these laws imply.

stoichiometry_matrix(*, sparse=False)[source]¶

The net stoichiometry matrix S, shape (n_species, n_reactions).

S[i, r] is the net change in species i when reaction r fires: products positive, reactants negative, so A + A -> B puts -2 in A’s row. 0-based, in species_names and reaction-index order. For a mass-action .net model rhs(y) == S @ v(y) for the vector of ODE reaction rates, and every conservation law satisfies L @ S == 0 — this is the matrix conservation_laws was row-reduced from. (An SBML model with a reaction between compartments of different size is the exception: its laws are found for the amounts, and satisfy L @ rhs(y) == 0 and not L @ S == 0.)

A $-prefixed boundary species has an all-zero row: the RHS zeroes its derivative and an SSA firing never updates it, so the matrix that describes what the engines do has no entry there, whatever the reaction line says.

Parameters:: sparse (bool) – False (default) returns a dense float64 ndarray. True returns a scipy.sparse.csc_array of the same shape, which is the form to ask for on a large network — the dense matrix of a 10,000-species, 100,000-reaction model is 8 GB. Needs scipy.

Notes

The 1-based C++ entry list is converted once per model and cached; each call assembles a fresh matrix from it, so the caller may modify what it gets back.

rhs(y, t=0.0)[source]¶

dy/dt at state y and time t, under the live parameters.

The same interpreted RHS the CVODE callback evaluates, with the same observable-and-function refresh it runs — observable totals and function-bound rate parameters are recomputed at y for a model that has functions and skipped for a pure mass-action one (the rhs_evaluates_observables gate) — against the model’s current get_param() values. A $-fixed species reports a zero derivative, as it does in a solve.

The model’s stored concentrations (get_state()) are not read and not written; only the caller’s y enters. What does change is the model’s cached observable and function values, exactly as one integrator step would change them.

Parameters:

y (array_like, shape (n_species,)) – The state to evaluate at, in species_names order.
t (float) – The time rateOf, time-dependent functions and table functions read. 0.0 by default.

Returns:

f(t, y), float64.

Return type:

JacobianMatrix or scipy.sparse.csc_array

Raises:

ValueError – If y is not 1-D or its length is not n_species.
ModelError – If rates of change are written in terms of each other’s rateOf in a ring (an algebraic loop, which has no right-hand side).

See also

jacobian: ∂f/∂y at the same point.
propensities: the SSA-convention reaction rates at y.

jacobian(y, t=0.0, *, sparse=False)[source]¶

∂f/∂y at state y and time t — the ODE Jacobian.

J[i, j] = ∂f_i/∂x_j, both indices in species_names order. Built from the model’s closed-form terms whenever they cover every reaction, and from a one-sided difference quotient of rhs() otherwise; the result says which. That is the rule the ODE integrator and the steady-state solver apply under jacobian="auto", and the difference quotient is stepped by the steady-state solver’s own rule, so a caller gets the matrix a solve would have factored at that point rather than a third approximation. A model whose analytical Jacobian is partial (a rate law the differentiator declined) is never handed the partial matrix: it is differenced, and reported as such.

The first call on a model with functional rate laws derives the closed-form terms (prepare_analytical_jacobian()), which an ODE solve would do at its own setup; later calls reuse them.

Where there is no closed form, the matrix is exact all the same, by dual numbers through the interpreter (issue #1330): each rate law and function as the interpreter evaluates it, differentiated branch by branch, with no step. Only a model the dual evaluation cannot follow (a ring of functions that read each other, rateOf) is differenced.

Parameters:

y (array_like, shape (n_species,)) – The state to linearize at.
t (float) – The time, for time-dependent rate laws. 0.0 by default.
sparse (bool) – False (default) returns a JacobianMatrix, a dense (n_species, n_species) float64 ndarray with a source attribute. True returns a scipy.sparse.csc_array over the model’s structural sparsity pattern — the same pattern from either source, so a consumer can allocate for it once — with the same source attribute set on it. Needs scipy.

Returns:

.source is "analytical", "automatic-differentiation" or "finite-difference", the spellings bngsim.SteadyStateResult.sens_jacobian_source uses for the same three cases.

Return type:

Raises:

ValueError – If y is not 1-D or its length is not n_species.
ModelError – If rates of change are written in terms of each other’s rateOf in a ring (an algebraic loop, which has no right-hand side).
SimulationError – If the matrix is exact by dual numbers and y sits where a rate law changes branch (if(X > c, …) at X = c, abs(X) or max(X, 0) at X = 0, a step of floor()) such that an entry is not the same from the two sides (issue #1330): the Jacobian has no value there. Or if the matrix is a difference quotient whose step decides an entry (issue #1002): the entry moves by more than 1e-4 of itself (1e-10 of its column’s largest, at least) when its column is taken again at the species’ own step or at half the first.

propensities(y, t=0.0)[source]¶

Every reaction’s SSA propensity at state y and time t.

One entry per reaction, in reaction-index order: the amount/time propensity a stochastic step samples from, in the SSA volume convention rather than the ODE one. Two things differ from the ODE rate of the same reaction. A repeated reactant takes the falling factorial over its molecule count n when bngsim reads the reaction as mass action — A + A -> B fires at k·n(n−1)/2, where the ODE rate is k·n²/2, with n = x·V for a species stored as a concentration in a compartment of size V (issue #692). An SBML law it does not read as mass action is evaluated as written, still with the falling factorial when the law is a product of its species; one that writes its own combinatorics (k*X*(X-1)/2) is left alone — and a reaction in a compartment of volume V is multiplied by V (Reaction::ssa_volume_factor, or the live value of the compartment-size parameter it names), converting the ODE’s storage-units rate to a per-event rate. For a .net model with volume 1 and no repeated reactant the two conventions agree.

Observable totals and function-bound parameters are refreshed at y first, as the SSA loop refreshes them before its propensity pass, under the same gate rhs() takes. The stored concentrations are not touched. Values are read from y as given — pass molecule counts for a count-valued answer.

Parameters:

y (array_like, shape (n_species,)) – The state to evaluate at, in species_names order.
t (float) – The time, for time-dependent rate laws. 0.0 by default.

Returns:

The propensity vector a(y), float64.

Return type:

Raises:

ValueError – If y is not 1-D or its length is not n_species.

See also

stoichiometry_matrix: the state change each of these reactions makes.

property param_names: list[str]¶: List of all parameter names.

property param_is_expression: list[bool]¶

Per-parameter is_expression flag, parallel to param_names.

True for derived ConstantExpression parameters such as the _rateLaw{N} symbols BNG2.pl emits when a BNGL rate law is a compound expression (e.g. chi*kon). These are not independent knobs — their values are computed from primary parameters and are re-evaluated automatically by set_param().

Derived means the value expression references another of the model’s symbols, which is the condition that makes the chain rule necessary and the line BNG2.pl itself draws between # ConstantExpression and # Constant. A constant written as arithmetic — gamma 1/7, pi 2*asin(1), c6 ln(2)/120 — names nothing, so it differentiates as a leaf exactly like a literal and reports False here (issue #227). It used to report True, on the narrower reading that the value text does not parse as a float, and the recovery rate of SIR.net was consequently missing from primary_param_names.

This is live state, not a property of the declaration: writing a derived parameter a value its expression does not currently produce overrides the expression and flips its entry to False until the override is lifted (issue #188 — see set_param()). param_expressions keeps the defining expression throughout, so an entry that is False here while its expression there names another of the model’s symbols is exactly an overridden derived parameter, as distinct from a genuine primary, whose text is a number (2e7, 1/7) or empty. That pairing is also the whole of what the generated C source depends on: an override moves the chain rule, and nothing else about a parameter write does.

property param_expressions: list[str]¶

Per-parameter defining text, parallel to param_names.

For a derived parameter (param_is_expression), the expression it is re-evaluated from, such as 1/S0. It stays on record while a set_param() override pins the parameter to a literal, which is how an overridden derived parameter is told from a genuine primary.

Anything else has a number or nothing here: a primary read from a .net file keeps its value as the file wrote it (2e7, 1/7), and most SBML primaries and every synthesized slot (param_is_internal) have "". This is the declaration, which a write does not change, so read a parameter’s value with get_param().

property param_is_internal: list[bool]¶

Per-parameter synthesized-slot flag, parallel to param_names.

True for a slot bngsim created for its own bookkeeping rather than one the model declared. Two kinds, and neither is a knob — each is left out of primary_param_names and refused a value-changing set_param():

_V0_<comp> (issue #170) — an SBML compartment’s size as it was at load, which the rate constants in that compartment are normalised against. Moving it rescales those rates while the volume stays put.
a function’s backing slot (issue #227) — where the engine stores what the function last evaluated to. Every function that does not already name a declared parameter has one, so every model with a functions block carries these.

Distinct from param_is_expression: a derived parameter is recomputed from primaries, these are not primaries in the first place. The two flags are disjoint in practice, and primary_param_names is param_names minus both.

property primary_param_names: list[str]¶

List of parameter names that are the model’s independent knobs.

These are the genuine knobs of the model — primary rate constants, initial-condition parameters, etc. Use this when you want to expose the model to an external optimizer or sampler that should treat each parameter as an independent variable; varying a primary via set_param() automatically propagates to derived parameters. Leave out frozen_params where a model has any: those are listed here, and a write that changes one is refused.

Two kinds are left out, one per flag. A derived ConstantExpression (param_is_expression, e.g. _rateLaw{N}) is recomputed from its primaries, so it is not independent of them. And a synthesized slot (param_is_internal) is not a parameter at all: _V0_<comp>, an SBML compartment’s size as it was at load, which the rate constants in that compartment are normalised against — moving it would rescale those rates without moving the volume, so set the compartment size itself, an ordinary writable parameter since issue #170 — or the slot holding a function’s evaluated value, which the engine overwrites before every derivative evaluation (issue #227).

Both exclusions are about the same failure: a coordinate an optimizer cannot move. Issue #227 is where this list stopped being either. It listed every function name — 327 of 327 function-carrying .net models leaked them — so jax.grad spent a coordinate on one and got exactly 0.0, forever. And it dropped a constant written as arithmetic (gamma 1/7), which is a working knob by every other measure, so a fit over this list held the recovery rate of SIR.net fixed with no warning.

property species_names: list[str]¶: List of all species names.

property observable_names: list[str]¶: List of all observable group names.

pure_sink_species()[source]¶

Names of the write-only accumulator species in this network.

A pure sink is a species this network only ever writes to — the degraded / produced / secreted pool a BNGL model carries to count cumulative flux. Its derivative is a non-zero constant for as long as its producing reactions fire, so ||f(y)||_2 / n_species has a floor above tol and Simulator.steady_state() reports failure however long it integrates — even when every other species has settled. These are the species to hand (negated) to steady_state(mask=...).

Detection is purely structural — nothing is measured, no annotation is needed. A species qualifies when all of the following hold:

it is a product of at least one reaction,
it is a reactant of none,
no species’ derivative depends on it, and
it is not a $-prefixed (fixed / boundary-condition) species.

Clause 3 is not implied by clauses 1-2 and is what makes excluding the species provably harmless to the rest of the system: an Elementary rate law reads only its reactants, but a Functional one reads observables, so a product-only species named in an observable that a rate law consumes still feeds back into the dynamics. Clause 4 drops the boundary conditions, whose derivative is zeroed anyway, so they never hold the residual up.

What this does not claim is that the returned species are the only obstacle to convergence, or that a species not returned is settled. It answers exactly one question: which species can be dropped from the convergence test without changing the problem the remaining ones solve.

Returns:: Species names, in species order. Empty for the common case.
Return type:: list[str]

Examples

>>> model.pure_sink_species()
['Ad()']
>>> ss = sim.steady_state(mask=~model.is_pure_sink())
>>> ss.converged
True

is_pure_sink()[source]¶

Boolean mask of the pure-sink species, shape (n_species,).

The array form of pure_sink_species(), so the convergence-test mask is a negation:

ss = sim.steady_state(mask=~model.is_pure_sink())

Returns:: True where the species is a write-only accumulator.
Return type:: ndarray of bool

add_table_function(name, *, file=None, times=None, values=None, index='time', method=None)[source]¶

Give a table function of the model new data, or add a table function.

Under the name of a table function the model has (table_function_names: a tfun the .net declares, or a tfun written inside a function, registered as <function>__tfun<k>), the table takes the new data in place, and every engine reads it from then on: the interpreter, the compiled code, the SSA and the steady-state solvers, and a Simulator made before the call. The new data has to be read over the same index as the table’s, since that is compiled into the model’s code. It keeps the table’s interpolation unless method names another. A clone keeps the data it was cloned with.

Under a new name the table is added, but nothing in the model reads it: a loaded model’s expressions are fixed, and none of them can name a table the model’s source did not declare. A UserWarning says so, and the table changes no run.

Parameters:

name (str) – Table function name (e.g., "cumNcases").
file (str | Path | None) – Path to a .tfun file. Mutually exclusive with times/values.
times (list[float] | None) – X (index) values. Must be used with values.
values (list[float] | None) – Y (function) values. Must be used with times.
index (str) – Index variable name. Default "time". Can also be the name of a parameter, an observable or a function, for a table not read on the clock. New data for a table has to name the table’s own index.
method (str | None) – Interpolation method, "linear" or "step". None (the default) keeps the table’s own for new data, as a .tfun file carries none, and is "linear" for a new table.

Raises:

ModelError – If the file cannot be read or data is invalid; if name is the name of a function of the model that is not a table, whose body a table cannot replace; or if the model’s table of that name is read over another index (issue #772). The model is left as it was.
ValueError – If arguments are inconsistent (e.g., both file and times).

Return type:

Examples

>>> model.add_table_function("cumNcases", file="case_data.tfun")
>>> model.add_table_function("response", file="dose.tfun", index="drug_conc")
>>> model.add_table_function("drive", times=[0, 1, 2], values=[0, 5, 10])

property n_table_functions: int¶: Number of registered table functions.

property table_function_names: list[str]¶: Names of all registered table functions.

class bngsim.Simulator(model, method='ode', *, xml_path='', poplevel=None, gml=None, connectivity=None, nfsim_v1143_compat=False, block_same_complex_binding=True, traversal_limit=None, jacobian='auto', force_dense_linear_solver=False, force_sparse_linear_solver=False, codegen=None, net_path='', sensitivity_params=None, sensitivity_ic=None, sensitivity_method='staggered', strict_ssa=True, reaction_stats=False)[source]¶

Bases: object

Unified simulation interface for ODE, SSA, PSA, and network-free methods.

Parameters:

model (Model) – The model to simulate.
method (str) –
Simulation method:

Deterministic / network-based:
- "ode" — CVODE adaptive BDF integrator (deterministic)
- "ssa" — Variant of Gillespie’s direct method (exact stochastic)
- "psa" — Partial Scaling Algorithm (approximate stochastic). Lin, Feng, Hlavacek, J. Chem. Phys. 150, 244101 (2019). Requires poplevel keyword argument.
Network-free (canonical tokens):
- "nf" — Network-free simulation (default policy; currently routes to nf_reject).
- "nf_reject" — Rejection/null-event algorithm (NFsim-style). Requires xml_path keyword argument.
- "nf_exact" — Exact non-local network-free algorithm (RuleMonkey). Requires xml_path keyword argument.
- "nf_fixed" — Legacy compatibility token. Recognized but unavailable in this build.
Legacy aliases (accepted for compatibility):
- "nfsim" → "nf_reject"
- "rulemonkey" / "rm" → "nf_exact"
- "dynstoc" / "ds" → unavailable compatibility alias
xml_path (str)
poplevel (float | None) – Critical population size N_c for PSA. Required when method="psa". Must be > 1. Larger values are more conservative (less acceleration, less approximation error). Typical values: 100–1000.
gml (int | None)
connectivity (bool | None) – Only used for method="nf" / "nf_reject" / "nfsim". Controls NFsim’s reaction-connectivity optimization at XML initialization. False uses the conservative full membership update path; True enables the inferred dependency-graph path. If omitted, the underlying NFsim wrapper default is used (currently False).
nfsim_v1143_compat (bool) – Only used for method="nf" / "nf_reject" / "nfsim". When true, preserve NFsim v1.14.3’s extra selector draw for same-seed trajectory compatibility with the standalone CLI.
block_same_complex_binding (bool) – Used for method="nf" / "nf_reject" / "nfsim" and method="rm" / "rulemonkey". When True (the default), a bimolecular rule cannot bind two molecules of one complex: BNGL’s semantics, as BNG2.pl’s network applies them. False allows it. With the NFsim methods, a rule whose products are written apart (A.B -> A + B) and that deletes a bond or a molecule does not fire where the molecules it acts on would stay in one complex across products, whatever this setting (see NfsimSession for the cases not yet tested), so with False a gelling model can run far slower: it fills with ring bonds such a rule may never break, and each attempt on one is refused. With method="rm", False also stops RuleMonkey making that test, so on a model with rings such a rule can break a bond inside a ring, which gives wrong counts with no error. RuleMonkey 3.13.0 deprecates False with method="rm": it prints a warning to stderr the first time in a process, and RuleMonkey 4.0 removes the mode. Complex bookkeeping for Species-typed observables is enabled whenever the model declares one, whatever this setting.
traversal_limit (int | None) – Only used for method="nf" / "nf_reject" / "nfsim". NFsim -utl N: universal traversal limit. None (default) lets NFsim auto-compute a suggested limit from the XML.
jacobian (str)
codegen (bool | None) –
For method="ode": when true, use a compiled C RHS, generated from the model bngsim built – whatever it was loaded from (.net, BNGL, SBML, Antimony or the builder). None (default) compiles automatically at or above BNGSIM_CODEGEN_THRESHOLD species (256), where native code beats the interpreter. False runs the interpreter: it never compiles, and it ignores an artifact that an earlier Simulator or Model.clone() left on the model.

For method="ssa": unless False, a small pure mass-action model with no events gets a compiled propensity vector and the recompute-all loop, reported as ssa_diagnostics["propensity_backend"]. False runs the interpreted propensities, on every engine the Simulator builds, run_replicates and run_batch included (issue #1213). Other methods ignore it.
net_path (str) – Deprecated and ignored (issue #803). Every backend reads the model it is given – codegen since step 3, jacobian="jax" since step 4 – so passing it warns and changes nothing.
sensitivity_params (list[str] | None) –
Parameter names to integrate forward sensitivities for, alongside the state ODEs. The result then carries a (n_times, n_species, n_params) sensitivities tensor whose (t, i, k) entry is ∂x_i(t) / ∂p_k evaluated at the baseline parameter values. Only valid for method="ode".

Each name must be a coordinate some write can move. A slot a model function owns — a functions block entry, or an SBML <parameter> an <assignmentRule> defines — is not one: the engine rewrites it from the function’s own expression before every derivative evaluation, so its column can only be identically zero. Naming one raises ValueError rather than returning that zero (issue #329); differentiate the parameters the function’s expression reads instead. Model.primary_param_names omits them.
sensitivity_ic (list[str] | None) – Species names to integrate forward initial-condition sensitivities for. The result carries a (n_times, n_species, n_ic) sensitivities_ic tensor whose (t, i, k) entry is ∂x_i(t) / ∂x_k(0). Useful when fitting IC parameters via chain rule from a Python-side reparameterization (e.g., model.set_concentration("Epo", 10**theta)) without a corresponding model parameter to differentiate against. Requires the codegen sensitivity RHS path; codegen is auto-enabled for any sensitivity workflow. Only valid for method="ode".
reaction_stats (bool) – method="ssa" only. Record, at every output time, each reaction’s cumulative firing count and integrated propensity — the two accumulators a likelihood-ratio (Girsanov) parameter gradient of an ensemble is built from — on Result.reaction_firing_counts and Result.reaction_propensity_integrals (bngsim.girsanov turns them into scores and gradients). Default False; enabling it leaves every trajectory byte-identical. Refused for PSA, whose scaled channels are not the exact process, and for ODE and network-free runs, which have no reaction channels (GH #616).
strict_ssa (bool) –
Only used for method="ssa" / "psa". Default True.

SBML loader records SSA-compatibility issues at load time (e.g. reversible_non_mass_action, assignment_rule_on_reactant). When True the Simulator raises SsaValidationError on any error-severity issue — this is the safe default that prevents fitting workflows from silently consuming wrong dynamics on broken-under-SSA constructs.

Pass False to downgrade most error-severity issues to warnings (logged via logging) and let the Simulator construct anyway. This mirrors roadrunner’s warn-and-run behavior under gillespie integration. Useful when comparing bngsim against roadrunner on the same model, or when the user understands that the dynamics under SSA will be approximate for these constructs.

Two issue codes remain non-overridable even with strict_ssa=False: non_integer_stoichiometry (SSA requires ±1 fire deltas) and fast_reaction (no fast-equilibrium constraint solver).
sensitivity_method (str) –
CVODES corrector strategy for the coupled state + sensitivity system. Both modes integrate state and all sensitivity ODEs as one extended ODE in a single CVODES pass; they differ only in how each integration step’s nonlinear solve is structured:
- "staggered" (default, CV_STAGGERED): advance the state first, then — with the new state in hand — advance the sensitivity ODEs as a separate solve. Two smaller nonlinear solves per step instead of one big one. Often more robust on stiff or large systems; this is CVODES’ / BNGsim’s default.
- "simultaneous" (CV_SIMULTANEOUS): solve state and all sensitivity variables together as one coupled nonlinear system at every step. Often a touch faster per step on small / well-conditioned problems; the per-step solve is larger so it can struggle on stiff or large systems. This is AMICI’s default, so this is the value to use when you want apples-to-apples timing against AMICI.
CVODES has a third mode (CV_STAGGERED1, one parameter at a time) that BNGsim does not currently expose.
force_dense_linear_solver (bool) –
Only used for method="ode". Default False. Force CVODE’s dense direct linear solver even for large, low-density models that would otherwise auto-select sparse KLU. This is orthogonal to jacobian (which selects how the Jacobian is computed) — it overrides only the linear-solver kind. Intended for benchmarking the dense path against KLU on the same model; it has no effect in a build compiled without KLU (already always dense).

Since issue #128 this reaches steady_state() and steady_state_batch() too — specifically their CVODE march, which routes by this same rule and reports the outcome as SteadyStateResult.linear_solver. The KINSOL polish and the dY_ss/dp solve factor a reduced matrix, whose sparsity pattern is a different object from the model’s, and stay dense either way.
force_sparse_linear_solver (bool) –
Only used for method="ode". Default False. The mirror image of force_dense_linear_solver: force sparse KLU even on a model the auto rule would send to the dense solver for being too small (n_species < 50), too dense (Jacobian density >= 10%), or without a structural nonzero in its Jacobian (issue #1296). Those three gates and the LU-fill test below are bypassed — KLU still needs a sparsity pattern and a non-JAX Jacobian — so it is likewise a no-op in a build without KLU. Passing both force flags raises ValueError.

It also keeps KLU where the LU-fill test would leave it. A model of 256 to 5,000 species with a sparse Jacobian but an LU factor estimated to fill 10% or more of n² — typical of rule-derived networks, where a few hub species couple to most complexes — is factored with the BLAS dense solver by default, on a build that has one (bngsim.HAS_LAPACK_DENSE) and with the analytical Jacobian in use. force_sparse_linear_solver=True is the way back to KLU for a model that rule misjudges; model._core.jacobian_sparsity["lu_fill_estimate"] reports the estimate.

Intended for measuring the auto-selection rule against its own alternative: forced-dense shows what KLU buys on large sparse networks, forced-sparse shows KLU’s setup and indexing overhead on the small dense ones. A model that is both too dense to have been graph-colored and without a usable analytical Jacobian has no way to fill a sparse matrix at all; run() raises there rather than quietly reverting to dense.

Applies to steady_state() as well (issue #128) — with one difference: a model whose Jacobian has no structural nonzero is left on the dense solver there even when forced, where run() takes KLU with the analytical Jacobian and raises with jacobian="fd". f(y) ≡ 0 makes such a model a steady state that solves immediately, as it did before #128.

Examples

>>> model = bngsim.Model.from_net("model.net")
>>> sim = bngsim.Simulator(model, method="ode")
>>> result = sim.run(t_span=(0, 100), n_points=101)
>>> result.time.shape
(101,)

>>> ssa = bngsim.Simulator(model, method="ssa")
>>> result = ssa.run(t_span=(0, 100), n_points=101, seed=42)

>>> psa = bngsim.Simulator(model, method="psa", poplevel=300)
>>> result = psa.run(t_span=(0, 100), n_points=101, seed=42)

>>> # Network-free (all equivalent):
>>> nf1 = bngsim.Simulator(model, method="nf", xml_path="m.xml")
>>> nf2 = bngsim.Simulator(model, method="nf_reject", xml_path="m.xml")
>>> nf3 = bngsim.Simulator(model, method="nfsim", xml_path="m.xml")
>>> rm = bngsim.Simulator(model, method="rm", xml_path="m.xml")

METHODS = {'ds', 'dynstoc', 'nf', 'nf_exact', 'nf_fixed', 'nf_reject', 'nfsim', 'ode', 'psa', 'rm', 'rulemonkey', 'ssa'}¶

auto_atol(*, rtol=None, floor=None)[source]¶

Derive a per-species absolute tolerance from this model’s state.

atol[i] = rtol * max(|y_i|, floor), where y is the model’s current species vector — the state the next run() would start from. Every species is therefore resolved to rtol of its own magnitude, which is the statement a scalar atol cannot make about a model whose species span decades.

This is what atol="auto" runs internally. It is public so a caller can see the numbers before committing to them, and adjust: the return value is an ordinary array, and passing a modified copy back as atol= is the supported way to override any single species.

If you need the tolerance to be a constant of the model rather than of the point, do not use this — call bngsim.derive_atol() with the state you chose. Reading the live state is right for a one-off run and wrong for a parameter fit that moves initial conditions: there the vector would be re-derived at every evaluation, and the objective would step wherever the derivation crossed a rounding boundary. Derive once from the nominal state, then pass that same array as atol= for the whole fit.

Parameters:

rtol (float | None) – Relative tolerance to scale by. Defaults to this Simulator’s (set_tolerances, else 1e-8).
floor (float | None) – Magnitude substituted for a species whose current value is smaller — the case where there is no magnitude to scale. Defaults to the smallest strictly positive species value in the model, but not above 1.0 (see bngsim.derive_atol()).

Returns:

n_species absolute tolerances, ordered like Model.species_names.

Return type:

See also

bngsim.derive_atol: the same rule against a state you supply.

Notes

Derived from initial values, this cannot see a species that starts at order one and decays to something tiny — a within-species, over-time mismatch that needs an error-weight function, not a vector. What it removes is the cross-species compromise. The vector returned here is exactly what bngsim.TrackingAtol takes as its ceiling and re-evaluates against the trajectory (CVodeWFtolerances, issue #213).

Examples

>>> atol = sim.auto_atol()
>>> atol[sim.model.species_names.index("IRp")] *= 10   # loosen one
>>> result = sim.run(t_span=(0, 100), atol=atol)

run(t_span=(0.0, 100.0), n_points=101, *, sample_times=None, seed=None, rtol=None, atol=None, max_steps=None, max_step=None, timeout=None, steady_state=False, steady_state_tol=None, carry_sensitivities=False)[source]¶

Run a simulation.

Parameters:

t_span (tuple[float, float]) – (t_start, t_end) time interval.
n_points (int) – Number of output time points (including t_start).
sample_times (list[float] | None) –
Explicit output time points. When provided, overrides t_span and n_points. Must contain at least 2 values, all finite. Values are sorted automatically. A time may repeat, as replicate measurements do ([0, 0, 0, 5, 5, 5]): each copy is a row of its own, and the copies are equal.

The override covers the interval, not just the output grid: the run starts at the first sample time and ends at the last, so sample_times=[5.0, 10.0] integrates over [5, 10] from the model’s initial conditions and never sees [0, 5]. Pass the start time as the first entry whenever the pre-history matters.
seed (int | None) – Random seed for stochastic methods. When omitted (or None), bngsim draws a fresh seed from system entropy so consecutive run() calls produce independent trajectories. Pass an explicit integer for reproducibility, from 0 to 2**64 - 1, or to 2**32 - 1 under NFsim, whose generators read 32 bits (issue #1515); a seed outside its range raises ValueError. The actual seed used is exposed via Result.seed. For method="ode" it orders only simultaneous equal-priority events (GH #242).
rtol (float | None) – Relative tolerance for ODE solver. Default 1e-8.
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) –
Absolute tolerance for the ODE solver. A float, or the default 1e-8 where none is given, is a ceiling: each species takes min(atol, rtol * max(|x_i|, floor)) from the state the run starts from (a run_until() step keeps the one its trajectory started with), where floor is the smallest value of a species the solver integrates, but not under 1e-6 of the largest (issue #1420). A species below the floor keeps its own value wherever rtol of it is more than double precision resolves beside the largest (issue #1439). That is the number itself wherever the number is a tolerance for the species, and in a model of counts, whose species start at 0 or at least 1, it is the scalar on the scalar’s path. In a model in molar units it is what keeps a species at 1e-9 from running under the noise floor. Where the default still cannot scale a species, as in a model whose species all start at 0, bngsim.AtolScaleWarning names it.

A bngsim.ExactAtol holds every species to its value exactly, the scalar tolerance CVODE has always taken (CVodeSStolerances), as a float was read before issue #1420. An array of n_species values, ordered like Model.species_names, is a per-species tolerance (CVodeSVtolerances, issue #196) — entry i is the absolute tolerance for species i. A wrong-length array raises ValueError; it is never broadcast or truncated, because the consequence of guessing is species i held to the number written for species j.

"auto" derives the vector from the model’s own current state as rtol * max(|y_i|, floor) — see auto_atol(), which returns the same array so it can be inspected or adjusted first. Note the current in that sentence: if you need a tolerance that is a constant of the model rather than of the point being integrated — a fit that moves initial conditions is the case — derive it once with bngsim.derive_atol() from the state you chose and pass that array here instead of "auto".

A per-species tolerance is what a model spanning decades of species magnitude needs: one scalar cannot be a tolerance for a species at 1e-09 and one at 1e+01 at the same time, and picking either end makes the model unintegrable or leaves the small species under the noise floor.

"tracking" (or bngsim.TrackingAtol, which is the same thing with the depth and ceiling spelled out) is the third mode, issue #213. It hands CVODE an error-weight function (CVodeWFtolerances) that re-evaluates the vector above against the state being integrated, so a species that starts at order one and decays to something tiny keeps a tolerance that means something for it. A vector says which species; this says when. It costs steps — roughly 4x on a pure decay resolved twelve decades down — and it is the difference between an answer and noise for the species it is about.

Passing a float leaves the run on exactly the code path it was on before, bit for bit. Neither of the other two is the same thing: a constant vector is the same tolerance reaching CVODE through CVodeSVtolerances, and a depth-0 TrackingAtol is the same tolerance again through CVodeWFtolerances. All three compute the same error weights by differently contractable expressions, so they agree to about one ulp rather than exactly.
max_steps (int | None) – Max internal solver steps per output point. Default 10000.
max_step (float | None) – ODE-only. Upper bound on a single internal integrator step (time units). None (default) leaves the step unconstrained, except that a model loaded from SBML with a periodic floor()/modulo dosing schedule auto-applies a bound that keeps the integrator from stepping over a narrow dose pulse (GH #88). Pass an explicit value to override that (or to bound any model); <= 0 disables the bound.
timeout (float | None) – Wall-clock budget in seconds. When set (and positive), the simulator raises bngsim.SimulationTimeout if elapsed wall-clock time exceeds this limit. None or <= 0 disables the budget. Supported on every backend (ODE/SSA/PSA/NFsim/RuleMonkey); RuleMonkey polls every ~1024 SSA events via its upstream cancellation hook. Partial results are not attached to the timeout exception.
steady_state (bool) – ODE-only. When True, the integrator checks ||f(t,y)||_2 / n_species after recording each output point and stops once it falls below steady_state_tol. The returned Result is truncated to only the rows actually integrated (BNG2.pl simulate({steady_state=>1}) parity, i.e. run_network -c). Default False. Result.solver_stats["steady_state_reached"] reports whether the criterion fired before t_end. A model that reads the time, or that has an event, is integrated to the end of its span whatever this says, and is never marked steady (issue #710): the criterion is of the right-hand side at one instant, which for such a model says nothing of the next.
steady_state_tol (float | None) – Tolerance for the steady_state check above. None or <= 0 falls back to atol (matching BNG2.pl, which reuses the integrator atol as the steady-state cutoff) — the scalar atol, also when a per-species vector is in force, since ||f(t,y)||_2 / n_species is one norm over every species and has no per-species reading. Say what “steady” means explicitly when running with a vector.
carry_sensitivities (bool) – ODE-only, pre-equilibration (GH #210, ADR-0052). When True and this run continues a carried-over species state from a prior run() on the same persistent Simulator (a two-phase equilibrate-then-measure protocol with no reset between phases), the forward-sensitivity initial conditions yS(0) are seeded from the prior phase’s final steady-state sensitivity dx_ss/dθ instead of a fresh start. This makes output_sensitivities() correct across the pre-equilibration boundary: the measurement phase’s IC is x_ss(θ), so ∂x(0)/∂θ is the equilibration sensitivity, not zero. Requires the equilibration phase to have been run on the same Simulator with the same sensitivity_params (and no reset). Requesting sensitivities on a carried-over state without this flag raises (no silent wrong derivatives); a fresh single run is unaffected. Default False.

Returns:

Simulation results with time, species, observables.

Return type:

Raises:

SimulationError – If the solver fails.
SimulationTimeout – If timeout is set and the wall-clock budget is exceeded.
StopConditionMet – If a stop condition triggers (partial result attached).
ValueError – If t_span, n_points or sample_times are invalid (a time that is not finite, t_end <= t_start, n_points < 2, fewer than two sample times), or if output sensitivities were requested (sensitivity_params / sensitivity_ic, including the carry_sensitivities path) on a model that contains events. Events reinitialise the CVODE state discontinuously without a matching forward-sensitivity reinitialisation, so derivatives go silently stale at and after an event; bngsim refuses rather than return wrong numbers (GH #205). Discontinuity triggers (forcing pulses / piecewise-time dosing) do not jump state and are unaffected.

run_batch(t_span=(0.0, 100.0), n_points=101, *, params=None, seed=None, rtol=None, atol=None, max_steps=None, max_step=None, num_processors=None, squeeze=False, timeout=None, steady_state=False, steady_state_tol=None)[source]¶

Run a batch of simulations over parameter sets.

For each parameter set: 1. Clone the model (independent copy) 2. Apply parameters via set_params 3. Reset species to initial conditions 4. Run the simulation (GIL released during each run) 5. Collect the result

Step 3 starts each point from the initial conditions (Model.reset()), not from the model’s current state, and warns where the two differ (issue #1362).

Parameters:

t_span (tuple[float, float]) – (t_start, t_end) time interval for each simulation.
n_points (int) – Number of output time points per simulation.
params (Sequence[dict[str, float]] | None) – Parameter sets. Each dict maps parameter names to values.
seed (int | None) – Base random seed for stochastic methods. Simulation 0 uses base_seed itself and simulation i > 0 a seed drawn from numpy’s SeedSequence(base_seed), so batches from nearby bases share nothing (issue #1365); the seeds of one batch are distinct. When omitted (or None), base_seed is drawn fresh from system entropy on each call so consecutive batches produce independent trajectories. The actual per-sim seed is exposed via Result.seed on each result; simulation 0’s is the base, so seed=results[0].seed reruns the whole batch.
rtol (float | None) – Relative tolerance for ODE solver.
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) – Absolute tolerance for the ODE solver — scalar, or one value per species ordered like Model.species_names (issue #196). See run() for the full contract. A vector, "auto" or a tracking tolerance is resolved once for the whole batch: every point integrates at the same tolerance, so the results are comparable across the sweep. "auto", and a tracking ceiling it gives, take each species’ smallest scale over the rows’ starts, so no row is held above a species’ size where it starts. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from. Each row reads it from the model’s initial conditions (where a row that writes none starts) and from its own start, and each species takes the tighter of the two: a row that starts below the model is the run() from there, and no row is held looser than the model’s start holds it (issue #1660). The model’s live state is no row’s start, and after a run, read from there, it held every row to where that run ended.
max_steps (int | None) – Maximum internal solver steps per output point.
max_step (float | None)
num_processors (int | None) – Number of threads for parallel execution. Default None (sequential). The GIL is released during each simulation, so threads parallelize effectively.
squeeze (bool) – If True, return a single Result with 3D arrays (n_sims, n_times, n_cols) instead of a list. Refused while a stop condition is registered (add_stop_condition()): a row one cuts is shorter than the rest, and a stack does not say which stopped (issue #1369).
timeout (float | None)
steady_state (bool) – ODE-only. When True, every simulation in the batch stops early once ||f(t,y)||_2 / n_species falls below steady_state_tol and its Result is truncated to the rows actually integrated (BNG2.pl simulate({steady_state=>1}) / run_network -c parity, applied per parameter point). Default False. Because each point truncates independently, the per-Result row counts may differ; use squeeze=False (the default) when mixing steady-state early-stop with heterogeneous equilibration times. A model that reads the time, or that has an event, is integrated to the end of its span whatever this says, and is never marked steady (issue #710): the criterion is of the right-hand side at one instant, which for such a model says nothing of the next.
steady_state_tol (float | None) – Tolerance for the steady_state check above. None or <= 0 falls back to atol (matching BNG2.pl).

Returns:

One Result per parameter set (list), or a single squeezed Result with 3D arrays if squeeze=True. Each row is checked against the registered stop conditions as run() checks its result: a row that meets one comes back cut at the first sample it holds at, with Result.stop_condition naming it, and nothing is raised (issue #1369).

Return type:

Raises:

SimulationError – If any simulation fails.
ValueError – If params is empty or t_span/n_points are invalid, or squeeze=True while a stop condition is registered.

Examples

>>> param_sets = [{"k1": v} for v in [0.1, 1.0, 10.0]]
>>> results = sim.run_batch(
...     t_span=(0, 100), n_points=101,
...     params=param_sets, num_processors=4,
... )
>>> len(results)
3

>>> batch = sim.run_batch(
...     t_span=(0, 100), n_points=101,
...     params=param_sets, squeeze=True,
... )
>>> batch.species.shape
(3, 101, n_species)

run_replicates(n_replicates, t_span=(0.0, 100.0), n_points=101, *, seed=None, timeout=None, num_processors=None, squeeze=False)[source]¶

Run n_replicates stochastic replicates of the same model.

Unlike run_batch() — a parameter scan that clones the model per point — replicates share identical parameters and differ only in RNG seed, so this reuses a single simulator and calls reset() between replicates instead of cloning + reconstructing one each time. Reusing the simulator also reuses its cached SSA dependency graph (built once), so on low-activity models the per-replicate cost collapses to the actual trajectory work rather than the fixed clone + graph-rebuild overhead.

Every replicate starts from the model’s state at this call, as run() does: a Model.set_concentration() made first is the start of each (issue #1362), and after a run() so is that run’s end state, as it is for a second run(). Call Model.reset() first to start each from the initial conditions. The model is left in the state it was in at the call.

Replicate 0 uses seed_base itself and replicate i > 0 a seed drawn from numpy’s SeedSequence(seed_base) (seed_base resolved once from seed, or from system entropy when seed is None; each value is exposed via the corresponding Result.seed, and run() from it repeats that replicate), matching the seed schedule run_batch() uses across parameter points. Replicate 0’s seed is the base, so seed=reps[0].seed reruns the whole ensemble. It was seed_base + i, so an ensemble from seed 1 was the one from seed 0 shifted by one replicate, 49 of 50 trajectories the same (issue #1365); the derived seeds of nearby bases are unrelated, and the seeds of one ensemble are distinct.

Sequential execution (num_processors None or 1) reuses this simulator directly. Parallel execution clones the model once per worker thread (not per replicate) for thread-safety, each thread reusing its clone across the replicates it handles.

SSA/PSA only — ODE has no replicate concept; use run_batch() for ODE parameter scans.

Parameters:

n_replicates (int) – Number of replicate trajectories (>= 1).
t_span (tuple[float, float]) – As in run() / run_batch(); squeeze=True is refused while a stop condition is registered.
n_points (int) – As in run() / run_batch(); squeeze=True is refused while a stop condition is registered.
seed (int | None) – As in run() / run_batch(); squeeze=True is refused while a stop condition is registered.
timeout (float | None) – As in run() / run_batch(); squeeze=True is refused while a stop condition is registered.
num_processors (int | None) – As in run() / run_batch(); squeeze=True is refused while a stop condition is registered.
squeeze (bool) – As in run() / run_batch(); squeeze=True is refused while a stop condition is registered.

Returns:

One Result per replicate, or a single squeezed Result with 3D arrays (n_replicates, n_times, n_cols) when squeeze=True. A replicate that meets a registered stop condition comes back cut there, as a run_batch() row does (issue #1369).

Return type:

Examples

>>> ssa = bngsim.Simulator(model, method="ssa")
>>> reps = ssa.run_replicates(30, t_span=(0, 100), n_points=101, seed=0)
>>> len(reps)
30

parameter_scan(parameter, par_scan_vals=None, *, par_min=None, par_max=None, n_scan_pts=None, log_scale=False, t_span=(0.0, 100.0), n_points=101, reset_conc=True, reset_to=None, on_point=None, seed=None, rtol=None, atol=None, max_steps=None, max_step=None, timeout=None, steady_state=False, steady_state_tol=None, squeeze=False)[source]¶

Sweep one parameter, running a simulation per value (BNG parameter_scan).

This is the native scan primitive whose reset_conc semantics match BNG2.pl (issue #11). Unlike a hand-rolled loop that re-derives every point’s species from the .net seed initializers, each point here resets to the state at scan invocation (or to a named snapshot) — so a pre-equilibrate → intervene → scan protocol carries its post-intervention state into the sweep faithfully, instead of discarding it.

For each scanned value the model is: reset to the snapshot (when reset_conc), assigned the scanned parameter, passed through the optional on_point hook (for coupled setConcentration overrides that track the scanned parameter), then integrated over t_span.

Supported on the stateful model-backed backends (ODE / SSA / PSA) only; the XML network-free backends have no in-process scan path.

Parameters:

parameter (str) – Name of the parameter to scan. Must exist in the model.
par_scan_vals (Sequence[float] | None) – Explicit values to scan. Mutually exclusive with the par_min / par_max / n_scan_pts triple.
par_min (float | None) – Endpoints of a generated scan range (inclusive).
par_max (float | None) – Endpoints of a generated scan range (inclusive).
n_scan_pts (int | None) – Number of points in the generated range (>= 1).
log_scale (bool) – When generating a range, space points geometrically (requires positive endpoints) rather than linearly. Default False.
t_span (tuple[float, float]) – (t_start, t_end) for each per-point simulation.
n_points (int) – Output time points per simulation (including t_start).
reset_conc (bool) – BNG reset_conc. When True (default), every point resets to the snapshot before applying the scanned parameter — points are independent. When False, points are not reset between values; each continues from the previous point’s end-state (a continuation scan — see bifurcate()).
reset_to (str | None) – Name of a saved concentration snapshot (Model.save_concentrations(label=...)) to reset each point to. When None (default), the reset target is the model’s live state captured at the moment this method is called. Only consulted when reset_conc=True.
on_point (Callable[[Model, float], None] | None) – on_point(model, value) invoked after the reset + scanned-parameter assignment and before integration, for each point. Use it to apply coupled setConcentration overrides whose value tracks the scanned parameter (e.g. a ligand species whose count is value * NA * Vecf) — the model-specific bookkeeping the primitive cannot infer on its own.
seed (int | None) – Base seed for stochastic methods; point 0 uses seed_base itself and point i > 0 a seed drawn from numpy’s SeedSequence(seed_base), as run_batch() does (issue #1365; drawn fresh from entropy when None, and point 0’s Result.seed reruns the whole scan). Ignored for ODE.
rtol (float | None) – Per-simulation solver options, forwarded to run(). atol may be a per-species vector or "auto" (issue #196); that is resolved once, so every point of the scan is held to the same tolerance rather than to one derived from its own initial conditions. With reset_conc=True, "auto" (and a tracking ceiling it gives) takes each species’ smallest scale over the points’ starts, the reset with the scanned value written; a point whose on_point hook writes a species lower holds that species tighter still. A continuation scan’s later points start where the last ended, and its "auto" is read from its first point’s start, the state at the call with the first value written. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from; each point reads it from the invocation state and from its own start, and each species takes the tighter of the two (issue #1660).
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) – Per-simulation solver options, forwarded to run(). atol may be a per-species vector or "auto" (issue #196); that is resolved once, so every point of the scan is held to the same tolerance rather than to one derived from its own initial conditions. With reset_conc=True, "auto" (and a tracking ceiling it gives) takes each species’ smallest scale over the points’ starts, the reset with the scanned value written; a point whose on_point hook writes a species lower holds that species tighter still. A continuation scan’s later points start where the last ended, and its "auto" is read from its first point’s start, the state at the call with the first value written. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from; each point reads it from the invocation state and from its own start, and each species takes the tighter of the two (issue #1660).
max_steps (int | None) – Per-simulation solver options, forwarded to run(). atol may be a per-species vector or "auto" (issue #196); that is resolved once, so every point of the scan is held to the same tolerance rather than to one derived from its own initial conditions. With reset_conc=True, "auto" (and a tracking ceiling it gives) takes each species’ smallest scale over the points’ starts, the reset with the scanned value written; a point whose on_point hook writes a species lower holds that species tighter still. A continuation scan’s later points start where the last ended, and its "auto" is read from its first point’s start, the state at the call with the first value written. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from; each point reads it from the invocation state and from its own start, and each species takes the tighter of the two (issue #1660).
max_step (float | None) – Per-simulation solver options, forwarded to run(). atol may be a per-species vector or "auto" (issue #196); that is resolved once, so every point of the scan is held to the same tolerance rather than to one derived from its own initial conditions. With reset_conc=True, "auto" (and a tracking ceiling it gives) takes each species’ smallest scale over the points’ starts, the reset with the scanned value written; a point whose on_point hook writes a species lower holds that species tighter still. A continuation scan’s later points start where the last ended, and its "auto" is read from its first point’s start, the state at the call with the first value written. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from; each point reads it from the invocation state and from its own start, and each species takes the tighter of the two (issue #1660).
timeout (float | None) – Per-simulation solver options, forwarded to run(). atol may be a per-species vector or "auto" (issue #196); that is resolved once, so every point of the scan is held to the same tolerance rather than to one derived from its own initial conditions. With reset_conc=True, "auto" (and a tracking ceiling it gives) takes each species’ smallest scale over the points’ starts, the reset with the scanned value written; a point whose on_point hook writes a species lower holds that species tighter still. A continuation scan’s later points start where the last ended, and its "auto" is read from its first point’s start, the state at the call with the first value written. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from; each point reads it from the invocation state and from its own start, and each species takes the tighter of the two (issue #1660).
steady_state (bool) – Per-simulation solver options, forwarded to run(). atol may be a per-species vector or "auto" (issue #196); that is resolved once, so every point of the scan is held to the same tolerance rather than to one derived from its own initial conditions. With reset_conc=True, "auto" (and a tracking ceiling it gives) takes each species’ smallest scale over the points’ starts, the reset with the scanned value written; a point whose on_point hook writes a species lower holds that species tighter still. A continuation scan’s later points start where the last ended, and its "auto" is read from its first point’s start, the state at the call with the first value written. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from; each point reads it from the invocation state and from its own start, and each species takes the tighter of the two (issue #1660).
steady_state_tol (float | None) – Per-simulation solver options, forwarded to run(). atol may be a per-species vector or "auto" (issue #196); that is resolved once, so every point of the scan is held to the same tolerance rather than to one derived from its own initial conditions. With reset_conc=True, "auto" (and a tracking ceiling it gives) takes each species’ smallest scale over the points’ starts, the reset with the scanned value written; a point whose on_point hook writes a species lower holds that species tighter still. A continuation scan’s later points start where the last ended, and its "auto" is read from its first point’s start, the state at the call with the first value written. A scalar, the default or the caller’s, is a ceiling that run() reads from the state it starts from; each point reads it from the invocation state and from its own start, and each species takes the tighter of the two (issue #1660).
squeeze (bool) – When True, stack the per-point results into a single Result with 3-D arrays (like run_batch()); otherwise return a list of per-point results. Refused while a stop condition is registered, as run_batch() refuses it.

Returns:

One Result per scanned value (in order), each carrying custom_attrs["scan_parameter"] and custom_attrs["scan_value"]. A squeezed Result when squeeze=True. A point that meets a registered stop condition comes back cut there, with Result.stop_condition naming it, rather than raising StopConditionMet out of the scan (issue #1369).

Return type:

Notes

The persistent model + backend simulator are left as they were before the call: the scanned parameter, the reset-target concentrations and the carried sensitivity state are restored afterward, so a Simulator can be scanned repeatedly (and the returned trajectories, not the live model, are the product).

Output sensitivities across the scan boundary (issue #81). On a Simulator built with sensitivity_params, each point’s forward-sensitivity seed ∂x(0)/∂θ is the carried dx/dθ of the state it starts from — the pre-equilibration’s steady-state sensitivity — because that state, not the model’s seed initial conditions, is the point’s initial condition. So a scan is run with sensitivities only when that state carries a matching dx/dθ:

reset_conc=True — every point restores the reset target’s state and its dx/dθ, then integrates with carry_sensitivities=True. The target must carry one: the live state at invocation does after an equilibration run on this same Simulator, and a reset_to snapshot does when it was saved while that dx/dθ was pending.
reset_conc=False (bifurcate()) — each point continues the previous point’s state and its dx/dθ, so the whole continuation is one differentiable protocol.

An on_point hook assigns the initial condition its point starts from, so for the species it writes the seed is its derivative, not the carried one (issue #111). Row by row, the most specific thing available wins: a row the hook installed wholesale, then one declared with Model.declare_ic_sensitivity(), then one measured through the hook (bngsim calls the hook at perturbed inputs, so a literal dose measures 0, a dose computed from a differentiated parameter measures its true derivative, and an increment of the carried pool measures the carried row plus the dose’s), and otherwise the carried row unchanged. Measuring calls the hook several extra times per point, so it must be a deterministic function of (model, value); declaring a row skips its measurement.

The scanned parameter must not be a sensitivity_params entry (each point overwrites it, which cannot be composed with the derivative carried into the point), sensitivity_ic is not supported across the boundary, an on_point hook may not move a differentiated parameter, and a hook whose assigned initial condition is not differentiable in a parameter (a dose rounded to whole molecules) must declare that row rather than have a difference quotient of the jump reported as a derivative. Every one of those is a raise, never a silently-reseeded gradient. For a sweep whose points start from the model’s own seed initial conditions, use run_batch().

bifurcate(parameter, par_scan_vals=None, *, par_min=None, par_max=None, n_scan_pts=None, log_scale=False, t_span=(0.0, 100.0), n_points=101, seed=None, rtol=None, atol=None, max_steps=None, max_step=None, timeout=None, steady_state=False, steady_state_tol=None, squeeze=False)[source]¶

Continuation scan of one parameter (BNG bifurcate, reset_conc=0).

A parameter_scan() sibling that does not reset concentrations between points: each point continues from the previous point’s end-state, so the sweep traces a branch of steady states as the parameter is stepped. Sweep par_scan_vals up then down (two calls) to expose hysteresis. The first point starts from the model’s live state at invocation.

Accepts the same arguments as parameter_scan() except reset_conc (pinned to False), reset_to, and on_point (which pertain to the per-point reset that continuation omits). See parameter_scan() for the shared parameters.

Parameters:

parameter (str)
par_scan_vals (Sequence[float] | None)
par_min (float | None)
par_max (float | None)
n_scan_pts (int | None)
log_scale (bool)
t_span (tuple[float, float])
n_points (int)
seed (int | None)
rtol (float | None)
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None])
max_steps (int | None)
max_step (float | None)
timeout (float | None)
steady_state (bool)
steady_state_tol (float | None)
squeeze (bool)

Return type:

compute_all_sensitivities(t_span=(0.0, 100.0), n_points=101, *, params=None, chunk_size=2, n_workers=None, rtol=None, atol=None, max_steps=None)[source]¶

Compute full sensitivity tensor via parallel chunked CVODES jobs.

Splits Np parameters into ⌈Np/chunk_size⌉ independent CVODES forward-sensitivity jobs, runs them in parallel via model.clone() + ThreadPoolExecutor (GIL released during C++ CVODE integration), and stitches the partial sensitivity arrays into a complete (n_times, n_species, n_params) tensor.

Benchmarks showed that 2-parameter sensitivity chunks add only ~1.2× overhead for large models (593–1281 species). With ⌈Np/2⌉ parallel jobs, the full sensitivity tensor costs ~1.2× wall-clock of a plain ODE solve — making gradients nearly free with cores.

Parameters:

t_span (tuple[float, float]) – (t_start, t_end) time interval.
n_points (int) – Number of output time points (including t_start).
params (list[str] | None) – Parameter names to compute sensitivities for. Default (None): Model.primary_param_names — the model’s independent knobs, which is a subset of Model.param_names. A derived (expression-backed) parameter is left out because its column is not a coordinate of its own: the primaries underneath it already carry its effect by the chain rule, so the two columns are not independent and Result.gradient() would count the same physical effect twice (issue #203). A synthesized _V0_<comp> is left out because it is bngsim’s record of a compartment’s load-time size rather than a knob at all (issue #170). Both are dropped with a warning naming them, exactly as an unwritable compartment size is. Naming either in params keeps its column — an explicit ask is a statement that you want that derivative on its own terms, which is how bngsim.jax.differentiable_solve(..., flat=True) obtains one.
chunk_size (int) – Number of sensitivity parameters per CVODES job. Default 2, which benchmarking found to work well for large models because 2-parameter chunks add only ~1.2× overhead. Chiefly a throughput knob, but not only one: a chunk’s columns share one CVODES error test, so the grouping can decide whether a marginal column’s step is accepted (issue #243). A chunk that fails as a group is therefore retried column by column before anything is reported — see _recover_failed_chunk().
n_workers (int | None) – Number of parallel threads. Default: min(⌈Np/chunk_size⌉, os.cpu_count()). Set to 1 for serial execution (debugging/profiling).
rtol (float | None) – Relative tolerance for ODE solver.
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) – Absolute tolerance for the ODE solver — scalar, or one value per species ordered like Model.species_names (issue #196). See run(). Resolved once and shared by every chunk, so the columns being stitched back together were error-controlled identically. Each column’s own tolerance is then built from the per-species value of the state it differentiates.
max_steps (int | None) – Max internal solver steps per output point.

Returns:

Simulation result with full sensitivities tensor of shape (n_times, n_species, n_params). The sensitivity_params attribute lists all parameter names in the order they appear in the tensor, and is the authority on the column set — under params=None it is narrower than Model.param_names, as described above. Species trajectories are from the first chunk’s ODE solve (all chunks produce identical trajectories since they share the same model and parameters).

Return type:

Raises:

ValueError – If method is not ‘ode’, or t_span / n_points are ones run() refuses (a time that is not finite, t_end <= t_start, n_points < 2; issue #1523), or params contains unknown names, or the analytical sensitivity RHS every chunk needs cannot be built — codegen=False / BNGSIM_NO_CODEGEN, no codegen backend, or rate laws that do not differentiate to closed form. That is the same hard requirement Simulator(..., sensitivity_params=...) carries (GH #214): the interpreted finite-difference sensitivity path is not reliable at tight tolerances, so it is refused rather than used. Also if params names a slot a model function owns (issue #329) — under params=None those are dropped with the warning below, and naming one gets the same answer as a hard refusal rather than the identically-zero column it would otherwise produce.
SimulationError – If a chunk simulation fails and the column-by-column retry shows at least one of its columns cannot be integrated on its own either (issue #243). A chunk that fails only as a group is recovered from those single-column solves, with a warning, rather than raising.

Examples

>>> model = bngsim.Model.from_net("model.net")
>>> sim = bngsim.Simulator(model, method="ode")
>>> result = sim.compute_all_sensitivities(
...     t_span=(0, 100), n_points=101,
...     chunk_size=2, n_workers=8,
... )
>>> result.sensitivities.shape  # (101, n_species, n_params)
(101, 149, 40)
>>> fim = result.fisher_information(sigma=0.1)
>>> grad = result.gradient(
...     lambda species, time: np.sum((species - data)**2)
... )

Notes

Architecture: Each chunk clones the model (deep copy, thread-safe), creates a fresh CvodeSimulator, runs CVODES with chunk_size sensitivity parameters, and returns its partial (n_times, n_species, chunk_size) sensitivity tensor. The main thread stitches these along axis 2 (the parameter axis).

Why ThreadPoolExecutor works: The GIL is released during C++ CVODE integration (py::call_guard<py::gil_scoped_release>), so threads achieve true parallelism for the compute-intensive portion. Python overhead is negligible (model clone + setup).

Optimal chunk_size: Benchmarks show that chunk_size=2 minimizes per-chunk overhead while keeping thread count manageable. chunk_size=1 works but doubles the number of threads needed.

steady_state(*, tol=1e-09, max_time=1000000.0, method='integration', rtol=None, atol=None, max_steps=None, sensitivity_params=None, mask=None, criterion='scaled')[source]¶

Find the steady state of the ODE system f(y) = 0.

A model with an event, or one that reads the time (time() or a table function indexed by time, in a rate law, a rule or a reported function), is refused with SimulationError (issue #710): the solvers read the right-hand side at t = 0 and fire no event. Integrate such a model with run().

Solver methods:

"integration" (default): CVODE BDF integrated until criterion is met (by default "scaled"; "residual" is BNG2.pl’s ||f(y)||_2 / n_species < tol, run_network -c).
"newton": two-tier integrate-first solver. The same CVODE burst carries the state into the physical root’s basin, then KINSOL polishes. The polish is accepted only once it is seed-stable (agrees across two successively tighter bursts) and carries no eigenvalue in the right half-plane (issue #78, see SteadyStateResult.root_stability); otherwise integration continues. Seeding Newton at the raw initial condition instead can converge to a spurious root of f(y)=0 the trajectory never reaches, or walk a species negative into NaN (GH #27) — hence the burst.
"kinsol": accepted alias for "newton".

Both tiers factor the model’s closed-form Jacobian when it has one and jacobian= asks for it (issue #127): CVodeSetJacFn on the march, KINSetJacFn on the polish, the latter projected onto the polish’s own unknown set. jacobian="fd" pins the difference quotient instead — one RHS evaluation per unknown per Jacobian setup, which is what both tiers paid unconditionally before #127. ss.solver_jacobian_source reports which matrix was factored. The same option still selects how the stability certificate and dY_ss/dp build theirs.

Because "newton" integrates first and only then polishes, it is "integration" plus extra work: across six published dose-response models it cost 1.4-3.9x more wall clock (GH #28), which is why "integration" is the default. Two things still argue for "newton":

A much tighter root. The polish lands near a residual of ~1e-13 where integration stops the moment it crosses tol (~1e-9) — worth having when the steady state feeds a stiff downstream solve.
A tight ``max_time`` budget. Newton reaches tol from a looser burst than integration needs on its own, so when max_time is cut well below the default it can converge where integration runs out of time. At the default max_time=1e6 this does not happen on any model in the benchmark corpus, but at max_time=1e3 several models flip.

Parameters:

tol (float) – Convergence tolerance on ||f(y)||_2 / n_species. Default 1e-9. Both criteria ask it.
criterion (str) –
What converged asks of a state (issue #1001).
- "scaled" (default): the parity test, ||f(y)||_2 / n_species < tol, and the run test. A run taken on from the state for max_time (at most 1e6) moves no species by more than rtol·|y_i| + atol_i: what tells a state the system rests at from one under tol that a run leaves (a model in small units, or a slow one). An unstable root that a run from it does not leave passes (the cubic’s 0.3); on the Newton path SteadyStateResult.root_stability says so. atol_i is the march’s own, never more than rtol times the largest species the solve starts with or reaches. The Newton step δ = −J⁻¹f from the state to the root, within the same tolerances, decides where the solve stops (up to 512 unknowns), which takes the issue’s cubic at a scale of 1e-9 on from 0.894 to its root at 1, and the state it stops at is polished by KINSOL; the step decides nothing about converged. mask says which species are measured, not which move: the step moves every species but a $-fixed one and a masked one that nothing reads or that reads nothing (a counter), which is held; the run moves every species but a $-fixed one. Where neither atol nor set_tolerances() gives one, each species is marched at the smaller of what atol="auto" gives it and the scalar 1e-8: the scalar is no tolerance for a model in small units, nor "auto" for a species that starts large and settles small.
- "residual": ||f(y)||_2 / n_species < tol, BNG2.pl’s test alone, kept for comparisons with run_network. The norm is in the model’s units of concentration per unit time, so a model in molar units, or a slow one, passes it far from its steady state: the cubic X' = -k·x·(x - 0.3)·(x - 1) at a scale of 1e-9 came back converged=True at x = 0.894, 1.003 and 0.131 after one step, where every run ends at 1.
Under "scaled" the run test and the step are reported as SteadyStateResult.continued_run_distance and SteadyStateResult.newton_distance.
max_time (float) – Max integration time for the integration path, and how long the run of criterion="scaled"’s run test is taken on for, up to 1e6. Default 1e6.
method (str) – "integration" (default), "newton", or "kinsol" (alias for "newton").
rtol (float, optional) – Relative tolerance for the CVODE march, and under criterion="scaled" for the run test. Default 1e-8.
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) – Absolute tolerance for the CVODE march — scalar, or one value per species ordered like Model.species_names (issue #196). See run() for the contract. This is the accuracy the march is held to on the way to the root, and under criterion="scaled" the run test’s; tol above is the separate test for having arrived, a single norm over all species.
max_steps (int, optional) – Max CVODE internal steps per output point.
sensitivity_params (list[str], optional) – Parameter names for dY_ss/dp sensitivity. Requires code generation, exactly as run() and compute_all_sensitivities() do — see the sensitivity note below.
mask (array-like of bool, or sequence of str, optional) – Which species the convergence test covers (issue #74). A boolean array of length n_species (True = keep), or the species names to keep. None (default) tests every species, which is the BNG2.pl parity criterion. See the accumulator note below.

Return type:

SteadyStateResult

Notes

Write-only accumulators (issue #74). A species some reaction produces and none consumes — a degraded / produced / secreted pool counting cumulative flux, a common BNGL idiom — has a constant non-zero derivative forever. ||f(y)||_2 / n_species then has a floor above tol and this method reports converged=False however far it integrates, even when every other species is settled to 1e-10. mask restricts the test to the subspace that does have a steady state; Model.is_pure_sink() finds the accumulators structurally, so the recipe needs no hand-listed species:

ss = sim.steady_state(method="newton", mask=~model.is_pure_sink())

Everything still integrates — the excluded species’ equations stay in the RHS and their trajectories come back in ss.concentrations. What the mask restricts is the residual norm (over n_included, so tol keeps its meaning), the KINSOL polish’s unknown set, and the dY_ss/dp linear system; those last two have to follow, because an accumulator contributes a structurally zero Jacobian column and makes both systems singular at every seed. Excluded species come back with a NaN dY_ss/dp row: no steady value, no steady-state gradient.

When a solve fails, ss.unconverged_pure_sinks names any accumulator that was in the test and is carrying flux — that is a structural floor, not a slow tail, and no max_time will move it. Detection is not a convergence verdict, though: A -> B with nothing feeding A has a textbook pure sink and converges anyway, because the flux dies out. The time-course early stop, run(steady_state=True), keeps the unrestricted criterion and takes no mask.

Codegen (issue #63). The solve runs whatever RHS this Simulator’s codegen_backend reports — the compiled .so, the MIR JIT, or the ExprTk interpreter. Until #63 the steady-state path read no codegen option at all, so a Simulator built with codegen=True still solved interpreted; ss.rhs_backend now reports which backend actually ran.

Sensitivity requires codegen. dY_ss/dp = -J⁻¹·(∂f/∂p) used to build both factors from finite differences. It now prefers the model’s analytical Jacobian and the analytical ∂f/∂p the codegen sensitivity RHS emits, and — like run() and compute_all_sensitivities() since GH #214 — refuses rather than silently degrading when codegen is unavailable. ss.sens_jacobian_source / ss.sens_dfdp_source report which path each factor took. Since GH #67 a Functional model whose rate laws are smooth algebra has an analytical ∂f/∂p here too — it is the same emitted bngsim_codegen_sens_rhs, read at yS = 0. What still differences that factor, with a warning, is Michaelis-Menten and the Functional laws carrying a condition or a non-smooth builtin.

The expression block d(func)/dp is compiled code too (issue #75): it is evaluated by bngsim_codegen_output_sens, the GH #198 chain rule a CVODES forward-sensitivity run() already uses, fed the solved dY_ss/dp columns — so a steady-state gradient and a long-run one come from the same evaluator.

Where there is no closed form, the derivatives are exact all the same (issue #1330): the Jacobian, ∂f/∂p and a function’s row are taken by dual numbers through the interpreter, each law and function as it is evaluated, differentiated branch by branch with no step. That covers jacobian="fd", a model whose closed form declines a law, a rate law the sensitivity emitter declines (Michaelis-Menten, a condition, a non-smooth builtin), and a function the compiled chain rule declines (a table function, if(), abs(), max()). At a steady state that sits exactly where a law changes branch, the answer is taken from each side, and refused where the two differ. Only a model the dual evaluation cannot follow (a ring of functions that read each other, rateOf) falls back: its Jacobian and ∂f/∂p to difference quotients checked at more than one step (issue #1002), and its function rows to NaN, which ss.output_sensitivities refuses by name, with the reason. ss.sens_jacobian_source, ss.sens_dfdp_source and ss.sens_output_source say which path each took. The observable block is an exact linear projection through the group factors and never differences anything.

steady_state_batch(params=None, *, states=None, tol=1e-09, max_time=1000000.0, method='integration', rtol=None, atol=None, max_steps=None, n_workers=None, mask=None, newton_burst=True, criterion='scaled', certify=True)[source]¶

Compute steady states for multiple parameter sets, starting states, or both.

A model with an event, or one that reads the time, is refused, as steady_state() refuses it (issue #710).

Parameters:

params (sequence of dict[str, float], optional) – Parameter sets. Omitted with states: every entry runs at this model’s parameters. Without states, each entry starts from the model’s initial conditions (Model.reset()), as run_batch()’s points do, not from its current state, and the batch warns where the two differ (issue #1521): pass states=[model.get_state()] * len(params), or call Model.save_concentrations() first, to start from it. Either start drops a dose an entry writes to an initial-condition parameter, which only the reset applies; for a dose scan after a run, call Model.reset() first.
states (array-like, optional) – Starting states, shape (n_entries, n_species) ordered like Model.species_names (issue #524). Entry i resets its clone and then starts from states[i]; with params, the two have the same length and pair up row by row. A conservation law’s total is the one the starting state holds.
tol (float)
max_time (float)
newton_burst (bool) – With method="newton", False solves from each starting state directly, with no integration burst first, and reports the root’s stability (root_stability, eigenvalues) instead of rejecting an unstable one: what steady_states() enumerates roots with. A $-fixed species is then left out of the unknowns, whose Newton system it would make singular. Default True, the two-tier solver of steady_state().
certify (Union[bool, Literal['start']]) – With newton_burst=False, False skips the stability certificate, a dense eigendecomposition of the reduced Jacobian, and a converged result then has root_stability == "" and no eigenvalues. "start" takes no Newton step: each starting state is the result, converged where it passes tol, and certified there. steady_states() solves its starts with False and certifies each distinct root once with "start" (issue #1191). Default True. Ignored otherwise.
method (str) – "integration" (default), "newton", or "kinsol" (alias for "newton"). See steady_state().
criterion (str) – "scaled" (default) or "residual", as in steady_state() (issue #1001).
rtol (optional) – Solver options for each march, as in steady_state(). atol may be a per-species vector or "auto" (issue #196), resolved once and applied identically to every entry. The "auto" criterion="scaled" takes when no atol is given is each entry’s own, from its starting state.
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) – Solver options for each march, as in steady_state(). atol may be a per-species vector or "auto" (issue #196), resolved once and applied identically to every entry. The "auto" criterion="scaled" takes when no atol is given is each entry’s own, from its starting state.
max_steps (optional) – Solver options for each march, as in steady_state(). atol may be a per-species vector or "auto" (issue #196), resolved once and applied identically to every entry. The "auto" criterion="scaled" takes when no atol is given is each entry’s own, from its starting state.
n_workers (int, optional) – Number of parallel threads.
mask (array-like of bool, or sequence of str, optional) – Which species the convergence test covers, applied identically to every entry (issue #74). See steady_state(). The species set is structural, so it does not move with the parameter set — which is what makes one mask correct for a whole dose scan.

Return type:

list[SteadyStateResult]

Notes

Every entry runs this Simulator’s codegen_backend, the same as steady_state() (issue #63). Each entry resolves the artifact for itself — one dlopen of the shared path, or one JIT compile of the shared source, per solve — so a dose-response sweep pays that per entry rather than once for the batch. r.rhs_backend reports what ran.

steady_states(*, n_starts=200, seed=0, upper=None, decades=4.0, dedupe_rtol=1e-06, tol=1e-09, n_workers=None, mask=None)[source]¶

Enumerate the steady states of this parameter set by multistart Newton (issue #524).

n_starts starting states are sampled on the model’s conservation manifold, KINSOL is run from each with no integration burst (steady_state_batch(states=..., method="newton", newton_burst=False)), and the roots it reaches are clustered. Each distinct root comes with the stability certificate’s verdict, so the saddles a run of the dynamics never rests on are reported too. A start that KINSOL stops short of its root, under tol but at a Newton distance above 1, is solved once more from where it stopped before the roots are grouped (issue #1518).

The root one ordinary steady_state() solve reaches (with the same mask) is solved as one more start, with each held species set as the draws hold it, so the call returns that root too wherever Newton holds it at tol (issue #1189). It is no draw: it adds nothing to multiplicity or basins, and a root only it found has multiplicity 0. The seed is then put on the conservation laws at the model’s totals, as every draw is. Two cases return less. Where Newton cannot hold the root at tol, as at the large rates below, it is missed. And where the solve had moved a held species, it solved another problem: a masked sink at the end of a chain, as the pure-sink mask below makes it, moves its law’s total, and a held species the rates read moves the rates. Newton from the seed then reaches a root of the problem the draws solve, or none, and it need not be the solve’s.

This is a search, not a proof: a root with a small basin, or a saddle KINSOL never lands on, can be missed, and nothing here claims the list is complete. Completeness for polynomial mass-action systems is what homotopy continuation offers (#502). Nor is a continuum of steady states a list of roots: where one root’s spectrum holds a zero eigenvalue (SIR’s line of disease-free states), each start lands somewhere on it and comes back as a root of its own. A write-only accumulator makes Newton’s system singular at every start; exclude it with mask=~model.is_pure_sink().

Sampling: each species a start moves is drawn log-uniformly from decades decades below a bound up to it. For a species in no conservation law the bound is upper, which defaults to ten times the largest concentration of the model’s state, law total, and such species of the steady state an ordinary steady_state() solve reaches (one solve, for the scale), or 10 when all are 0. The species in a law are drawn the same way below the largest total and then moved onto the laws together, at the totals the model’s state holds: to the point that satisfies them nearest the draw in relative entropy, which keeps every species positive however many a law holds. Members of a law therefore start within about decades decades of each other, and a root with one further below its law’s total is reached only if Newton walks there. A $-fixed species, and one mask leaves out, keeps its value in every start and is no unknown of the solve, though a law’s dependent species is rebuilt from the law.

tol is absolute, as in steady_state(). A model at so small a scale that most starts are already under it is refused. Each solve rebuilds, from each conservation law, the member that carries the most of the law’s total where the solve is, so a large total’s round-off lands in a species at its own scale rather than in a small one (issue #1190). Where a root’s rates are themselves large, their own rounding can still exceed tol, and the root is then reported as a failed start: A <-> B at rates 1 and 2 and a total of 1e8 ends at a residual of 5e-9 when Newton starts at the root itself. Below concentrations of about 1e-9, KINSOL’s step limits are absolute too, and most starts fail even with a tol scaled to match.

Parameters:

n_starts (int) – How many starting states to solve from.
seed (int) – Seed for the sampling, so a call is reproducible.
upper (float or array-like, optional) – The sampling bound for a species in no conservation law: a scalar or one value per species. Unused when every species a start moves is in a law.
decades (float) – How many decades below its bound a species is drawn from. Wider reaches roots near zero and spends more starts where Newton diverges: on the Gardner toggle, 6 decades lost 95 of 200 starts and 2 found one of its three roots; 4 found all three and lost 52.
dedupe_rtol (float) – Two roots are the same when every species agrees to within dedupe_rtol of the larger magnitude, or to within a floor for a species at or near 0, where a root at 0 is reached at 0 or a hair beside it (-5.9e-24 on #1001’s cubic at a nanomolar scale): 1e-9 of the largest value the species takes at any converged start, or double precision’s resolution beside the largest magnitude the solves saw, whichever is larger. A species in a conservation law agrees, too, within 1e-9 of its smallest law’s scale, Σ|L·y|: the solve places it no closer than that. The residual tol is no distance in concentration and plays no part: at 1e-9 it merged the 0.2 nM and 0.5 nM roots of a nanomolar switch (issue #1518).
tol (float) – As in steady_state_batch().
n_workers – As in steady_state_batch().
mask – As in steady_state_batch().

Return type:

SteadyStates

add_stop_condition(condition, *, label='')[source]¶

Add a stop condition checked after each simulation.

Parameters:

condition (str | Callable) –
- str: A Python expression evaluated at each sample. It reads the observables by name, a species by its name with ( and ) written as _ (A() reads as A__), and the time as t or time; nothing else is bound, Python builtins included. It is met at the first sample where it is true, or, for a number, > 0. Arithmetic numpy warns of (0/0, x/0, an overflow) is read at the value numpy gives, so a NaN is not met, and one RuntimeWarning names it (issue #1556). Example: "A_tot < 10"
- callable: A function f(result) -> bool called after the simulation. If it returns a true value (in Python’s sense), the stop condition is triggered.
label (str) – Human-readable label for the condition.

Raises:

ValueError – Here, for an expression that reads a name it cannot (a misspelled observable, a parameter, a model function, a builtin such as abs), that uses ^ (bitwise XOR in Python; write ** for a power), or that is a tuple. From the run, for an expression that raises at a sample or gives something other than a truth value or a real number. Each used to be accepted: a name it could not read made it a condition never met, and a tuple one met at the first sample (issue #1509). An exception a callable raises propagates.

Return type:

Notes

A result that meets a condition is cut at the first sample it holds at, and Result.stop_condition names the condition (its label, else its expression). run() raises StopConditionMet with that result. run_batch(), run_replicates() and parameter_scan() (and so bifurcate()) check every row the same way and return the cut rows in their list, beside the rest, raising nothing; they refuse squeeze=True while a condition is registered (issue #1369). Conditions are checked in the order they were added, and the first one met cuts the row where it first holds, even where a condition added later held earlier. A callable is met at the last sample, so the row it cuts keeps every sample. A cut result carries the reported blocks only, as the one run() raises with always has: Result.state refuses on it.

Examples

>>> sim.add_stop_condition("A_tot < 10", label="low_A")
>>> sim.add_stop_condition(
...     lambda r: r.species[-1, 0] < 5,
...     label="very_low_A",
... )

clear_stop_conditions()[source]¶

Remove all stop conditions.

Return type:: None

run_until(t, *, n_points=None, seed=None, rtol=None, atol=None, max_steps=None)[source]¶

Run simulation from current time to t.

This enables interactive simulation: run to a time point, inspect or modify the model, then continue. Supported for stateful model-backed solvers (ODE / SSA / PSA) only.

Parameters:

t (float) – Target time. Must be finite and > current time.
n_points (int | None) – Number of output points. Default: max(2, int(dt)+1).
seed (int | None) – Random seed for stochastic methods. None (default) draws a fresh seed; pass an integer for reproducibility. See Simulator.run for the full contract.
rtol (float | None) – Solver options (ODE only). atol takes the same scalar, per-species vector, "auto" or "tracking" as run() (issues #196, #213); an "auto" vector — including a tracking ceiling — is derived from the state this leg starts from, which for an interactive protocol is where the previous leg ended.
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) – Solver options (ODE only). atol takes the same scalar, per-species vector, "auto" or "tracking" as run() (issues #196, #213); an "auto" vector — including a tracking ceiling — is derived from the state this leg starts from, which for an interactive protocol is where the previous leg ended.
max_steps (int | None) – Solver options (ODE only). atol takes the same scalar, per-species vector, "auto" or "tracking" as run() (issues #196, #213); an "auto" vector — including a tracking ceiling — is derived from the state this leg starts from, which for an interactive protocol is where the previous leg ended.

Returns:

Simulation result for the [current_time, t] interval.

Return type:

Examples

>>> sim.run_until(t=50)        # simulate to t=50
>>> sim.intervene({"k1": 0.0}) # knock out a reaction
>>> result = sim.run_until(t=100)  # continue to t=100

Notes

A stop condition does not stop this leg early. Conditions are evaluated against the completed result, so the backend integrates the whole [current_time, t] interval and leaves the model at t; the StopConditionMet carries the trajectory truncated at the trigger point, and current_time reads t rather than the trigger time (issue #553). To continue from the trigger point instead, take a snapshot() first and restore() it — the partial result cannot be used to rewind the model, because on a model with a promoted state entry it carries only the reported species block.

intervene(params)[source]¶

Apply a perturbation (parameter change) mid-simulation.

Use between run_until() calls to modify the model during an interactive simulation session.

Parameters:: params (dict[str, float]) – Parameter name → value mapping.
Return type:: None

Examples

>>> sim.run_until(t=50)
>>> sim.intervene({"k1": 0.0})  # knock out reaction
>>> sim.run_until(t=100)        # continue

save_concentrations(label=None)[source]¶

Snapshot the model’s current concentrations (BNG saveConcentrations).

Thin delegator to Model.save_concentrations() on the underlying model. With no label this overwrites the default slot (a later Model.reset() returns here); with a label it stores a named snapshot that parameter_scan() can reset each point to via reset_to=label. Use it to capture a post-intervention state between run_until phases and a following scan (issue #11).

Parameters:: label (str | None)
Return type:: None

restore_concentrations(label=None)[source]¶

Restore the model’s concentrations from a snapshot (BNG resetConcentrations).

Thin delegator to Model.restore_concentrations(). With no label this restores the default slot (identical to Model.reset()); with a label it restores that named snapshot. The backend simulator is rebuilt so a subsequent run() / run_until() seeds from the restored state.

Parameters:: label (str | None)
Return type:: None

snapshot()[source]¶

Capture the current simulation state.

Returns a dict that can be passed to restore() to return to this point.

Returns:: Opaque snapshot of model + simulator state.
Return type:: dict

Examples

>>> sim.run_until(t=50)
>>> snap = sim.snapshot()
>>> sim.run_until(t=100)
>>> sim.restore(snap)  # back to t=50

restore(snapshot=None)[source]¶

Restore simulation state from a snapshot.

A snapshot of a state a run advanced restores the model’s carried-over mark (issue #1528): a sensitivity run from it needs a carried seed, and a later set_param of an initial-condition parameter leaves the restored species as they are, as after a run; dose them with Model.set_concentration.

Parameters:: snapshot (dict | None) – A snapshot from snapshot(). If None, restores the most recent snapshot from the internal stack.
Return type:: None

Examples

>>> snap = sim.snapshot()
>>> sim.run_until(t=100)
>>> sim.restore(snap)  # back to snapshot point
>>> sim.restore()      # same (uses internal stack)

get_state()[source]¶

Bulk-copy the live species-concentration vector (GH #102).

Thin delegator to Model.get_state() on the underlying model. After a stateful run_until/run the model holds the final state (the ODE/SSA backends write it back), so this returns the post-step state. It is the per-step get half of driving bngsim as a reaction kernel from an external orchestrator; pair with set_state().

Return type:: ndarray

set_state(state, *, time=None)[source]¶

Bulk-assign the live species-concentration vector (GH #102).

Thin delegator to Model.set_state(). The C++ simulator reads the model’s current concentrations as the initial condition at the start of the next run_until/run, so a bulk set_state between steps is the per-step set half of the kernel exchange (e.g. injecting the SSA-subset coupling species before advancing the ODE subset).

time, when given, also sets current_time (see set_time()), so the state and the clock it belongs to move together — what a predictor-corrector loop needs to redo a step. Both are validated before either is assigned, so a rejected call changes neither.

Parameters:

state (ndarray)
time (float | None)

Return type:

set_time(t)[source]¶

Set current_time, the time the stored state is at.

The next run_until() integrates from here, and everything that reads the clock sees it: time-indexed table functions, time() in rate laws and functions. Paired with set_state() (or set_state(state, time=t)) this rolls a stepped simulation back to a saved point — a predictor-corrector coupling loop advances a step, corrects the coupling input, and redoes the step from the saved state and time. Unlike restore(), it does not rebuild the backend, so the ODE solver’s warm state (and its linear-solver setup) is kept.

Parameters:: t (float) – The new time; must be finite.
Raises:: ValueError – If t is not finite.
Return type:: None

set_tolerances(rtol=1e-08, atol=None)[source]¶

Set ODE solver tolerances for every subsequent call.

Parameters:

rtol (float) – Relative tolerance.
atol (Union[float, Sequence[float], NDArray[double], Literal['auto'], Literal['tracking'], TrackingAtol, ExactAtol, None]) –
Absolute tolerance; omitted, each run takes its own per species, never looser than 1e-8 (see run()), and a steady_state(criterion="scaled") march its own (see steady_state()). A float is read by every run as a ceiling, per species, the same way (see run()); a bngsim.ExactAtol is held exactly (CVodeSStolerances). A sequence of n_species values, ordered like Model.species_names, is a per-species tolerance (CVodeSVtolerances, issue #196); "auto" derives one now from the model’s current state (see auto_atol()) and pins it, so a later state change does not silently re-derive it. A bngsim.TrackingAtol (or "tracking") installs the error-weight function of issue #213, which re-evaluates that vector against the state being integrated; its ceiling is pinned here for the same reason "auto" is.

Each of the three clears the others rather than sitting behind them — the alternative is a call that appears to set the tolerance and changes nothing.

Raises:

ValueError – If a per-species vector is not n_species long, or holds an entry that is not finite and >= 0 (strictly > 0 when it is a tracking ceiling).

Return type:

set_max_steps(max_steps)[source]¶

Set maximum internal solver steps per output point.

Parameters:: max_steps (int) – Maximum steps.
Return type:: None

property method: str¶: Internal dispatch method (‘ode’, ‘ssa’, ‘psa’, ‘nfsim’, or ‘rulemonkey’).

property requested_method: str¶

Original method token as provided by the user.

Useful for logging and reproducibility metadata. For example, if the user passed method="nf", this returns "nf" while method returns the backend dispatch key ("nfsim" for nf/nf_reject).

property model: Model¶: The model being simulated.

property codegen_backend: str¶

The RHS codegen backend this Simulator hands the ODE engine.

Returns one of:

"mir" — in-process MIR micro-JIT (GH #78): the generated C source is JIT-compiled inside C++ (reached only when BNGSIM_CODEGEN_JIT=mir on a MIR-enabled build prepared codegen for this model).
"cc" — native C compiled to a .so by cc and dlopen’d (auto-selected at/above BNGSIM_CODEGEN_THRESHOLD species — 256 by default — or when codegen=True was passed explicitly).
"exprtk" — the ExprTk bytecode interpreter, no native code (the default below the codegen threshold).

This is the backend that actually runs, not a request: it reflects exactly what run() passes the engine — a non-empty JIT source selects MIR, else a non-empty .so path selects cc, else ExprTk (mirroring the opts.codegen_* dispatch). Only meaningful for method="ode"; other backends never codegen and report "exprtk".

property jacobian_strategy: str¶

"analytical", "fd", or "jax".

This is the post-resolution strategy, not the requested jacobian= mode. With jacobian="auto" (the default) the engine uses the analytical Jacobian when the model has one — every Elementary mass-action law, plus Functional laws whose derivatives were symbolically derived within the build-time budget (GH #76/#95) — and finite differences otherwise. So this reports "analytical" only when the analytical Jacobian is genuinely complete and not overridden:

jacobian="fd" → always "fd".
BNGSIM_ANALYTICAL_FUNCTIONAL_JAC=0 (a Functional model loaded with this set never attaches its analytical terms) → "fd".
a derivation that blew the budget and fell back → "fd".
jacobian="auto" whose analytical attempt failed to integrate and fell back to FD at run time (GH #176) → "fd".
jacobian="jax" → "jax".

Mirrors the engine’s callback selection exactly (cvode_simulator.cpp): analytical iff the requested mode is auto/analytical and analytical_jacobian_complete. Only meaningful for method="ode".

Type:: The Jacobian strategy the ODE engine actually uses

property last_codegen_sec: float¶

Wall seconds spent generating/compiling this model’s codegen RHS.

≈0.0 for an ExprTk model (no codegen runs) or a codegen cache hit; the cc compile time on a cold "cc" model; the source-generation time on a "mir" model. Recorded once at setup by the bngsim._codegen.prepare_* entry points (whether codegen ran at model load or in this Simulator), so a single run() exposes the codegen cost directly — no run-twice-and-subtract needed. Purely a setup-time figure; the per-step integration hot path is never instrumented.

property last_libsbml_parse_sec: float¶: Wall seconds the SBML loader spent in the libSBML parse phase (readSBML* + document-level error check) — the shared C++ core both engines use. Recorded once at load; setup-time only, never the hot path. 0.0 for a model not loaded from SBML (e.g. a .net model).

property last_interpret_sec: float¶: Wall seconds spent interpreting the parsed libSBML document into the internal _core model (bngsim’s Python interpretation layer, including the builder.build() core construction; excludes libSBML parse, the analytical-Jacobian derivation, and codegen, which are timed separately). Recorded once at load; setup-time only.

property last_jacobian_sec: float¶: Wall seconds spent symbolically deriving this model’s analytical Functional Jacobian (GH #76, sympy sp.diff), with SymPy already imported — the one-time SymPy import is process-warmup, measured separately, not here. ≈0.0 for an all-Elementary model, an FD fallback, an over-budget derivation (GH #95), or BNGSIM_ANALYTICAL_FUNCTIONAL_JAC=0. A bngsim-only per-model cost — RoadRunner uses a difference-quotient Jacobian and has no analog. Setup-time only; the per-step integration hot path is never instrumented. Recorded wherever the derivation runs (at load today; at first ODE-solve setup after the GH #145 lazy deferral), so this accessor is stable across that change.

property codegen_cache_hit: bool | None¶

Whether this model’s compiled .so was reused from the on-disk codegen cache (~/.cache/bngsim/codegen/).

True — the .so was found in the cache and loaded without recompiling (the cc compile was skipped).
False — no cached .so matched, so it was compiled fresh.
None — no .so was involved at all: an ExprTk model (no codegen) or a MIR model (in-process JIT, which has no on-disk .so cache).

This is the definitive cache signal recorded by the codegen pipeline at the get_cached_so / memo branch — not inferred from last_codegen_sec (a model-based cache hit still spends time on source generation, so a small wall time does not imply a cache hit). Only meaningful for method="ode" with the "cc" backend.

property has_analytic_sens_rhs: bool¶

Whether a forward-sensitivity run on this Simulator uses bngsim’s analytic ∂f/∂p rather than CVODES’ internal difference quotient (issue #438).

This is a property of the (build, model) pair, not of the build. CVodeSensInit1 takes ONE sensitivity-RHS callback for every column, so a single rate law bngsim cannot differentiate declines the analytic ∂f/∂p for the whole model — there is no per-reaction fallback to mix in — and the difference quotient that replaces it costs one extra RHS evaluation per column per step. An N-parameter fit therefore pays roughly N times the sensitivity cost, which on a fit measured in hours is what ends runs: on Smith_BMCSystBiol2013 all 25 columns fell back and every gradient start timed out. Two models on one bngsim get different answers here, and the same model can flip on a derivation-budget timeout (GH #90), so no bngsim.capabilities() key can answer it — the question is per run.

Read off the artifact this run installs: it either exports bngsim_codegen_sens_rhs, the exact symbol the C++ resolves with try_symbol to choose the analytic RHS, or it does not. That is the ground truth whatever the codegen cache did, and it is the same answer _apply_switch_time_sens() and _raise_if_uncompensated_crossing_sensitivities() already act on.

False on a Simulator built for a plain solve is not a decline. Since issues #209/#217 the sensitivity RHS is emitted only when a sensitivity is actually requested, so an artifact built for Simulator(model) carries no analytic ∂f/∂p simply because nobody asked for one; ask a Simulator built with sensitivity_params= (or after compute_all_sensitivities(), which rebuilds), and read sens_rhs_decline_reason for the half that says why.

A run with no codegen artifact at all — codegen=False, an interpreted ExprTk RHS — reports False too, and correctly: the analytic ∂f/∂p is a compiled symbol, so a run that compiles nothing has none and CVODES’ difference quotient carries every column there as well.

property sens_rhs_decline_reason: str | None¶

Why this model is on CVODES’ difference quotient, or None (issue #438).

None means one of two things, and has_analytic_sens_rhs is what separates them: with the analytic ∂f/∂p in use there is nothing to report, and without it there is nothing recorded — a build older than issue #438, a cache directory the note could not be written into, or a plain run that never asked for a sensitivity RHS at all. Absence is never evidence that a model is on the analytic path.

The reason is what makes the verdict actionable. “uses unsupported construct: abs()” tells a modeller precisely what to re-encode, while “no analytic RHS, expect roughly 25x the sensitivity cost” tells them only to wait or give up.

It also carries the one distinction that is about correctness rather than cost. Where the decline is an UncompensatedCrossingReason the difference quotient does not answer the same question more slowly: it integrates the variational equation smoothly through a branch crossing whose time moves, dropping the saltation jump (f⁻−f⁺)·dt*/dθ, so every column is wrong at and after that crossing (GH #146). The returned string is an instance of that class exactly when it applies, so an isinstance check tells “wrong” from “slow”. bngsim itself refuses such a run since GH #414 (_raise_if_uncompensated_crossing_sensitivities()), so on this build the distinction is mostly of interest to a consumer reporting on artifacts built elsewhere.

Resolved from the artifact first and the model second. A decline is derived once, while source is being generated, and since issue #174 the codegen cache key is structural — so a warm cache resolves the .so with no source generated and no decline derived. The reason is written into a small note beside the artifact for that reason, and a run that hits the cache reads it back; a MIR JIT build, which compiles no .so, has the reason its own codegen recorded on the model instead.

property current_time: float¶

The time the model’s stored state is at.

Every solve advances it to where that solve left off, run and run_until alike, because the backend writes the end-of-span concentrations back into the model and this names the state the model then holds (issue #553). It is the point a following run_until integrates from, so it moves even when the solve ends by raising — notably the StopConditionMet raised after the interval was integrated, which truncates the reported result without rewinding the model. parameter_scan rewinds the model when it finishes and puts this back with it.

class bngsim.ReactionKernel(model, *, method='ode', **simulator_kwargs)[source]¶

Bases: object

Drive a bngsim model per step as a pluggable reaction kernel (GH #102).

Parameters:

model (Model) – The model to drive. Its live species concentrations are the kernel’s state vector (ordered like state_names).
method (str) – Simulation method passed to bngsim.Simulator. Default "ode". Stepping (advance()) requires a stateful backend (ode / ssa / psa); nfsim / rulemonkey may be wrapped for introspection but cannot be advanced.
**simulator_kwargs – Forwarded verbatim to bngsim.Simulator (e.g. codegen, jacobian, poplevel, xml_path, sensitivity_params).

See also

ReactionKernel.from_simulator: wrap an already-configured Simulator.
bngsim.Model.get_state: the bulk state primitive the kernel exchanges.

Notes

Not thread-safe — it owns mutable model + simulator state. For parallel workers, build one kernel per bngsim.Model.clone().

classmethod from_simulator(simulator)[source]¶

Wrap an already-constructed bngsim.Simulator.

Use this when the caller has configured solver options the kernel does not surface directly (custom tolerances via set_tolerances, a prepared codegen .so, sensitivity parameters, …). The kernel adopts the simulator as-is, including its current interactive time.

Parameters:: simulator (Simulator) – The simulator to drive.
Return type:: ReactionKernel

get_state()[source]¶

Bulk-copy the live species-concentration vector (GH #102).

Returns a fresh float64 array of length n_species, ordered like state_names. After an advance() this is the post-step state. One O(n_species) Python call — the get half of the per-step kernel exchange.

Return type:: NDArray[double]

set_state(state, *, time=None)[source]¶

Bulk-assign the live species-concentration vector (GH #102).

Parameters:

state (NDArray[double]) – 1-D array of length n_species, ordered like state_names. Copied into the model’s live concentrations; the next advance() reads them as its initial condition. The set half of the per-step kernel exchange.
time (float | None) –
Also set time to this value (Simulator.set_time). With a state saved by get_state() and its time, this redoes a step — the rollback a predictor-corrector coupling loop needs, and the only way to get it right when the model reads the clock (a time-indexed table function, time() in a rate law):
```
>>> saved, t0 = kernel.get_state(), kernel.time
>>> kernel.advance(dt)                    # predictor
>>> kernel.set_state(corrected, time=t0)  # roll back, correct
>>> kernel.advance(dt)                    # corrector
```

Raises:

ValueError – If state has the wrong length or time is not finite; the call then changes neither the state nor the clock.

Return type:

advance(dt, *, n_points=2, seed=None, **run_kwargs)[source]¶

Advance the simulation by dt and return the post-step state.

Integrates (or steps, for stochastic backends) from the current time to time + dt, then returns get_state(). The model is left holding the post-step state, so the typical per-step loop is set_state(...) → advance(dt) → get_state() (or just use the returned array).

Parameters:

dt (float) – Coupling step. Must be > 0.
n_points (int) – Output points recorded over the step. Default 2 (endpoints only), which is all an orchestrator needs and keeps recording cheap; the full sub-step trajectory is available via last_result.
seed (int | None) – Random seed for stochastic backends (ssa / psa). Ignored by ode. None draws a fresh seed each step.
**run_kwargs – Forwarded to Simulator.run_until (e.g. rtol, atol, max_steps).

Returns:

The post-step state vector (a copy; ordered like state_names).

Return type:

Raises:

ValueError – If dt <= 0.
NotImplementedError – If the backend is network-free (nfsim / rulemonkey), which has no stateful per-step API.

reset()[source]¶

Reset to the model’s initial concentrations and time = 0.

Restores species to their initial values (bngsim.Model.reset()) and rewinds the interactive clock, so the kernel can be re-driven from scratch. Clears any cached step result.

Return type:: None

observables()[source]¶

Observable values at the current simulation state.

Returns a float64 array of length n_observables, ordered like observable_names. Right after an advance() these are the post-step observables (read straight from the step result, no recomputation). Otherwise — before the first advance, after a set_state() (with or without time=), or after an advance that raised having moved the state (a stop condition) — they are computed once, side-effect-free, from the current state via a throwaway model clone. last_result is left as it was either way.

Return type:: NDArray[double]

property state_names: list[str]¶: Species names, in the order of the get_state() vector.

property species_names: list[str]¶: Alias for state_names.

property observable_names: list[str]¶: Observable names, in the order of observables().

property n_species: int¶: Number of species (length of the state vector).

property n_observables: int¶: Number of observables.

property time: float¶: Current interactive simulation time.

property method: str¶: Internal backend dispatch key ('ode' / 'ssa' / …).

property model: Model¶: The wrapped model (its live concentrations are the kernel state).

property simulator: Simulator¶: The underlying bngsim.Simulator driving the steps.

property last_result: Result | None¶

The bngsim.Result from the most recent advance().

None before the first advance. Carries the full sub-step trajectory of the last step (per n_points), beyond the endpoint state get_state() returns. It keeps describing that step after a set_state(), when it no longer describes the live state.

class bngsim.UnitConverter(volume_factors, *, names=None)[source]¶

Bases: object

Bulk count ↔ concentration ↔ amount conversion over the state vector.

Wraps the per-species storage→amount factor vector volume_factor (V_c; Species.volume_factor, baked at load time) and converts a whole get_state / set_state storage vector to and from the units an external orchestrator works in. See the module docstring for the unit algebra; the short version is amount = storage * V_c (volume-invariant) and concentration = amount / V (needs a volume).

Build one with from_model() (or from_kernel()); it is immutable and cheap to keep alongside a kernel.

Parameters:

volume_factors (ArrayLike) – Per-species V_c, ordered like the model’s state vector. Must be finite and > 0.
names (Sequence[str] | None) – Species names parallel to volume_factors (for error messages and for_species()). Not required for the array conversions.

classmethod from_model(model)[source]¶

Gather the per-species volume_factor vector from a model/kernel.

Parameters:: model (Model | ReactionKernel)
Return type:: UnitConverter

classmethod from_kernel(model)¶

Alias — a kernel exposes .model, so from_model() already accepts it; kept for symmetry with the rest of the API.

Parameters:: model (Model | ReactionKernel)
Return type:: UnitConverter

property volume_factors: NDArray[float64]¶: Per-species V_c (a copy).

property n_species: int¶: Length of the state vector this converter expects.

property names: list[str] | None¶: Species names parallel to the state vector, if supplied.

to_amounts(storage)[source]¶

Storage vector → molecule amounts (storage * V_c).

Parameters:: storage (ArrayLike)
Return type:: NDArray[double]

from_amounts(amounts)[source]¶

Molecule amounts → storage vector (amounts / V_c).

Parameters:: amounts (ArrayLike)
Return type:: NDArray[double]

to_counts(storage, *, policy='nearest', rng=None)[source]¶

Storage vector → integer molecule counts (rounded amounts).

Parameters:

storage (ArrayLike)
policy (Literal['nearest', 'floor', 'ceil', 'stochastic'])
rng (Generator | None)

Return type:

from_counts(counts)[source]¶

Integer molecule counts → storage vector (counts / V_c).

Parameters:: counts (ArrayLike)
Return type:: NDArray[double]

to_concentrations(storage, *, volume=None)[source]¶

Storage vector → concentrations (storage * V_c / volume).

volume=None returns the at-load concentration (== storage). Pass the current compartment volume to report concentrations against a framework-grown compartment whose baked V_c is now stale.

Parameters:

storage (ArrayLike)
volume (ArrayLike | None)

Return type:

from_concentrations(concentrations, *, volume=None)[source]¶

Concentrations → storage vector (concentration * volume / V_c).

Inverse of to_concentrations(); pass the same volume you would read the concentration against.

Parameters:

concentrations (ArrayLike)
volume (ArrayLike | None)

Return type:

for_species(names)[source]¶

A sub-converter over names (must have been supplied at build).

Parameters:: names (Sequence[str])
Return type:: UnitConverter

class bngsim.CouplingMap(all_names, shared_names)[source]¶

Bases: object

Name ↔ index addressing of the shared coupling subset within a state vector.

Two coupled subsets share species by name but order their state vectors independently. A CouplingMap pins the index of each shared name in one subset’s vector, so the orchestrator gathers / scatters only the coupling subset and exchanges it in a single canonical order across subsets.

Parameters:

all_names (Sequence[str]) – The full state-vector ordering of one subset (its species_names).
shared_names (Sequence[str]) – The coupling subset, in the canonical exchange order. Every name must appear in all_names; duplicates are rejected.

classmethod from_model(model, shared_names)[source]¶

Build from a model/kernel’s species_names.

Parameters:

model (Model | ReactionKernel)
shared_names (Sequence[str])

Return type:

CouplingMap

classmethod from_kernel(model, shared_names)¶

Build from a model/kernel’s species_names.

Parameters:

model (Model | ReactionKernel)
shared_names (Sequence[str])

Return type:

CouplingMap

property names: list[str]¶: The shared species, in exchange order.

property indices: NDArray[int64]¶: Indices of the shared species into the full state vector (a copy).

property n_shared: int¶: Number of shared (coupling) species.

property n_full: int¶: Length of the full state vector this map addresses.

gather(state)[source]¶

Extract the shared subset from a full state vector, in exchange order.

Parameters:: state (ArrayLike)
Return type:: NDArray[double]

scatter(state, values, *, copy=True)[source]¶

Write values (exchange order) into the shared slots of state.

Returns the updated full vector. With copy=True (default) state is not mutated; copy=False writes in place and returns the same array.

Parameters:

state (ArrayLike)
values (ArrayLike)
copy (bool)

Return type:

read(source)[source]¶

Pull the shared subset straight from a live model/kernel’s state.

Parameters:: source (Model | ReactionKernel)
Return type:: NDArray[double]

write(target, values)[source]¶

Inject values into the shared slots of a live model/kernel’s state.

Parameters:

target (Model | ReactionKernel)
values (ArrayLike)

Return type:

class bngsim.DiscreteExchange(n_species, *, policy='nearest', dither=True, nonneg=True, rng=None)[source]¶

Bases: object

Explicit, inspectable rounding at the SSA/NFsim hand-off with leak accounting.

The SSA engine already rounds whatever continuous storage it is handed at the next advance (round_initial_population_to_storage in src/ssa_simulator.cpp), but that entry rounding is implicit and silent — repeatedly injecting continuous amounts and reading integer counts back sheds the fractional remainder every step, which leaks mass over a long hybrid run. DiscreteExchange makes the boundary explicit and conserving: it carries the per-species fractional residual forward (error feedback / dithering), so the discrete counts track the continuous amounts with bounded, accounted error rather than a silent downward drift.

Parameters:

n_species (int) – Length of the amount vectors handed across the boundary.
policy (Literal['nearest', 'floor', 'ceil', 'stochastic']) – Per-step rounding rule (RoundingPolicy).
dither (bool) – Carry the fractional residual into the next step (default True). With dithering the cumulative leak stays bounded by the carry; without it, each step rounds independently and leak can accumulate — set False only when you explicitly want memoryless rounding.
nonneg (bool) – Clamp counts at 0 (default True) so the SSA boundary never receives a negative population. When dithering, a clamp keeps the unrepresentable deficit in the carry, so conservation accounting is undisturbed.
rng (Generator | None) – Used by policy="stochastic".

Examples

>>> dx = DiscreteExchange(3, policy="nearest")
>>> counts = dx.discretize([0.4, 0.4, 0.4])   # rounds down this step
>>> later = dx.discretize([0.4, 0.4, 0.4])    # carry has built up to 0.8 → rounds up
>>> dx.leak                                    # net mass injected vs the continuous input

discretize(amounts)[source]¶

Round amounts to integer counts, carrying the residual forward.

Returns integer-valued float64 counts of length n_species. Updates carry (dithering), last_residual (this step’s counts - amounts) and leak (cumulative net mass discretization has added or removed since construction / reset()).

Parameters:: amounts (ArrayLike)
Return type:: NDArray[double]

property carry: NDArray[float64]¶: Current per-species fractional residual buffer (a copy).

property last_residual: NDArray[float64]¶: counts - amounts from the most recent discretize() (a copy).

property leak: float¶

Cumulative net molecules added (>0) or removed (<0) by rounding.

With dither=True this stays bounded (it equals -carry.sum()); a growing magnitude under dither=False is exactly the silent SSA-entry leak this boundary is meant to surface.

reset()[source]¶

Clear the carry, leak ledger, and last residual.

Return type:: None

bngsim.round_to_counts(amounts, policy='nearest', *, rng=None)[source]¶

Discretize continuous amounts to integer molecule counts (stateless).

Parameters:

amounts (ArrayLike) – Continuous molecule amounts (the volume-invariant exchange currency; use UnitConverter.to_amounts() to get here from storage).
policy (Literal['nearest', 'floor', 'ceil', 'stochastic']) – See RoundingPolicy. "nearest" rounds half away from zero, matching the SSA engine’s entry rounding so an explicit hand-off and the implicit one agree. "stochastic" rounds up with probability equal to the fractional part (unbiased), and needs rng.
rng (Generator | None) – Required for policy="stochastic"; ignored otherwise.

Returns:

Integer-valued float64 counts (same dtype as the state vector, so it feeds straight back through UnitConverter.from_counts()).

Return type:

bngsim.round_half_up(value)[source]¶

Round a fractional molecule count to the nearest integer, ties away from zero.

bngsim’s repo-wide stochastic initial-count policy (GH #51): floor(x + 0.5) for x >= 0 and ceil(x - 0.5) for x < 0. The two branches differ only at negative half-integers, which never arise for molecule populations. The rule is idempotent on integers (round_half_up(n) == n), so warm-restart paths that re-round already-integer state are no-ops.

Examples

>>> [round_half_up(x) for x in (0.389, 0.10001, 5.7, 155.6747, 466.98)]
[0, 0, 6, 156, 467]
>>> [round_half_up(x) for x in (0.5, 1.5, 2.5)]  # ties away from zero, not bankers'
[1, 2, 3]

Parameters:: value (float) – Continuous molecule count, already in molecule-number space. Any volume / concentration factors must be applied before calling.
Returns:: Nearest integer count.
Return type:: int
Raises:: ValueError – If value is not finite.

class bngsim.ConservationLedger(weights=None, *, atol=1e-09, rtol=1e-09, name='total')[source]¶

Bases: object

Track a conserved moiety across exchange-boundary round-trips (GH #102).

Records the moiety total at a baseline state and on every subsequent state, so an operator-split loop can assert the shared moiety is preserved across get_state → orchestrator → set_state exchanges and the discrete rounding.

Parameters:

weights (ArrayLike | None) – Moiety weight vector (see moiety_total()). Default: total count.
atol (float) – Absolute / relative tolerance for check() and assert_conserved(). The drift bound is atol + rtol * |baseline|.
rtol (float) – Absolute / relative tolerance for check() and assert_conserved(). The drift bound is atol + rtol * |baseline|.
name (str) – Label used in error messages.

record(state)[source]¶

Record the moiety total of state; set the baseline on first call.

Returns the moiety total.

Parameters:: state (ArrayLike)
Return type:: float

check(state)[source]¶

Record state and return (within_tolerance, signed_drift).

Parameters:: state (ArrayLike)
Return type:: tuple[bool, float]

assert_conserved(state)[source]¶

Record state; raise ConservationError if drift exceeds tol.

Returns the signed drift from baseline.

Parameters:: state (ArrayLike)
Return type:: float

property baseline: float | None¶: The first recorded moiety total (None before any record).

property last: float | None¶: The most recently recorded moiety total.

property max_abs_drift: float¶: Largest absolute drift from baseline seen so far.

property n_records: int¶: Number of states recorded.

exception bngsim.ConservationError[source]¶

Bases: BngsimError

A conserved moiety drifted beyond tolerance across the exchange boundary.

bngsim.moiety_total(state, weights=None)[source]¶

Total of a conserved moiety: weights · state (or state.sum()).

A conserved moiety of a reaction network is a left-null-space vector w of the stoichiometry matrix; w · n is invariant under the reactions. The default weights=None is the all-ones vector — the total molecule count of a closed transfer network, the moiety an operator split must not leak.

Parameters:

state (ArrayLike)
weights (ArrayLike | None)

Return type:

class bngsim.Divider(*, method='binomial', rng=None)[source]¶

Bases: object

Partition molecule counts across daughter cells at division (GH #102).

Cell division is a pure count-space operation with no analogue elsewhere in bngsim: the molecules of each species are dealt out among the daughters, exactly conserving the parent total (sum(daughters) == parent for every partitioned species). It composes on top of UnitConverter — convert storage → counts, divide, then convert each daughter’s counts back to storage via UnitConverter.from_counts() and inject with set_state. Volume is the orchestrator’s to halve (see set_compartment_volume()).

Parameters:

method (Literal['binomial', 'multinomial', 'deterministic']) – "binomial" / "multinomial" (synonyms) deal each molecule to a uniformly-random daughter — the physically faithful stochastic split, exactly conserving. "deterministic" splits as evenly as possible (floor share + largest-remainder distribution of the leftover), also exactly conserving; reproducible without an RNG.
rng (Generator | None) – Required for the stochastic methods.

divide(counts, n_daughters=2, *, partition_mask=None)[source]¶

Partition integer counts into n_daughters daughter vectors.

Parameters:

counts (ArrayLike) – Non-negative, integer-valued molecule counts (e.g. from UnitConverter.to_counts()).
n_daughters (int) – Number of daughters (default 2).
partition_mask (ArrayLike | None) – True where a species is a partitioned molecule pool (split and conserved); False where it is shared environment (copied identically to every daughter, not split). Default: partition all.

Returns:

n_daughters integer-valued float64 count vectors. For every partitioned species the daughters sum exactly to the parent count.

Return type:

list[NDArray[double]]

bngsim.make_subset_model(model, *, keep_reactions=None, fixed_species=None, compute_conservation_laws=False)[source]¶

Build an operator-split subset as its own model (GH #102 Stage 1, #6).

bngsim runs a whole network with one method — there is no native reaction-subset integration. For a static operator split the continuous (ODE) subset is supplied as its own model: the same species namespace, only the subset’s reactions, and the other operator’s species marked fixed=True so the integrator holds them at the boundary values the orchestrator writes each step via set_state. The fixed mechanism already clamps a species in both the ODE (zeroed derivative) and SSA (skipped fire) backends — this helper just constructs the partitioned model; no engine change is involved.

The full species set is kept (so the two subsets share one addressing space for CouplingMap); only the reaction set is subset and the boundary species fixed. A species touched by no kept reaction is already inert, so fixed_species is needed only for boundary species whose value a kept rate law reads but should not evolve.

Parameters:

model (Model | ReactionKernel) – The full network to partition.
keep_reactions (Sequence[int] | None) – 0-based indices of the reactions this subset integrates. None keeps all reactions (a pure re-fix with no reaction split).
fixed_species (Sequence[str] | None) – Species names to additionally mark fixed (the other operator’s species). Species already fixed in the source stay fixed.
compute_conservation_laws (bool) – Forwarded to the builder (default False — the dense O(n³) detector is unused by stepping and costly at scale; see GH #102 MVP).

Returns:

A freshly built subset model at the source’s initial concentrations (taken from a reset clone, so it is immune to a prior simulation having left the source mid-trajectory). Inject a custom starting state with set_state on the subset afterwards.

Return type:

Raises:

NotImplementedError – If the source uses features this reconstruction cannot reproduce faithfully (events, table functions, discontinuity triggers, amount_valued species, or per_species_volume_scaling reactions). The supported class — mass action, functional, and Michaelis–Menten rate laws, parameters, observables, functions — covers the issue’s ~100K first-order target.

bngsim.get_compartment_volume(model, name)[source]¶

Read a compartment volume (a model parameter) by name.

There is no Compartment object in bngsim — a compartment volume is a plain parameter, so this is Model.get_param(), named for the coupling use it serves: feeding the current volume into UnitConverter.to_concentrations() / UnitConverter.from_concentrations() as the live-volume override.

Parameters:

model (Model | ReactionKernel)
name (str)

Return type:

bngsim.set_compartment_volume(model, name, volume)[source]¶

Couple a framework’s volume growth into bngsim’s compartment volume.

Sets the compartment-volume parameter. On the ODE subset this flows through bngsim’s variable-volume machinery natively (the integrator dilutes by the live compartment symbol; GH #74/#85), so growing or halving (at division) a compartment is just this call. On the SSA subset the baked per-reaction ssa_volume_factor does not track the new value — the exchange layer compensates by reading concentrations against the live volume (the UnitConverter override), but the stochastic propensities using the live volume are Stage 2.

Parameters:

model (Model | ReactionKernel)
name (str)
volume (float)

Return type:

class bngsim.Result(core=None, *, custom_attrs=None, _time=None, _species=None, _observables=None, _expressions=None, _species_names=None, _observable_names=None, _expression_names=None, _solver_stats=None, _species_volume_factors=None, _seed=None, _ic_sensitivity_seed=None, _sensitivities=None, _sensitivity_params=None, _sensitivities_ic=None, _sensitivity_ic_species=None, _observable_sensitivities=None, _expression_sensitivities=None, _observable_sensitivities_ic=None, _expression_sensitivities_ic=None, _reaction_firing_counts=None, _reaction_propensity_integrals=None, _reaction_labels=None, _ssa_diagnostics=None)[source]¶

Bases: object

Simulation result container.

Variables:

time (ndarray, shape (n_times,)) – Time points.
species (ndarray, shape (n_times, n_species)) – Species concentrations at each time point.
observables (ndarray, shape (n_times, n_observables)) – Observable values at each time point.
solver_stats (dict) – Solver diagnostics (n_steps, n_rhs_evals, etc.).

Parameters:

core (ResultCore | None)
custom_attrs (dict[str, Any] | None)
_time (NDArray | None)
_species (NDArray | None)
_observables (NDArray | None)
_expressions (NDArray | None)
_species_names (list[str] | None)
_observable_names (list[str] | None)
_expression_names (list[str] | None)
_solver_stats (dict[str, int] | None)
_species_volume_factors (list[float] | list[list[float]] | None)
_seed (int | None)
_ic_sensitivity_seed (dict[str, dict[str, float]] | None)
_sensitivities (NDArray | None)
_sensitivity_params (list[str] | None)
_sensitivities_ic (NDArray | None)
_sensitivity_ic_species (list[str] | None)
_observable_sensitivities (NDArray | None)
_expression_sensitivities (NDArray | None)
_observable_sensitivities_ic (NDArray | None)
_expression_sensitivities_ic (NDArray | None)
_reaction_firing_counts (NDArray | None)
_reaction_propensity_integrals (NDArray | None)
_reaction_labels (list[str] | None)
_ssa_diagnostics (dict[str, Any] | None)

Examples

>>> result = sim.run(t_span=(0, 100), n_points=101)
>>> result.time.shape
(101,)
>>> result.observables.shape
(101, 5)
>>> result.observables["Atot"]  # named access
array([100., 99.5, ...])

Notes

A Result holds a handle into bngsim’s compiled extension, so it cannot be pickled or deep-copied; save() and Result.load(), or one of the conversions (as_roadrunner(), to_xarray()), carry it across a process boundary instead. copy.copy(result) does work, and the copy shares that handle and the result arrays with the original. That is left working, where the same sharing makes copy.copy refuse on a Model or a Simulator (issue #643), because a Result has no setter and no mutating method: no call you can make on the copy writes into the original. Writing into its arrays in place does — copy.copy(result).species[0] = ... changes result — but so does result.species[0] = ..., because every access returns the same live view rather than a fresh array. Take .copy() of a block you mean to modify, whichever name you reach it through.

Parameters:

core (Optional[ResultCore])
custom_attrs (dict[str, Any] | None)
_time (NDArray | None)
_species (NDArray | None)
_observables (NDArray | None)
_expressions (NDArray | None)
_species_names (list[str] | None)
_observable_names (list[str] | None)
_expression_names (list[str] | None)
_solver_stats (dict[str, int] | None)
_species_volume_factors (list[float] | list[list[float]] | None)
_seed (int | None)
_ic_sensitivity_seed (dict[str, dict[str, float]] | None)
_sensitivities (NDArray | None)
_sensitivity_params (list[str] | None)
_sensitivities_ic (NDArray | None)
_sensitivity_ic_species (list[str] | None)
_observable_sensitivities (NDArray | None)
_expression_sensitivities (NDArray | None)
_observable_sensitivities_ic (NDArray | None)
_expression_sensitivities_ic (NDArray | None)
_reaction_firing_counts (NDArray | None)
_reaction_propensity_integrals (NDArray | None)
_reaction_labels (list[str] | None)
_ssa_diagnostics (dict[str, Any] | None)

custom_attrs: dict[str, Any]¶

property seed: int | None¶

The integer RNG seed actually used for this simulation.

None for deterministic (ODE) results and for results loaded from HDF5 files predating seed exposure. For stochastic methods (SSA, PSA, NFsim, RuleMonkey), this is the integer that was passed down to the backend — equal to the seed= keyword when the caller supplied one, or the freshly drawn integer when the caller passed seed=None (or omitted it). A seed is from 0 to 2**64 - 1, and under NFsim from 0 to 2**32 - 1, the width its generators read, so the seed passed is the seed used (issue #1515).

For squeezed batch results, this is the single seed if every underlying sim used the same one, otherwise None; the per-sim seeds remain accessible on the individual Result objects from Simulator.run_batch(..., squeeze=False).

property stop_condition: str | None¶

The stop condition this run met, by its label, or None.

A condition registered with Simulator.add_stop_condition() cuts the result at the first sample it holds at: Simulator.run() raises StopConditionMet with that result, and Simulator.run_batch() and Simulator.run_replicates() return it in their list, with this set, beside the rows that ran to the end (issue #1369). A condition registered without a label reads as its expression, or a callable as its repr(). save() and to_xarray() carry it.

A cut result holds the reported blocks only, so state refuses on it, as it does on the result Simulator.run() raises with.

property ic_sensitivity_seed: dict[str, dict[str, float]] | None¶

∂x(0)/∂θ — the initial-condition seeding this run’s gradients carry.

The per-run record of the matrix Model.effective_ic_sensitivity() reports. Prefer that one to build gradient routing: it answers at setup, from model structure alone, with no simulation. Reach for this to see what a particular run used — the only correct answer for a batch or scan over a nonlinear derived initial condition (Rtot = R0*scale), whose coefficients move from point to point, and the one that reports None for a carried-seed phase (below).

{species_name: {param_name: ∂x_k(0)/∂θ}}, the forward-sensitivity seed CVODES was actually initialized with for this simulation: the parameter-graph derivative of each species’ initial condition (issues #43/#147), after the issue #113 retirement of a species moved off its declared initial condition and the issue #111 Model.declare_ic_sensitivity() overlay.

Read this to compose your own initial-condition chain rule without double-counting. output_sensitivities() with axis="parameter" is a total derivative:

d_param[θ] = (right-hand-side path) + Σ_k (∂x_k(0)/∂θ)·d_ic[x_k]

so for any θ appearing here, the seeding term is already inside d_param[θ] and adding an ic-axis term for it counts it twice. A θ absent from this mapping has an RHS-only d_param column, and an initial-condition term your own code supplies is the missing piece (or hand it to Model.declare_ic_sensitivity() and let the parameter axis carry it instead).

Only parameters in sensitivity_params appear — a row the engine computed for a parameter you did not request has no column to describe. Coefficients are evaluated at this run’s parameter values, so a batch or scan over a nonlinear derived initial condition (Rtot = R0*scale) reports a different matrix per point.

Returns:: {} means recorded, and nothing seeds — every d_param column is right-hand-side only. None means not recorded: no parameter sensitivities were requested; the result came from a non-ODE method or an HDF5 file predating this field; a squeezed batch whose rows disagreed (the per-row Result objects from run_batch(..., squeeze=False) still carry theirs); or the run used carry_sensitivities=True, where the seed is the prior phase’s dx/dθ and these rows are discarded by the engine (GH #210/#81 — and no double-count can arise there, because an ic axis across a carry boundary is refused outright). The two are deliberately distinct — an empty mapping is a positive statement.
Return type:: dict or None

Examples

>>> sim = Simulator(model, method="ode", sensitivity_params=["R0", "kf"])
>>> r = sim.run(t_span=(0, 10), n_points=11)
>>> r.ic_sensitivity_seed
{'R(r)': {'R0': 1.0}}

property time: NDArray[float64]¶: Time points, shape (n_times,).

property species: NDArray[float64]¶: Species concentrations, shape (n_times, n_species).

property state: NDArray[float64]¶

The full integrator state trajectory, shape (n_times, n_state).

Every core state entry, exactly as the integrator held it and in bngsim.Model.species_names order, so a row is a state bngsim.Model.rhs() and bngsim.Model.jacobian() accept as is (issue #508). It differs from species on two kinds of model. An SBML model whose event assigns a parameter or a compartment promotes the target to a state entry that species and species_names omit (GH #71), so this has more columns. A model whose reported species column is remapped after the solve — an assignment-rule target reported at the rule’s live value, or a species in an event-resized or rate-ruled compartment reported at amount / V_live(t) (GH #85, #131) — reports the remapped value in species and the stored one here. Elsewhere the two are identical.

Only a result returned by a solve carries it: a result loaded from disk or stacked from a batch holds the reported species block alone.

Raises:: ValueError – If this result was not produced directly by a solve.

property state_names: list[str]¶

the producing model’s bngsim.Model.species_names, promoted entries included.

Type:: Names of the state columns

property observables: _ObservableAccessor¶

Observable values with named access.

Can be used as: - result.observables → full array (n_times, n_obs) - result.observables["name"] → single column (n_times,)

property expressions: _ObservableAccessor¶: Expression (function) values with named access.

property n_times: int¶: Number of time points.

property n_species: int¶: Number of species.

property n_observables: int¶: Number of observables.

property n_expressions: int¶: Number of expressions.

property has_simulation_data: bool¶

True iff the simulator produced at least one numeric column.

Returns False for runs that complete without any species, observable, or expression columns — e.g. a BNGL file with a simulate action but an empty model body. Used by batch harnesses to distinguish “ran and produced data” from “ran and produced nothing meaningful” without an exception path.

property species_names: list[str]¶: Species names.

property observable_names: list[str]¶: Observable names.

property expression_names: list[str]¶

Expression (function) names (bare, in-memory column keys).

Auto-generated _rateLawN rate-law intermediates are filtered out; recover them via raw_expression_names.

property gdat_expression_names: list[str]¶

Function-column labels for .gdat/.scan headers.

bngsim emits bare function headers (no () suffix, issue #58), so this is identical to expression_names. It is retained for consumers that assemble a BNG-native file header (e.g. a consumer-built .scan) and want an intent-revealing name.

property raw_expression_names: list[str]¶

Unfiltered function names, including internal _rateLawN.

On a freshly-simulated result these include the auto-generated _rateLawN rate-law intermediates that expression_names filters out (the #58 recoverability fix). On a result loaded from HDF5 or assembled from raw arrays, only the filtered (bare) set was persisted, so this returns the same columns as expression_names.

property raw_expressions: _ObservableAccessor¶

Unfiltered function values, including internal _rateLawN.

Columns correspond to raw_expression_names. See that property for the freshly-simulated vs. loaded distinction.

property raw_n_expressions: int¶: Number of unfiltered function columns (includes _rateLawN).

resolve_outputs(selectors)[source]¶

Resolve typed output selectors to structured metadata.

A selector names one output column of this result. Fitting frontends compare against named outputs (species, observables, global functions) rather than raw column indices; this is the unified, typed lookup that maps a selector string onto the column it refers to. No sensitivity or value computation happens here — it is pure name resolution.

Parameters:

selectors (str | Iterable[str]) –

One selector or a sequence of them. A single string is treated as a one-element list (so resolve_outputs("observable:Atot") and resolve_outputs(["observable:Atot"]) are equivalent).

Accepted forms:

"species:<name>" — a species column.
"observable:<name>" — an observable column.
"expression:<name>" — a global-function / expression column. A trailing () is stripped ("expression:foo()" → "expression:foo"), matching the BNG .gdat/.scan header convention while in-memory column keys stay bare (issue #58).
Aliases: "state:" → "species:", "function:" → "expression:".
A bare name (no prefix) resolves only if it is unique across species, observables, and expressions. A bare "foo()" that does not match any column verbatim is retried as the expression foo (the function-call convention).

Returns:

One dict per input selector, in input order, each with keys:

"selector" — the canonical typed selector ("<kind>:<name>"), suitable as a stable downstream key.
"kind" — "species", "observable" or "expression".
"name" — the bare in-memory column name (the key used by the named accessors, e.g. result.observables[name]).
"index" — column index within that kind’s array.
"column_label" — the label as written to a .gdat/.cdat header. bngsim emits bare headers (issue #58), so this is currently identical to "name".

Return type:

list[dict[str, Any]]

Raises:

ValueError – If a selector uses an unknown kind prefix, is empty, is unresolved, or is a bare name that matches more than one column. Ambiguity errors list every matching typed selector; unresolved errors list the available candidates.
TypeError – If a selector is not a string.

Examples

>>> result.resolve_outputs("observable:Atot")
[{'selector': 'observable:Atot', 'kind': 'observable',
  'name': 'Atot', 'index': 0, 'column_label': 'Atot'}]
>>> [m["selector"] for m in result.resolve_outputs(
...     ["species:A()", "function:scaled()", "Btot"])]
['species:A()', 'expression:scaled', 'observable:Btot']

outputs(selectors)[source]¶

Return the value columns named by selectors.

Thin value accessor on top of resolve_outputs(): resolves each selector and stacks the named columns into a single array.

Parameters:: selectors (str | Iterable[str]) – Selectors accepted by resolve_outputs().
Returns:: Shape (n_times, n_outputs) for a single-simulation result, with one column per selector in input order. For a stacked batch result (3-D arrays) the shape is (n_sims, n_times, n_outputs). An empty selector list yields a (n_times, 0) array.
Return type:: NDArray[double]
Raises:: ValueError, TypeError – Propagated from resolve_outputs().

output_sensitivities(selectors, *, axis='parameter')[source]¶

Return d(named output)/dθ for each selector, stacked.

Selector-addressed companion to outputs(): resolves each selector and stacks the matching sensitivity slice — the chain-rule derivative of that output column with respect to the sensitivity parameters (default) or the differentiated initial conditions.

Observable sensitivities are computed at simulation time from the CVODES species sensitivities via the linear chain rule d obs_j/dθ = Σ_i c_ji·dx_i/dθ (GH #197); expression: selectors carry the codegen output-sensitivity chain rule for global functions (GH #198). species: selectors read the species sensitivities directly — except for an SBML AssignmentRule-target species, whose reported value is the rule’s live value (its state slot is frozen), so its output sensitivity follows the assignment expression: the sensitivity of the rule’s observable (linear-on-species rule) or function (everything else), GH #205. Since GH #221 the sensitivities tensor carries that same redirected row, so this method and the tensor agree for those species; the redirect is kept here because it raises by name where the tensor can only hold NaN (an expression whose output sensitivity a construct declines). A rule target in a compartment whose size changes is read from the tensor, whose row has the size’s share in it (issue #787), and is refused by name where that row is NaN.

The two axes are not independent (GH #155)

axis="parameter" is the total derivative dy(t)/dθ: it carries every path by which θ reaches the trajectory, the right-hand side and the initial-condition seeding ∂x(0)/∂θ (issues #43/#147). axis="ic" is ∂y(t)/∂x_k(0) with the initial value held as an independent variable. So for a parameter that seeds an initial condition the two overlap:

d_param[θ] = (right-hand-side path) + Σ_k (∂x_k(0)/∂θ)·d_ic[x_k]

Do not sum the axes. A consumer that routes a fitted parameter to several columns and adds them must add an ic term only for the part of ∂x(0)/∂θ this engine does not already carry — which ic_sensitivity_seed reports exactly, per run. The common case that makes the overlap easy to miss: a parameter appearing only in the initial condition, with coefficient 1, gives a parameter column bit-identical to the seeded species’ ic column (same seed vector, same variational equation), so summing reports it at exactly 2×.

On SteadyStateResult the question does not arise: that parameter axis is the implicit-function derivative of the algebraic system, there is no seeding term in it, and axis="ic" raises.

Parameters:

selectors (str | Iterable[str]) – Selectors accepted by resolve_outputs().
axis (str) – Which sensitivity axis to return. "parameter" (default) gives d output/dp over sensitivity_params — the total derivative, seeding included, see above; "ic" gives d output/dY(0) over sensitivity_ic_species.

Returns:

Shape (n_times, n_outputs, n_axis) for a single-simulation result ((n_sims, n_times, n_outputs, n_axis) for a stacked batch), one slice per selector in input order. n_axis is the number of sensitivity parameters (or IC species). An empty selector list yields a (n_times, 0, n_axis) array.

Return type:

Raises:

ValueError – If axis is invalid; if the requested sensitivity axis was not computed for this result; or if a selector names a kind whose sensitivities are unavailable (e.g. an expression: selector before GH #198).
TypeError – Propagated from resolve_outputs().

Examples

>>> sim = Simulator(model, method="ode", sensitivity_params=["k1"])
>>> result = sim.run(t_span=(0, 10), n_points=11)
>>> result.output_sensitivities("observable:Atot").shape
(11, 1, 1)

fisher_information(sigma=1.0, *, outputs=None, axis='parameter')[source]¶

Compute the Fisher Information Matrix from sensitivity data.

The FIM quantifies how much information the simulated output trajectories carry about each parameter. It is a pure function of the model + output-sensitivity tensor — no measurement data, no residuals, no objective — so it measures the model’s practical identifiability (sloppiness), independent of any fit. It is computed as:

\[\text{FIM} = \sum_t \left(\frac{\partial Y}{\partial \theta}\right)^T \Sigma^{-1} \left(\frac{\partial Y}{\partial \theta}\right)\]

where $\partial Y / \partial \theta$ is the (n_outputs, n_params) sensitivity matrix at time t, and $\Sigma$ is the diagonal output-noise covariance.

By default the FIM is built over species (back-compatible with the original species-only method). Pass outputs to build it over named outputs — any mix of species:/observable:/expression: selectors (GH #202) — using the output-sensitivity tensor (GH #197/#198). For the richer eigenvalue / identifiability readout, see identifiability().

Unlike the gradient methods, the FIM weights every output it is built over, so there is no unweighted row to drop: a model with a refused assignment-rule row (ar_sensitivity_refused, GH #221) returns a NaN matrix by default. Name the outputs you want in outputs to build it over the rows that are known.

Parameters:

sigma (float | NDArray[double]) –
Output-noise standard deviation used only to scale the FIM (it is not measurement data). Defaults to 1.0 (unscaled).
- scalar — homogeneous: every output shares the same σ.
- 1-D array — per-output σ. Length must match the number of outputs: n_species when outputs is None, otherwise the number of selectors.
outputs (str | Iterable[str] | None) – Output selectors (see resolve_outputs()) to build the FIM over. None (default) uses every species — the original species-only behaviour.
axis (str) – Sensitivity axis: "parameter" (default) builds the FIM over sensitivity_params; "ic" over sensitivity_ic_species.

Returns:

Symmetric positive-semidefinite FIM over n parameters (or IC species, for axis="ic").

Return type:

Raises:

ValueError – If the requested sensitivity axis was not computed (run with sensitivity_params / sensitivity_ic); if sigma’s shape is wrong; if axis is invalid; or if this is a batch result (the FIM is defined per single simulation — iterate replicates).

Notes

The FIM is the Cramér–Rao lower bound on the inverse of the parameter covariance: $\text{Cov}(\hat\theta) \geq \text{FIM}^{-1}$. Large diagonal entries indicate identifiable parameters; small entries indicate practical non-identifiability. np.linalg.cond(fim) (or IdentifiabilityReport.condition_number) gauges the overall identifiability of the parameter set.

Examples

>>> sim = Simulator(model, method="ode",
...                 sensitivity_params=["k1", "k2"])
>>> result = sim.run(t_span=(0, 100), n_points=101)
>>> fim = result.fisher_information(sigma=0.1)           # species
>>> fim.shape
(2, 2)
>>> result.fisher_information(outputs=["observable:Atot"]).shape
(2, 2)

identifiability(sigma=1.0, *, outputs=None, axis='parameter', rtol=None)[source]¶

Model-identifiability readout from the Fisher Information Matrix.

Builds the FIM (see fisher_information()) and returns its eigen-decomposition with practical-identifiability flags: which parameter combinations are well-constrained by the model and which are “sloppy” (near-zero eigenvalues / practically non-identifiable). This is sensitivity analysis of the model — no data, no residuals, no objective.

Parameters:

sigma (float | NDArray[double]) – Output-noise σ; see fisher_information(). Defaults to 1.0.
outputs (str | Iterable[str] | None) – Output selectors to build the FIM over; None (default) uses every species.
axis (str) – Sensitivity axis; see fisher_information().
rtol (float | None) – Relative eigenvalue cutoff for flagging non-identifiable / sloppy directions: a direction is flagged when its eigenvalue is below rtol * lambda_max. Defaults to n * eps (NumPy’s numerical-rank tolerance), which flags only numerically singular directions; pass a larger value (e.g. 1e-6) to surface sloppy-but-nonzero directions.

Returns:

Eigenvalues/eigenvectors, numerical rank, condition number, per-direction identifiability flags, and the Cramér–Rao bound FIM⁻¹ (all-NaN + a warning when the FIM is rank-deficient).

Return type:

IdentifiabilityReport

Examples

>>> report = result.identifiability(sigma=0.1)
>>> report.rank, report.condition_number
(2, 3.4)
>>> report.non_identifiable_directions
[]

gradient(loss_fn)[source]¶

Compute parameter gradient from sensitivity data and a loss function.

Given forward sensitivities $\partial Y / \partial p$ and a user-supplied loss function that returns the per-species loss gradient $\partial L / \partial Y$ at each time point, computes the parameter gradient:

\[\nabla_p L = \sum_t \left(\frac{\partial Y}{\partial p}\right)^T \frac{\partial L}{\partial Y}(t)\]

This enables gradient-based optimization via scipy.optimize.minimize(method='L-BFGS-B').

Parameters:

loss_fn (Callable[[NDArray[double], NDArray[double]], NDArray[double]]) –

A function with signature loss_fn(species, time) -> dL_dY where:

species — species, the (n_times, n_species) array of species values
time — (n_times,) array of time points
returns dL_dY — (n_times, n_species) array of $\partial L / \partial Y_{i,t}$.

dL_dY is contracted against sensitivities, whose species axis is the same one species is indexed on (GH #202), so a dL_dY shaped like the species argument always fits.

Common example (sum-of-squares vs data):

def loss_fn(species, time):
    return 2 * (species - data)  # dL/dY for SSE

Returns:

Parameter gradient $\nabla_p L$, one entry per column of sensitivities, ordered like sensitivity_params. That list — not Model.param_names — is what the returned vector is aligned with, and it is what an optimizer’s parameter vector must be built from (issue #203; see Notes).

Return type:

Raises:

ValueError – If sensitivity data is not available.
TypeError – If loss_fn is not callable.

Examples

>>> # Sum-of-squares loss vs observed data
>>> data = np.random.randn(101, 2)  # observed data
>>> result = sim.compute_all_sensitivities(
...     t_span=(0, 10), n_points=101
... )
>>> grad = result.gradient(
...     lambda species, time: 2 * (species - data)
... )
>>> grad.shape
(n_params,)

>>> # Use with scipy L-BFGS-B. The fitted vector is the *default*
>>> # sensitivity column set, so `grad` lines up with `p_vec`.
>>> from scipy.optimize import minimize
>>> names = model.primary_param_names
>>> def objective(p_vec):
...     model.set_params(dict(zip(names, p_vec)))
...     model.reset()
...     result = sim.compute_all_sensitivities(...)
...     assert result.sensitivity_params == names
...     loss = np.sum((result.species - data)**2)
...     grad = result.gradient(
...         lambda sp, t: 2 * (sp - data)
...     )
...     return loss, grad
>>> x0 = [model.get_param(n) for n in names]
>>> minimize(objective, x0, method='L-BFGS-B', jac=True)

Notes

The columns must be independent coordinates, and under the ``compute_all_sensitivities()`` default they are (issue #203). This method contracts the whole parameter axis into one (n_params,) vector, so it is only a gradient if no two columns describe the same physical effect. A derived parameter’s column does describe an effect the primaries underneath it already carry — _rateLaw1 = chi*kon reaches the trajectory only through chi and kon, whose own columns are total derivatives through it — so summing all three counts that effect twice. compute_all_sensitivities(params=None) drops the derived columns for this reason and warns; a sensitivity_params= list you wrote yourself is your own to keep independent, and if it names both a derived parameter and a primary underneath it, do not hand this vector to an optimizer over both. (Same shape as issue #155 one level down: the parameter axis of an initial-condition sensitivity is a total derivative and must not be summed with the ic axis.)

The gradient computation is O(n_times × n_species × n_params) — a single matrix multiply per time point. With parallel compute_all_sensitivities(), the total cost of computing both the loss value and its gradient is ~1.2× a plain ODE solve (with sufficient cores).

sse_gradient(data, *, target='species', species_indices=None, observable_indices=None)[source]¶

Compute SSE loss and parameter gradient in one call.

Sum-of-squared-errors: $L = \sum_{t,i} (Y_{t,i} - D_{t,i})^2$.

This is the most common objective in parameter estimation. The analytical derivative is $\partial L / \partial Y = 2(Y - D)$, so no user-supplied loss_fn is needed.

Parameters:

data (NDArray[double]) – Observed data, of the shape of what target names (or of the columns species_indices / observable_indices select).
target (str) – What the data is compared with: species, through sensitivities, or observables, through sensitivities_observables (issue #1363). Named, not read from the data’s shape, which is the same for both when a model has as many observables as species.
species_indices (list[int] | None) – If given, only these species columns are used for the loss (target="species"). Sensitivities for all parameters are still included.
observable_indices (list[int] | None) – The same, over observables (target="observables").

Returns:

(loss, gradient) where loss is the scalar SSE and gradient is shape (n_params,).

Return type:

tuple[float, NDArray[double]]

Raises:

ValueError – If sensitivity data is missing, target is neither, or shapes don’t match.

Examples

>>> result = sim.compute_all_sensitivities(...)
>>> loss, grad = result.sse_gradient(observed_data)
>>> loss, grad = result.sse_gradient(observed_obs, target="observables")
>>> # Use with scipy L-BFGS-B:
>>> minimize(objective, x0, method='L-BFGS-B', jac=True)

chi2_gradient(data, sigma, *, target='species', species_indices=None, observable_indices=None)[source]¶

Compute chi-squared loss and parameter gradient.

$\chi^2 = \sum_{t,i} \left(\frac{Y_{t,i} - D_{t,i}}{\sigma_i}\right)^2$

Equivalent to weighted SSE with weights $1/\sigma_i^2$. The derivative is $\partial L / \partial Y_{t,i} = 2(Y - D)/\sigma_i^2$.

Parameters:

data (NDArray[double]) – Observed data, of the shape of what target names.
sigma (float | NDArray[double]) –
Measurement noise standard deviation.
- scalar — same σ for all columns.
- 1-D array, shape ``(n_columns,)`` — per-column σ.
- 2-D array, shape ``(n_times, n_columns)`` — per-point σ.
target (str) – What the data is compared with, as in sse_gradient().
species_indices (list[int] | None) – Subset of species columns to include.
observable_indices (list[int] | None) – Subset of observable columns to include.

Returns:

(chi2, gradient) where chi2 is the scalar value and gradient is shape (n_params,).

Return type:

tuple[float, NDArray[double]]

Examples

>>> loss, grad = result.chi2_gradient(data, sigma=0.1)
>>> loss, grad = result.chi2_gradient(data, sigma=per_species_sigma)

neg_log_likelihood_gradient(data, sigma, *, target='species', species_indices=None, observable_indices=None)[source]¶

Compute negative Gaussian log-likelihood and gradient.

$-\ln \mathcal{L} = \frac{1}{2} \sum_{t,i} \left[\left(\frac{Y_{t,i} - D_{t,i}}{\sigma_i}\right)^2 + \ln(2\pi\sigma_i^2)\right]$

The constant term $\ln(2\pi\sigma^2)$ doesn’t affect the gradient but is included in the loss value for correct log-likelihood.

Parameters:

data (NDArray[double]) – Observed data, of the shape of what target names.
sigma (float | NDArray[double]) – Measurement noise standard deviation (scalar, per-column, or per-point).
target (str) – What the data is compared with, as in sse_gradient().
species_indices (list[int] | None) – Subset of species columns.
observable_indices (list[int] | None) – Subset of observable columns.

Returns:

(nll, gradient) — negative log-likelihood and gradient.

Return type:

tuple[float, NDArray[double]]

Examples

>>> nll, grad = result.neg_log_likelihood_gradient(data, sigma=0.1)

property sensitivities: NDArray[float64]¶

Forward sensitivity data dY/dp.

Shape (n_times, n_species, n_params) when sensitivities are available, or empty (0, 0, 0) otherwise.

The species axis is the same axis as species / species_names: row i is d species_names[i]/dp, so a dL/dY built from species contracts against this tensor row-for-row (gradient(), sse_gradient(), …). On a model where an SBML parameter or compartment is an event-assignment target, bngsim promotes that symbol to full integrator state but does not report it as a trajectory column (GH #71); its row is likewise absent here, and its derivative is output_sensitivities("observable:<name>") — the loader gives every promoted symbol a same-named observable (GH #202).

Row i is the derivative of the value in column *i*, which for an SBML AssignmentRule-target species is the rule’s live value, not the frozen state slot the integrator carries (GH #221): those rows are the chain rule through the assignment, the same quantity output_sensitivities("species:<name>") returns. Where that chain rule is unavailable the row is NaN — never a structural zero, which no consumer could tell from a measured one — and the run warns naming the species (ar_sensitivity_refused).

Example

>>> sim = Simulator(model, method="ode",
...                 sensitivity_params=["k"])
>>> result = sim.run(t_span=(0, 10), n_points=11)
>>> result.sensitivities.shape
(11, 2, 1)
>>> result.sensitivities[-1, 0, 0]  # dA/dk at t=10

property has_sensitivities: bool¶: Whether sensitivity data is available.

property sensitivity_params: list[str]¶: Parameter names for sensitivity analysis.

property sensitivities_ic: NDArray[float64]¶

Forward IC sensitivity data dY/dY(0).

Shape (n_times, n_species, n_ic_species) when IC sensitivities are available, or empty (0, 0, 0) otherwise. Column k is ∂Y(t)/∂Y_k(0) where the ordering of k matches sensitivity_ic_species. The row axis is species_names, exactly as in sensitivities (GH #202).

property has_sensitivities_ic: bool¶: Whether IC sensitivity data is available.

property sensitivity_ic_species: list[str]¶: Species names whose ICs were differentiated against.

property reaction_firing_counts: NDArray[float64]¶

Cumulative firings of each reaction at each output time, N_r(t).

Shape (n_times, n_reactions) for a single run and (n_sims, n_times, n_reactions) for a stacked batch, or empty (0, 0) unless the run was made with Simulator(model, method="ssa", reaction_stats=True). Column r is labelled by reaction_labels. Integer-valued, held as float64 so it subtracts directly from reaction_propensity_integrals: N_r(t) − ∫₀ᵗ a_r ds is the compensated counting process of reaction r, whose mean over replicates is zero, and which bngsim.girsanov turns into likelihood-ratio scores for the rate constants.

property reaction_propensity_integrals: NDArray[float64]¶

Integrated propensity of each reaction at each output time, ∫₀ᵗ |a_r| ds.

Same shape and labelling as reaction_firing_counts. The integrand is the propensity magnitude the sampler used to draw the reaction — a reaction whose rate law goes negative fires in reverse at |a_r| (GH #110) — so the block describes the process as sampled.

property reaction_labels: list[str]¶

Labels for the reaction axis of the per-reaction statistics.

"R<index> (A + B -> C)": the .net file’s 1-based reaction index, then the reaction written in species names, 0 standing for no reactant or no product. Empty unless the statistics were recorded.

property has_reaction_stats: bool¶: Whether per-reaction firing counts and propensity integrals were recorded.

property sensitivities_species: NDArray[float64]¶

Alias for sensitivities (species forward sensitivities).

Provided so the four output-sensitivity blocks have a parallel, self-describing name on the same object: sensitivities_species, sensitivities_observables, sensitivities_expressions (+ the _ic variants). Identical array to sensitivities.

Every row is the derivative of the value the matching species column reports. For an SBML AssignmentRule-target species that is the rule’s live value: GH #205 routed output_sensitivities("species:…") through the rule, and GH #221 put the same row in this tensor, because the raw integrated yS of the frozen slot it used to hold is identically zero — not a low-level escape hatch, just a silent zero that gradient() read as a flat objective.

property ar_sensitivity_refused: frozenset[str]¶

assignment-rule targets, and amounts in a compartment whose size’s derivative is not to be had.

An AR-target species’ sensitivity row is the chain rule through its assignment (GH #221). Where that is unavailable the row is NaN rather than left at the structural zero no consumer could distinguish from a measurement, and the run warns naming these species. Two causes: codegen declined the rule’s own output sensitivity (GH #198 refuses an if() conditional or a table function rather than guess — the common one), or the target is an amount in a compartment whose size changes and codegen declined the output sensitivity of the size’s own rule (issue #787).

Also listed: an amount-valued species (hasOnlySubstanceUnits) that is no rule’s target, in a compartment an assignment rule resizes, where the rule’s own output sensitivity is not available. Its row is the quotient rule through the compartment’s size (issue #742), and without that derivative it is NaN. output_sensitivities() returns that row, where for an assignment-rule target it raises.

Membership means some entry of the row is NaN, not necessarily all of them. gradient() and the built-in objective gradients drop these rows wherever dL/dY is exactly zero, so a fit that never scores these species still gets a number; an entry the loss does weight stays NaN. output_sensitivities() raises for an assignment-rule target among them, with the specific reason.

Empty on every model that resolves: 590 of the 639 AR rows reachable across the 215 rr_parity corpus models that run under sensitivities.

Type:: Species whose sensitivity row carries NaN

property sensitivities_observables: NDArray[float64]¶

Observable parameter sensitivities d observable / dp.

Shape (n_times, n_observables, n_params) when computed, or empty (0, 0, 0) otherwise. The row axis matches observable_names (the loader’s internal network-rewrite observables are dropped here as well, GH #202) and the parameter axis matches sensitivity_params.

property has_sensitivities_observables: bool¶: Whether observable parameter sensitivities are available.

property sensitivities_expressions: NDArray[float64]¶

Expression (function) parameter sensitivities d expression / dp.

Shape (n_times, n_expressions, n_params) when computed, or empty (0, 0, 0) otherwise. The parameter axis matches sensitivity_params.

property has_sensitivities_expressions: bool¶: Whether expression parameter sensitivities are available.

property sensitivities_observables_ic: NDArray[float64]¶

Observable IC sensitivities d observable / dY(0).

Shape (n_times, n_observables, n_ic_species) when computed, or empty (0, 0, 0) otherwise. The row axis matches observable_names (GH #202) and the IC axis matches sensitivity_ic_species.

property has_sensitivities_observables_ic: bool¶: Whether observable IC sensitivities are available.

property sensitivities_expressions_ic: NDArray[float64]¶

Expression (function) IC sensitivities d expression / dY(0).

Shape (n_times, n_expressions, n_ic_species) when computed, or empty (0, 0, 0) otherwise. The IC axis matches sensitivity_ic_species.

property has_sensitivities_expressions_ic: bool¶: Whether expression IC sensitivities are available.

property solver_stats: dict[str, int]¶

Solver diagnostics dict.

Keys: n_steps, n_rhs_evals, n_jac_evals, n_jac_rhs_evals, n_jac_formations, n_err_test_fails, n_nonlin_iters, n_nonlin_conv_fails, n_sens_err_test_fails, n_sens_nonlin_conv_fails.

n_rhs_evals is CVODE’s count of the integration’s own evaluations. n_jac_rhs_evals counts those spent forming Jacobians by difference quotients (issue #1333); it is 0 where every Jacobian is closed-form. n_jac_evals counts the linear solver’s setups, which can reuse a Jacobian; n_jac_formations counts the Jacobians the solver asked for, one that a failed right-hand-side call cut short included.

property ssa_diagnostics: dict[str, Any]¶

SSA boundary diagnostics dict (GH #110).

The exact SSA evaluates rate laws literally: it does not floor species at zero, and a reaction whose rate law evaluates negative fires in reverse (propensity |rate|) so the SSA mean tracks the ODE. Both events are surfaced here instead of being silent.

Keys:

n_negative_crossings (int): times a species count crossed from >= 0 to < 0 (no floor applied).
first_negative_species (str): name of the first species to go negative; "" if none.
n_reverse_fires (int): reactions fired in reverse because their rate law was negative.
first_reverse_reaction (str): label of the first reaction reversed; "" if none.
n_rounded_populations (int): molecule counts that were not whole numbers when the run (or run_until leg) started and were rounded to one. Event-assigned parameters and compartments and rate-rule targets hold continuous values and are never rounded (issue #718).
first_rounded_species (str): name of the first rounded species; "" if none.

All zero/empty on every non-SSA backend.

property psa_diagnostics: dict[str, Any]¶

PSA partial-scaling diagnostics dict (GH #15).

A cheap audit layer for a partial-scaling run: per-reaction, time- integrated scaling statistics plus a measured event-rate speedup. These are measured accept/reject signals for a fixed PSA run, not a certification of accuracy. Inactive (active=False, all zero/empty) on exact SSA and every non-SSA backend.

For reaction r, the integer leap multiplier is m_r(N) = max(1, floor(N'_min / N_c)) (N'_min = smallest participating population, N_c = poplevel), a_r its unscaled propensity, and T the simulated horizon.

Keys:

active (bool): whether this run used partial scaling.
reaction_index (list[int]): 1-based reaction indices, parallel to every per-reaction list below (matches .net reaction numbering).
time (float): T, the simulated horizon (integral denominator).
mbar (list[float]): dwell-time-averaged multiplier m̄_r = (1/T)∫ m_r dt. A non-scaling reaction reports 1.0.
qexc (list[float]): propensity-weighted excess q̄_r = (1/T)∫ (m_r−1)·a_r dt — the correct per-reaction statistic for attributing variance inflation (E[m·a] ≠ E[m]·a(E[N])), and the source term for an affine/monomolecular covariance correction.
mbar_integral / qexc_integral (list[float]): the raw ∫ m_r dt and ∫ (m_r−1)·a_r dt (divide by time for the averages; kept so a transient correction can use the un-normalized integral).
exact_event_integral (float): ∫ Σ_r a_r dt (would-be exact event count).
scaled_event_integral (float): ∫ Σ_r a_r/m_r dt (actual scaled event count).
speedup (float): measured per-path event-rate ratio ŝ = exact/scaled ≥ 1, computable from the scaled run alone. It is not a wall-clock speedup and ignores correction-estimator variance.
peak_population (float): running max over the populations (counts) participating in scaling-eligible reactions.
activation_crossed (bool): whether peak_population ≥ 2·poplevel. Scaling can only engage past 2·poplevel; if this is False no channel ever scaled and the run was pathwise identical to exact SSA.

To turn speedup into a cost/precision go/no-go, pair it with the observable’s variance-inflation ratio ρ (from a calibration/anchor run) via bngsim.psa_cost_decision().

On a stacked batch (squeeze(), squeeze=True) the per-run statistics are not combined into one: active is whether any run used partial scaling, peak_population the largest and activation_crossed whether any run crossed, and the per-reaction lists, the integrals and speedup are left empty (NaN). Each run’s own report is in psa_diagnostics_per_run.

property psa_diagnostics_per_run: list[dict[str, Any]]¶

Each run’s psa_diagnostics, in the order the runs were stacked.

One report for a single run; one per row for a stacked batch (squeeze()), each exactly as that run reported it.

property dataframe¶

Return a pandas DataFrame with time + observables.

Requires pandas (pip install bngsim[pandas]).

Returns:: Columns: time, then one per observable.
Return type:: pandas.DataFrame
Raises:: ImportError – If pandas is not installed.

property xr: _XarrayAccessor¶

Per-field xarray DataArray accessor (AMICI-style).

Each attribute access (result.xr.species, result.xr.observables, etc.) builds a fresh xarray.DataArray with labeled coords. Mirrors AMICI’s rdata.xr.x / .y / .sx ergonomics so downstream code accustomed to that ecosystem can slice by name:

>>> result.xr.species.sel(state="A").values
>>> result.xr.observables.sel(observable="A_tot")
>>> result.xr.sensitivities.sel(parameter="k1", state="A")

Available fields (empty when the corresponding block has zero width):

species — dims (time, state)
observables — dims (time, observable)
expressions — dims (time, expression)
sensitivities — dims (time, state, parameter) (only when sensitivities were computed)
sensitivities_ic — dims (time, state, ic_state) (only when IC sensitivities were computed)

Requires the optional xarray dependency (pip install xarray). For a one-shot xarray.Dataset of every field with a shared time coord, see to_xarray().

Raises:

ImportError – If xarray is not installed (raised on first attribute access).
AttributeError – If the requested field is not a known xarray-exposable block, or if it has zero width.

to_xarray()[source]¶

Bundle every field into a single xarray.Dataset.

One-call complement to xr. The returned Dataset has a shared time coord and exposes species, observables, expressions, and (when present) sensitivities / sensitivities_ic as data variables, each with the same labeled coords as the per-field DataArrays on xr.

custom_attrs is mirrored onto ds.attrs; the stochastic seed (when present) is stored as ds.attrs["seed"], and the stop_condition the result was cut at as ds.attrs["stop_condition"].

Returns:

A Dataset with all populated trajectory fields as data vars.

Return type:

xarray.Dataset

Raises:

ImportError – If xarray is not installed.
RuntimeError – If the result holds 3-D batch arrays. Iterate per-replicate and call to_xarray on each Result.

Examples

>>> ds = result.to_xarray()
>>> ds["species"].sel(species="A").plot()
>>> ds.to_netcdf("result.nc")  # archive via xarray's writer

as_roadrunner(selections=None)[source]¶

Return a bngsim.NamedArray mimicking rr.simulate.

Drop-in shape replacement for libRoadRunner’s simulate(...) return value: a 2-D array whose columns are labeled per selections. Used to swap RR for BNGsim in PyBNF stochastic- fitting workflows without changing call sites.

Parameters:

selections (list[str] | None) –

Column selectors to include, in order. Each entry is one of:

"time" — the time column.
"[X]" — the concentration of species X (BNGsim’s stored species value, which equals amount/V_c uniformly per Phase 2).
"X" — the amount of species X (concentration × V_c); equals the stored value when the species’s compartment volume is 1 (the default and the .net case).

When None, defaults to RoadRunner’s ["time"] + [f"[{s}]" for s in species_names].

Returns:

Shape (n_times, len(selections)); columns ordered to match selections. The result subclasses numpy.ndarray, so all numpy operations work. Indexing carries the column names over — arr[:, 1:] is labeled with the columns the slice kept — and any other transform (arithmetic, a transpose, a reduction) returns an array with no names at all, so a lookup by name on it raises rather than answering with a column that has moved (issue #561). The names also survive a pickle round trip, so an array collected from a worker process arrives labeled (issue #629).

Return type:

Any

Raises:

ValueError – If a selector is unrecognized or names a species not in the model, or if it asks for an amount ("X") from a result loaded from a file that does not record compartment sizes (one written before bngsim saved them, issue #787).
RuntimeError – If the result is a stacked batch (3-D species, from squeeze() or run_batch(..., squeeze=True)); call this on each row’s Result from run_batch(..., squeeze=False).

Warns:

UserWarning – When a "[X]" column is asked of a result loaded from such a file: it is the stored value, which is not the concentration where a compartment’s size changes during a stochastic run, or for an amount-valued species an event resizes.

Examples

Replace rr.simulate(0, 100, 101) (RoadRunner) with the BNGsim equivalent:

sim = bngsim.Simulator(model, method="ssa")
arr = sim.run(t_span=(0, 100), n_points=101).as_roadrunner()
# arr has columns ["time", "[X]", "[Y]", ...]
x_trace = arr["[X]"]

Custom selections:

arr = result.as_roadrunner(
    selections=["time", "[X]", "X", "[Y]"]
)

summary()[source]¶

Return a compact, JSON-serializable description of this result.

The HPC scheduler-free contract (GH #203) pairs a serializable EvaluationSpec (what was run) with a lightweight summary (what came back), so a cluster driver can index/log/checkpoint thousands of evaluations without re-reading every full HDF5 payload (use save()/load() for the complete arrays). The summary carries shapes, output names, sensitivity availability, solver diagnostics, and the seed — every value a built-in JSON-encodable type.

Returns:: Keys: version (bngsim version), is_batch, n_sims (None unless batch), n_times, t_start/t_end (None for an empty time grid), shapes (per-array shape lists), species_names/observable_names/expression_names, the has_sensitivities* flags, sensitivity_params, sensitivity_ic_species, solver_stats, and seed.
Return type:: dict[str, Any]

Examples

>>> import json
>>> s = result.summary()
>>> json.dumps(s)  # always succeeds
'...'

save(path)[source]¶

Save result to HDF5 file.

Parameters:: path (str | Path) – Output file path (e.g. "results.h5").
Raises:: ImportError – If h5py is not installed.
Return type:: None

classmethod load(path)[source]¶

Load result from HDF5 file.

Parameters:

path (str | Path) – Input file path.

Returns:

Loaded result.

Return type:

Raises:

ImportError – If h5py is not installed.
FileNotFoundError – If the file does not exist.

Notes

A file written before bngsim saved compartment sizes (issue #787) loads, but does not say what the sizes were. as_roadrunner() on it refuses a bare-id amount selector and warns on a "[X]" one; the arrays themselves are as they were saved.

to_gdat(path, *, print_functions=False, print_rate_laws=False)[source]¶

Export observables in BNG .gdat format.

Function-column headers are always bare (no () suffix) and identical across every simulation method (issue #58).

Parameters:

path (str | Path) – Output file path.
print_functions (bool) – When True, append the user-named function columns after the observables (auto-generated _rateLawN rate-law columns are still omitted). Default False keeps the observables-only output, matching BNG2.pl’s behavior when print_functions=>1 is not set.
print_rate_laws (bool) – When True, additionally append the auto-generated _rateLawN rate-law columns (bare names). These are internal rate-law intermediates, omitted by default. On a result loaded from HDF5 or assembled from raw arrays they are unavailable (only the filtered set is persisted), so this flag is a no-op there. Default False.

Return type:

to_cdat(path)[source]¶

Export species in BNG .cdat format.

The columns are species, the values this Result reports.

Parameters:

path (str | Path) – Output file path.

Raises:

ValueError – For a squeezed run_batch or run_replicates Result, which holds one trajectory per row: write each row’s .cdat from the unsqueezed results (squeeze=False).
RuntimeError – If the file cannot be opened for writing.

Return type:

to_csv(path, *, kind='observables', sep=',', include_time=True, header=True)[source]¶

Export the trajectory as a plain delimited text file.

Unlike to_gdat() / to_cdat() (BNG-native, #-prefixed, fixed-width space-padded), this writer produces a format SBML/ RoadRunner/Tellurium users expect: a header row of column names followed by data rows, with a user-chosen separator (, for CSV by default; "\t" for TSV; any single character).

The first column is time (unless include_time=False); the remaining columns are the observable or species columns named exactly as they appear on the in-memory Result. No # comment prefix is emitted, so the file can be loaded directly by pandas.read_csv or numpy.loadtxt.

Parameters:

path (str | Path) – Output file path. The caller chooses the path; nothing is written to the current working directory implicitly.
kind (str) – Which trajectory block to write. "observables" (default) matches what .gdat would contain; "species" matches .cdat. Result.expressions is not currently exported via to_csv — use Result.save() (HDF5) for a lossless capture.
sep (str) – Column separator. Default ",". Use "\t" for TSV. Must be a single character.
include_time (bool) – Whether to prepend a time column. Default True.
header (bool) – Whether to write the header row of column names. Default True.

Raises:

ValueError – If kind is not "observables" or "species", if sep is not a single character, or if the trajectory is not a 2-D single-sim array (call this per-replicate on batch results).
OSError – If the file cannot be opened for writing.

Return type:

Notes

Backend coverage: works on every method (ODE / SSA / PSA / NFsim / RuleMonkey) and on results loaded from HDF5. What is not exported by this writer: expressions, sensitivities, solver_stats, custom_attrs, seed. Use HDF5 (Result.save()) when any of those need to round-trip.

Examples

Plain CSV of observables, ready for pandas.read_csv:

sim = bngsim.Simulator(model, method="ode")
result = sim.run(t_span=(0, 100), n_points=101)
result.to_csv("out.csv")

Tab-separated species amounts:

result.to_csv("out.tsv", kind="species", sep="\t")

Headerless data block (e.g. to append to an existing file):

result.to_csv("data.csv", header=False)

static squeeze(results)[source]¶

Stack a list of single-sim Results into one batch Result.

The returned Result has 3D arrays: - species.shape == (n_sims, n_times, n_species) - observables.shape == (n_sims, n_times, n_observables)

Any populated sensitivity block (species, observable, expression, and their IC variants — GH #196) is carried through with the sim axis prepended, e.g. sensitivities.shape == (n_sims, n_times, n_rows, n_cols); a block that is empty on the inputs stays empty. Each row keeps its own compartment sizes and per-sample volume factors.

Parameters:: results (list[Result]) – Results from individual simulations. All must have the same time points, the same species, observable and expression names in the same order, the same sensitivity parameters and IC species, the same reaction labels, and the same blocks present: the stack has one time axis and one set of column names, which every row is read by (issue #1513).
Returns:: A batch result with 3D arrays.
Return type:: Result
Raises:: ValueError – If results is empty, any of them met a stop condition (stop_condition, issue #1369), or a row’s time points, names or blocks differ from the first row’s (issue #1513).

bngsim.run_bngl(path, *, seed=None, rtol=None, atol=None, allow_lossy=False, bng2_pl=None, net_out=None, timeout=600, cache=True)[source]¶

Load a .bngl file and run its actions, as BNG2.pl runs the file.

The network comes from Model.from_bngl() (bng2_pl, net_out, timeout and cache are its arguments), the actions from parse_bngl_protocol(), and they run through run_protocol() (seed, rtol, atol, allow_lossy).

Examples

>>> out = bngsim.run_bngl("egfr_path.bngl")
>>> [run.experiment.t_span for run in out]
[(0.0, 100000.0), (0.0, 120.0)]
>>> out[1].result.observables["Sos_act"][-1]

Parameters:

path (str | Path)
seed (int | None)
rtol (float | None)
atol (float | None)
allow_lossy (bool)
bng2_pl (str | Path | None)
net_out (str | Path | None)
timeout (int)
cache (bool)

Return type:

ProtocolResult

bngsim.run_protocol(model, protocol, *, seed=None, rtol=None, atol=None, allow_lossy=False)[source]¶

Run a BNGL protocol’s steps on model, in order, as BNG2.pl runs them.

Parameters:

model (Model) – The protocol’s network, as Model.from_bngl() loads it from the same file: the state at the file’s generate_network, which already holds the actions before it. It is not changed; the protocol runs on a clone.
protocol (ProtocolSpec) – The steps, from Model.from_bngl(path, protocol=True) or parse_bngl_protocol(). A spec with no generation point (generation_step None: SED-ML, SBML events) runs from its first step.
seed (int | None) – Seeds the stochastic runs that do not name their own seed: each draws its seed from this one, in order. None draws fresh seeds.
rtol (float | None) – Override every ODE run’s tolerances. By default a run uses its own rtol/atol, or run_network’s 1e-8; atol is held exactly.
atol (float | None) – Override every ODE run’s tolerances. By default a run uses its own rtol/atol, or run_network’s 1e-8; atol is held exactly.
allow_lossy (bool) – Run a protocol that has an action bngsim cannot apply, or a run argument it holds only approximately, as nearly as it can, with a ConversionWarning, instead of refusing it.

Returns:

One ProtocolRun per experiment, and the model as the protocol left it.

Return type:

ProtocolResult

Raises:

ConversionError – Before anything runs, if a step cannot be run here: a network-free or PLA run, a run argument such as stop_if, an action that changes results but carries no step (unless allow_lossy), or an argument the protocol holds only approximately (unless allow_lossy). While running, if a step names a parameter or species the model does not have, or resets to a state nothing saved, which BNG2.pl refuses too.

class bngsim.ProtocolResult(runs, model, protocol)[source]¶

Bases: Sequence[ProtocolRun]

What run_protocol() ran: one ProtocolRun per experiment.

Index it by position, or by a run’s suffix. model is the model in the state the protocol left it; protocol is the spec that ran.

Parameters:

runs (Sequence[ProtocolRun])
model (Model)
protocol (ProtocolSpec)

property runs: tuple[ProtocolRun, ...]¶

property labels: list[str | None]¶: Each run’s suffix, in order (None where it has none).

class bngsim.ProtocolRun(experiment, result=None, scan_values=(), scan_results=(), backward_values=(), backward_results=())[source]¶

Bases: object

One experiment of a protocol, as it ran.

A simulate has its time course in result. A parameter_scan has one Result per value in scan_results, in the order they ran, each carrying custom_attrs["scan_parameter"] and custom_attrs["scan_value"] as Simulator.parameter_scan()’s do. A bifurcate has its upward sweep there and the sweep back down in backward_results.

Parameters:

experiment (Experiment)
result (Result | None)
scan_values (tuple[float, ...])
scan_results (tuple[Result, ...])
backward_values (tuple[float, ...])
backward_results (tuple[Result, ...])

experiment: Experiment¶

result: Result | None = None¶

scan_values: tuple[float, ...] = ()¶

scan_results: tuple[Result, ...] = ()¶

backward_values: tuple[float, ...] = ()¶

backward_results: tuple[Result, ...] = ()¶

property kind: str¶: "simulate", "scan" or "bifurcate".

property label: str | None¶: The run’s suffix, if it has one.

property method: str¶

class bngsim.IdentifiabilityReport(fim, parameters, eigenvalues, eigenvectors, rank, condition_number, identifiable, non_identifiable_directions, cramer_rao_bound, threshold)[source]¶

Bases: object

Model-identifiability readout of a Fisher Information Matrix (GH #202).

A pure function of the model + output-sensitivity tensor (no data, no residuals, no objective): the eigen-decomposition of the FIM and the practical-identifiability flags derived from it. This is sloppiness / practical-identifiability analysis of the model, independent of any fit.

Produced by Result.identifiability().

Variables:

fim (ndarray, shape (n, n)) – The (symmetric, positive-semidefinite) Fisher Information Matrix the readout was computed from, over n parameters (or IC species).
parameters (list of str) – Names labelling the FIM axes — the Result.sensitivity_params (axis="parameter") or Result.sensitivity_ic_species (axis="ic") in matrix order.
eigenvalues (ndarray, shape (n,)) – FIM eigenvalues in ascending order. Tiny negative values from round-off are clipped to 0. eigenvalues[0] is the smallest — the sloppiest (least-constrained) direction.
eigenvectors (ndarray, shape (n, n)) – Orthonormal eigenvectors as columns, aligned with eigenvalues: eigenvectors[:, i] is the parameter combination for eigenvalues[i]. A near-zero eigenvalue marks its eigenvector as a practically non-identifiable parameter combination.
rank (int) – Numerical rank — the number of eigenvalues above threshold.
condition_number (float) – $\lambda_\max / \lambda_\min$; inf when rank-deficient. A large value signals an ill-conditioned (sloppy) parameter set.
identifiable (ndarray of bool, shape (n,)) – Per-eigen-direction flag (aligned with eigenvalues): True where the eigenvalue exceeds threshold, False for a practically non-identifiable / sloppy direction.
non_identifiable_directions (list of int) – Indices into eigenvalues / eigenvectors columns that are flagged non-identifiable (the False entries of identifiable), smallest eigenvalue first.
cramer_rao_bound (ndarray, shape (n, n)) – $\text{FIM}^{-1}$ — the Cramér–Rao lower bound on parameter variance, an identifiability aid only (a property of the FIM). This is not a data/noise-weighted fit covariance. All-NaN when the FIM is rank-deficient (the inverse is undefined; a warning is emitted).
threshold (float) – Eigenvalue cutoff (rtol * lambda_max) below which a direction is flagged non-identifiable.

Parameters:

fim (NDArray[double])
parameters (list[str])
eigenvalues (NDArray[double])
eigenvectors (NDArray[double])
rank (int)
condition_number (float)
identifiable (NDArray[bool])
non_identifiable_directions (list[int])
cramer_rao_bound (NDArray[double])
threshold (float)

fim: NDArray[float64]¶

parameters: list[str]¶

eigenvalues: NDArray[float64]¶

eigenvectors: NDArray[float64]¶

rank: int¶

condition_number: float¶

identifiable: NDArray[bool]¶

non_identifiable_directions: list[int]¶

cramer_rao_bound: NDArray[float64]¶

threshold: float¶

property is_identifiable: bool¶: True when every direction is identifiable (full-rank FIM).

class bngsim.EvaluationSpec(model_source, model_format='net', method='ode', t_span=(0.0, 100.0), n_points=101, params=dict(), sensitivity_params=(), sensitivity_ic=(), sensitivity_method='staggered', outputs=(), rtol=None, atol=None, max_steps=None, model_sha256=None, output_start=None, seed=None)[source]¶

Bases: object

A serializable, stateless description of one bngsim evaluation.

Parameters:

model_source (str) – For a path format (net/sbml/antimony), the model file path. For a *_string format, the inline model text.
model_format (str) – One of "net", "sbml", "antimony", "sbml_string", "antimony_string".
method (str) – Simulation method. Default "ode" (the only method with sensitivities).
t_span (tuple[float, float]) – (t_start, t_end) integration interval.
n_points (int) – Number of output time points (including t_start, or output_start when one is set).
output_start (float | None) – Where the output starts, when that is after t_start: the run integrates from t_start and returns n_points even times over (output_start, t_end). This is a SED-ML uniform time course whose outputStartTime is past its initialTime (issue #1493). None (the default) outputs from t_start. It takes sensitivity_params across the unreported part as carry_sensitivities does, and refuses sensitivity_ic, which does not cross it.
params (Mapping[str, float]) – Parameter overrides applied to the model before running (the θ vector).
sensitivity_params (Sequence[str]) – Parameter names to compute forward output sensitivities for. Empty disables parameter sensitivities.
sensitivity_ic (Sequence[str]) – Species names whose initial conditions to differentiate against.
sensitivity_method (str) – CVODES sensitivity method ("staggered" / "simultaneous").
outputs (Sequence[str]) – Output selectors (species:/observable:/expression: …) that the caller intends to read off the result. Recorded for provenance; not applied during integration.
rtol (float | None) – Relative solver tolerance. None uses the Simulator default.
atol (float | Sequence[float] | str | ExactAtol | None) – Absolute solver tolerance. None uses the Simulator default. A sequence is the per-species tolerance of issue #196, ordered like Model.species_names and normalized to a tuple so the spec stays hashable and byte-stable; "auto" asks the worker to derive one from the model it materializes. Both survive the JSON round trip, so a decade-spanning model can be distributed at a tolerance that means something for every species. "tracking" is the over-time mode of issue #213 at its default depth — a token, so it round-trips too. A bngsim.TrackingAtol with a non-default depth or an explicit ceiling does not: it is not JSON, and a spec is a wire format. A float is read as a ceiling, as bngsim.Simulator.run() reads one (issue #1420); a bngsim.ExactAtol holds every species to it exactly, and goes out as {"exact": value}. Distribute the token, or the ceiling as a plain sequence.
max_steps (int | None) – Max internal solver steps per output point. None uses the default.
model_sha256 (str | None) – Hex SHA-256 of the model source (file bytes for a path format, the UTF-8 text for a *_string format). When set, build_model() verifies the live source matches and raises on mismatch — a cluster integrity guard against a stale/edited model file on a shared filesystem.
seed (int | None) – The run’s random seed, from 0 to 2**64 - 1, passed to bngsim.Simulator.run() (issue #1517). A stochastic spec with a seed gives one trajectory on every node; without one, each evaluate() draws a fresh seed and reports it on Result.seed, so replace(spec, seed=result.seed) reruns that trajectory. An ODE spec reads it only to order simultaneous equal-priority events (GH #242). With output_start, the output leg runs at seed, which Result.seed reports, and the leg up to the output start at the second seed drawn from it, as bngsim.Simulator.run_batch() draws its row seeds (issue #1365), so the two legs do not replay one stream. An unseeded stochastic spec draws that seed before the first leg, so its Result.seed reruns both. Emitted by to_dict() only when set, so a spec without one keeps its bytes. SED-ML carries it as KISAO:0000488.

model_source: str¶

model_format: str = 'net'¶

method: str = 'ode'¶

t_span: tuple[float, float] = (0.0, 100.0)¶

n_points: int = 101¶

params: Mapping[str, float]¶

sensitivity_params: Sequence[str] = ()¶

sensitivity_ic: Sequence[str] = ()¶

sensitivity_method: str = 'staggered'¶

outputs: Sequence[str] = ()¶

rtol: float | None = None¶

atol: float | Sequence[float] | str | ExactAtol | None = None¶

max_steps: int | None = None¶

model_sha256: str | None = None¶

output_start: float | None = None¶

seed: int | None = None¶

to_dict()[source]¶

Return a deterministic, JSON-encodable dict of this spec.

params is emitted sorted by name and sequences as lists, so the dict (and to_json()) is byte-stable for a fixed spec — usable directly as a content key for caching/deduplication. output_start and seed are emitted only when set, so a spec without them keeps the bytes it always had. A stochastic spec without a seed is one key for many trajectories: set seed for a key that names one run.

Return type:: dict[str, Any]

classmethod from_dict(data)[source]¶

Reconstruct a spec from a to_dict() mapping (extra keys rejected).

Parameters:: data (Mapping[str, Any])
Return type:: EvaluationSpec

to_json(*, indent=None)[source]¶

Serialize to a JSON string (keys sorted for byte-stability).

Parameters:: indent (int | None)
Return type:: str

classmethod from_json(text)[source]¶

Deserialize from a to_json() string.

Parameters:: text (str)
Return type:: EvaluationSpec

with_params(params, *, merge=False)[source]¶

Return a copy with params replaced (or merged when merge=True).

The ergonomic path for a parameter sweep / multistart: serialize one base spec, then stamp each θ row through with_params on the worker. Returns a new frozen instance; the original is untouched.

Parameters:

params (Mapping[str, float])
merge (bool)

Return type:

EvaluationSpec

compute_source_sha256()[source]¶

SHA-256 of the model source (file bytes for a path, UTF-8 text inline).

Return type:: str

build_model()[source]¶

Load the model from model_source per model_format.

When model_sha256 is set, the live source is hashed and compared first; a mismatch raises ValueError (cluster integrity guard).

Return type:: Model

build_simulator()[source]¶

Build a Simulator with θ applied and sensitivities wired.

Return type:: Simulator

evaluate()[source]¶

Materialize and run this evaluation, returning the Result.

Deterministic for fixed inputs (model, θ, sensitivity set, solver options, and seed for a stochastic method). Read named outputs via Result.outputs() and gradients via Result.output_sensitivities() with outputs as selectors.

Return type:: Result

class bngsim.SteadyStateResult(core)[source]¶

Bases: object

Result of a steady-state computation.

Variables:

concentrations (ndarray) – Species concentrations at steady state, shape (n_species,).
species_names (list[str]) – Species names.
residual (float) – max|f(y)| at convergence.
newton_distance (float) – The Newton step from the state returned to the root, in each species’ own tolerance, which under criterion="scaled" decided where the solve stopped (issue #1001); at most 1 is within it. NaN under criterion="residual", which does not take the step, and above 512 unknowns; inf where the Jacobian is singular at a state that is not a root, or a differenced one straddles a jump.
continued_run_distance (float) – The run test of criterion="scaled": the largest move of a measured species, over its tolerance, across the run taken on from the state returned; at most 1 passes (issue #1001). inf where the run could not be taken to its time; NaN where none was taken (under "residual", where the parity test already refused, and for a burst-free Newton solve).
criterion (str) – The test converged was decided by, "scaled" or "residual".
method_used (str) – "integration" or "newton".
converged (bool) – Whether steady state was found.
n_steps (int) – Number of solver steps.
n_rhs_evals (int) – Number of RHS evaluations.
n_state_switch_roots (int) – The crossings of a rate-law condition over model state (or of an SBML model’s discontinuity trigger) the march stopped at, so a window the approach passes through is integrated over (issue #904). It restarts on each one up to 1000; the 1001st, where it stops rooting them, is counted too.
state_switch_roots_dropped (bool) – Whether the march stopped rooting those conditions after 1000 crossings, as a state that slides along its threshold makes, and marched over the switch from there as it did before #904. If that march ends further from a steady state than where it stopped rooting, the state there is returned.
n_residual_species (int) – How many species entered ||f||_2 / n — n_species unless steady_state(mask=...) restricted it (issue #74). Assert on this rather than inferring from a residual that moved.
excluded_species (list[int]) – 0-based indices the mask excluded, ascending; empty without a mask. Their concentrations are returned but their settling was never tested, and their sensitivity rows are NaN.
unconverged_pure_sinks (list[str]) – On a failed solve, the write-only accumulator species that were in the convergence test and are carrying flux at the returned state — the structural reason the residual has a floor (issue #74). Empty otherwise, including on every converged solve. See Model.pure_sink_species().
root_stability (str) – Whether the system can rest on the returned root (issue #78), from the eigenvalues of the Jacobian restricted to the species the Newton polish solved for: "stable" (every eigenvalue in the closed left half-plane), "undetermined" (the certificate declined — more than 512 unknowns, or a Jacobian whose spectrum it could not compute), or "unstable". "" when the result came from integration, which needs no certificate: a trajectory cannot come to rest on an unstable equilibrium. "" too from a burst-free solve told to skip it (steady_state_batch(newton_burst=False, certify=False), issue #1191). "unstable" is returned only when the initial condition itself was that root; a polish that lands on one is discarded instead.
n_unstable_roots_rejected (int) – How many candidate Newton roots this solve discarded as unstable (issue #78). Non-zero means the polish landed on a saddle the trajectory was merely passing near — on a bistable model the seed-stability guard cannot catch that, because near a separatrix the trajectory slows down and two successively tighter bursts hand Newton the same seed. Explains a method_used of "integration" from a method="newton" solve.
eigenvalues (ndarray, complex128) – The spectrum root_stability was read off (issue #523): the eigenvalues of the Jacobian restricted to the species the Newton polish solved for, sorted by descending real part, a conjugate pair adjacent with its +imaginary member first. One per unknown — n_species less one per conservation law less the masked species — so the full system’s spectrum is this plus one zero per conservation law. Filled whenever the certificate computed a spectrum, whatever it concluded; empty (shape (0,)) when it did not: an integration result, or a root it declined on before the eigensolver — more than 512 unknowns, a solver failure — which is every "undetermined" except an all-zero spectrum, reported as-is. Compare against numpy.linalg.eigvals(model.jacobian(ss.concentrations)) to see the conservation-law zeros the restriction removes.
sensitivity (ndarray or None) – Species dY_ss/dp matrix, shape (n_species, n_params). None if no sensitivity was requested.
rhs_backend (str) – Which RHS actually evaluated f(y): "exprtk" (interpreted), "codegen-so" (compiled + dlopen’d), or "codegen-jit" (in-process MIR). Issue #63 — before it, a steady-state solve ignored the Simulator’s codegen artifact, so this was always effectively "exprtk" no matter what Simulator.codegen_backend reported.
solver_jacobian_source (str) – Which Newton matrix the solver tiers factored (issue #127): "codegen" (the compiled analytical Jacobian), "analytical" (the interpreted fill), or "finite-difference" (CVODE’s and KINSOL’s own difference quotients, one RHS evaluation per unknown per Jacobian setup). Reported on every solve, and it covers both tiers — the CVODE march and the KINSOL polish take the same gate, so they cannot disagree. Before #127 neither installed a callback and this was always the difference quotient, on models whose closed form was already loaded.
solver_jacobian_retried (bool) – Whether the closed-form Jacobian had to be called off: CVODE gave up on the march with it installed, and the solve was retried on difference quotients (issue #127, the steady-state half of GH #176). The usual cause is a rate law that is discontinuous in a state variable, whose exact derivative omits the jump. The returned answer is the retry’s, so solver_jacobian_source reads "finite-difference" alongside this. Only jacobian="auto" retries; "analytical" surfaces the failure.
linear_solver (str) – Which direct linear solver the CVODE march factored its Newton matrix with (issue #128): "klu" (sparse CSC), "dense" (the built-in dense LU) or "lapack-dense" (the GH #84 BLAS factor). Chosen by the same size/density/force-flag rule Simulator.run() uses, so the two agree on a model. Before #128 the steady-state paths had no sparse route and this was always dense — 1.9x to 3.1x of wall clock on the 1000+ species corpus models, for the same answer. The KINSOL polish and the dY_ss/dp solve factor a reduced matrix and are always dense; this field does not describe them.
sens_dfdp_source (sens_jacobian_source,) – How each factor of dY_ss/dp = -J⁻¹·(∂f/∂p) was built. sens_jacobian_source is "codegen" (compiled analytical Jacobian), "analytical" (interpreted analytical), "automatic-differentiation" (exact, by dual numbers through the interpreter, where there is no closed form or jacobian="fd" asked for none; issue #1330), or "finite-difference" (a model the dual evaluation cannot follow: a ring of functions that read each other, rateOf); sens_dfdp_source is "codegen", "automatic-differentiation" or "finite-difference". Both are "" when no sensitivity was requested. Issue #63 — both factors were unconditionally finite-differenced before it, and the result gave no way to tell.
sens_output_source (str) – What answered the expression block d(func)/dp (issues #75, #1002, #1330): "codegen" (the compiled bngsim_codegen_output_sens chain rule, the same evaluator a CVODES forward-sensitivity run() uses, for every function), "automatic-differentiation" (dual numbers through the interpreter, for every function: those the chain rule declines, or every one of a model with none), "codegen+automatic-differentiation" (every function, by the two), "partial" (some), or "none" (not one). A row that has no answer is NaN, and output_sensitivities() refuses it by name, with the reason: a function’s partials are not differenced at a steady state (issue #1002). "" when no sensitivity was requested or the model has no global functions. The observable block is an exact linear projection either way and is not covered by this field.
sens_jacobian_rcond (float) – min|U_jj| / max|U_jj| from the LU of the (reduced) Jacobian that was inverted — how close to singular the sensitivity system was. dY_ss/dp exists only when that Jacobian has full rank; a steady state that is a continuum rather than an isolated point makes it rank-deficient. Well-posed corpus models measure 1e-4 to 1e-1; rank-deficient ones 1e-12 to 1e-9. 0.0 when no sensitivity was requested. It depends on the units of the species and is not what decides whether the columns are returned: the sens_root_* measures are (issue #995), and a result that carries sensitivities has passed them.
sens_root_determinant_ratio (float) – The determinant of the system dY_ss/dp is solved on at the state one Newton step on, over the one at the returned state (issue #995). 1, to the accuracy of the solve, at an isolated root; next to nothing where the Jacobian is singular at the steady state, a continuum of steady states; 1/2 at a double root; negative or large where the returned state is far from its root or a rate law is discontinuous between the two. steady_state raises outside 0.6 to 1.67. 1.0 when no sensitivity was requested.
sens_root_pivot_share (float) – The least a pivot of that system’s factorization is of the terms it was computed from: 1 for a pivot nothing cancelled into, 1e-16 for one that is what rounding left of a zero, where the Jacobian is singular whatever the state. Read at the state the solver stopped at, as the condition number is. steady_state raises below 1e-13. 1.0 when no sensitivity was requested.
sens_root_condition (float) – The componentwise condition number of that system, the Perron root of |A⁻¹|·|A|: how many times a relative error in each entry of the Jacobian is magnified in the columns. The same in any units. steady_state raises above 1e13. 1.0 when no sensitivity was requested.
sens_root_column_shift (float) – The largest move of a column of dY_ss/dp when it is solved again a Newton step on, as a fraction of the column’s largest entry at the two states and of no less than 1e-3/|p|, each entry over its species’ scale (sens_species_scale): how far the columns are from those of the root. Where a column moves by more than 0.01 of itself the solve steps the state on, up to ten times, and this is the move at the last step; steady_state raises where it is still above 0.01. 0.0 when no sensitivity was requested.
sens_root_state_shift (float)
sens_root_state_species (str or None) – The largest move of a species under that Newton step, over its scale, and the species. A state that moves by more than 0.01 is stepped on as one whose columns move is, and steady_state raises where it still does after ten steps.
sens_root_newton_steps (int) – The Newton steps the returned state and columns are from the state the solver stopped at: 0 where neither it nor its columns moved, and the result is the solver’s own. Otherwise concentrations, residual and the sensitivities are those of the stepped state, which is nearer the root: tol bounds the residual and not the distance to it, and a model whose concentrations are small passes it a long way off. steady_state raises where a state that was stepped to is more than 1% in any species from where the run that is taken on ends.
sens_root_hold_drift (sens_root_hold_shift,) – A run is taken on for max_time from the state the solver stopped at, with every concentration moved by up to a millionth of itself (a species that is at nothing is not moved). sens_root_hold_shift is the largest move of a column when it is solved again where that run ends, measured as the column shift is, and steady_state raises above 0.01, or where the run could not be made (not a number): the state is one a run is passing, or a root the system leaves. sens_root_hold_drift is the largest move of a species over the run, against what its entries are taken over there: the larger of its two values, or for a species at a zero what stands beside it. 0.0 when no sensitivity was requested.
sens_root_hold_steps (int)
sens_root_hold_time (float) – The steps that run took and the time it reached. The run has max_steps steps, and steady_state raises where they are used up short of max_time: the run was not seen to stay.
sens_root_spectral_radius (sens_root_growth_rate,) – The largest real part among the eigenvalues of that system at the returned state, and the largest eigenvalue in size, for systems of up to 512 unknowns (above that, those of the species the state has at a zero). A species that is at nothing is among them, with what it would do if there were any, also where nothing makes it. steady_state raises where the first, times max_time, is above 0.01: a state beside this one is more than 1% further off after a run of that length. Not a number, and 0.0, where the eigenvalues were not taken and when no sensitivity was requested.
sens_root_relaxation (float) – The most of a column of dY_ss/dp that a run of max_time would leave unestablished, as a fraction of the column’s largest entry: the bound A⁻¹·column / max_time where that is under 1e-3, and what eight implicit steps of max_time/8 leave of the column where it is not. steady_state raises above 0.01. 0.0 when no sensitivity was requested.
sens_species_scale (numpy.ndarray) – What an entry for each species is small against: its own concentration, the larger of the returned one and the corrected one. A species at a zero has none. It is one the corrected state has at nothing and that nothing left there makes (its rate is zero with every such species set to zero, and it does not grow from next to nothing), and it is taken over where it has been, its starting value among them, and the largest such among the species the Jacobian couples it to, no more than a conserved total it belongs to allows. A species that is small and has a steady value, 1e-12 beside another at 1, keeps its own. Empty when no sensitivity was requested.
sens_mask_reader_species (sens_mask_held_species,) – A species mask= left out that is held where the solve left it, and a kept species whose rate reads it. steady_state raises where there is one: the columns are not those of the steady state then.
sens_root_pivot_species (sens_root_determinant_species,)
sens_root_hold_species (sens_root_condition_species,) – The species each of those points at: for the determinant, the one whose pivot moved furthest.
sens_root_hold_param (sens_root_column_param,)
sens_root_relaxation_param (str or None) – The parameter whose column each of those is read at.
sens_fd_own_step_species (list of str) – Where sens_jacobian_source is "finite-difference": the species far below the state’s scale whose columns of the Jacobian were taken again at their own step, for the columns to be solved a second time with (issue #1002). Empty where there is none, or the Jacobian is in closed form.
sens_fd_half_step_species (list of str) – The same for the columns taken again whole at half the floored step, for a third solve: those of the species far below the state’s scale, and of species at 0. A secant across the step moves with it and a slope does not, also in a row that a much larger species is in, which the own step cannot read.
sens_fd_disagree_param (sens_fd_disagree_species,)
sens_fd_half_step (sens_fd_floored, sens_fd_own_step,) – The entry of dY_ss/dp furthest from the floored answer, by more than 1e-4 of itself (1e-10 of its column’s largest, at least), and its three values: the one with every species stepped against the state’s scale, which is what sensitivity holds, the one with those species at their own step, and the one with them at half the floored step (each the floored one where that matrix has no column taken). steady_state raises where there is one. None where they agree.
sens_tie_param (sens_tie_species,)
sens_tie_below (sens_tie_value, sens_tie_above,) – Where the Jacobian or ∂f/∂p was taken by dual numbers ("automatic-differentiation") and the steady state sits where a rate law changes branch — if(X > c, ...) at X = c, abs(X) or max(X, 0) at X = 0, a step of floor() — so that the derivative of the branch it is on is one of two one-sided ones: the entry of dY_ss/dp furthest apart when the columns are solved with the matrices from the two sides, by more than 1e-8 of itself (1e-10 of its column’s largest, at least), its value in sensitivity and its values from above and below (issue #1330). steady_state raises where there is one: the steady state has no derivative there. None where the side does not matter.
sens_tie_jump (bool) – Whether that entry was refused because a rate law’s value jumps along its species there and the parameter moves the species, so that the steady state crosses the jump (the two sides may agree then).
sens_held_species (str or None) – Where the state is held on a jump in a rate law’s value within the solver’s tolerance of it, with no root of its branch beside it: a Newton step to that branch’s root changed some law’s branch and left a larger residual than it started from. The species whose residual the step left largest (issue #1330); steady_state raises where there is one. None where the step changed no branch or took the residual down.
sens_held_after (sens_held_before,) – The largest residual before and after that step, each over its species’ terms (at least 1e-9 of the largest species’).

Notes

When sensitivity_params is passed to Simulator.steady_state(), the result also carries the observable- and expression-level steady-state forward sensitivities (GH #12): read them by name with output_sensitivities(), mirroring bngsim.Result.output_sensitivities() on a CVODE run. These are the chain-rule projection of the species dY_ss/dp onto the model’s observables and global functions, so a gradient consumer gets ∂(observable)/∂θ directly without re-deriving the output Jacobian.

Examples

>>> ss = sim.steady_state()
>>> ss.converged
True
>>> ss["A"]  # species A at steady state
50.0
>>> ss.concentrations
array([50., 25., ...])
>>> ss = sim.steady_state(sensitivity_params=["k_deg"])
>>> ss.output_sensitivities(["observable:Stot"])  # (n_sel, n_params)
array([[-1.25]])

residual¶

newton_distance¶

continued_run_distance¶

criterion¶

method_used¶

converged¶

n_steps¶

n_rhs_evals¶

n_state_switch_roots¶

state_switch_roots_dropped¶

n_residual_species¶

excluded_species¶

unconverged_pure_sinks¶

root_stability¶

n_unstable_roots_rejected¶

eigenvalues¶

rhs_backend¶

solver_jacobian_source¶

solver_jacobian_retried¶

linear_solver¶

sens_jacobian_source¶

sens_dfdp_source¶

sens_output_source¶

sens_jacobian_rcond¶

sens_root_determinant_ratio¶

sens_root_condition¶

sens_root_column_shift¶

sens_root_relaxation¶

sens_root_pivot_share¶

sens_root_hold_drift¶

sens_root_hold_shift¶

sens_root_hold_steps¶

sens_species_scale¶

sens_root_hold_time¶

sens_root_newton_steps¶

sens_root_state_shift¶

sens_root_growth_rate¶

sens_root_spectral_radius¶

sens_root_pivot_species¶

sens_root_hold_species¶

sens_root_state_species¶

sens_mask_held_species¶

sens_mask_reader_species¶

sens_root_determinant_species¶

sens_root_condition_species¶

sens_root_hold_param¶

sens_root_column_param¶

sens_root_relaxation_param¶

sens_fd_own_step_species¶

sens_fd_half_step_species¶

sens_fd_disagree_species¶

sens_fd_disagree_param¶

sens_fd_floored¶

sens_fd_own_step¶

sens_fd_half_step¶

sens_tie_species¶

sens_tie_param¶

sens_tie_value¶

sens_tie_above¶

sens_tie_below¶

sens_tie_jump¶

sens_held_species¶

sens_held_before¶

sens_held_after¶

property concentrations¶: Steady-state species concentrations.

property species_names: list[str]¶: Species names.

property sensitivity¶

Sensitivity matrix dY_ss/dp, shape (n_species, n_params).

None if no sensitivity was requested.

property sensitivity_params: list[str]¶: Parameter names for sensitivity.

property observable_names: list[str]¶

Observable names labelling the observable output-sensitivity rows.

Populated on a sensitivity run (sensitivity_params=[...]); empty otherwise. Provided for parity with bngsim.Result.observable_names.

property expression_names: list[str]¶

Expression (global-function) names for the expression rows.

Bare, user-facing names (the auto-generated _rateLawN intermediates are filtered out, matching bngsim.Result.expression_names). Populated on a sensitivity run; empty otherwise.

property sensitivities_observables: ndarray¶

Observable steady-state sensitivities d(observable)/dp.

Shape (n_observables, n_params) on a sensitivity run, aligned with observable_names (rows) and sensitivity_params (columns); empty (0, 0) otherwise. This is the exact linear projection of the species sensitivity through each observable’s group factors — including the amount-valued volume factor for an SBML hasOnlySubstanceUnits="true" species, whose observable denotes an amount rather than the stored concentration (issue #119). So for such a species this block is not simply sensitivity re-weighted by the group coefficients; it is the derivative of the observable’s value.

property sensitivities_expressions: ndarray¶

Expression (global-function) steady-state sensitivities d(func)/dp.

Shape (n_expressions, n_params) on a sensitivity run, aligned with expression_names (rows) and sensitivity_params (columns); empty (0, 0) otherwise. Carries the function’s full total derivative — the state-chain term (∂func/∂x)·dY_ss/dp plus the function’s explicit parameter dependence ∂func/∂p.

resolve_outputs(selectors)[source]¶

Resolve typed output selectors to structured metadata.

Same selector grammar as bngsim.Result.resolve_outputs() (species:/observable:/expression: with state:/function: aliases, () handling, and bare-name uniqueness). Observable and expression names resolve only on a sensitivity run.

Parameters:: selectors (str | Iterable[str])
Return type:: list[dict[str, Any]]

output_sensitivities(selectors, *, axis='parameter')[source]¶

Return steady-state d(named output)/dp for each selector, stacked.

The steady-state analogue of bngsim.Result.output_sensitivities(): resolves each selector and stacks the matching steady-state sensitivity row, so a gradient consumer reads ∂(observable)/∂θ / ∂(expression)/∂θ directly instead of re-deriving the output Jacobian.

species: selectors read the species dY_ss/dp rows; observable: selectors the exact linear group projection; expression: selectors the compiled chain rule’s total derivative of the global function (state chain + explicit parameter dependence), and refuse by name a function it has none for (issue #1002).

Parameters:

selectors (str | Iterable[str]) – Selectors accepted by resolve_outputs().
axis (str) – Only "parameter" (the default) is computed here. A stable steady state keeps nothing of its initial conditions but the conserved totals, so the "ic" axis is zero in a model with no conservation law and the projection onto the totals in one with any; it is not computed, and requesting it raises ValueError. A parameter that sets an initial amount has its share of the totals in its own "parameter" column (issue #704).

Returns:

Shape (n_selectors, n_params), one row per selector in input order (no time axis — a steady state is a single point). An empty selector list yields a (0, n_params) array.

Return type:

Raises:

ValueError – If axis is not "parameter"; if no parameter sensitivities were computed (run with sensitivity_params=[...]); or if a selector names a kind whose sensitivities are unavailable.
TypeError – Propagated from resolve_outputs().

Examples

>>> ss = sim.steady_state(sensitivity_params=["k_deg"])
>>> ss.output_sensitivities(["observable:Stot", "expression:foo"]).shape
(2, 1)

to_dict()[source]¶

Return species concentrations as a dict.

Return type:: dict[str, float]

class bngsim.SteadyStates[source]¶

Bases: object

The distinct steady states Simulator.steady_states() found (issue #524).

One row per distinct root, in the order first reached.

Variables:

states (ndarray, shape (n_roots, n_species)) – Each root’s concentrations, ordered like species_names: the start that reached it with the smallest residual, the extra start from Simulator.steady_state()’s root among them.
stability (list[str]) – "stable", "unstable" or "undetermined" per root, from the stability certificate.
eigenvalues (list[ndarray]) – The reduced Jacobian’s spectrum at each root (empty where the certificate did not compute it).
multiplicity (ndarray of int) – How many starts reached each root: a rough measure of its basin. The root of Simulator.steady_state(), solved as one more start, is no draw and counts for nothing, so a root only it reached reads 0 (issue #1189).
basins (ndarray of int, shape (n_starts,)) – Root index per start, -1 where the solve did not converge.
residuals (ndarray) – ||f||_2 / n at each root.
n_failed (int) – Starts whose solve did not converge.
species_names (list[str])

states: ndarray¶

stability: list[str]¶

eigenvalues: list[ndarray]¶

multiplicity: ndarray¶

basins: ndarray¶

residuals: ndarray¶

n_failed: int¶

species_names: list[str]¶

ROUNDOFF: float = 1e-09¶: What a solve leaves in a species it places at a root, as a share of the scale that species is computed against: its smallest law’s for a law member, and for every species the largest value it takes at any converged start (a root at 0 is reached at 0 or a hair beside it: -5.9e-24 on #1001’s cubic at a nanomolar scale).

classmethod cluster(results, species_names, *, dedupe_rtol, atol=None, laws=None, scale=0.0)[source]¶

Group converged results into distinct roots.

Two results are one root when every species agrees to dedupe_rtol of the larger magnitude, or to an absolute floor in concentration. atol is that floor, a scalar or one value per species. None derives it per species as the larger of two terms: ROUNDOFF of the largest value the species takes in any converged result, the round-off a solve leaves in it at a root where it is zero; or double precision’s resolution (ε) beside the largest magnitude the solves saw, the larger of scale and the largest value any result holds. scale is the largest magnitude of the starts, 0 to read it from the results alone. The second term is for a species at 0 at every root, where the first is round-off itself: where every species is at 0 at every root, the results hold nothing else, and without scale a model whose only root is 0 came back as 108 roots. The first is the species’ own scale, not the model’s: a switch in a pool of 1e-3 beside a pool of 1e8 (#1265) is told apart at its own 1e-3. And it is a distance in concentration, not the solve’s residual tolerance: Simulator.steady_states() once passed its tol, and at 1e-9 in molar units it grouped the 0.2 nM and 0.5 nM roots of a nanomolar switch as one (issue #1518).

laws, shape (n_laws, n_species): a species in a conservation law is told apart only to ROUNDOFF of the smallest of its laws’ scales, Σ_j |L_kj·y_j|, the round-off a solve leaves in it. egfr_path’s EGF, 5e-8 in a law over 1.8e5, came back from 40 starts at 34 values, and each counted as a root of its own. The floor is not dedupe_rtol’s: at 1e-6 of a pool of 1e8 it merged the toggle’s saddle into an attractor. It still matters beside the floor above: a member far below its law’s total, as EGF is, carries the law’s round-off rather than its own, and a coefficient above 1 multiplies it.

A result joins a root by agreeing with the first result that reached it. The state reported is the member with the smallest residual, and comparing with that instead made the grouping depend on the order of the residuals.

Parameters:

dedupe_rtol (float)
atol (float | ndarray | None)
scale (float)

class bngsim.NfsimSession(xml_path, *, molecule_limit=None, connectivity=None, v1143_compat=False, block_same_complex_binding=True, traversal_limit=None)[source]¶

Bases: object

Stateful NFsim simulation session with per-step control.

Use this class when you need fine-grained control over an NFsim session: setting parameters before initialization, running multiple simulation segments, and mutating molecule or exact species counts between segments.

For simple one-shot NFsim simulations, bngsim.Simulator(model, method="nf", xml_path=...) is easier.

Parameters:

xml_path (str | Path) – Path to a BNG XML file (generated by BNG2.pl writeXML).
molecule_limit (int | None) – Global molecule limit (NFsim -gml flag).
connectivity (bool | None) – Enable NFsim’s reaction-connectivity optimization.
v1143_compat (bool) – Enable NFsim v1.14.3 selector compatibility mode for same-seed trajectory reproducibility with the standalone CLI.
block_same_complex_binding (bool) – When True (the default), a bimolecular rule cannot bind two molecules of one complex: BNGL’s semantics, as BNG2.pl’s network applies them. NFsim blocks same-complex binding under -cb (or -bscb, the same flag), which BNG2.pl’s simulate_nf passes by default (complex=>1). False allows it, so a binding rule can close a ring inside one complex. Either way, a rule whose products are written apart (A.B -> A + B) and that deletes a bond or a molecule does not fire where the molecules it acts on would stay in one complex across products, as BNG2.pl’s network requires. NFsim v1.14.3 (BioNetGen 2.9.3) never makes that test, and current NFsim makes it only under -cb, for one-reactant rules. Not tested yet: a product held apart only by a molecule the rule deletes, or tied to the rest of its pattern only by .; and a product written joined (A.B -> A.B) that breaking the bond splits. With False, a bimolecular rule that deletes nothing (A + B -> A + B changing states) may match two molecules of one complex and fire. With False a gelling model can run far slower: it fills with ring bonds such a rule may never break, and each attempt on one is refused. Complex bookkeeping for Species-typed observables is enabled whenever the model declares one, whatever this setting.
traversal_limit (int | None) – NFsim -utl N: universal traversal limit. None (default) lets NFsim auto-compute a suggested limit from the XML.

Raises:

RuntimeError – If NFsim support is not available in this bngsim build.

Examples

Basic one-shot simulation:

with bngsim.NfsimSession("model.xml") as nf:
    nf.initialize(seed=42)
    result = nf.simulate(0, 100, n_points=101)

Multi-step protocol:

with bngsim.NfsimSession("model.xml", molecule_limit=100000) as nf:
    nf.set_param("kf", 0.5)
    nf.initialize(seed=42)
    r1 = nf.simulate(0, 50, n_points=51)
    nf.add_species("L(r)", 500)
    r2 = nf.simulate(50, 100, n_points=51)

Parameter scan (new session per point):

for kf_val in [0.1, 1.0, 10.0]:
    with bngsim.NfsimSession("model.xml") as nf:
        nf.set_param("kf", kf_val)
        nf.initialize(seed=42)
        result = nf.simulate(0, 100, n_points=101)

set_molecule_limit(limit)[source]¶

Set the global molecule limit.

Can be called before or after initialize() to change the NFsim -gml cap.

Parameters:: limit (int) – Maximum number of molecules per molecule type.
Return type:: None

set_param(name, value)[source]¶

Set a parameter value and propagate to dependents.

May be called before or after initialize(). Pre-init writes are stored and applied on session start; post-init writes update the live NFsim System in place. In both cases bngsim re-evaluates every parameter whose BNG XML expr= transitively depends on name, so derived values like LT = LT_conc_M*NA*V_sim and _rateLaw* follow automatically.

Pre-init writes are baked into the XML that NFsim parses, so any <Species concentration="X"> seed-species expression that references name (directly or through a derived parameter) is evaluated against the override-resolved namespace when the agent population is created. This matches BNG2.pl’s behavior when setParameter is invoked between generate_network and a method=>"nf" simulate.

Post-init writes do not alter the live agent population — they only update reaction rates from the new parameter namespace. To change live counts after initialize(), use add_species() / remove_species().

NfsimSession refuses a compartment model (GH #1067). Once that refusal is lifted, note that neither kind of write reaches a compartment’s size: BioNetGen writes the size into the XML as a number (issue #1061).

Parameters:

name (str) – Parameter name.
value (float) – Parameter value.

Raises:

ParameterError – If the parameter name is not recognized or if expression propagation fails for some dependent parameter.

Return type:

clear_param_overrides()[source]¶

Clear all parameter overrides set via set_param.

When called on a live session, parameters revert to their XML-time values and dependent expressions are re-evaluated so the System matches the original BNG XML namespace.

Return type:: None

evaluate(expr, overrides=None)[source]¶

Evaluate a BNG expression in the current parameter namespace.

Uses the same ExprTk evaluator that powers set_param propagation, so expressions can reference any BNG parameter, BNG built-in (_pi, _NA, …), or BNG-style construct (if(cond, then, else), &&, ||, mratio).

Parameters:

expr (str) – Expression string. May reference any BNG parameter name.
overrides (dict[str, float] | None) – Extra name → value bindings layered on top of the simulator’s persistent set_param overrides for this single evaluation. Useful for “what if X were Y?” probes without mutating session state.

Returns:

Numeric value of the expression.

Return type:

Raises:

ParameterError – If the expression fails to compile.

initialize(seed=None)[source]¶

Initialize the NFsim session.

Parses the XML, seeds the RNG, and prepares the simulation system. Must be called before simulate().

Parameters:

seed (int | None) – Random number generator seed. When omitted (or None), bngsim draws a fresh seed from system entropy so consecutive session lifetimes produce independent trajectories. Pass an explicit integer for reproducibility, from 0 to 2**32 - 1: NFsim seeds its generators from 32 bits (issue #1515). The actual seed used is exposed via the seed property and stamped onto every Result returned by simulate().

Raises:

ValueError – If seed is negative or past 2**32 - 1. A live session is left as it was.
SimulationError – If initialization fails.

Return type:

destroy()[source]¶

Destroy the session and release C++ resources.

Safe to call multiple times. Called automatically by __exit__ when using the context manager.

Return type:: None

simulate(t_start=None, t_end=None, n_points=None, *, timeout=None, relative_time=False, sample_times=None)[source]¶

Run a simulation segment and return a Result.

Session state is preserved across calls, enabling multi-step protocols (simulate, mutate, simulate again).

Parameters:

t_start (float | None) – Segment start time. Required unless sample_times is given.
t_end (float | None) – Segment end time. Required unless sample_times is given.
n_points (int | None) – Number of output time points (including endpoints). Required unless sample_times is given.
sample_times (list[float] | None) – Explicit output instants (GH #184). When provided, overrides t_start / t_end / n_points: the returned Result.time array equals these times exactly (sorted ascending, treated as absolute), instead of a uniform grid. The live NFsim clock advances from the current session time by the span of the list, so this works mid-protocol (after setConcentration / a prior segment), continuing the same trajectory. Must contain at least two finite points. Mutually exclusive with relative_time (sample_times are absolute).
timeout (float | None) – Wall-clock budget in seconds. None or 0 disables the budget. Positive values arm a check between each stepTo() output point; on overrun this method raises bngsim.SimulationTimeout and the session must be destroyed (or the context manager exited) before reuse — the live NFsim System has advanced an indeterminate distance into the segment and the session clock is not updated.
relative_time (bool) – When True, the returned Result.time array is offset so it begins at 0.0 and ends at t_end - t_start — the convention used by BNG2.pl’s simulate({method=>"nf", ...}) action, which prints “NFsim timepoints are reported as time elapsed since t_start=$t_start” at runtime. The internal NFsim clock and session_logical_time are unaffected, so multi-segment threading (simulate, mutate, simulate again) keeps the same absolute physical time it would have under the default. Recommended for hosts that thread NF protocols across multiple bngsim simulate() calls and want their stitched time axis to match BNG2.pl’s. Default False preserves the absolute-time labelling used by every other bngsim backend.

Returns:

Simulation result for this segment.

Return type:

Raises:

SimulationError – If the session is not initialized or simulation fails.
SimulationTimeout – If timeout is set and the wall-clock budget is exceeded.

step_to(time)[source]¶

Advance the simulation to time without recording output.

Useful for fast-forwarding before a protocol action.

Parameters:: time (float) – Target time to advance to.
Raises:: SimulationError – If the session is not initialized.
Return type:: None

get_molecule_count(mol_type)[source]¶

Get the current count of a molecule type.

Parameters:: mol_type (str) – Molecule type name (e.g. "L", "R").
Returns:: Current molecule count.
Return type:: int

add_molecules(mol_type, count)[source]¶

Add molecules of a given type in their default (unbound) state.

Parameters:

mol_type (str) – Molecule type name.
count (int) – Number of molecules to add. Must be > 0. A fractional value is rounded to the nearest integer (round-half-up, GH #51) so warm / dosing / scan re-seeds match bngsim SSA and the cold-start seed.

Raises:

ValueError – If count is not positive.

Return type:

get_species_count(pattern)[source]¶

Get the current count of an exact BNGL species pattern.

The current implementation supports exact single-molecule species patterns such as "X(p~0,y)" or "TNF()". Mutation-oriented species patterns must list every component and specify every stateful component state so that the live-session operation is unambiguous. Multi-molecule complex patterns fail clearly.

Symmetric components (e.g. L(r~u~c~g, r~u~c~g) declares two sites named r) are matched by sorted multiset of state/bond constraints across class members, so heterogeneous patterns like L(r~u,r~c) and L(r~c,r~u) resolve to the same species count.

Parameters:: pattern (str) – BNGL species pattern.
Returns:: Current exact species count.
Return type:: int

add_species(pattern, count)[source]¶

Add exact BNGL species instances to the live NFsim session.

Parameters:

pattern (str) – Exact single-molecule BNGL species pattern.
count (int) – Number of instances to add. Must be > 0. A fractional value is rounded to the nearest integer (round-half-up, GH #51).

Return type:

remove_species(pattern, count)[source]¶

Remove exact BNGL species instances from the live NFsim session.

Parameters:

pattern (str) – Exact single-molecule BNGL species pattern.
count (int) – Number of instances to remove. Must be > 0. A fractional value is rounded to the nearest integer (round-half-up, GH #51).

Return type:

set_species_count(pattern, count)[source]¶

Set the live count for an exact BNGL species pattern.

This is implemented as a diff against get_species_count(pattern) and may add or remove exact species instances. A fractional count is rounded to the nearest integer (round-half-up, GH #51).

Parameters:

pattern (str)
count (int)

Return type:

get_parameter(name)[source]¶

Evaluate and return a parameter value from the live session.

Can evaluate expressions, not just named parameters.

Parameters:: name (str) – Parameter or expression name.
Returns:: Evaluated value.
Return type:: float

get_observable_names()[source]¶

Return observable names from the live session.

Returns:: Observable names defined in the BNG XML.
Return type:: list[str]

get_observable_values()[source]¶

Return current observable values from the live session.

Returns:: Current observable values.
Return type:: list[float]

save_species(path)[source]¶

Write the live System’s molecular species to a BNG .species file.

Mirrors the artifact BNG2.pl emits when running simulate({method=>"nf", get_final_state=>1, ...}). The file is a plain-text listing of species patterns and counts, formatted exactly as NFsim’s standalone CLI writes it. A host like PyBioNetGen that drives multi-segment NF protocols can read this file to thread state across segments without needing in-process snapshot support (see issue #52 for the in-process path).

Parameters:: path (str | Path) – Output file path. Overwritten if it exists.
Raises:: SimulationError – If the session is not initialized or the file cannot be written.
Return type:: None

save_concentrations(label=None)[source]¶

Snapshot the live session’s full molecular state for later restore.

Captures every molecule’s count, component states, and bonds into an in-memory snapshot held by the session. A later restore_concentrations() call rewinds the session to exactly this state. Mirrors the BNG saveConcentrations() / saveConcentrations("name") action.

Unlike save_species() (which writes a .species file for out-of-process state threading), this keeps the state in process — useful for equilibrate → snapshot → perturb → restore workflows without touching disk.

Parameters:

label (str | None) –

Name for the snapshot, so a multi-phase protocol can refer to it by name in restore_concentrations() (BNG saveConcentrations("name")). Mirrors bngsim.Model.save_concentrations(). When omitted (or None), the snapshot goes to the default (unlabeled) slot.

Each distinct label owns its own snapshot, so multiple named NFsim states coexist and round-trip faithfully (issue #11): a later save_concentrations("other") does not disturb a snapshot saved under a different name, and restore_concentrations(label) rewinds to exactly the state saved under that name. Saving to the same label again overwrites just that slot.

Raises:

SimulationError – If the session is not initialized.

Return type:

restore_concentrations(label=None)[source]¶

Restore the molecular state captured by save_concentrations().

Rewinds the live session to the snapshotted state (counts, component states, and bonds), discarding any simulation progress since the snapshot. Mirrors the BNG resetConcentrations() / resetConcentrations("name") action.

Parameters:: label (str | None) – Name of the snapshot to restore. When omitted (or None), restores the default (unlabeled) slot saved by save_concentrations() with no label. When given, restores the state saved under that exact name by save_concentrations(label); named and default slots are independent (issue #11).
Raises:: SimulationError – If the session is not initialized, or no snapshot has been saved under the requested label (or the default slot when label is omitted).
Return type:: None

has_saved_concentrations(label=None)[source]¶

Whether a save_concentrations() snapshot is available to restore.

Parameters:: label (str | None) – When given, reports whether a snapshot saved under that exact name exists. When omitted, reports whether any snapshot is held (a named slot or the default/unlabeled one).
Return type:: bool

property saved_concentration_labels: list[str]¶

Sorted names of the currently saved named concentration snapshots.

Does not include the default (unlabeled) slot, which is saved by save_concentrations() with no label and restored via restore_concentrations() with no label. Mirrors bngsim.Model.saved_concentration_labels. Empty before initialize().

property initialized: bool¶: Whether initialize() has been called.

property seed: int | None¶: The integer RNG seed used by initialize(), or None if the session has not been initialized yet.

property destroyed: bool¶: Whether destroy() has been called.

property xml_path: str¶: Path to the BNG XML file.

class bngsim.RuleMonkeySession(xml_path, *, molecule_limit=None, block_same_complex_binding=True)[source]¶

Bases: object

Stateful RuleMonkey simulation session with per-step control.

RuleMonkey reads BNG XML and implements the nf_exact network-free backend. Use this class for multi-action workflows that need to initialize, simulate, add molecules, and simulate again.

The session clock is the model’s clock: a rate law that reads time(), or a table function indexed by time, reads absolute model time, as method="ode" and method="ssa" do (issue #1022). See initialize() for where it starts and simulate() for a segment that starts elsewhere.

Events are placed on time elapsed since the session’s start wherever the model allows, which keeps their waiting times from rounding away at a large start (issue #1063): RuleMonkey reads a copy of the XML whose clock is time() plus the start, as Simulator.run does, and every time this class takes or reports is model time. A model with a table function indexed by time runs on model time itself (see bngsim/_rulemonkey_clock.py).

block_same_complex_binding is as for Simulator. RuleMonkey 3.13.0 deprecates False: it prints a warning to stderr the first time in a process, and RuleMonkey 4.0 removes the mode.

Parameters:

xml_path (str | Path)
molecule_limit (int | None)
block_same_complex_binding (bool)

set_molecule_limit(limit)[source]¶

Set the global molecule limit before initialization.

Parameters:: limit (int)
Return type:: None

set_param(name, value)[source]¶

Set a parameter value before initialization.

Parameters:

name (str)
value (float)

Return type:

clear_param_overrides()[source]¶

Clear all parameter overrides set via set_param.

Return type:: None

initialize(seed=None, t_start=None)[source]¶

Initialize the RuleMonkey session.

Parameters:

seed (int | None) – RNG seed. None (default) draws a fresh seed; pass an integer from 0 to 2**64 - 1 for reproducibility (a seed outside that range raises ValueError). The actual seed used is exposed via the seed property and stamped onto every Result returned by simulate().
t_start (float | None) – For a model that reads the clock (time(), or a table function indexed by time), where the seed state sits on the model’s clock: current_time reads it, and a rate law that reads time() starts from its value there. None (default) starts the clock at 0, and lets the first simulate() move the seed state to its own start if nothing has moved the session first (step_to(), or a change to its molecules), so that initialize(seed) then simulate(5, 7, ...) runs the model from its initial state at t = 5, as Simulator.run(t_span=(5, 7)) does (issue #1022). A model that reads no clock cannot tell, and its session always starts at 0, which keeps event times small; a segment from t_start is labelled with its own times.

Return type:

destroy()[source]¶

Destroy the session and release C++ resources.

Return type:: None

simulate(t_start=None, t_end=None, n_points=None, *, timeout=None, sample_times=None)[source]¶

Run a simulation segment and return a Result.

A segment that starts at current_time runs on the session clock, the model’s. One that starts elsewhere continues from current_time by the same span and is labelled with the times asked for, which a model that does not read the clock cannot tell apart. For a model that reads the clock (time(), or a table function indexed by time) the labels would not be the times it read, so such a segment raises bngsim.SimulationError, unless initialize() left the start unset and nothing has moved the session since: the seed state then moves to the segment’s start.

The result’s solver_stats["n_steps"] counts the events since initialize(), not this segment’s alone, as an bngsim.NfsimSession’s does.

Parameters:

t_start (float | None) – Segment endpoints and number of output samples. Required unless sample_times is given.
t_end (float | None) – Segment endpoints and number of output samples. Required unless sample_times is given.
n_points (int | None) – Segment endpoints and number of output samples. Required unless sample_times is given.
sample_times (list[float] | None) – Explicit output instants (GH #184). When provided, overrides t_start / t_end / n_points: the returned Result.time array equals these times exactly (sorted ascending, treated as absolute), instead of a uniform grid. The live session advances from its current time by the span of the list, so this works mid-protocol (after a prior segment or setConcentration), continuing the same trajectory. Must contain at least two finite points.
timeout (float | None) – Wall-clock budget in seconds. None or 0 disables the budget. Positive values arm a check polled by upstream every ~1024 SSA events; on overrun this method raises bngsim.SimulationTimeout and the session should be destroyed (or the context manager exited) before reuse — the live RuleMonkey session sits at the last completed event and subsequent state is undefined for bngsim purposes.

Return type:

step_to(time, *, timeout=None)[source]¶

Advance the simulation to time without recording output.

Parameters:

time (float) – Target time to advance to.
timeout (float | None) – Wall-clock budget; see simulate() for semantics.

Return type:

get_molecule_count(mol_type)[source]¶

Get the current count of a molecule type.

Parameters:: mol_type (str)
Return type:: int

add_molecules(mol_type, count)[source]¶

Add molecules of a given type in their default unbound state.

A fractional count is rounded to the nearest integer (round-half-up, GH #51), matching bngsim SSA and the cold-start seed.

Parameters:

mol_type (str)
count (int)

Return type:

get_parameter(name)[source]¶

Return a parameter value from the loaded XML plus overrides.

Parameters:: name (str)
Return type:: float

get_observable_names()[source]¶

Return observable names from the loaded XML.

Return type:: list[str]

get_observable_values()[source]¶

Return current observable values from the live session.

Return type:: list[float]

get_species_count(pattern)[source]¶

Get the current count of an exact BNGL species pattern.

Peer of bngsim.NfsimSession.get_species_count(). pattern is parsed and canonicalized by RuleMonkey’s runtime species-pattern parser, so it need not be byte-identical to a label RuleMonkey emits. Scope is exact, fully-specified, connected species: every molecule lists every component, every stateful component carries a concrete ~state, and bonds are numeric labels. Wildcards (!+, !?) or omitted components raise.

Parameters:: pattern (str) – Exact BNGL species pattern.
Returns:: Current exact species count.
Return type:: int

add_species(pattern, count)[source]¶

Add exact BNGL species instances to the live RuleMonkey session.

Parameters:

pattern (str) – Exact, fully-specified, connected BNGL species pattern.
count (int) – Number of instances to add. Must be > 0. A fractional value is rounded to the nearest integer (round-half-up, GH #51).

Return type:

remove_species(pattern, count)[source]¶

Remove exact BNGL species instances from the live RuleMonkey session.

Parameters:

pattern (str) – Exact, fully-specified, connected BNGL species pattern.
count (int) – Number of instances to remove. Must be > 0. A fractional value is rounded to the nearest integer (round-half-up, GH #51). Raises if fewer than count copies are live.

Return type:

set_species_count(pattern, count)[source]¶

Set the live count for an exact BNGL species pattern.

Drives the live count of pattern to exactly count, adding or removing the difference. A fractional count is rounded to the nearest integer (round-half-up, GH #51). Mirrors bngsim.NfsimSession.set_species_count(); PyBioNetGen’s _apply_nfsim_concentration_changes probes this method via getattr to propagate setConcentration to rm runs.

Parameters:

pattern (str)
count (int)

Return type:

evaluate(expr, overrides=None)[source]¶

Evaluate a BNG expression against the live session’s current state.

Resolvable symbols: every declared parameter, the clock t / time(), every observable, and every global function — all settled against the current pool, exactly as a rate law would see them between events. overrides supplies extra name → value bindings layered on top for a single evaluation; an override name shadows a model symbol on a clash.

Note

Unlike bngsim.NfsimSession.evaluate() (which only requires a live simulator), the RuleMonkey engine resolves expressions against the active session’s pool, so this method requires an initialized session — call initialize() first.

Parameters:

expr (str) – Expression string.
overrides (dict[str, float] | None) – Extra name → value bindings for this evaluation only.

Returns:

Numeric value of the expression.

Return type:

Raises:

ParameterError – If the expression fails to compile or reference an unknown symbol.

save_species(path)[source]¶

Write the live session’s molecular species to a BNG .species file.

Exact peer of bngsim.NfsimSession.save_species(). Enumerates the live complex pool, deduplicates by graph isomorphism, and writes a BNG-format .species file (# comment header followed by one <pattern> <count> line per species) readable by BNG2.pl’s readNFspecies. This is the artifact PyBioNetGen’s simulate_nf hook reads to thread get_final_state continuation across saveConcentrations / resetConcentrations segments — binding it makes multi-segment method=>"rm" runs reproduce the subprocess trajectory with no PyBioNetGen change.

Parameters:: path (str | Path) – Output file path. Overwritten if it exists.
Raises:: SimulationError – If the session is not initialized or the file cannot be written.
Return type:: None

save_state(path)[source]¶

Save the active session state to a binary snapshot file.

Serializes the full molecular pool plus the RNG state, keyed by a fingerprint of the molecule-type schema. A later load_state() (on a session built from the same XML) resumes the trajectory exactly. This is RuleMonkey’s in-process state-continuation path; NFsim has no equivalent (here RuleMonkey is ahead).

For threading state into a BNG2.pl-driven host like PyBioNetGen, prefer save_species() — that goes through BNG’s text .species + readNFspecies mechanism, which is what the backend hook actually calls.

Caveats¶

The schema fingerprint covers molecule-type schema only (molecule types, components, allowed states) — not ReactionRule or Observable patterns. Loading into a simulator built from a different XML that keeps the molecule schema but changes rules / observables succeeds silently and continues with the new rules.
The RNG state is serialized via the C++ stdlib’s mt19937_64 stream operator, which is not specified to be byte-identical across stdlib implementations. Save/load is reliable only between binaries built against the same toolchain.

param path:: Output snapshot path. Overwritten if it exists.
type path:: str | Path
raises SimulationError:: If the session is not initialized or the file cannot be written.

Return type:: None
Parameters:: path (str | Path)

load_state(path)[source]¶

Load session state from a snapshot written by save_state().

Creates a new live session from the snapshot, replacing any existing one. The snapshot’s schema fingerprint must match this session’s XML (see save_state() caveats). On success the session is initialized and current_time reflects the snapshot’s logical time — feed it to simulate(t_start, ...) to resume sampling.

The RNG seed is not recoverable from a snapshot, so seed reads back None after load_state() even though the underlying RNG stream is restored.

Parameters:: path (str | Path) – Snapshot path written by save_state().
Raises:: SimulationError – If the snapshot cannot be read or its schema fingerprint does not match this session’s model.
Return type:: None

property initialized: bool¶

property seed: int | None¶: The integer RNG seed used by initialize(), or None if the session has not been initialized yet.

property destroyed: bool¶

property xml_path: str¶

property current_time: float¶

class bngsim.NamedArray(data, colnames)[source]¶

Bases: ndarray

A 2-D ndarray with named columns (RoadRunner-compatible).

Constructed via bngsim.Result.as_roadrunner(). End users rarely instantiate it directly; subclass primarily exists so that callers can identify the array shape (arr.colnames) and look up columns by name (arr["[X]"], arr[:, "[X]"]).

Variables:

colnames (list[str]) –

One name per column, in column order — the invariant __new__() checks, and one that holds for the lifetime of every array of this class (issue #561). Names are positional, so they are carried exactly where this class can see the whole operation and knows what became of the columns: construction; indexing, which relabels whatever survives the key (arr[:, 1:], arr[:, ["time", "[X]"]], arr[:, ::-1], a row slice); and a pickle round trip, so an array that crosses a process boundary arrives labeled (issue #629).

Any other array numpy hands back — a product, a transpose, a reduction, arr.copy(), copy.deepcopy(), np.copy(arr, subok=True), astype, view — has colnames == [] and raises KeyError on a lookup by name, rather than answering with whichever column now sits at the name’s old position. A copy is in that list although its columns are plainly the parent’s: .copy() is one of numpy’s own primitives (np.sort(a, axis=1) is a copy followed by an in-place sort), so a copy that carried labels would hand them to operations that permute the columns underneath them. The comment on the method that is not overridden has the measurement.

libroadrunner keeps the labels through pickle as this class does, and drops them on every derived array, copies and slices included. Where this class keeps more of them, it is additive — code written against RoadRunner cannot depend on a label being absent — and it keeps one only where the column is established to be the same column.

Parameters:

data (NDArray[np.float64])
colnames (list[str])

Return type:

NamedArray

Examples

>>> result = sim.run(t_span=(0, 10), n_points=11)
>>> arr = result.as_roadrunner()
>>> arr.colnames
['time', '[X]', '[Y]']
>>> arr["[X]"]                      # 1-D column view
array([ ... ])
>>> arr[:, "[X]"]                   # equivalent
array([ ... ])
>>> arr[:, 1:].colnames             # the columns the slice kept
['[X]', '[Y]']

Parameters:

data (NDArray[double])
colnames (list[str])

Return type:

NamedArray

class bngsim.JacobianMatrix(data, source)[source]¶

Bases: ndarray

A dense (n_species, n_species) Jacobian that knows how it was built.

Returned by bngsim.Model.jacobian() with sparse=False. An ordinary float64 numpy.ndarray in every respect — J @ v, numpy.linalg.eigvals(J), slicing — with one extra attribute, below.

J[i, j] is ∂f_i/∂x_j — row i is the species whose derivative is being differentiated, column j the species it is differentiated with respect to, both in bngsim.Model.species_names order.

Variables:

source (str) – "analytical" when every reaction’s derivative came from the model’s closed form; "automatic-differentiation" when it has none and the matrix is exact by dual numbers through the interpreter (issue #1330); "finite-difference" when the matrix is the one-sided difference quotient of bngsim.Model.rhs(), for a model with neither. A model whose analytical Jacobian is incomplete never hands back the partial matrix under the first label.

Parameters:

data (NDArray[double])
source (str)

Return type:

JacobianMatrix

source: str¶

exception bngsim.BngsimError[source]¶

Bases: RuntimeError

Base exception for all bngsim errors.

exception bngsim.CodegenDeclinedError[source]¶

Bases: BngsimError

codegen=True on a model codegen declines, such as one whose functions read each other (issue #621): a statement about the model, which runs on the interpreted RHS with codegen=None or codegen=False (issue #1313). The deprecated bngsim.prepare_codegen() raises it for a decline too.

Distinct from a codegen build failure (a compiler error, a timeout), which says nothing about the model and stays a plain RuntimeError. A RuntimeError (through BngsimError), as it was raised before it had a type, so an existing except RuntimeError still catches it.

Unrelated to bngsim._codegen.CodegenDeclined, the emitter’s private signal to decline, which the prepare calls turn into a None return; this is what a caller who asked for compiled code is told.

exception bngsim.CodegenToolchainError[source]¶

Bases: BngsimError

The C toolchain cannot build a codegen library for any model (issue #1544).

Raised when a codegen compile has failed and the same command fails on a trivial source too, or when the compiler is not there at all: the failure is the environment’s (a linker that cannot read the SDK’s libSystem.tbd, a missing compiler), not the model’s. The message names the compiler, the linker and, on macOS, the SDK, and suggests an installed SDK that links when one does.

The verdict holds while the compiler and its environment stay as they are (it is probed again every 5 minutes), so a later build is refused without compiling, and before its source is generated (a model codegen would have declined is refused with this too); a library already in the cache is still used. codegen=True, forward sensitivities and, under BNGSIM_CODEGEN_STRICT, automatic codegen raise it, each in its own words. A RuntimeError (through BngsimError), as a build failure was raised before it had a type.

exception bngsim.ConversionError[source]¶

Bases: BngsimError

A format conversion (e.g. SBML→.net) cannot be completed faithfully.

Raised by bngsim.convert.sbml_to_net() when the source model uses a construct that the target format cannot represent without changing the model’s meaning — for example an hasOnlySubstanceUnits species in a compartment whose volume is not 1, a cross-compartment reaction that needs per-species volume scaling, a rate-rule ODE, a Michaelis–Menten rate-law type, or a table (interpolation) function. The message names the offending construct and the model so the cause is actionable.

Pass strict=False (API) / --allow-lossy (CLI) to downgrade these to a ConversionWarning and emit a best-effort .net anyway.

exception bngsim.ConversionWarning[source]¶

Bases: UserWarning

A format conversion dropped or approximated part of the source model.

Emitted when the network channel is preserved but something outside it was left behind — most commonly SBML <event> elements, which belong to the simulation-protocol channel (a SED-ML sidecar), not the .net network. Also used when strict=False downgrades an otherwise-fatal ConversionError. Filter or escalate it like any UserWarning:

import warnings, bngsim
warnings.simplefilter("error", bngsim.ConversionWarning)   # promote
warnings.simplefilter("ignore", bngsim.ConversionWarning)  # silence

exception bngsim.ModelError[source]¶

Bases: BngsimError

Error loading or manipulating a model.

Raised when: - A .net file cannot be parsed - A model is in an invalid state for simulation - A reserved name conflict is detected

exception bngsim.ParameterError[source]¶

Bases: BngsimError

Error with parameter access.

Raised when: - A parameter name is not found in the model - A parameter value is invalid (NaN, Inf, wrong type)

exception bngsim.SimulationError[source]¶

Bases: BngsimError

Error during simulation.

Raised when: - CVODE fails to converge - NaN/Inf detected in species concentrations - Maximum number of steps exceeded

exception bngsim.SimulationTimeout(message, *, timeout=0.0, elapsed=0.0, partial_result=None)[source]¶

Bases: BngsimError

Raised when a simulation exceeds its wall-clock budget.

Distinct from SimulationError so callers (e.g. PyBNF’s wall_time_sim) can classify wall-clock terminations separately from solver/convergence failures. Inherits from BngsimError (and therefore RuntimeError).

Variables:

timeout (float) – The configured wall-clock budget, in seconds.
elapsed (float) – Actual elapsed wall-clock time at the point the timeout fired, in seconds. elapsed >= timeout always holds.
partial_result (Result | None) – Reserved for future use. Currently always None — bngsim does not yet salvage a partial Result from a timed-out run.

Parameters:

message (str)
timeout (float)
elapsed (float)
partial_result (object)

Return type:

None

exception bngsim.SparseJacobianUnavailableError[source]¶

Bases: SimulationError

force_sparse_linear_solver=True on a model that cannot supply a sparse Jacobian for KLU to factorize: its analytical Jacobian is unavailable (or suppressed by jacobian="fd") and its sparsity pattern has no structural nonzero to colour for finite differences (issue #1313).

A refusal by name, not a solver failure: drop the flag to use the dense solver. A SimulationError, as it was raised before it had a type, so an existing except SimulationError still catches it.

exception bngsim.SsaBoundaryWarning[source]¶

Bases: UserWarning

An SSA run hit a literal-rate-law boundary condition (GH #110).

Emitted once per run when the exact SSA either drove a species count negative (no non-negativity floor is applied — non-negativity is the rate law’s responsibility, matching the CVODE path) or fired a reaction in reverse because its rate law evaluated negative. Both keep the SSA mean consistent with the ODE; the warning surfaces what would otherwise be a silent boundary behavior. Filter or escalate it like any UserWarning:

import warnings, bngsim
warnings.simplefilter("error", bngsim.SsaBoundaryWarning)  # promote
warnings.simplefilter("ignore", bngsim.SsaBoundaryWarning)  # silence

The structured counts are always available on result.ssa_diagnostics regardless of warning filters.

exception bngsim.SsaRoundingWarning[source]¶

Bases: UserWarning

An SSA run rounded a fractional molecule count to a whole one (issue #718).

The exact SSA fires whole molecules, so a species whose count is not a whole number when a run starts is rounded to the nearest one, as run_network does. This is emitted for every run and every run_until leg that rounds, not only when the Simulator is built: a count set between legs, or left fractional by an event assignment, is rounded at the next leg. A parameter or compartment that an event assigns and a rate-rule target hold continuous values and are never rounded. result.ssa_diagnostics carries the count (n_rounded_populations) regardless of filters:

import warnings, bngsim
warnings.simplefilter("error", bngsim.SsaRoundingWarning)  # promote

exception bngsim.DenseSolverFallbackWarning[source]¶

Bases: UserWarning

A large ODE model is running on the dense solver for lack of KLU (GH #209).

Emitted once per process when Simulator.run() integrates a large model (n_species past the warn threshold) on the dense linear solver only because this bngsim install was built without SuiteSparse/KLU — not because the user asked for force_dense_linear_solver or jacobian="jax" (both legitimately dense). Dense LU factorizes the full N×N Jacobian at O(N³); for a sparse genome-scale network that is the difference between minutes and hours. The fix is to rebuild bngsim with KLU (see bngsim.capabilities() and GH #209). Filter or escalate it like any UserWarning:

import warnings, bngsim
warnings.simplefilter("ignore", bngsim.DenseSolverFallbackWarning)  # silence
warnings.simplefilter("error", bngsim.DenseSolverFallbackWarning)  # promote

exception bngsim.AtolScaleWarning[source]¶

Bases: UserWarning

An ODE run’s default absolute tolerance is above species that matter (issue #1361).

BNG2.pl’s scalar atol of 1e-8 was chosen for species held as molecule counts. A model in molar units (an SBML model of species at 1e-12 to 1e-6 M, say) holds species below it, and CVODE then controls none of their error: the binding example of the issue returned a negative concentration, and BIOMD0000000002, whose ligand is above it and whose complexes are below it, errors of 1e-2 against its reference. The default now takes each species’ scale from the state the run starts from (issue #1420, see Simulator.run()), which fixes both. What it cannot scale is a species that starts at 0, or below every scale the start gives, in a model where no species starts small: one whose species all start at 0 runs at the scalar.

Emitted by Simulator.run() (and so run_until), run_batch, parameter_scan and compute_all_sensitivities when the caller gave no atol and a species the solver integrates stays below its tolerance all run while reaching at least 1e-6 of the run’s largest such species. A $-fixed species or an assignment-rule target is not integrated and is never named. atol="tracking" (a tolerance that follows each species’ magnitude, issue #213) is the remedy it names. Filter or escalate it like any UserWarning:

import warnings, bngsim
warnings.simplefilter("ignore", bngsim.AtolScaleWarning)  # silence
warnings.simplefilter("error", bngsim.AtolScaleWarning)  # promote

exception bngsim.RuleMonkeyWarning[source]¶

Bases: UserWarning

RuleMonkey runs the model but reports how it reads part of it (issue #1022).

Emitted when method="rm" ("nf_exact") or a bngsim.RuleMonkeySession loads a BNG XML for which RuleMonkey reports a warning-level diagnostic: a rule whose products can never be what it writes (BioNetGen generates no reaction for it, and RuleMonkey never fires it), a Michaelis-Menten law over a pattern that matches several species (BNG2 expands it per species pair, RuleMonkey sums the counts), an observable read at local and global scope in one function, a rule priority, TotalRate, and the like. Each message names the rule or construct and says what RuleMonkey does with it. Errors are refused instead. Filter or escalate it like any UserWarning:

import warnings, bngsim
warnings.simplefilter("error", bngsim.RuleMonkeyWarning)  # promote
warnings.simplefilter("ignore", bngsim.RuleMonkeyWarning)  # silence

exception bngsim.SensitivityRefusedDuringRunError[source]¶

Bases: SensitivityUnsupportedError, SimulationError

A forward-sensitivity run the solver refused partway, by declaration.

The solver reached a construct whose derivative it cannot produce. Raised as this type itself, by a condition no choice of requested parameters changes: an event that fires as its trigger leaves its threshold with no first-order motion (issue #934), a compartment that jumps twice at one instant (issue #1181), a delay a run_batch row made positive (issue #144), and the like. The message names it and what to change. A SensitivityUnsupportedError, as the refusals raised before a run are, so a fitting driver or a parity suite tells it from a solver failure by its type (issue #1279). And a SimulationError, which all but one of these were raised as before they had a type, so an existing except SimulationError still catches it.

A refusal whose condition turns on a requested column (a crossing or a fire the column moves, two times it parts) is its subclass SensitivityRefusedForParametersError, because Simulator.compute_all_sensitivities() retries a failed chunk of those one column at a time and a column alone can run, where it raises this type as it is. The one exception is the refusal of a held equality the state leaves (issue #1279), which had this type first: where a restart lands on the stop, a column that moves it is refused and one that does not runs.

exception bngsim.SensitivityRefusedForParametersError[source]¶

Bases: SensitivityRefusedDuringRunError

A forward-sensitivity run the solver refused partway, for the parameters requested.

The condition can turn on which sensitivities were asked for: a crossing or a fire a requested parameter moves, two times it moves apart, or an edge it moves the run onto, where the result has a kink or a jump and no derivative (issues #49, #150, #760, #767, #945, #1025 and others). The message names it and, where it can, the parameters to drop. A smaller set may run, so Simulator.compute_all_sensitivities() retries a chunk refused this way one column at a time (GH #243) and raises this type, naming the parameters refused alone and the ones that ran, when some are refused alone too.

A SensitivityRefusedDuringRunError, so a SensitivityUnsupportedError and a SimulationError: a fitting driver or a parity suite tells it from a solver failure by its type, and an except SimulationError that caught these before they had one still does (issue #1435).

exception bngsim.SensitivityUnsupportedError[source]¶

Bases: BngsimError, ValueError

Forward sensitivity is unavailable for this model, by declaration.

A clean refusal, not a failure. bngsim inspected the model, found a construct whose derivative it cannot produce correctly, and declined rather than return silently-wrong numbers. Two constructs raise it today:

an event whose crossing time moves with a requested parameter and whose dt*/dp bngsim cannot compute, so every post-jump sensitivity would be stale (GH #205, remaining subclasses tracked in issue #144); and
a rate law codegen could not differentiate to closed form — a non-smooth min/max/abs/floor, rateOf() inside a rate law, or an unparseable expression (GH #214). bngsim refuses rather than fall back to finite differences.

The solver’s refusals partway through a run are its subclass SensitivityRefusedDuringRunError, which is a SimulationError too.

This is the forward-sensitivity peer of SsaValidationError: both say “this model exercises a feature this regime does not support”, and both are distinguishable by type, so a caller can separate a declared non-run from an actionable bug without pattern-matching on message text. A fitting driver can catch it to fall back on a derivative-free optimizer; the parity suites bucket it as UNSUPPORTED (non-scoring) rather than EXCEPTION.

Also inherits ValueError, which both sites raised before this class existed, so existing except ValueError handlers keep working.

Environment failures are deliberately NOT this exception: “no C compiler and no JIT available” is a fixable local problem, not a property of the model, and stays a plain RuntimeError.

exception bngsim.SsaValidationError(issues, *, override_attempted=False)[source]¶

Bases: BngsimError

An SBML model contains constructs that cannot be simulated under SSA.

Raised by bngsim.Simulator when method="ssa" is requested on a model whose loader-captured SSA issues include any with severity="error". The full issue list (errors and warnings) is available as issues.

The error message advertises strict_ssa=False as a workaround when at least one of the raised codes is overridable. override_attempted flips the message to call out that the user did pass strict_ssa=False but the listed codes are non-overridable.

Parameters:

issues (list[SsaIssue])
override_attempted (bool)

Return type:

None

NON_OVERRIDABLE_CODES = frozenset({'fast_reaction', 'non_integer_stoichiometry'})¶

exception bngsim.StopConditionMet(message, *, result=None, condition='')[source]¶

Bases: BngsimError

A stop condition was triggered during simulation.

The partial result up to the trigger point is attached as self.result.

Conditions are checked against a completed result, so this does not interrupt the solve: the backend integrates the whole requested interval and leaves the model at its end, and self.result is that trajectory truncated at the first row where the condition held. Simulator’s current_time therefore reads the end of the interval, not the trigger time — the two differ on purpose, and the clock is the one that names the state the model holds (issue #553). Take a snapshot() beforehand to continue from the trigger point instead.

Variables:

result (Result) – Simulation result truncated at the stop point. The reported species block only: Result.state is not available on it, so it cannot be used to restore the model’s integrator state.
condition (str) – Description of the condition that triggered.

Parameters:

message (str)
result (object)
condition (str)

Return type:

None

exception bngsim.UnderSpecifiedModelError[source]¶

Bases: ModelError

The model reads a symbol it never gives a value.

Raised at load when a <parameter> has no value attribute and no <initialAssignment>, yet is referenced by a kineticLaw / rule / event / initialAssignment expression, or when a kineticLaw’s local parameter has no value and the law reads it (issue #580). SBML leaves such a value undefined, so there is nothing to integrate; bngsim refuses rather than default it to 0.0, which would silently zero out every expression that consumes it and return a wrong-but-plausible trajectory. RoadRunner refuses the global case too (it runs a Level 2 local one at 0).

That failure mode is not hypothetical. MODEL1006230032 under BNGSIM_ALLOW_UNSET_PARAMS=1 integrates to all five species identically zero for the whole run, and every cross-engine check — state, sensitivity, both engines — agrees, because AMICI defaults the same symbol the same way. A model computing nothing, with every verdict green.

Typed (issue #323) so a caller can tell “this model is under-specified” from the other things ModelError covers — a .net parse failure, an invalid model state. Parity suites bucket it UNSUPPORTED: a clean, documented refusal, not an actionable bngsim bug. Set the parameter’s value (a global one can take an <initialAssignment> instead), or opt into the legacy lenient behavior with BNGSIM_ALLOW_UNSET_PARAMS=1 (which downgrades this to a warning).

class bngsim.SsaIssue(severity, code, message, location=None)[source]¶

Bases: object

A single SSA-compatibility finding from an SBML model.

Variables:

severity (Literal['error', 'warning']) – "error" (Simulator init will raise) or "warning" (logged).
code (str) – Stable machine-readable identifier for the issue category.
message (str) – Human-readable description, ideally pointing the user at a fix.
location (str | None) – Optional locator string (e.g. "reaction:R1", "compartment:cell", "species:S").

Parameters:

severity (Literal['error', 'warning'])
code (str)
message (str)
location (str | None)

severity: Literal['error', 'warning']¶

code: str¶

message: str¶

location: str | None = None¶

bngsim.validate_for_ssa(model)[source]¶

Return SSA-compatibility issues recorded for model.

The list is populated by the SBML loader at load time, with an additional model-agnostic warning when the current species state contains non-integer molecule populations. Fractional populations are rounded by the C++ SSA setup path before simulation.

PyBNF and other consumers can call this before constructing a Simulator to pre-flight a model — the same list will be re-checked by Simulator(model, method="ssa").

Parameters:: model (Model)
Return type:: list[SsaIssue]

bngsim.reserved_names()[source]¶

Return dict of reserved constant and function names.

Returns:: {"constants": [...], "functions": [...]}
Return type:: dict[str, list[str]]

Example

>>> import bngsim
>>> names = bngsim.reserved_names()
>>> "_pi" in names["constants"]
True
>>> "time" in names["functions"]
True

bngsim.configure_logging(level=20, *, handler=None, fmt='%(asctime)s [bngsim] %(levelname)s %(message)s')[source]¶

Configure the bngsim logger.

By default, bngsim is silent (no handler attached). Call this function to enable log output.

Parameters:

level (int) – Logging level (e.g. logging.DEBUG, logging.INFO).
handler (Handler | None) – Custom handler. Default: StreamHandler to stderr.
fmt (str) – Log message format string.

Returns:

The configured bngsim logger.

Return type:

Logger

Examples

>>> import bngsim, logging
>>> bngsim.configure_logging(logging.DEBUG)
>>> # Now all bngsim operations produce log output

bngsim.normalize_method(requested)[source]¶

Normalize a user-requested method token to its canonical form.

Parameters:: requested (str) – The method string as provided by the user.
Returns:: (canonical, dispatch) where canonical is the normalized algorithm name (e.g. "nf_reject") and dispatch is the internal backend key used for simulator creation (e.g. "nfsim").
Return type:: tuple[str, str]
Raises:: ValueError – If the method token is not recognized at all, or if it maps to an unavailable backend.

bngsim.derive_atol(state, rtol, *, floor=None)[source]¶

Build a per-species absolute tolerance from a state vector.

The rule is atol[i] = rtol * max(|y[i]|, floor): every species is resolved to rtol of its own magnitude, which is what a scalar cannot say and what every model spanning decades needs said.

This derives from the state you hand it, which is the difference that matters between this function and bngsim.Simulator.auto_atol() (and atol="auto", which is the same thing): those read the model’s live state, the one the next run() would start from. Reach for this one when the tolerance has to be a constant of the model rather than of the point being integrated — a parameter fit that moves initial conditions is the case. There, derive once from the nominal state, hold the vector for the whole fit, and pass it as atol=; "auto" would re-derive at every evaluation and put a step in the objective wherever the derivation crossed a rounding boundary, which is invisible in the usual way (the objective still looks right and only the search behaves oddly).

floor exists for exactly one case — a species whose value here is zero has no magnitude to scale, and atol[i] = 0 would put it under pure relative error control, which CVODE will not integrate the moment that species is genuinely zero. The default takes the smallest strictly positive entry of state, but not above 1: a species with no scale of its own is treated as living at the smallest scale the model actually exhibits, or at one unit (one molecule, in a model of counts) where every species that has a scale is larger. That is the tight choice, deliberately — it can cost step size on a zero-initialized species that later grows, but it never quietly drops one below the noise floor, and it is the direction a caller can measure and then override. Without the cap, a model whose every start is large held its species at 0 to that scale: MODEL2009230001’s, beside one at 1.1e10, to an absolute 110, and a species it grew from 0 came out 100% off.

Parameters:

state (Sequence[float] | NDArray[double]) – Species values to scale from, ordered like bngsim.Model.species_names. Normally the model’s initial concentrations.
rtol (float) – Relative tolerance the run will use. Must be > 0.
floor (float | None) – Magnitude to substitute for a species whose |y[i]| is below it. Must be > 0. Defaults to the smallest strictly positive |y[i]|, but not above 1.0, which is also the floor when every entry is zero.

Returns:

float64 array of len(state) absolute tolerances.

Return type: