bngsim SBML bridge reloads the model (re-derives the analytical Jacobian) on every objective evaluation

[wshlavacek]

Summary

The bngsim SBML bridge (pybnf/bngsim_sbml_model.py) fully reloads the bngsim
model — libSBML parse → interpretation → analytical-Jacobian derivation → codegen —
inside execute(), which runs once per objective evaluation. None of these
per-model costs amortize across a fit; in particular the symbolic Jacobian (SymPy
sp.diff) is re-derived from scratch on every evaluation, even though it depends only
on the rate-law structure and is identical across parameter sets.

The .net/BNGL bridge already avoids this (loads the model once, applies parameter
values in place via set_param), so the fix is to bring the SBML path in line.

Where (v1.4.0, 50bc026)

pybnf/bngsim_sbml_model.py, execute() — both the TimeCourse (L397) and
ParamScan (L422) branches:

doc = self._build_sbml_doc(mut=mut)
engine_model = self._load_bngsim_model_from_text(_sbml_doc_to_text(doc))  # full load, every eval
result = self._run_simulation(engine_model, ...)

_load_bngsim_model_from_text (L171) calls bngsim.Model.from_sbml_string /
from_sbml — the full loader (parse + interpret + analytical Functional Jacobian +
auto-codegen) — with no caching. execute() is invoked once per parameter set during
fitting/sampling, so the entire pipeline reruns each evaluation.

Cost

The dominant re-paid cost is the analytical-Jacobian derivation (SymPy). Measured with
bngsim 0.9.36:

model	analytical-Jacobian derivation (re-paid per eval)
BIOMD0000000001 (12 species)	~27 ms
BIOMD0000000051	~785 ms

It scales with model size. Parse + interpretation add a further ~10–40 ms re-paid per
eval, and the first evaluation in each worker process additionally pays the one-time
SymPy import (~0.3 s). Because the symbolic Jacobian is independent of parameter
values, every one of these re-derivations after the first is redundant work. As an
illustration, a 100k-evaluation fit of a BIOMD51-class model spends on the order of
~22 h re-deriving the identical Jacobian.

The .net bridge already does the right thing

pybnf/bngsim_model/net_model.py, execute():

model = self._engine_model          # loaded once (Model.from_net), cached on the instance
model.set_param(pname, value)       # apply new parameter values in place
model.reset()

Parameter values don't change the symbolic Jacobian, so set_param keeps the derived
Jacobian valid and the load cost is paid once for the whole fit.

Suggested fix

In the SBML path, load the bngsim model once (cache it on the instance, as the
.net path does) and apply per-evaluation parameter values via
model.set_param(...) + model.reset() instead of rebuilding the SBML document and
reloading. Reserve a reload for cases that genuinely change model structure (e.g.
mutants), not parameter-value changes.

Secondary: integration runs cold on every evaluation (both bridges)

Both bridges construct a fresh bngsim.Simulator(model, ...) per execute()
(bngsim_sbml_model.py:335; bngsim_model/net_model.py:267), so each integration
pays the one-time CVODE workspace / linear-solver setup (cold) instead of reusing a
warm solver. Measured on the same models, the cold→warm gap is ~2.5× on a fast model
(0.49 → 0.20 ms) and negligible on a slow/stiff one. Reusing the Simulator across
evaluations (with reset() to restore initial conditions) would let integration run
warm — a smaller win than the reload fix above, but the same idea of amortizing
per-model setup across a fit.

Note

This is a bridge issue, not a bngsim engine issue: bngsim reuses its solver workspace
correctly when a Simulator is reused, and supports set_param for in-place
parameter updates.

[wshlavacek]

Fixed the primary cost in aab848c (on main).

What changed

pybnf/bngsim_sbml_model.py now loads the bngsim engine model once per
worker process and clones it per evaluation, instead of reloading
(parse + analytical Jacobian + codegen) inside every execute().

