Solver configuration¶

sim.set_tolerances(rtol=1e-10, atol=1e-10)
sim.set_max_steps(50000)

# Per-run overrides
result = sim.run(
    t_span=(0, 100), n_points=101,
    rtol=1e-12, atol=1e-12, max_steps=100000,
)

# Solver diagnostics
print(result.solver_stats)
# {'n_steps': 1247, 'n_rhs_evals': 2891, 'n_jac_evals': 43, ...}

Jacobian strategy¶

The ODE solver (CVODE BDF/Newton) needs the Jacobian matrix ∂f/∂y at each Newton iteration. BNGsim supports three strategies, selectable via the jacobian keyword:

# Default: auto-select best available strategy
sim = bngsim.Simulator(model, method="ode")  # jacobian="auto"

# Force analytical Jacobian (error if model has Functional/MM rates)
sim = bngsim.Simulator(model, method="ode", jacobian="analytical")

# Force finite-difference Jacobian (baseline, for benchmarking)
sim = bngsim.Simulator(model, method="ode", jacobian="fd")

Strategy	Description	Cost per Jacobian	Availability
`"auto"` (default)	Analytical if available, else finite-difference; falls back to FD if the analytical attempt fails to integrate	—	All models
`"analytical"`	Exact derivatives from mass-action stoichiometry	O(nnz) ops, zero RHS evals	All-Elementary models only
`"fd"`	Finite-difference approximation (SUNDIALS DQ for dense, colored FD for sparse)	O(N) or O(n_colors) RHS evals	All models

After a run, sim.jacobian_strategy reports the strategy that actually produced the result ("analytical", "fd", or "jax") — including "fd" when an "auto" run fell back (see below).

When to use "fd": For A/B benchmarking, or if you suspect the analytical Jacobian is causing issues. The "fd" option uses SUNDIALS’ internal difference-quotient (DQ) approximation for dense models, and graph-colored finite differences for sparse (KLU) models.

Analytical Jacobian: Available automatically when all reactions use Elementary (mass-action) rate laws. For models with Functional or Michaelis-Menten rates, the solver falls back to finite-difference. The analytical Jacobian provides ~10-27% speedup on large models by eliminating O(N) RHS evaluations per Jacobian update.

Auto-fallback on a solver failure (GH #176): jacobian="auto" is a bet — an analytical Jacobian is a strict speedup where it integrates, but it is not guaranteed to. A rate law that is discontinuous in a state variable (e.g. an if() whose condition crosses a threshold the state sits at) has an exact derivative that omits the jump, which can de-stabilize CVODE’s implicit corrector even though the Jacobian is mathematically correct. Under "auto", such a CVODE failure is caught and the integration is transparently retried once with the finite-difference Jacobian (which straddles the step and regularizes the corrector), so the default config still integrates the model. An explicit jacobian="analytical" is not second-guessed — it surfaces the failure.

Build-time derivation budget (GH #95 / #187): for models with Functional or Michaelis-Menten rate laws, the analytical Jacobian terms are derived once (symbolically) at model load. That derivation is wall-clock-budgeted so a pathological small model cannot hang the load — it simply falls back to the FD Jacobian, which is just as fast to solve at small scale. The budget scales with species count and becomes unbounded on genome-scale models (≥ 20,000 species): there an FD Jacobian needs ~n_species RHS evaluations per Newton step and is not a viable solver path, so the analytical Jacobian is treated as mandatory and is always derived to completion. Override the budget with BNGSIM_JAC_DERIV_BUDGET_S (seconds, or inf/none/0 to disable it entirely — the manual genome-scale escape hatch).

Logging¶

import logging
bngsim.configure_logging(logging.DEBUG)
# Now all bngsim operations produce log output