API reference¶

`bngsim.Model`¶

Factory methods:

Model.from_net(path) — Load from BNG .net file
Model.from_antimony(path) — Load from Antimony .ant file (requires antimony, python-libsbml)
Model.from_antimony_string(text) — Load from Antimony string
Model.from_sbml(path) — Load from SBML .xml file (requires python-libsbml)
Model.from_sbml_string(text) — Load from SBML XML string

Properties:

n_species, n_reactions, n_observables, n_parameters, n_functions, n_events
species_names, param_names, observable_names

Methods:

set_param(name, value) — Set a parameter value
get_param(name) — Get a parameter value
set_params(dict) — Set multiple parameters (atomic)
reset() — Reset species to initial concentrations
clone() — Deep copy for parallel workers
save_concentrations() — Snapshot current state as new initial conditions
set_concentration(name, value) — Set a single species concentration
get_concentration(name) — Get a single species concentration
add_table_function(name, *, file, times, values, index) — Add a piecewise-linear table function
n_table_functions, table_function_names — TFUN introspection

`bngsim.Simulator`¶

Constructor:

Simulator(model, method="ode", *, jacobian, codegen, net_path, sensitivity_params, poplevel, xml_path, gml)

method: "ode", "ssa", "psa", "nf", "nf_reject", "nfsim", "nf_exact", "rulemonkey", "rm" (and aliases)

ODE-specific kwargs: jacobian ("auto", "analytical", "fd"), codegen (bool), net_path (BioNetGen .net path only), sensitivity_params (list[str])

For SBML and Antimony models, use codegen=True without net_path. net_path is not a generic model path and should not point to SBML XML.

PSA: poplevel (float, required). Accepts .net, SBML, and Antimony models (same dispatch as SSA, sharing the validate_for_ssa gate). NFsim/RuleMonkey: xml_path (str, required), gml (int).

Simulation:

run(t_span, n_points, *, seed, rtol, atol, max_steps, timeout, steady_state, steady_state_tol) → Result
run_batch(t_span, n_points, *, params, seed, num_processors, squeeze, timeout, steady_state, steady_state_tol) → list[Result] or Result
compute_all_sensitivities(t_span, n_points, *, params, chunk_size, n_workers, rtol, atol, max_steps) → Result with full sensitivity tensor

timeout is a wall-clock budget in seconds (or None to disable, the default). When the budget is exceeded, run / run_batch raise bngsim.SimulationTimeout (a typed BngsimError) carrying timeout (the configured limit) and elapsed (actual wall-clock time at the trip). Supported on every backend (ODE, SSA, PSA, NFsim, RuleMonkey). NFsim and RuleMonkey poll the budget at a coarser granularity than ODE/SSA/PSA (between stepTo() output points for NFsim; every ~1024 SSA events for RuleMonkey, via its upstream cancellation hook), so the surfaced elapsed time may overshoot the budget by one such interval.

Interactive:

run_until(t, *, n_points, seed) → Result
intervene(params) — Change parameters mid-simulation
snapshot() → dict — Save state
restore(snapshot) — Restore state

Stop conditions:

add_stop_condition(condition, *, label) — str expression or callable
clear_stop_conditions()

Steady-state (method ∈ "newton" (default), "integration", "kinsol" alias):

steady_state(*, tol, max_time, method, rtol, atol, max_steps, sensitivity_params) → SteadyStateResult
steady_state_batch(params, *, tol, max_time, method, rtol, atol, max_steps, n_workers) → list[SteadyStateResult]

Configuration:

set_tolerances(rtol, atol) — ODE solver tolerances
set_max_steps(max_steps) — Max internal solver steps

Properties:

method, model, current_time

`bngsim.NfsimSession`¶

Stateful NFsim session API for multi-action workflows where parameters or live particle counts are changed between simulation segments.

Constructor:

NfsimSession(xml_path, *, molecule_limit) — Create an NFsim session from a BNGL-generated NFsim XML file

Session control:

initialize(seed) — Initialize the live NFsim system
simulate(t_start, t_end, *, n_points, timeout) → Result — timeout is a wall-clock budget in seconds (None disables, the default). Checked between stepTo() output points; on overrun raises bngsim.SimulationTimeout and the session must be destroyed before reuse.
destroy() — Release the live NFsim session

Parameters:

set_param(name, value) — Set a parameter before initialization
get_parameter(name) — Evaluate a parameter in the live session

Live counts:

