[Montana]

Fixes #1505.

The problem

Result.to_cdat() wrote species from before bngsim's post-solve corrections, so the .cdat disagreed with Result.species.

Simulator._stamp corrects species columns in Python after the solve:

• a linear assignment-rule species (S := 2*A) gets the rule's live value;
• a species in a compartment an event resizes is rescaled to amount over the live volume (SSA reports the wrong concentration after an event resizes a compartment: Result.species divides by the load-time size (exactly 2x after c: 2 -> 4) #1368).

The C++ core's copy is never updated, and to_cdat handed off to Result::to_cdat, which wrote that copy. On main (0b2055a9), against each model's closed form:

case	closed form	Result.species, to_csv(kind="species"), Result.load(save).to_cdat	to_cdat on main
S := 2*A, ODE	S = 20 e^(−t/2): 12.1306 at t = 1	right	20 at every sample
event resize c: 1 → 4, ODE	[Y] = 10 e^(−0.2t)/c: 1.372 at t = 3	right	5.4881 (amount/V_static)
event resize c: 2 → 5, SSA	[X] = 20/c: 4 after the resize	right	10

So one run wrote two different .cdat files depending on whether it went through HDF5. to_csv(kind="species"), which its docstring says matches .cdat, did not.

The fix

One .cdat writer, handed the species the Result reports. The formatter in Result::to_cdat moves into a free function, bngsim::write_cdat, that reads pointers to time, species and names. Result::to_cdat calls it with its own stored data, so its output is unchanged. A new internal binding, _bngsim_core._write_cdat, passes numpy's buffers straight through, without copies.

Result.to_cdat now always calls it with Result.species and its names: the same arrays to_csv, save and to_xarray read. This gives:

• Corrected species are written as corrected.
• Unchanged files where nothing was corrected: a Result whose species were not corrected writes the bytes it wrote before.
• One layout for a loaded Result too. A core-less Result (a loaded file, a stop-condition-cut run_batch row) took a second, Python writer with another layout (# time …, 20-wide fields), so its values were right but not its bytes. It now writes the run's bytes.

A squeezed run_batch or run_replicates Result is refused by name. It holds one trajectory per row, and a .cdat file holds one. It used to fail inside string formatting with a TypeError that named neither the method nor the reason. It now raises a ValueError that says both. That is the one exception type this PR changes.

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, to_cdat now equals the closed form, Result.species, to_csv(kind="species") and the save/load .cdat on all three cases.

3. Corpus .cdat bytes, every .net in the repository plus 300 BioModels (seeded sample), 485 models, written from an ODE and an SSA run and hashed on main and here:

   • Byte-identical: 448 ODE and 354 SSA files.
   • Changed: 47 files, every one from a model whose species the solve corrects: assignment-rule totals (totp53, Myc_total, Total_lymphocyte_count, …) and resized compartments. In each of these, 46 of 47, the new file is Result.species to 1e-11 and main's was the core's uncorrected copy.
   • The 47th, Eftimie2017_2 SSA, differed only in the 6-way parallel run. Alone it is byte-identical on main and here, and deterministic, as were the two such files a first run flagged (Cortes2019, BIOMD0000000434).
   • An observation, not this PR's: an SSA run's .cdat differed under load, a different model each time, and not alone. That suggests something wall-clock dependent in SSA, like rr_ssa: a wall-time replicate cap lets runner speed decide MODEL2003200001's verdict #1603's replicate cap, and is not changed here.

4. Cost: to_cdat of a 100,001 × 60 result (116 MB) takes 1.35–1.43 s on both main and here, best of 5, alternating. A first version copied the arrays on the way in and cost 5%; it now passes numpy's buffers through.

5. New tests fail on main for the right reason. test_cdat_writes_reported_species.py has 10 cases:

   • On main, 8 fail: the three cases against the closed form and Result.species; the three save/load byte checks, whose layouts differ; every route with a core (run_batch, parameter_scan, run_replicates rows and copy.copy); and the squeezed refusal.
   • 2 controls pass on main: an uncorrected Result writes the core's bytes, and so does a GH fix(codegen): a zero-arg observable call emits uncompilable C, and differentiates to zero #71-projected one (an event-mutated parameter, not reported).

6. Mutation tests: 5 mutants, all killed:

   • Python: the old core path restored; the squeezed check removed.
   • C++, each rebuilt into the core: the binding dropping the last column; the binding reversing the columns; Result::to_cdat ignoring the fix(codegen): a zero-arg observable call emits uncompilable C, and differentiates to zero #71 projection.

   The projection mutant survived the first set of tests, since Python no longer calls Result::to_cdat. The file now pins it through the core's own writer.

7. Existing tests nearest the change: the 16 files that name to_cdat, .cdat, _write_bng_file or Result.load, with the new file and the stub-coverage test, 18 files. 520 passed, the same under forced cc, forced MIR and BNGSIM_SKIP_AUDIT=strict.

8. The stub: python/bngsim/_bngsim_core.pyi was regenerated with pybind11_stubgen 2.5.5 and scripts/rebuild_editable.py's own normalization. It differs from main by the one new function. test_core_stub_covers_bindings.py passes.

9. 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 skips nothing.

10. C++ unit tests: ctest 8 of 8 passed.

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

    • The full Python suite: 18,749 passed, 594 skipped. Four failures fail identically on main built the same way in the same image: two tests of a compiler that cannot be started, plus test_run_tests_sh and test_vendor_script_git_cwd, which need a git work tree the tarball does not have.
    • The save/load case: it failed 3 times in that image, which has no h5py, because Result.save needs it. The case now calls pytest.importorskip("h5py"), as the suite's other save tests do. CI installs --extra dev, which carries hdf5, so it runs and does not skip there.
    • With [antimony,hdf5]: the new file passes, 10 of 10. A probe of the run and its reloaded copy shows the same names, dtypes and values, and byte-identical .cdat files.

12. pre-commit --all-files passes on the final tree, clang-format included. ci/changelog.py check and required pass, and local sphinx-build -W exits 0.

Checked against past review findings

• Sibling sites. Every route the issue says reaches the stale writer is pinned:

  • Simulator.run;
  • run_batch(squeeze=False) rows and run_replicates rows;
  • parameter_scan rows;
  • copy.copy(r), which shares the core;
  • a loaded Result, whose bytes now match the run's.

  to_gdat hands off to the core the same way. It is right today, because no pass rewrites _observables or _expressions, and it is left as is. A Python-side fix for SSA observables and expressions read count/V_static in a resized compartment: B := 2*X reports 64 for 32 after c: 2 -> 4 #1464 would make it a sibling, as the issue says.

• One writer, not two. The .cdat layout now has a single definition, which a run, a loaded Result and the C++ API all reach. A Python copy of the format would have had to match std::scientific edge cases (NaN, infinities, three-digit exponents) by hand.

• Bytes where nothing was corrected: measured across the corpus (step 3), and pinned for a plain and a fix(codegen): a zero-arg observable call emits uncompilable C, and differentiates to zero #71-projected model (step 5).

• Cost at scale (cvodes: restore nominal parameters after each difference-quotient probe, so nothing outside the RHS reads p − Δ (#690) #832): measured, and a 5% copy cost removed before opening (step 4).

• Exception types (net: read block markers as BNG2.pl reads them, and refuse what it refuses (#1498) #1539): one changed, deliberately, and stated: a squeezed Result's to_cdat raises ValueError where it raised an accidental TypeError from string formatting. No test pinned the old one.

• Docstrings: to_csv says its "species" output matches .cdat, and it does now. docs/user-guide/results.md describes .cdat as species export, which is unchanged.

• Count bounds, both directions: step 9. The one gate the new file adds is h5py, present on every CI leg (step 11).

• Claims are measured: each number here comes from a run above.

The diff (abridged)

@@ src/result.cpp — the formatter, over pointers
+void write_cdat(const std::string &path, const double *time, std::size_t n_times,
+                const double *species, std::size_t n_species,
+                const std::vector<std::string> &species_names,
+                const std::vector<std::size_t> &cols) {
+    ... // the header and data lines, as Result::to_cdat wrote them
+}
+
 void Result::to_cdat(const std::string &path) const {
     // GH #71: project to the reported species subset ...
     std::vector<std::size_t> cols = reported_species_indices_;
     ...
-    ... // header and data lines, inline
+    write_cdat(path, time_.data(), static_cast<std::size_t>(n_times_), species_.data(),
+               static_cast<std::size_t>(n_species_), species_names_, cols);
 }

The formatting lines move out unchanged, so Result::to_cdat writes what it wrote.

@@ src/_bngsim_core.cpp — the internal binding
+    m.def(
+        "_write_cdat",
+        [](const std::string &path, py::array_t<double, ...> time,
+           py::array_t<double, ...> species, std::vector<std::string> names) {
+            // 1-D time, 2-D species, and the shapes agree, or ValueError
+            ...
+            bngsim::write_cdat(path, time.data(), n_rows, species.data(), n_cols, names, cols);
+        }, ...);

species is exactly the columns to write: Result.species is already the reported projection.

@@ python/bngsim/_result.py — Result.to_cdat
-        if self._core is not None:
-            self._core.to_cdat(str(path))
-        else:
-            _write_bng_file(path, self._time, self._species, self._species_names)
+        from bngsim._bngsim_core import _write_cdat
+
+        if np.ndim(self._species) != 2:
+            raise ValueError("to_cdat writes one trajectory; this Result holds N ...")
+        _write_cdat(str(path), np.ascontiguousarray(self._time, dtype=float),
+                    np.ascontiguousarray(self._species, dtype=float),
+                    list(self._species_names))

_write_bng_file is still used by to_gdat for a Result without a core, which this PR leaves as it was.

Tests

• python/tests/test_cdat_writes_reported_species.py (new, 10 cases):
  • The issue's three cases against the closed form, Result.species and to_csv(kind="species").
  • Save/load bytes: a loaded Result writes the run's bytes.
  • The core's bytes for an uncorrected and a fix(codegen): a zero-arg observable call emits uncompilable C, and differentiates to zero #71-projected Result.
  • Every route with a core.
  • The squeezed refusal.
• No existing test changed.

[elyfmiller]

Montana, thanks. The fix reproduces here. On the PR head (6eb2a7f, arm64 with MIR and h5py), the new file's 10 cases pass; on its base (0b2055a), 8 fail and the two controls pass, as the body says. The 18 nearest files pass: 496 passed and 7 skipped, the same under forced MIR with BNGSIM_SKIP_AUDIT=strict and under forced cc. Those are every test that names to_cdat, .cdat, _write_bng_file, _write_cdat or Result.load, and the stub-coverage test. The branch merges cleanly with today's main (74977c7), and onto it with #1640, which also touches _result.py.

What else we checked:

• Result kinds the tests don't write, head against base:
  • a plain ODE run and a plain SSA run: the same bytes as before;
  • NFsim and RuleMonkey on tests/data/nfsim/context_symmetry.xml: they report no species, and write the same time-only .cdat as before;
  • a loaded plain Result: now the run's bytes, where the base wrote # time x;
  • the second leg of run_until on S := 2*A: S = 7.357590082380 at t = 2, which is 2·A; the base wrote 20;
  • a squeezed run_batch of 3 rows: the ValueError says 3, for species of shape (3, 6, 1).
• More Result kinds, each file equal to Result.species under the reported header:
  • a steady-state run, truncated and reached;
  • run_replicates rows with squeeze=False;
  • run_batch rows a stop condition cut, and StopConditionMet.result (no core);
  • both legs of run_until;
  • a Fortran-ordered and a float32 _species;
  • a single species.
• Every .net file in the repository (185), base against head: no hash changes, and on the head each file equals the core writer's bytes. None of the 185 has a corrected species, so the corrected cases rest on the tests and your BioModels run.
• The full Python suite on our commit: 19,370 passed, 145 skipped, 1 xfailed, 0 failed (19,516 collected). The C++ unit suite, configured as native-tests.yml configures it, passes 6 of 6.
• The docs build: sphinx-build -W passes, and the new Raises section renders.
• simulator: leave assignment-rule species at the value a run reports (#1510) #1637 (Model.get_concentration / get_state and ReactionKernel.get_state return an assignment-rule species' t=0 value after a run: S := 2*A reads 20 for 2.706706 #1510): it writes assignment-rule targets back into the model after a run. It touches neither Result.species nor the writer, so the two compose.
• benchmarks/perf_ab.py --quick, base against head: no slowdown over 21 models, with the solver's work counters the same on every one. The geometric-mean warm-run ratio is 0.994 interpreted and 0.993 compiled.
  • A first run, started right after the head was reinstalled, flagged loads and cold runs about ten times slower on every model.
  • Timed by hand in fresh processes once the environments were warm, both builds load in 1.1 to 1.2 ms and run cold in 4.1 to 4.6 ms. The next perf_ab run found nothing.
  • So that was the reinstall, not the change.

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

1. The file's byte checks couldn't see a change to the format. Each compares two files from the one writer: Result.to_cdat against r._core.to_cdat, which now calls the same write_cdat, and a run's file against its reloaded copy's. A mutant writing setw(19) for the time column changes every file bngsim writes, and all 10 cases pass on it.

   • The new test: test_the_layout_is_the_one_main_wrote pins the header and the t = 0 row that main (0b2055a) wrote. Those values are the initial conditions exactly, and they include a three-digit exponent (1e-100) that overruns its 18-wide field.
   • Its checks: it passes on main and here, and fails on that mutant. It is a change detector, and its docstring says so.

2. A second exception type changes. A loaded Result's to_cdat to a path it can't open raised FileNotFoundError, from its own writer. It now raises RuntimeError("Cannot open file for writing: …"), as a run's Result always has.

   • Why keep it: unifying the two on the run's type moves fewer callers than the other way round.
   • What we added: the changelog now says so, and test_a_path_that_cannot_be_opened_is_refused_alike pins both kinds of Result. It fails on main.

3. Each tolerance says why, and to_cdat's docstring has a Raises section for the squeezed refusal and for a path it can't open. The 1e-6 against the closed forms is room: the ODE cases run at rtol = 1e-10, and the SSA case is exact. The 1e-11 against Result.species and to_csv is the .cdat's 13 significant digits (12 after the point).

Not changed:

4. A loaded Result's .gdat keeps the second layout. to_gdat for a Result without a core still takes _write_bng_file (# time S_obs P_obs), so a run's .gdat and its reloaded copy's differ in bytes: the case this PR closes for .cdat.

   • The values are the same: we checked a .net model whose loader adds observables of its own, and none leak into either file.
   • The same split in errors: to_gdat also still raises FileNotFoundError for a loaded Result's bad path, where a run's raises RuntimeError.
   • Why a note: the core's .gdat writer has the same layout, so routing the loaded case through a shared writer would close both gaps. But it's outside Result.to_cdat() writes species from before the post-solve corrections: S := 2*A stays at 20 for 12.1306 at t = 1, and an SSA resize c: 2 -> 5 leaves [X] at 10 for 4 #1505, and you left to_gdat alone on purpose.

5. Two limits a loaded Result now shares with a run's, through std::ofstream:

   • the writer doesn't check the stream after writing, so a full disk leaves a truncated file with no error, where Python's writer raised OSError;
   • on Windows, a non-ASCII path goes through the narrow-string constructor.

   Both already held for a run's Result, and no file writer in src/ widens its path, so these are notes.