[Montana]

Fixes #1510.

The problem

After a run, an SBML assignment-rule species read its t=0 value through the model, while Result.species reported its live value. On main (0b2055a9), with S := 2*A (linear in species) and N := A^2/(5 + A) (nonlinear), A(0) = 10, A' = −0.5 A, at t = 4:

	rule value	Result.species	get_concentration / get_state on main
ODE, S	2.706706	2.706706	20.000000
ODE, N	0.288283	0.288283	0.288283
SSA, S	4.000000 (2·A at the model's A)	4.000000	20.000000
SSA, N	0.571429	0.571429	6.666667
ReactionKernel, 4 × advance(1.0), S	2.706706		20.000000

So did Simulator.get_state, ReactionKernel.advance's return value, CouplingMap.read (a co-simulation orchestrator reads S through it), and Simulator.snapshot. The same happened on the .net written by sbml_to_net and reloaded with Model.from_net, and with codegen=True.

The cause: a rule target is a fixed slot that the integrator does not move, and only the report (Simulator._apply_ar_report_map) reads the rule.

• A rule linear in species is read as an observable, which the C++ copy-back (it matches functions only) never writes into the slot.
• An SSA run copies back neither kind.

The fix

Simulator._stamp now writes each rule target back into the model at the value the run reports for its last row (_write_ar_targets_back), right after _apply_ar_report_map.

• One value: it is Result.species' own number, in the model's storage convention. It runs before the variable-volume passes, so the issue's separate amount-over-V_live question for resized compartments is left as it is.
• One place: run, run_until (so ReactionKernel), run_replicates, run_batch rows and the chunked sensitivity runs all pass through _stamp, under every method and for the .net path's report map.

The carry-over state goes back around the write. The write uses the core's set_concentration, which also drops a pending carry-over sensitivity seed (GH #210), since a set value is a fresh initial condition. So the write is bracketed by the Simulator's own _capture_carryover_state / _restore_carryover_state, the mechanism the GH #176 retry uses. A write the model did not ask for must not end a pre-equilibration. It calls the core directly rather than Model.set_concentration, which would also log an initial-condition write.

Nothing that follows moves. The slot is read by no rate law and no event. A fixed component's local error is 0, so its value does not enter CVODE's step control either, and a continuation run takes the steps it took (step 4).

How this was checked

1. Reproduced on main with the issue's script: every row of the table above.

2. Fixed, on the same script: every getter equals Result.species and the rule, under ODE and SSA and through the kernel.

3. The sibling routes the issue names, each against the rule's value:

   • codegen=True;
   • CouplingMap.read through a kernel;
   • snapshot() after run_until;
   • an SSA set_state(get_state()) round trip;
   • an hOSU=true target in a compartment of size 2, where the getter equals Result.species;
   • the .net route: Antimony-made SBML through sbml_to_net and Model.from_net.

4. Bit-identity of what follows a run: a corpus A/B over every .net in the repository plus 300 BioModels (seeded sample), 485 models. Each runs 0 → 5, then continues 5 → 10 from the model's own state, on main and here.

   • ODE: 472 of 472 continuation legs are bit-identical.
   • SSA: 373 of 373. Four SSA models timed out on main's side under the 6-way load (Tsur2019, Solis-perez2019, Hu2018, BIOMD0000000434); run alone, both sides finish and their second legs are bit-identical.
   • The model's final state: it differs from main's only on models with rule targets (24 ODE, 22 SSA), as it should.

   Two more cases are pinned by tests: a carry_sensitivities=True pre-equilibration (0 → 3, carried 3 → 6) gives the single 0 → 6 run's dA/dk, and a continuation reports the rule at its new end.

5. New tests fail on main for the right reason. test_rule_species_getters_after_a_run.py has 11 cases:

   • On main, 9 fail: every getter route above, and a target in a growing compartment.
   • 2 controls pass on main: a model without rules, and the carried pre-equilibration, which main never disturbs.

6. Mutation tests: 5 mutants, all killed:

   • the write-back not called (9 failing);
   • the write moved before the report map corrects the value (9 failing);
   • the carry-over state not restored (the pre-equilibration case);
   • the first row written instead of the last (9 failing);
   • the write moved after the variable-volume pass (the growing-compartment case). This one passed every test until that case was added (37d9421).
   • Each tolerance names its reason and the error it was measured at. Four closed-form checks had held a run at rtol 1e-10 to 1e-8, with 2 to 5x room (measured up to 4.2e-09). They now share CLOSED_FORM = 1e-7, 20x room, which still tells the rule's 2.71 from main's 20.

7. Existing tests nearest the change: every file that names get_concentration, get_state(, ReactionKernel, CouplingMap, .snapshot(, carry_sensitivities or _ar_report_map, 80 files and the new one. 3,062 passed, 52 skipped, the same under BNGSIM_SKIP_AUDIT=strict; 3,056 passed and 58 skipped under forced cc and forced MIR.

8. Full Python suite, macOS arm64, no corpus (as CI runs it), on 0b2055a9: 19,372 passed, 141 skipped, 1 xfailed, 0 failed, so 19,514 collected.

   • Collected floor: 63 above PYTEST_COLLECTED_FLOOR (19,451), inside the 500 slack.
   • Skips: the new file's one gate is the C compiler, present on every CI leg.

9. Docker, Linux x86-64 (python:3.12, --platform linux/amd64, gcc 14, built from this tree with the [antimony,hdf5] extras and scipy):

   • The new file: 10 passed.
   • Full suite: 18,809 passed, 537 skipped, 1 xfailed, 4 failed. The 4 are this container's standing failures (test_codegen_compile_reaping, test_codegen_toolchain_verdict, test_run_tests_sh, test_vendor_script_git_cwd), which fail the same way on main in it; none touches the changed code.

10. pre-commit --all-files passes on the final tree. ci/changelog.py check and required pass, and local sphinx-build -W exits 0. No C++ changed.

Checked against past review findings

• Sibling sites (the top review finding): every getter route the issue lists is pinned in a test (step 3). One the issue names as related, steady_state() (steady_state() reports an SBML assignmentRule species at its frozen initial value, and that species' zero Jacobian row refuses the whole model's sensitivity solve #247), is left alone deliberately. steady_state does not write the model's state for any species: A also reads its initial 10 after it, where its steady value is 4, so S reading 20 is consistent with an unmoved model rather than a stale slot.
• A change that could reach other callers: the write's one side effect, dropping the carry-over seed, is undone and pinned by a test that a mutant without the undo fails. Continuation runs are bit-identical across the corpus (step 4).
• A test that can still fail: each mutant fails at least one case (step 6).
• Exception types: none changed.
• Count bounds, both directions: step 8.
• Claims are measured: each number here comes from a run above.

The diff (abridged)

@@ python/bngsim/_simulator.py — Simulator._stamp
         ar_map, ar_blocked = self._ar_sensitivity_metadata(model)
         self._apply_ar_report_map(result, ar_map)
+        self._write_ar_targets_back(result, ar_map, model)
         self._apply_ar_sensitivity_map(result, ar_map)
         self._apply_varvol_conc_map(result, model)

The write sits after the report map, so it writes the number Result.species holds. It comes before the variable-volume passes, so the slot keeps the model's storage convention.

@@ python/bngsim/_simulator.py — the write-back
+    def _write_ar_targets_back(self, result, amap, model):
+        """Leave each AssignmentRule-target species in the model at the value
+        the run reports for it at its last row (#1510). ..."""
+        if not amap or result._species.ndim != 2 or result._species.shape[0] == 0:
+            return
+        core = (self._model if model is None else model)._core
+        last = result._species[-1]
+        idx = {n: i for i, n in enumerate(result._species_names)}
+        writes = [(name, float(last[idx[name]])) for name in amap if name in idx]
+        if not writes:
+            return
+        snapshot = self._capture_carryover_state(core)
+        try:
+            for name, value in writes:
+                core.set_concentration(name, value)
+        finally:
+            self._restore_carryover_state(snapshot, core)

A run on a clone (a run_batch row) writes the clone. A model with no rule targets has an empty map, and nothing is written.

Tests

• python/tests/test_rule_species_getters_after_a_run.py (new, 11 cases):
  • The getters after an ODE run, interpreted and compiled, and after an SSA run, for the linear and the nonlinear rule.
  • The other routes: the kernel, its advance return and CouplingMap.read; snapshot() after run_until; an hOSU=true target in a sized compartment; and the sbml_to_net → from_net route.
  • The state round trip and a continuation.
  • The carried pre-equilibration.
  • A target in a growing compartment: an amount S := 2*A in a compartment grown from 1 to 2 keeps the storage convention, 20 e^-1, where Result.species reports it over the live size.
  • A control without rules.
• No existing test changed. Model.get_concentration's Returns says what a rule species reads after a run, in which convention, and that steady_state() leaves it.

[elyfmiller]

Montana, thanks. The fix reproduces here. On the PR head (37d9421, arm64 with MIR), the new file's 11 cases pass; on its base (0b2055a), 9 fail and the two controls pass, as the body says, and v0.17.1 fails the same 9. The branch merges cleanly with today's main (23ac85b), and onto it with #1638, #1640, #1642 and #1643, which also touch _model.py or _simulator.py.

What else we checked:

• Calls that leave the model where it was: threaded run_replicates and run_batch run on clones, and sequential run_replicates and parameter_scan put the state back. Each leaves the caller's rule targets as they were. compute_all_sensitivities did not (item 1).
• A rule a reaction reads (J1: S => B; kb*S, S a boundary species): sbml_to_net writes the rate against the rule's group, not the slot. B's continuation is the same bits on main and here, on both routes, as the body says.
• A rule that ends non-finite (S := 1/A once SSA takes A to 0): the model now holds inf, the rule's value, and a new ODE or SSA Simulator on it runs.
• A compartment resized by set_param: the rule target and the other species agree (T/A = 3 for T := 3*A), both over the size set.
• ReactionKernel.advance: 33 µs a call here, 35 on main.
• The full Python suite on our commit: 18,953 passed, 402 skipped, 1 xfailed and 1 failed. The one, test_run_tests_sh, wants the optional extras (jax, diffrax, roadrunner, vivarium-core) this machine lacks, and fails the same on v0.17.1. The C++ unit suite, configured as native-tests.yml configures it, passes 6 of 6.
• benchmarks/perf_ab.py --quick, base against our commit: no slowdown over 21 models, with the solver's work counters the same on every one. The geometric-mean warm-run ratio is 0.997 interpreted and 0.998 compiled.
  • A first run flagged loads and cold runs about ten times slower on every model. This machine writes no bytecode, so the newly installed side compiled bngsim from source in every process.
  • With bytecode on both sides they time the same, and the run above found nothing.

We've pushed one commit on top of yours, 0951155, for items 1 to 4:

1. compute_all_sensitivities wrote the end of its run into a model it never advances. Each chunk runs on a clone, but _run_sensitivity_chunk stamped its Result without model=clone, so the write-back went to the caller's model.

   • What it did: after the call A read 10, the start, and S 2.706706, the end; main left S at 20.
   • With n_workers > 1: the chunk threads also wrote the one core, and captured and restored its carry-over state, at the same time, while later chunks cloned it.
   • The fix: each chunk stamps its own clone. Every output is the same bits, hashed over three models with one and two workers.
   • The new test: test_a_call_that_runs_on_clones_leaves_the_model_as_it_was runs compute_all_sensitivities, threaded run_replicates and run_batch, and checks the model's state is the bits it was. The first case fails without the fix.

2. On the .net route, every SSA leg after the first rounded the rule target and warned. A .net has no line that marks a slot a value rather than a molecule count (Species::continuous is false for every .net species), so SSA and PSA round it at the start of a run.

   • What it did: on main the slot kept its t = 0 value, so it was rounded once at most. With the rule's value written back, each later leg rounds it again. S := A/3 with A(0) = 12, so that S(0) = 4 is whole, through sbml_to_net and Model.from_net, warned SsaRoundingWarning on run_until legs 2 to 4, and 20 ReactionKernel.advance calls warned 13 times; main warned on none. Under -W error the second leg raised.
   • The fix: NetworkModel.set_species_continuous, the loaded model's twin of ModelBuilder.set_species_continuous. Model.from_net and build_model_from_parsed mark each target _ar_report_map_from_net finds, as the SBML loader marks them.
   • What else it changes: the first leg no longer rounds a fractional t = 0 value either, so main's first-leg warning and its non_integer_initial_population notice for the slot go too. Every species and expression column is the same bits as main's, under SSA, PSA and ODE, over three models and three legs. The slot's identity observable, _obs_S, reports the slot: under SSA it now holds the rule's value at the leg's start, exactly, where main held the rounded t = 0 value on every leg. Under ODE it already held the leg's start.
   • The new test: test_a_net_rule_target_is_a_value_under_ssa runs four legs through each loader, with the warning as an error, and checks the getter against the rule. Each loader's case fails without its mark.
   • It adds C++ (a method and its binding), so native-tests runs on this PR too.

3. get_concentration's Returns said "amount over the compartment's size at load". After set_param on a compartment, the core holds amount over the size set, so it now says "static size".

4. The changelog says what the .net route now does.