get_molecule_count(molecule_type) — Count all live molecules of a type
add_molecules(molecule_type, count) — Add default unbound molecules
get_species_count(pattern) — Count an exact single-molecule BNGL species
add_species(pattern, count) — Add exact single-molecule species instances
remove_species(pattern, count) — Remove exact single-molecule species instances
set_species_count(pattern, count) — Set the exact species count by adding or removing instances

Species count mutation currently supports exact, unbound, single-molecule BNGL patterns such as "X(p~0,y)", "L(r)", and "TNF()". Patterns must list every component and specify every stateful component state. Multi-molecule complex patterns fail with a clear SimulationError.

`bngsim.RuleMonkeySession`¶

Stateful RuleMonkey session API for exact network-free workflows where parameters are configured before initialization and live particle counts can be changed between simulation segments.

Constructor:

RuleMonkeySession(xml_path, *, molecule_limit, block_same_complex_binding) — Create a RuleMonkey session from a BNGL-generated NFsim XML file

Session control:

initialize(seed) — Initialize the live RuleMonkey system
simulate(t_start, t_end, *, n_points, timeout) → Result — timeout is a wall-clock budget in seconds (None disables, the default). Polled by upstream every ~1024 SSA events; on overrun raises bngsim.SimulationTimeout and the session must be destroyed before reuse.
step_to(time, *, timeout) — Advance without recording output; honors the same timeout semantics as simulate(...).
destroy() — Release the live RuleMonkey session

Parameters and counts:

set_param(name, value) — Set a parameter before initialization
clear_param_overrides() — Clear parameter overrides before initialization
get_parameter(name) — Evaluate a parameter from XML plus overrides
get_molecule_count(molecule_type) — Count all live molecules of a type
add_molecules(molecule_type, count) — Add default unbound molecules
get_observable_names() — Return observable names from the XML
get_observable_values() — Return current observable values

`bngsim.Result`¶

Properties:

time — ndarray (n_times,)
species — ndarray (n_times, n_species) with named access
observables — ndarray (n_times, n_obs) with named access (e.g. result.observables["A_tot"])
expressions — ndarray (n_times, n_expr) with named access
sensitivities — ndarray (n_times, n_species, n_params) forward sensitivity tensor
has_sensitivities — bool
sensitivity_params — list[str] parameter names for sensitivity
dataframe — pandas.DataFrame (requires pandas)
xr — AMICI-style xarray accessor; result.xr.species, .observables, .expressions, .sensitivities, .sensitivities_ic return labeled xarray.DataArrays with time/state/observable/expression/parameter/ic_state coords (requires xarray)
solver_stats — dict with solver diagnostics
custom_attrs — dict for user metadata
n_times, n_species, n_observables, n_expressions
species_names, observable_names, expression_names

Methods:

resolve_outputs(selectors) → list[dict] — resolve typed output selectors (species:/observable:/expression:, aliases state:→species:, function:→expression:, with expression:foo()→expression:foo) to per-output metadata {"selector","kind","name","index","column_label"}. Bare names resolve only if unique across all three kinds; ambiguous/unresolved selectors raise with the candidates listed.
outputs(selectors) → ndarray (n_times, n_outputs) — value columns for the named selectors (one column per selector, in order)
fisher_information(sigma) → ndarray (n_params, n_params) — Fisher Information Matrix from sensitivity data
gradient(loss_fn) → ndarray (n_params,) — parameter gradient ∇_p L from sensitivity tensor and loss function
to_gdat(path) — Export observables as BNG .gdat
to_cdat(path) — Export species as BNG .cdat
to_csv(path, *, kind="observables"|"species", sep=",", include_time=True, header=True) — Plain delimited-text export (CSV / TSV); no # prefix; SBML/RoadRunner-friendly
to_xarray() → xarray.Dataset bundling species/observables/expressions/sensitivities/sensitivities_ic with shared time coord (requires xarray); custom_attrs + seed propagate to ds.attrs
save(path) — Save to HDF5 (requires h5py)
Result.load(path) — Load from HDF5 (class method)
Result.squeeze(results) — Stack list into 3D batch result

Exceptions¶

All inherit from bngsim.BngsimError (which inherits RuntimeError):

ModelError — .net parse failures, invalid model state
SimulationError — Solver failures (convergence, NaN)
SimulationTimeout — Wall-clock budget exceeded; .timeout and .elapsed attributes
ParameterError — Unknown parameter name, invalid value
StopConditionMet — Stop condition triggered; .result has partial data

Universal `.net` reader¶

bngsim.parse_net_file(path) → dict — Parse a .net file into an engine-agnostic Python dict with keys: parameters, species, observables, functions, reactions. Pure Python — no C++ extension needed for parsing.
bngsim.build_model_from_parsed(parsed) → Model — Build a BNGsim Model from the dict returned by parse_net_file(). Routes through ModelBuilder for full optimization (analytical Jacobian, conservation laws, etc.).