• Process-level template cache keyed by a hash of the model's base SBML
  text (_ENGINE_TEMPLATE_CACHE). I verified empirically that the base model
  is unpickled fresh in each dask worker on every job (a 30-job /
  3-worker LocalCluster probe → exactly one __setstate__ per job on the
  worker side), so an instance attribute reloaded via __setstate__ would
  still re-derive the Jacobian per evaluation. A module-scope cache is what
  actually amortizes across the whole fit. (Worth noting: the .net bridge's
  instance-cache + reload-in-__setstate__ reloads from_net per job too —
  it's just cheap there, no SymPy Jacobian.)
• Per evaluation: clone() (~0.4 ms on raf) + set_param / set_concentration
  • reset(). Verified bit-identical (max |Δ| = 0) to the old reload for
    params, species initials, and combined changes.

Correctness guard for parameter-driven initials

One subtlety: if a species' initial value is set by an initialAssignment
(or assignmentRule) that references a fitted/mutated/scanned parameter,
in-place set_param + reset does not recompute the baked-in initial
(confirmed: S0(0)=2·k_init stayed at the base value). So a one-time libSBML
dependency analysis at load (_compute_initial_dependency_names, transitive
through initialAssignments/assignmentRules; algebraicRule → conservative)
decides per-evaluation whether the change can touch a species initial. If it
can, it falls back to the exact previous reload path — zero regression.
The common case (fitting rate constants, or fitting species directly) takes
the fast path.

Both the SBML and Antimony backends benefit (shared execute()).

Result

• ~3.8× on tiny raf.xml (6 species; ~19 → ~5 ms/eval); the win scales with
  model size — for a BIOMD51-class model the per-eval Jacobian re-derivation
  (~785 ms) is eliminated after the first eval in each worker.
• ruff clean; full fast tier 1465 passed / 8 skipped. New tests pin:
  engine loaded exactly once across 4 evaluations, fast-path vs forced-reload
  numerical parity, and the initialAssignment fallback (both directions).

Secondary (warm Simulator reuse) — not in this commit

Each evaluation still constructs a fresh bngsim.Simulator on the cloned
model (cold solver). Since I clone per eval for isolation, warming it would
mean reusing one cached model+simulator with reset() across evals, which
trades the clone's correctness/thread-safety guarantees for a sub-millisecond
gain (the issue measured 0.49 → 0.20 ms on a fast model, negligible on stiff)
— dominated by the ~0.4 ms clone and the now-eliminated reload. Left as an
optional follow-up; happy to take it on if you want it.

[wshlavacek]

Follow-up in 1504f7b (on main): eliminated the remaining reload.

Previously the one fallback case — fitting a parameter that sets a species'
initial via an initialAssignment — did a full reload each evaluation, which
re-derives the Jacobian purely as a side effect of recomputing one
initial-condition number (the Jacobian doesn't depend on parameter values or
ICs, so that was wasted work).

Now those initials are recomputed in place, reusing the cached engine
model: build the doc with the evaluation's changes, evaluate its
initialAssignments via libsbml.SBMLTransforms.expandInitialAssignments
(resolves assignmentRule intermediates), read the resulting initials, and
set_concentration them on the clone. Verified bit-identical to a full reload
for direct (S0 = 2*k_init) and chained (S0 = P, P := 3*k_init) cases.

Two things fell out, each now pinned by a test:

• assignmentRule-governed species are evaluated dynamically by bngsim
  under set_param (bit-identical to reload), so they need no special handling
  — only initialAssignment-on-species drive the recompute.
• Species value application is now unit-aware: bngsim works in
  concentrations, so an amount-based species (hasOnlySubstanceUnits) value is
  divided by its compartment size. This also fixes a latent unit error in the
  first commit's in-place path for amount-based species in non-unit
  compartments.

Full reload is now reserved only for what we deliberately don't analyze in
place: algebraicRules, and amount-based species in non-constant-volume
compartments (the unit factor would itself be parameter-dependent) — plus a
safety fallback if expandInitialAssignments is unavailable.

So to your question "why reload if the IC changes" — we no longer do. ruff +
full fast tier (1469 passed) green; 30 SBML-bridge tests (11 new for this
issue) cover the three paths and their numerical parity against reload.

[wshlavacek]

Closing: the headline cost is resolved and verified.

• Engine model (libSBML parse + analytical Jacobian + codegen) is now loaded
  once per worker process and cloned per evaluation (aab848c).
• Parameter-driven species initials (initialAssignment) are recomputed in
  place via libSBML, reusing the cached model — no reload (1504f7b).
• Full reload survives only for the cases we don't analyze in place
  (algebraicRules, variable-volume compartments).

So both per-eval re-payments called out in the issue — the Jacobian
re-derivation and the libSBML reparse — are gone for the normal fitting/sampling
path. ruff + full fast tier (1469 passed) green; 30 SBML-bridge tests (11 new)
pin the three paths and their bit-identical parity against a full reload.

Deferred (not blocking): the secondary "warm Simulator reuse" note —
each evaluation still constructs a fresh Simulator (cold CVODE setup). That's
~0.3 ms on a fast model and negligible on stiff ones, dominated by the ~0.4 ms
clone and the now-eliminated reload; warming it would mean reusing one
model+simulator across evaluations, trading per-eval isolation for a
sub-millisecond gain. Can spin this out into its own small issue if it's worth
tracking — say the word.