Utility functions¶

bngsim.reserved_names() → dict with "constants" and "functions" lists
bngsim.configure_logging(level) — Enable log output

Generated API (autodoc)¶

The curated summaries above cover the common surface. The full, always-in-sync reference below is generated directly from the package docstrings.

bngsim — Embeddable simulation engine for BioNetGen reaction networks.

Usage:

import bngsim

model = bngsim.Model.from_net("model.net")
sim = bngsim.Simulator(model, method="ode")

model.set_param("kf", 0.5)
result = sim.run(t_span=(0, 1000), n_points=1001)

result.time         # (1001,) ndarray
result.observables  # (1001, n_obs) ndarray
result.species      # (1001, n_species) ndarray

See the package README for installation, usage, and API overview.

class bngsim.Model(_core)[source]

Bases: object

A BioNetGen reaction network model.

A Model holds species, reactions, observables, parameters, and functions. It can be loaded from .net files and, via the factory methods below, from Antimony and SBML inputs.

Models are not thread-safe. For parallel workers, use clone() to create independent copies.

Parameters:: _core (NetworkModel) – Internal C++ model object. Users should not construct this directly; use the factory methods instead.

Examples

>>> model = bngsim.Model.from_net("model.net")
>>> model.n_species
5
>>> model.set_param("kf", 0.5)
>>> model.get_param("kf")
0.5
>>> model.set_params({"kf": 1.0, "kr": 0.1})

classmethod from_antimony(path)[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.

Requires: pip install antimony python-libsbml

Parameters:

path (str | Path) – Path to the .ant file.

Returns:

The loaded model.

Return type:

Model

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)[source]

Load a model from an Antimony string.

Parameters:: text (str) – Antimony model text.
Returns:: The loaded model.
Return type:: Model

classmethod from_sbml(path, *, defer_jacobian=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.

Returns:

The loaded model.

Return type:

Model

Raises:

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

classmethod from_sbml_string(text, *, defer_jacobian=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.

Returns:

The loaded model.

Return type:

Model

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:

Model

Raises:

ModelError – If the file cannot be parsed.
FileNotFoundError – If the file does not exist.

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 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.

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

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

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 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().

property primary_param_names: list[str]

List of parameter names that are not derived constant expressions.

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.

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

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

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

Add a table function (piecewise-linear interpolation of data).

The function is registered with the expression evaluator and can be referenced by name in rate law expressions.

Parameters:

name (str) – 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 a parameter or observable name for non-time-indexed table functions.
method (str) – Interpolation method: "linear" (default) or "step".

Raises:

ModelError – If the file cannot be read or data is invalid.
ValueError – If arguments are inconsistent (e.g., both file and times).

Return type:

None

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, codegen=None, net_path='', sensitivity_params=None, sensitivity_ic=None, sensitivity_method='staggered', strict_ssa=True)[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
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.
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) – Only used for method="nf" / "nf_reject" / "nfsim". NFsim -bscb: when True, two reactant patterns in a bimolecular rule cannot match molecules in the same complex. Default: True in bngsim — NFsim CLI defaults it off, but bngsim defaults it on for correctness on BLBR/aggregation models. Pass False to allow same-complex binding (BNG2.pl complex=>1). This governs only the binding policy; complex bookkeeping for Species-typed observable counting is enabled automatically when the model declares such an observable, independent of this flag.
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.
codegen (bool | None) – Only used for method="ode". When true, use a compiled C RHS. Models loaded from BioNetGen .net files use the .net codegen path. SBML, Antimony, and other already-built models use model-based codegen. For SBML/Antimony models, pass codegen=True without net_path.
net_path (str) – BioNetGen .net path for the .net codegen path. This is not a generic model path and should not point to SBML XML. Models loaded with Model.from_net() remember their source path, so most callers do not need to pass this manually.
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".
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".
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).
xml_path (str)
gml (int | None)
jacobian (str)

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")

Parameters:

xml_path (str)
gml (int | None)
jacobian (str)
model (Model)
method (str)
poplevel (float | None)
connectivity (bool | None)
nfsim_v1143_compat (bool)
block_same_complex_binding (bool)
traversal_limit (int | None)
force_dense_linear_solver (bool)
codegen (bool | None)
net_path (str)
sensitivity_params (list[str] | None)
sensitivity_ic (list[str] | None)
sensitivity_method (str)
strict_ssa (bool)

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

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 3 values. Values are sorted automatically.
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. The actual seed used is exposed via Result.seed. Ignored for method="ode".
rtol (float | None) – Relative tolerance for ODE solver. Default 1e-8.
atol (float | None) – Absolute tolerance for ODE solver. Default 1e-8.
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.
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).
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:

Result

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 or n_points are invalid, 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

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 i uses base_seed + i. 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.
rtol (float | None) – Relative tolerance for ODE solver.
atol (float | None) – Absolute tolerance for ODE solver.
max_steps (int | None) – Maximum internal solver steps per output point.
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.
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.
steady_state_tol (float | None) – Tolerance for the steady_state check above. None or <= 0 falls back to atol (matching BNG2.pl).
max_step (float | None)
timeout (float | None)

Returns:

One Result per parameter set (list), or a single squeezed Result with 3D arrays if squeeze=True.

Return type:

list[Result] | Result

Raises:

SimulationError – If any simulation fails.
ValueError – If params is empty or t_span/n_points are invalid.

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)

Parameters:

max_step (float | None)
timeout (float | None)
t_span (tuple[float, float])
n_points (int)
params (Sequence[dict[str, float]] | None)
seed (int | None)
rtol (float | None)
atol (float | None)
max_steps (int | None)
num_processors (int | None)
squeeze (bool)
steady_state (bool)
steady_state_tol (float | None)

Return type:

list[Result] | Result

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.

Replicate i uses seed_base + i (seed_base resolved once from seed, or from system entropy when seed is None; each value is exposed via the corresponding Result.seed), matching the seed schedule run_batch() uses across parameter points.

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().
n_points (int) – As in run() / run_batch().
seed (int | None) – As in run() / run_batch().
timeout (float | None) – As in run() / run_batch().
num_processors (int | None) – As in run() / run_batch().
squeeze (bool) – As in run() / run_batch().

Returns:

One Result per replicate, or a single squeezed Result with 3D arrays (n_replicates, n_times, n_cols) when squeeze=True.

Return type:

list[Result] | Result

Examples

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

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: all model parameters.
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.
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 (float | None) – Absolute tolerance for ODE solver.
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. 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:

Result

Raises:

ValueError – If method is not ‘ode’, or params contains unknown names.
SimulationError – If any chunk simulation fails.

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='newton', rtol=None, atol=None, max_steps=None, sensitivity_params=None)[source]

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

Solver methods:

"newton" (default): KINSOL Newton solver with analytical Jacobian; on non-convergence it falls back EXPLICITLY to the parity integration path.
"integration": CVODE BDF integrated until the BNG2.pl parity criterion ||f(y)||_2 / n_species < tol (run_network -c).
"kinsol": accepted alias for "newton".

Parameters:

tol (float) – Convergence tolerance on ||f(y)||_2 / n_species. Default 1e-9.
max_time (float) – Max integration time for the integration path. Default 1e6.
method (str) – "newton" (default), "integration", or "kinsol" (alias for "newton").
sensitivity_params (list[str], optional) – Parameter names for dY_ss/dp sensitivity.

Return type:

SteadyStateResult

steady_state_batch(params, *, tol=1e-09, max_time=1000000.0, method='newton', rtol=None, atol=None, max_steps=None, n_workers=None)[source]

Compute steady states for multiple parameter sets.

Parameters:

params (sequence of dict[str, float]) – Parameter sets.
method (str) – "newton" (default), "integration", or "kinsol" (alias for "newton"). See steady_state().
n_workers (int, optional) – Number of parallel threads.
tol (float)
max_time (float)

Return type:

list[SteadyStateResult]

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

Add a stop condition checked after each simulation.

Parameters:

condition (str | Callable) –
- str: An expression string evaluated at each time point using observable names as variables. Simulation stops when the expression becomes true (> 0). Example: "A_tot < 10"
- callable: A function f(result) -> bool called after the simulation. If it returns True, the stop condition is triggered.
label (str) – Human-readable label for the condition.

Return type:

None

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 > 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 (float | None) – Solver options (ODE only).
max_steps (int | None) – Solver options (ODE only).

Returns:

Simulation result for the [current_time, t] interval.

Return type:

Result

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

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

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.

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)[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).

Parameters:: state (ndarray)
Return type:: None

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

Set ODE solver tolerances.

Parameters:

rtol (float) – Relative tolerance.
atol (float) – Absolute tolerance.

Return type:

None

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 current_time: float: Current time in interactive simulation.

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:: GenericAlias[double]

set_state(state)[source]

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

Parameters:: state (GenericAlias[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.
Return type:: None

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:

GenericAlias[double]

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. After an advance() these are the post-step observables (read straight from the step result, no recomputation). Before the first advance — or after a set_state() with no advance since — they are computed once, side-effect-free, from the current state via a throwaway model clone.

Return type:: GenericAlias[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.

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:: GenericAlias[double]

from_amounts(amounts)[source]

Molecule amounts → storage vector (amounts / V_c).

Parameters:: amounts (ArrayLike)
Return type:: GenericAlias[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:

GenericAlias[double]

from_counts(counts)[source]

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

Parameters:: counts (ArrayLike)
Return type:: GenericAlias[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:

GenericAlias[double]

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:

GenericAlias[double]

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:: GenericAlias[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:

GenericAlias[double]

read(source)[source]

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

Parameters:: source (Model | ReactionKernel)
Return type:: GenericAlias[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:

None

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()), 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:: GenericAlias[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:

GenericAlias[double]

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:

float

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[GenericAlias[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:

Model

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() / from_concentrations() as the live-volume override.

Parameters:

model (Model | ReactionKernel)
name (str)

Return type:

float

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:

None

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, _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)[source]

Bases: object

Simulation result container.

time

Time points.

Type:: ndarray, shape (n_times,)

species

Species concentrations at each time point.

Type:: ndarray, shape (n_times, n_species)

observables

Observable values at each time point.

Type:: ndarray, shape (n_times, n_observables)

solver_stats

Solver diagnostics (n_steps, n_rhs_evals, etc.).

Type:: dict

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, ...])

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] | None)
_seed (int | 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)

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).

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 time: NDArray[float64]: Time points, shape (n_times,).

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

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:: GenericAlias[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. The raw integrated-state sensitivity for such a species stays available as the low-level sensitivities_species tensor.

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; "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:

GenericAlias[double]

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().

Parameters:

sigma (float | GenericAlias[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:

GenericAlias[double]

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 | GenericAlias[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[[GenericAlias[double], GenericAlias[double]], GenericAlias[double]]) –

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

species — (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}$.

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$.

Return type:

GenericAlias[double]

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
>>> from scipy.optimize import minimize
>>> def objective(p_vec):
...     model.set_params(dict(zip(param_names, p_vec)))
...     model.reset()
...     result = sim.compute_all_sensitivities(...)
...     loss = np.sum((result.species - data)**2)
...     grad = result.gradient(
...         lambda sp, t: 2 * (sp - data)
...     )
...     return loss, grad
>>> minimize(objective, x0, method='L-BFGS-B', jac=True)

Notes

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, *, species_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 (GenericAlias[double]) – Observed data to compare against species trajectories. Must match the shape of self.species (or the selected subset if species_indices is given).
species_indices (list[int] | None) – If given, only these species columns are used for the loss. Sensitivities for all parameters are still included.

Returns:

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

Return type:

tuple[float, GenericAlias[double]]

Raises:

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

Examples

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

chi2_gradient(data, sigma, *, species_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 (GenericAlias[double]) – Observed data.
sigma (float | GenericAlias[double]) –
Measurement noise standard deviation.
- scalar — same σ for all species.
- 1-D array, shape ``(n_species,)`` — per-species σ.
- 2-D array, shape ``(n_times, n_species)`` — per-point σ.
species_indices (list[int] | None) – Subset of species columns to include.

Returns:

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

Return type:

tuple[float, GenericAlias[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, *, species_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 (GenericAlias[double]) – Observed data.
sigma (float | GenericAlias[double]) – Measurement noise standard deviation (scalar, per-species, or per-point).
species_indices (list[int] | None) – Subset of species columns.

Returns:

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

Return type:

tuple[float, GenericAlias[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.

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.

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 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.

This is the raw integrated state tensor — for an SBML AssignmentRule-target species its slot is the frozen state’s sensitivity (~0), not the rule’s. The rule-following derivative is the one output_sensitivities() returns for species:<name> (GH #205); this attribute is the documented low-level escape hatch.

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 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 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_err_test_fails, n_nonlin_iters, n_nonlin_conv_fails.

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.

All zero/empty on every non-SSA backend.

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"].

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.

Return type:

Any

Raises:

ValueError – If a selector is unrecognized or names a species not in the model.
RuntimeError – Only if the result was loaded from disk and the array is 3-D batch-shaped (call squeeze() per-row first).

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:

Result

Raises:

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

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:

None

to_cdat(path)[source]

Export species in BNG .cdat format.

Parameters:: path (str | Path) – Output file path.
Return type:: None

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:

None

Notes

Return type:

None

Parameters:

path (str | Path)
kind (str)
sep (str)
include_time (bool)
header (bool)

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. The parameter / IC-species name lists are taken from the first result.

Parameters:: results (list[Result]) – Results from individual simulations (all must have the same n_times, n_species, n_observables).
Returns:: A batch result with 3D arrays.
Return type:: Result

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().

fim

The (symmetric, positive-semidefinite) Fisher Information Matrix the readout was computed from, over n parameters (or IC species).

Type:: ndarray, shape (n, n)

parameters

Names labelling the FIM axes — the Result.sensitivity_params (axis="parameter") or Result.sensitivity_ic_species (axis="ic") in matrix order.

Type:: list of str

eigenvalues

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.

Type:: ndarray, shape (n,)

eigenvectors

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.

Type:: ndarray, shape (n, n)

rank

Numerical rank — the number of eigenvalues above threshold.

Type:: int

condition_number

$\lambda_\max / \lambda_\min$; inf when rank-deficient. A large value signals an ill-conditioned (sloppy) parameter set.

Type:: float

identifiable

Per-eigen-direction flag (aligned with eigenvalues): True where the eigenvalue exceeds threshold, False for a practically non-identifiable / sloppy direction.

Type:: ndarray of bool, shape (n,)

non_identifiable_directions

Indices into eigenvalues / eigenvectors columns that are flagged non-identifiable (the False entries of identifiable), smallest eigenvalue first.

Type:: list of int

cramer_rao_bound

$\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).

Type:: ndarray, shape (n, n)

threshold

Eigenvalue cutoff (rtol * lambda_max) below which a direction is flagged non-identifiable.

Type:: float

Parameters:

fim (GenericAlias[double])
parameters (list[str])
eigenvalues (GenericAlias[double])
eigenvectors (GenericAlias[double])
rank (int)
condition_number (float)
identifiable (GenericAlias[bool])
non_identifiable_directions (list[int])
cramer_rao_bound (GenericAlias[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=<factory>, sensitivity_params=(), sensitivity_ic=(), sensitivity_method='staggered', outputs=(), rtol=None, atol=None, max_steps=None, model_sha256=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).
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) – Solver tolerances. None uses the Simulator defaults.
atol (float | None) – Solver tolerances. None uses the Simulator defaults.
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.

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 | None = None

max_steps: int | None = None

model_sha256: str | 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.

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). 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.

concentrations

Species concentrations at steady state.

Type:: ndarray, shape (n_species,)

species_names

Species names.

Type:: list[str]

residual

max|f(y)| at convergence.

Type:: float

method_used

"integration" or "newton".

Type:: str

converged

Whether steady state was found.

Type:: bool

n_steps

Number of solver steps.

Type:: int

n_rhs_evals

Number of RHS evaluations.

Type:: int

sensitivity

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

Type:: ndarray or None

Examples

>>> ss = sim.steady_state()
>>> ss.converged
True
>>> ss["A"]  # species A at steady state
50.0
>>> ss.concentrations
array([50., 25., ...])

residual

method_used

converged

n_steps

n_rhs_evals

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.

to_dict()[source]

Return species concentrations as a dict.

Return type:: dict[str, 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) – NFsim -bscb: when True, block bimolecular reactions whose two reactant patterns match molecules in the same complex. Default: True in bngsim — the NFsim CLI default is off, but bngsim defaults it on for correctness on BLBR/aggregation models that polymerize within a single complex otherwise. Pass False to allow same-complex binding (matching BNG2.pl’s complex=>1). Note: this flag governs only the binding policy; complex bookkeeping needed to count Species-typed observables is enabled automatically whenever the model declares one, regardless of 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().

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:

None

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:

float

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. The actual seed used is exposed via the seed property and stamped onto every Result returned by simulate().
Raises:: SimulationError – If initialization fails.
Return type:: None

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:

Result

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:

None

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:

None

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:

None

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:

None

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()[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. Overwrites any previous snapshot. Mirrors the BNG saveConcentrations() 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.

Raises:: SimulationError – If the session is not initialized.
Return type:: None

restore_concentrations()[source]

Restore the molecular state captured by save_concentrations().

Rewinds the live session to the most recently snapshotted state (counts, component states, and bonds), discarding any simulation progress since the snapshot. Mirrors the BNG resetConcentrations() action.

Raises:: SimulationError – If the session is not initialized, or no snapshot has been saved with save_concentrations().
Return type:: None

has_saved_concentrations()[source]

Whether a save_concentrations() snapshot is available to restore.

Return type:: bool

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.

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:

None

clear_param_overrides()[source]

Clear all parameter overrides set via set_param.

Return type:: None

initialize(seed=None)[source]

Initialize the RuleMonkey session.

Parameters:: seed (int | None) – RNG seed. None (default) draws a fresh seed; pass an integer for reproducibility. The actual seed used is exposed via the seed property and stamped onto every Result returned by simulate().
Return type:: None

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.

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:

Result

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:

None

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:

None

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:

None

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:

None

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:

None

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:

float

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.

Return type:: None
Parameters:: path (str | Path)

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.

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]"]).

Parameters:

data (NDArray[np.float64])
colnames (list[str])

Return type:

NamedArray

colnames

One entry per column, in column order. Always a fresh list — slicing / arithmetic that drops columns invalidates the original mapping, so the inherited list is left untouched on non-trivial transforms.

Type:: list[str]

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([ ... ])

exception bngsim.BngsimError[source]

Bases: RuntimeError

Base exception for all bngsim errors.

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).

timeout

The configured wall-clock budget, in seconds.

Type:: float

elapsed

Actual elapsed wall-clock time at the point the timeout fired, in seconds. elapsed >= timeout always holds.

Type:: float

partial_result

Reserved for future use. Currently always None — bngsim does not yet salvage a partial Result from a timed-out run.

Type:: Result | None

Parameters:

message (str)
timeout (float)
elapsed (float)
partial_result (object)

Return type:

None

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.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.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.

result

Partial simulation result truncated at the stop point.

Type:: Result

condition

Description of the condition that triggered.

Type:: str

Parameters:

message (str)
result (object)
condition (str)

Return type:

None

class bngsim.SsaIssue(severity, code, message, location=None)[source]

Bases: object

A single SSA-compatibility finding from an SBML model.

severity: "error" (Simulator init will raise) or "warning" (logged).

code: Stable machine-readable identifier for the issue category.

message: Human-readable description, ideally pointing the user at a fix.

location: 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.capabilities()[source]

Return a structured capability report for this bngsim install.

Returns:

A dict with three keys:

"version" — the bngsim package version string.
"features" — dict[str, bool] mapping each feature/backend name to its availability flag in this install.
"missing" — dict[str, str] mapping each unavailable feature to a human-readable explanation that distinguishes a missing compiled backend (rebuild flag) from a missing optional Python dependency (pip install ...).

"features" always contains the same keys regardless of build: nfsim, rulemonkey, klu, libsbml, antimony, vivarium, sbml_import, sbml_ssa, sbml_psa, antimony_import, codegen, output_sensitivities. "missing" is empty when every feature is available.

output_sensitivities reports whether this install can emit the (n_times, n_outputs, n_param) output-sensitivity tensor via Result.output_sensitivities() (species/observable/expression derivatives w.r.t. parameters and ICs). Like codegen it is always True — it is the capability handshake fitting frontends (e.g. PyBNF) gate their gradient path on (GH #207).

klu reports whether the SuiteSparse/KLU sparse linear solver was compiled in. When False the ODE backend has only the dense solver, so large/sparse models factorize the full N×N Jacobian at O(N³) — use this to detect a dense-only install before a slow genome-scale run (GH #209).

Feature names are stable across releases; new features may be added but existing names will not be renamed or removed.

Return type:

dict[str, Any]

Examples

>>> import bngsim
>>> caps = bngsim.capabilities()
>>> set(caps) == {"version", "features", "missing"}
True
>>> caps["features"]["nfsim"] == bngsim.HAS_NFSIM
True
>>> caps["features"]["sbml_ssa"] == bngsim.HAS_LIBSBML
True
>>> caps["features"]["klu"] == bngsim.HAS_KLU
True

bngsim.prepare_codegen(net_path, model=None, emit_jac=True)[source]

Generate C code, compile, and return .so path (with caching).

This is the main entry point for the codegen pipeline. It generates combined RHS + sensitivity RHS when possible. The sensitivity RHS is included for all-Elementary models (analytical df/dp + J*v). For Functional/MM models, only the RHS is generated and CVODES uses internal FD for sensitivity.

GH #162: when model (the built model for this .net) is supplied, emit_jac is set, and its analytical Jacobian is complete, the compiled analytical Jacobian (bngsim_codegen_jac dense / bngsim_codegen_jac_sparse CSC) is appended so a .net-loaded large sparse model gets a compiled per-step Jacobian rather than the interpreted fallback.

GH #163: when model is supplied and it has observables/functions and no rateOf csymbol, the compiled output evaluator (bngsim_codegen_outputs, GH #136) is appended so the warm recording loop fills the per-row observable + function buffers with one compiled call instead of the interpreted ExprTk pass. This is INDEPENDENT of emit_jac — outputs are emitted for every jacobian strategy (fd/jax record observables too).

The cache key gains a :codegen_jac, :codegen_outputs, and/or :codegen_output_sens (GH #198) suffix so a .so carrying any of these callbacks never collides with one without it. The suffixes key off cheap O(1) model flags (not the generated source — _codegen_emit_flags), so a cross-process .so cache hit still avoids regenerating the (large) RHS source; a Jacobian derivation that fails the GH #95 budget reports analytical_jacobian_complete == False and drops the :codegen_jac suffix — no cache poisoning.

Parameters:

net_path (str) – Path to the .net file.
model (optional) – The built model (Model or NetworkModel) for this .net. When an analytical Jacobian is wanted (emit_jac), the caller must have prepared it (prepare_analytical_jacobian). Pass None to keep RHS(+sens)-only.
emit_jac (bool) – Whether to append the analytical Jacobian (jacobian in auto/ analytical). Does not affect the output evaluator, which is emitted whenever model qualifies.

Returns:

Path to the compiled shared library.

Return type:

Path

bngsim.parse_net_file(path)[source]

Parse a BNG .net file into a structured dictionary.

Parameters:

path (str | Path) – Path to the .net file.

Returns:

Structured contents of the .net file, with keys:

parameters   : list of (name, value, expression, is_expression)
species      : list of (name, init_conc, is_fixed)
observables  : list of (name, entries), entries = [(sp_idx0, factor), ...]
functions    : list of (name, expression)
reactions    : list of dict with keys reactants, products (0-based
               species indices), type ("elementary"/"functional"),
               rate_law (parameter or function name), stat_factor

Return type:

dict[str, Any]

bngsim.build_model_from_parsed(parsed)[source]

Build a NetworkModel from parsed .net data via ModelBuilder.

Parameters:: parsed (dict[str, Any]) – Output of parse_net_file().
Returns:: The constructed model.
Return type:: bngsim.Model

bngsim.sbml_to_net(sbml_path, out_path=None, *, validate='L2', strict=True, sedml=None, t_span=(0.0, 100.0), n_points=101, rhs_tol=1e-06)[source]

Convert an SBML model to a BioNetGen .net network.

Parameters:

sbml_path (str | Path) – Path to the source SBML .xml.
out_path (str | Path | None) – Where to write the .net. If None, the text is returned in the report but no file is written.
validate (str | None) – "L2" (default) runs the lightweight structural check plus a direct ODE-RHS identity self-check — it reloads the .net and confirms it reproduces the source right-hand side (the project’s faithfulness measure) at the initial and a nonlinear-probe state, no integration. This is what closes the GH #223 silent-loss hole: a network whose constructs were emitted as constants the .net cannot drive (assignment/rate-rule forcing the writer can’t carry) diverges here and is flagged (strict=True raises; strict=False warns and records ConversionReport.rhs_faithful False). "L1" runs only the structural-equivalence check (counts + topology). "full" runs the complete L0–L4 conversion-validation ladder (GH #217) and gates on the hard levels L0–L3 — “convert and prove faithful”; the L4 symbolic check is recorded but never blocks. None skips validation. A "full" gate’s verdict is attached as ConversionReport.validation and decides ConversionReport.ok.
strict (bool) – Raise bngsim.ConversionError on constructs the .net text format cannot carry faithfully. False downgrades to a bngsim.ConversionWarning and emits a best-effort network.
sedml (str | Path | None) – The SED-ML simulation protocol that accompanies the SBML (the COMBINE sidecar). When given, its time course is parsed and the "full" gate’s L3 numerical comparison runs over the model’s own horizon — the exact mirror of net_to_sbml()’s bngl= (avoids the blanket-grid stiff-hang and exercises the trajectory the modeller actually ran). The parsed ProtocolSpec is attached to the report.
t_span (tuple[float, float]) – Fallback time grid for the "full" gate’s L3 comparison when no sedml protocol horizon is available (ignored for other validate values).
n_points (int) – Fallback time grid for the "full" gate’s L3 comparison when no sedml protocol horizon is available (ignored for other validate values).
rhs_tol (float) – Scale-relative tolerance for the "full" gate’s L2 ODE-RHS identity.

Return type:

ConversionReport

API reference¶

bngsim.Model¶

bngsim.Simulator¶

bngsim.NfsimSession¶

bngsim.RuleMonkeySession¶

bngsim.Result¶

Exceptions¶

Universal .net reader¶

Utility functions¶

Generated API (autodoc)¶

Caveats¶

`bngsim.Model`¶

`bngsim.Simulator`¶

`bngsim.NfsimSession`¶

`bngsim.RuleMonkeySession`¶

`bngsim.Result`¶

Universal `.net` reader¶