# Changelog All notable changes to PyBNF are documented below. This project adheres to [Keep a Changelog](https://keepachangelog.com) conventions. ## [Unreleased] ### Added - **Two opt-in settings for how scatter search pays for its noise handling's deferral, both measured as no better (#696, ADR-0145).** `ss_noise_optimistic = 1` gives the reference slot to whichever side of an undecided parent-versus-child contest leads on the mean, right away, so the next round's combinations are built from the leading point; the draws continue and hand the slot back if they settle the other way. `ss_noise_redraw_budget = N` caps the re-draws a round queues for the orderings it cannot settle, keeping open only the contests whose objective gap and rank in the reference set say a wrong call would cost most and deciding the rest now on their means. Both were built because #663 measured the deferral as the whole cost of the noise handling, and both do what they were built to do: on an instrumented fit they raise how often a reference slot's point changes from 69 to 113 and 97, back into the range of plain scatter search, the budget while spending a fifth fewer re-draws. Neither helps. Over twenty paired seeds on four problems, against the deferral of 3 that is the default, optimistic acceptance succeeded from 17 seeds of 80 against 24 and a budget from 15 (at 4) and 13 (at 2), the last two significantly worse on success (McNemar p = 0.035 and p = 0.007); against plain scatter search all three are indistinguishable, so each gives back exactly what the deferral had won. The checkpoint curves show why: both lead early and lose at the end, because a reference set that moves faster is following draws it has not earned. Both keys default to 0 and a configuration that does not set them runs exactly as before. - **The stochastic recovery benchmark scores a variant against the method it varies (#696).** `run_baseline.py compare --baseline ` pairs every (problem, seed) both methods were run from and reports a two-sided sign test on the final parameter error, McNemar's exact test on success, and the error of the reported best at checkpoints through the budget, pooled and per problem. This is the comparison the #660 and #663 studies did by hand, now in the runner where the next study can run it; it reproduces their published tables from the committed records. The pieces live in `protocol.py` (`sign_test`, `mcnemar`, `paired`, `best_so_far`, `compare`, `format_comparison`), which stays pure Python and is the copy shared with [stochbench](https://github.com/wshlavacek/stochbench). - **Differential evolution can cross a candidate with the member it will replace (#700, ADR-0144).** Wherever a candidate is not mutated it keeps the values of one parameter set. PyBNF keeps those of the parameter set it was built from, which under the `rand` and `best` strategies stays in the population when the candidate takes another member's place, so values can spread from member to member until a parameter every member shares can never move again. `de_cross_with_target = 1` keeps the values of the member being replaced instead, as published differential evolution does; the copy guarantee (#698) then moves one parameter to a value neither the candidate's base nor that member holds, and the learned mutation settings judge a success against that member, SHADE's own rule. Off by default, under every edition: on noise-free test objectives it removed the frozen parameters and reached the optimum from every seed, and it helped the learned settings on the stochastic recovery benchmark, but with fixed settings it did no better there over twenty seeds, and on two of the tutorial's ODE models it did worse. A configuration that does not set it runs exactly as before. - **Differential evolution can learn its mutation settings during a run (#667, ADR-0142).** `de_adapt_mutation = 1` makes `de` and `ade` draw each candidate's `mutation_rate` and `mutation_factor` from a short history of the pairs that recently produced a candidate better than the parameter set it was built from, the success-history adaptation of SHADE (Tanabe and Fukunaga 2013), so the two settings no one can choose in advance follow the search instead. The configured pair is where the learning starts, `de_adapt_memory` (default 6) is how many generations of successes the history remembers, and the record of which settings built a candidate travels with the candidate, so `ade`'s out-of-order results need no special case. A success is judged against the candidate's base rather than the population slot it competes for, because PyBNF crosses the mutant with the base: judged against the slot, copying a better member counts as a success and the learned rate drifts toward zero. On a ten-parameter Gaussian, five seeds and 4,500 evaluations each, the learned settings reached a median objective of 6e-05 against 0.13 for the fixed defaults; on a ten-parameter Rosenbrock valley 3.8 against 8.3. On the stochastic recovery benchmark at five seeds it turned 6 successes of 30 into 11 for `de` and 4 into 8 for `ade` (tight successes 0 into 5 and 1 into 4), most of them on Hlavacek_PNAS2001 where every seed now lands within a few percent of the truth, while the final error pooled over all seeds was not significantly different and the network-free problem got worse. A control with fixed settings and only the guarantee that one parameter is always mutated showed that guarantee alone accounts for part of the gain, and for `ade` most of it, since plain `ade` collapsed its population onto one parameter set and stopped early in 24 of 30 fits. Off by default: a configuration that does not set the key runs exactly as before. - **A stochastic parameter recovery benchmark, small version (#663, ADR-0140).** `benchmarks/stochastic_recovery/` holds six published stochastic models (Shahrezaei and Swain 2008, Lin and Doering 2016, McKane and Newman 2005, Hlavacek et al. 2001, Munsky et al. 2012, Yang et al. 2008), each with a frozen problem definition (true parameter values, bounds, sampling grid, observables, replicate count and seed offset, simulation budget), committed synthetic data, a scoring protocol (log-decade error, success within a factor of two, cost in simulations, success rate over seeds), a runner that scores the baseline methods (`de`, `ss` and `cmaes` each with and without their noise handling) at equal simulation budgets, and the baseline results. It is what makes the stochastic-fitting claims of #659, #660 and #661 measurable; the full suite the issue asks for is built on it. ### Changed - **Scatter search's noise handling defers an undecided contest for three draws, not five (#663, ADR-0141).** `ss_noise_max_draws` now defaults to 3. Measured on the stochastic recovery benchmark over twenty seeds and four problems, the deferral length was the cost of the noise handling: at five, scatter search accepted about half as many replacements into its reference set per run, spent a tenth of its budget re-drawing points it had seen, and ended no better than without the feature; at three it had the most successes of the settings tried and cost nothing measurable against plain scatter search. A conf that sets `ss_noise_max_draws` itself is unchanged. - **Scatter search gives idle processors another draw of a member whose rank is in doubt (#660 step 4, ADR-0139).** Scatter search waits for every simulation of a round before it builds the next, so toward the end of a round processors sit idle waiting for the slowest simulation, and for a stochastic model the spread in running times is wide. The run now records how many simulations it can execute at once, and when a stochastic scatter search has fewer in flight than that, the difference is filled with fresh draws of the members the noise cannot order, the same draws it would have spent at the round's end, each within `ss_noise_max_draws`. `ss_fill_idle = 0` turns just this off, for a fit that would rather end each round as soon as its own simulations do. A deterministic fit, a run whose processor count cannot be read, and a round with nothing idle are unchanged. - **Scatter search fills the diverse half of its first reference set by distance under `edition = 2` (#660 step 2, ADR-0137).** Glover's template builds the first reference set from the best half of the initial population and the most *diverse* half of the rest, chosen one at a time as the candidate farthest from the nearest member already in the set. PyBNF picked that half at random, which is diverse only on average and, with many parameters, far less so than choosing for it. Distance is measured in the parameter sampling space with each parameter divided by its spread over the initial population. The legacy edition keeps the random choice, and `ss_diverse_by_distance` sets it either way. - **PyBNF requires Python 3.12 or newer, and the `petab` extra requires petab 0.9 or newer (#591).** petab 0.9.0 is the first release that loads a `language: bngl` model natively, through the `BnglModel` loader PyBNF contributed upstream (PEtab-dev/libpetab-python#508), and it requires Python 3.12. Python 3.11 support is dropped with it. ### Removed - **The `pybnf.petab.bngl_model` module, its `BnglModel` adapter, and the `register_bngl()` shim that taught older petab releases to load BNGL models (#591).** petab now does this itself, so PEtab's own validator checks a BNGL-model problem with no PyBNF code involved; `register_bngl()` had already collapsed to a no-op on petab 0.9.0. The stand-alone BNGL reader and the parameter-expression evaluator (#666) stay: the importer and exporter use the reader, and the evaluator is the staging copy for an upstream port, since petab's native loader does not yet evaluate an expression-valued parameter. ### Added - **Scatter search gained the improvement method of Glover's template (#660 step 1, ADR-0138).** PyBNF's scatter search relied on recombination alone; the template refines the candidates combination produces with a local search, which Egea and colleagues' version also does and which the issue names the largest gap. Under `edition = 2` the best child a round accepts into the reference set now starts a Nelder-Mead simplex search, the same one `job_type = sim` runs, driven asynchronously alongside the rounds that follow so it costs wall clock only when processors are short. When it finishes, its best point replaces the member it started from if it improved on it, else the worst member if it beats that, else it is archived. Three filters keep it from running on every candidate: at most one start every `ss_local_every` rounds, at most `ss_local_max_running` at a time, and none within a small distance of a refinement already started or found. `ss_local_search` turns it on or off; `ss_local_max_iterations` bounds each search. It is off under noise handling, since a simplex over single draws of a stochastic model converges on the noise. The legacy edition is unchanged. - **An observable's scale and offset can be solved out of the search (`linear_profiling = 1`, #671, ADR-0132).** A free parameter that an `observable: , formula: ` line reads linearly, as a scale, an offset, or the coupled pair, has a closed-form optimum for any fixed dynamics: the weighted least-squares fit of the observable to its data. Searching it anyway spends evaluations ranking candidates by how wrong their scale happens to be, and a twenty-decade box on an offset becomes a search dimension. With the new switch every such coefficient leaves the search and is solved at each evaluation from the data, so every parameter set the fit scores is coefficient-optimal by construction. Coefficients that share an observable are solved jointly, and a coefficient tied across several experiments is one solve over all of them. The evaluation behind the switch (ADR-0130) found that on a linear-Gaussian fixture the profiled search reaches the optimum in about half the simulations, and that its ordering of candidates tracks the true rate constants where the searched ordering does not. The switch is all-or-nothing and is refused before the run starts, naming the parameter and the reason, for anything it cannot solve: a coefficient that is also a model parameter, one a noise source reads (so moving it would move sigma), one that enters an observable nonlinearly, an observable whose noise family is not Gaussian, a log-scale Gaussian read by anything but a single scale of the whole formula (that one case is solved in log space, the geometric-mean form of `normalization = scale`, ADR-0134), a formula affine in each coefficient but not jointly, a cumulative or analytically scaled observable, a prediction-dependent sigma, and, with `noise_profiling` also on, a group whose observables do not share one profiled sigma. It is also refused for the Bayesian samplers, and for `ms` and `design`, whose assembly is not the shared one. The gradient optimizers `lbfgs`, `gntr` and `trf` run with it: the gradient of the profiled objective is the partial derivative at the solved coefficients, and the residual Jacobian and Gauss-Newton curvature are projected off the span of the solved design, Kaufman's variable-projection form, so the curvature is the information about the dynamics with the coefficients estimated rather than known (ADR-0133). A profiled coefficient stays declared and estimated: it counts in `k`, and its value is reported in `Results/profiled_linear.txt`. Its declared bounds are respected rather than ignored, since the closed form can return a negative scale for a parameter declared positive; the solve stays inside the box and the file says when a bound held. - **A stochastic fit now confirms its best fit by running the top parameter sets again (#659).** When a model is stochastic, running it twice with the same parameter values gives two different answers, so the objective value PyBNF computes is a noisy measurement rather than a fixed number. A fit picked its answer by taking the best value it ever saw, and each of those values came from a single simulation. A long fit scores tens of thousands of parameter sets, so the winner of that comparison was very often the parameter set that happened to get a lucky simulation. Two things came out wrong and nothing said so: the reported objective value was the best of many noisy draws and so was optimistic by a wide margin, and the reported parameter values were not the best ones found, because a slightly worse parameter set with a lucky draw beats a better one with an average draw. A fit that uses at least one stochastic model now ends by running its top `best_fit_candidates` parameter sets `best_fit_replicates` more times each, ranking them by their average objective value, and reporting that winner as the best fit. Ten and ten by default, so a hundred simulations against a run that may have done a hundred thousand, all submitted at once on processors that would otherwise sit idle at the end of a run. `Results/best_fit_confirmation.txt` gives each candidate's average, its standard error, and the single value the search had recorded for it, so you can see how much of the search's answer was luck. Everything the run writes afterwards, including the saved simulations, the best-fit model file, the information criteria, a refine's starting point, and a bootstrap replicate's recorded answer, describes the confirmed parameter set. The stage is on under `edition = 2` and above. Under the legacy edition it is off, since an unchanged configuration file has to keep behaving as it always has, and a legacy fit with a stochastic model is told at startup what it is missing and which two keys turn it on. Nothing happens for a fit with no stochastic model. No fitting method searches any differently. This covers the answer the fit reports and not the search that produced it, which is a larger piece of work; the new table is a way to measure how much that matters for a given model. - **A multi-start fit now reports how every one of its starts did (#658).** Several fit types run more than one search, each from a different starting point, and report the best result. That one number cannot be checked. Twenty starts that all reached about the same objective value mean the fit has very likely found the best answer available, and running more starts would not help. Twenty starts that all landed somewhere different mean the reported answer is only the least bad of twenty poor ones, so more starts are needed or the model and the parameter bounds need another look. Both cases used to print one number and nothing else. A fit with more than one start now writes `Results/multistart_summary.txt`: one row per start, sorted by final objective value from best to worst, with the objective that start reached, the steps it took, the simulations it cost, and why it stopped. A short version, including how many starts came within a tenth of a percent of the best, is printed at the end of the run. Reading the objective column downward is the check the parameter fitting literature calls a waterfall plot. This covers `trf`, `lbfgs`, `gntr`, `powell`, `sim`, `ms`, the polishing phase of `profile_likelihood`, and the metaheuristics `de`, `ade`, `ss` and `pso`. A start that was still running when the fit ended, and a start the fit never reached, are both listed and labelled, so a run stopped early by `wall_time_fit` cannot be misread as a complete set of starts that agreed with each other. Nothing is written for a fit with a single start. No fitting method searches any differently. - **PyBNF now says what to measure next (`job_type = design`, #574, ADR-0129).** A profile-likelihood run ends by telling you a parameter is practically non-identifiable, which is a diagnosis with no prescription. The new design run answers the question that follows it. It reads the expected Fisher information PyBNF already assembles for the `gntr` optimizer, notices that the information is a plain sum over the measured points, and scores a planned measurement by the one term it would add. So every noise model, scale and transform the fit already supports comes along, and nothing is re-simulated: the sensitivities at every simulated time were computed when the best fit was scored. A recommendation is one observable, in one experiment, at one time. The observable has to be one that experiment already measures, so its precision is known rather than invented. `design_criterion` chooses what makes one design better than another: the average variance of the parameters (`a`, the default, aimed at the parameters named by `design_target`), the volume of the joint confidence region (`d`), or the worst-determined direction (`e`). Naming a single target is the classical c-criterion, which is what a profile-likelihood verdict about one parameter asks for. The same measurement can be recommended twice, meaning measure it twice. For a time course PyBNF simulates the times the data was measured at, so by default a design could only recommend repeating an existing measurement. `design_grid` adds extra simulated times to choose from and `design_t_end` moves the far end of that window past the last measurement. The measured times are always kept, so the scoring of the data is unchanged. The report in `Results/experimental_design.txt` has two halves: the measurements to make, and each parameter's confidence interval now and after, at the same threshold a profile-likelihood run quotes. If measuring every candidate at once still leaves a target parameter undetermined, the run says so instead of recommending measurements that cannot help. Set `profile_likelihood_design = 1` and an identifiability run ends by writing the same report, around the optimum it just found and aimed at the parameters it just flagged. Off by default. Both surfaces are documented under gradient-based fitting. ### Fixed - **A half-bounded prior no longer hands CMA-ES its shape parameter as a search width (#777, ADR-0118 amendment).** A prior truncated on one side (`lower: 1e-12, upper: inf`, the one-sided box ADR-0047 made first class) reports a bounded support, so the box-mode optimizers accept it -- deliberately. But a half-line has no width, and CMA-ES needs one per coordinate to scale its first population. The substitute was `abs(p2 - p1)`, which is not a length in the parameter's units: for a location-scale family it is `|scale - location|`, and for the shape-scale families (`gamma`, `inv_gamma`, `weibull`) and `beta` it subtracts a dimensionless *shape* from a scale. Measured on the same prior and the same floor, moving only the open side: `gamma, shape: 2, scale: 1e-9, lower: 1e-12, upper: inf` gave a width of **2.0** -- the shape parameter -- for a coordinate whose plausible range is 1e-12 to 1e-6, about two million times too large; with `upper: 1e-6` the same declaration gave `1e-6`, the box. CMA-ES squares the width into its initial covariance diagonal, so the first population was drawn that far out: at a population of 200 over two such coordinates, **none** of the first generation's 400 values landed between 1e-12 and 1e-6, and their median was 0.40 -- about 2e8 times the start point. CMA-ES adapts, so the run still converged; the generations spent walking the step back down were the whole cost, and nothing said so. An open side now takes the **truncated prior's central 80% interval**, `ppf(0.9) - ppf(0.1)` in the parameter's sampling space: a length in the coordinate's own units by construction, finite on a half-line for every family in the catalog (verified over all fifteen non-`uniform` families, both truncation directions and both scales), and derived from the distribution that was declared. The gamma case above becomes `3.36e-09`, and all 400 of that first generation then land in range; a `normal, parameter_scale: log10, mean: -9, sd: 0.5, lower: 1e-12, upper: inf` becomes `1.28` -- the prior's own 80% spread -- where it was `9.5`, `|sd - mean|`. A coordinate with a finite box on both sides is untouched, bit for bit. The substitute is the prior's spread, not a range you stated, so `docs/priors.rst` now says which number a half-bounded declaration gets and that writing both sides is how to state a range. `docs/config_keys.rst` stated the rule as "an unbounded prior is refused", which does not distinguish a half-bounded declaration from one with no bounds at all; it now names both cases and says what each gets. - **A simulation folder that cannot be created is no longer retried 1000 times under new names, and the reason reaches the log (#791).** `Job.run_simulation` creates its working folder in a retry loop whose recovery is to take a new name. That is exactly right for the case it was written for -- dask can run the same job twice, so the folder is already there -- and it cannot help with any other `OSError`. The loop caught bare `OSError` anyway, so a missing parent, a read-only or full filesystem, a permission denial or a stale mount each ran the full 1000 attempts and emitted 1001 warning lines, per job, before giving up; and the exception was discarded, so the message that did reach the user named neither the errno nor the strerror it was carrying. Measured against real conditions: `EEXIST` succeeded on attempt 1, while `ENOENT` and `EACCES` each took 1001 attempts to fail. `FileExistsError` is now caught on its own and keeps the rename and the 1000-attempt cap; every other `OSError` fails the job immediately with `Job failed because it was unable to write to the Simulations folder: could not create : (errno )`. Nothing else about the job changes, and a folder collision still recovers on the first retry. There is no simulator log to read for this failure, because no simulator ran; the troubleshooting docs now say so and point at the output directory rather than the model. - **MCMC checkpointing no longer re-serializes millions of small arrays every iteration (#789).** `chain_history` is appended to once per chain per iteration and never trimmed, `should_pickle` keeps it, and the checkpoint fires once per iteration by default (`get_backup_every()` is `backup_every * population_size * smoothing`, all 1). So `backup()` re-serialized the entire growing history every iteration, and the cost of a run was quadratic in `max_iterations`. Measured on a real algorithm at 20 chains, the whole-algorithm pickle went from 0.04 MB empty to 46.45 MB and 1.335 s at 40,000 iterations -- the history was essentially the whole checkpoint. Each chain is now stored as one growing array (`ChainRecord`) rather than a list of per-step arrays, so the pickle is a handful of buffers instead of millions of objects. Same algorithm, same draws: **25.64 MB and 0.0052 s at 40,000 iterations, about half the size and 258x faster.** Integrated over a 20-chain, 100,000-iteration run that is the difference between tens of hours of pickling and about ten minutes. Resident memory falls by about the same factor of two, because a `(n_dim,)` array costs 112 bytes of object overhead per recorded step regardless of `n_dim`. Nothing else changes. `ChainRecord` supports append, `len`, indexing, slicing, slice assignment, iteration and comparison against a plain list, so every reader -- `diagnostics.split_chains`, DREAM's `_update_preconditioner` and its outlier reset -- works on it unmodified, and a test that assigns a plain list in its place still does. Slicing returns a copy, as list slicing does, so no caller can alias a buffer a later append reallocates. R-hat, ESS and `split_chains` are byte-identical between the two representations, including through the `start_floor` path added for #787. An algorithm resumed from a backup written before this change unpickles its plain lists and keeps working, since every reader is duck-typed. - **DREAM's outlier reset no longer lets a chain be compared against a duplicate of itself (#787).** `detect_and_reset_outliers` resets an outlier chain by copying the donor's history over it rather than discarding the outlier's own. That keeps the archive `_update_preconditioner` pools well formed, but it files the donor's draws under the outlier's index, and the convergence diagnostics read the same structure. Two chains then agreed exactly over that window, so R-hat's between-chain variance was deflated and it read low -- toward declaring a convergence that had not happened, which is the direction that matters. ESS was inflated by the same duplication, counting the copied draws twice. Measured on four chains that have genuinely not converged (means 0.0, 0.8, 1.6, 2.4), with one chain's window replaced by another's exactly as the reset does: max R-hat fell from 1.3664, 1.3885, 1.4046 and 1.4079 to 1.2478, 1.2585, 1.2648 and 1.2639 -- about -0.14 on every seed, from a single duplicated chain in four over a third of the window. The copied window was in reach because resets are gated to burn-in but `chain_history` is never truncated, while `split_chains` reads the last 50% of the whole history. At iteration T the window is `[T/2, T]`, so it overlapped the copied burn-in ranges whenever `T < 2 * burn_in`, and `check_convergence` is live for `T > burn_in`. Every shipped `dream` and `p_dream` benchmark conf sets `max_iterations` to exactly twice `burn_in`, so that band was their entire sampling phase. The reset now records where the chain's own history resumes, and the diagnostics start no earlier than the latest such point across the replicas they compare. `split_chains`, `rhat` and `ess` take an optional `start_floor` that defaults to 0 and is byte-identical to the previous two-argument call, so every sampler that never overwrites a chain -- `mh`, `pt`, `am`, `hmc`, and any `dream` run with no outlier -- is unchanged. The floor only moves forward, counts only the replicas the diagnostics read (#782), and is created on demand so an algorithm resumed from a backup written before it existed keeps working. It is reported in `Results/diagnostics_meta.json` as `history_starts_at`. - **Parallel tempering's R-hat and ESS describe the replicas that carry the posterior, not every replica on the beta ladder (#782).** `should_sample` gates `sample_pset`, so under `pt` only the max-beta replicas reach `samples.txt`, `log_likelihood.txt`, the histograms and the credible intervals. `compute_rhat` and `compute_ess` ignored it and passed `self.chain_history` with `self.num_parallel` -- every replica, including the tempered ones. A replica at beta < 1 is sampling p(x)^beta, so the diagnostics were mixing chains with different targets into one statistic about the posterior. Both directions were wrong. R-hat picked up the spread of the beta ladder, which does not shrink as the run converges: a user tightening `rhat_threshold` toward 1.05 under `pt` was chasing a floor set by `beta_range`, not a convergence criterion. ESS counted draws that were never reported, overstating how much posterior information the run had -- and `ess_per_eval` inherited that. On a constructed ladder whose two max-beta replicas share one target and whose two hot replicas sample the flatter distribution their beta = 0.5 implies, PyBNF reported max R-hat 1.21 where the truth is 1.00, and bulk ESS 800 where the posterior carries 400 draws. `should_sample` now has a default on `BayesianAlgorithm` returning True, and the two diagnostics read it through a new `_posterior_chain_history`. `mh`, `am`, `dream`, `p_dream` and `hmc` run `population_size` copies of one target, take the default, and are byte-identical -- verified against the all-replica computation for each. Only `pt` changes. With the default `reps_per_beta = 1` there is exactly one max-beta replica, so R-hat is now a split-R-hat over that one chain. That is a real diagnostic -- halving the chain is what the split is for -- but it cannot see a chain that never left one mode, which is the failure pt is usually run to avoid. That is a property of the default configuration rather than a mistake a user made, so it is reported by naming the statistic, not by a warning every pt run would learn to scroll past: the reported line now carries the number of chains it compared, reading `Max R-hat: 1.0034 (split, 1 chain)` or `(4 chains)`. A warning is reserved for the one case where that statistic drives a decision -- `rhat_threshold > 0`, which makes the run stop itself on it. For machines, the same provenance is written to a new `Results/diagnostics_meta.json` beside the table: a schema version, how many chains were compared, which replicas they were, their betas (`null` for the untempered samplers) and `num_parallel`. `pybnf.inference_data.from_pybnf` surfaces it as the `pybnf_diagnostics_chains`, `pybnf_diagnostics_replicas` and `pybnf_diagnostics_betas` attrs. It is a sidecar rather than a column because `diagnostics.txt` is appended to across a `--resume` (only a fresh run clears `Results/`), so a header or column change would let a resumed run write rows that disagree with the header written above them by an older version -- silently. The sidecar is rewritten in full on each diagnostics write and so has no append protocol to get wrong, and `diagnostics.txt` stays byte-for-byte as it was. The remedy is not simply to raise `reps_per_beta`. The number of temperatures is `population_size // reps_per_beta`, so raising it alone shortens the ladder and weakens the exchange the method exists for: at `population_size = 8` over `beta_range = 0.01 1`, going from 1 to 2 takes the largest ratio between adjacent betas from 1.93 to 4.64, and 4 collapses the ladder to `[0.01, 1.0]`. Since exchanges are accepted with probability `min(1, exp(dbeta * dF))`, that suppresses them. A between-chain R-hat costs twice the replicas -- `reps_per_beta = 2` with `population_size = 16` -- and the `diagnostics_every` and `reps_per_beta` entries in the config-key docs now carry the numbers. - **Parallel tempering records each iteration of a chain once, instead of recording the iteration after a replica exchange twice (#710).** At an exchange barrier `replica_exchange` resumed each chain by calling `try_to_choose_new_pset`, then rewound the iteration counter by one as if that call had only bumped the counter. It does more than that: after advancing the counter it runs the whole per-iteration block -- record a statistical sample, refresh the histograms, write output, run the convergence diagnostics. So the block ran once at iteration E+1 on the post-exchange state, the rewind put the counter back to E, and the next result advanced to E+1 and ran it a second time on the post accept/reject state. One iteration of one chain contributed two draws to `samples.txt`, two rows to `log_likelihood.txt`, and double weight in `Results/Histograms/*.txt` and the credible intervals -- two draws that are not even the same point, since the accept/reject between them may have moved the chain. Nothing in the output said so; the fit ran to completion and reported a posterior built from more draws than iterations. A run hits this whenever `(k*exchange_every + 1) % sample_every == 0` for some k: every exchange under `sample_every = 1`, and periodically under pairings such as `exchange_every = 19` with `sample_every = 20`, or `exchange_every = 7` with `sample_every = 5`. The shipped defaults (`exchange_every = 20`, `sample_every = 100`) never satisfy it, which is why it went unnoticed. Output already written by an affected run cannot be repaired by de-duplicating it. `sample_pset` labels each row with `current_pset[index].name`, the name of the last *accepted* pset rather than the iteration being recorded, so the two rows for one iteration carry no label that distinguishes them from an ordinary pair of draws -- and they need not hold the same point anyway, since the accept/reject between them may have moved the chain. An affected `samples.txt`, the `log_likelihood.txt` that stays row-aligned with it, and every histogram and credible interval derived from the two have to be regenerated by re-running the fit. The per-iteration block moved into `_advance_iteration`, which opens an iteration and does its bookkeeping, and `try_to_choose_new_pset` grew an `advance` flag. The exchange now resumes each chain with `advance=False`: the chain stays on the barrier iteration it already counted and recorded, and the rewind is gone. The psets the exchange proposes carry the same names as before, the run performs the same number of simulations, and a chain that has to spend iterations stuck at the same point to find a move inside the box still counts and records each of those. A driven pt run now records iterations 1..N exactly once each; before the fix it recorded the iteration after every barrier twice. `exchange_every = 1` also runs now. Advancing to E+1 during the resume landed on a barrier under that setting, so every chain parked itself again without proposing a move and the run aborted at the first exchange with "I seem to have gone from one replica exchange to the next replica exchange without proposing a single valid move". Resuming on the barrier iteration does not re-enter that check -- its exchange has just happened -- and the run reaches `max_iterations`, exchanging once per iteration. - **A constraint on a normalized observable is differentiated on the column its penalty is actually scored on (#718).** `normalization` rescales a predicted column in place before scoring and leaves the forward-sensitivity tensor in raw units, which is why the objective's gradient threads the normalizer's own derivative through a quotient/chain rule (ADR-0053/0066). The constraint gradient did not. It read the raw tensor and returned the penalty slope times `d(raw q)/d theta`, while the penalty itself was read out of the rescaled column. A fit under `job_type = lbfgs`, `trf` or `gntr` with both a `normalization` on a measured observable and a `.con`/`.prop` constraint reading it -- one experiment may list `.exp` and `.prop` files together, and a constraint reaches any other suffix through `suffix.Observable` -- was stepping on a gradient whose constraint term was wrong, with nothing to say so: the fit ran to completion and reported a converged result. The error is not a uniform scale a line search absorbs. The chain rule's reference term was dropped entirely, so the gradient's columns came back wrong by *different* factors, and the sum with the objective gradient weighted the constraint by the normalizer instead of by its own weight. On an exactly-known decay column the rate component was too large by 20x under `init`, 28.9x under `peak`, 16.1x under `unit`, 5.25x under `zero`, 1.20x under `floor 0.3` and 37.6x under the `floor, peak` chain; and an initial-condition scale, whose true normalized-column derivative is exactly 0 because it cancels against its own divisor, came back at -0.878 under every divisive method. The Gauss-Newton constraint Hessian `gntr` adds shares the same accessor, so its curvature block carried that error squared. The constraint accessor now threads the same fold the objective's does, keyed by the `(model, suffix, observable)` the readout names -- so a cross-suffix readout reads the records of the `Data` it indexes, a chain folds stage by stage rather than being refused, and the z-score's reductions inherit the measured-row masking added in #726/#727. A column that was never normalized keeps the bare tensor accessor and is byte-identical. Central differences of the real normalize-then-score path now agree with the assembled gradient for every method and for a chain, and a test pins the two accessors to one value for one column. The configuration that reaches this -- one experiment listing both a `.exp` and a `.prop`, with a `normalization` on the measured observable the constraint reads -- is now asserted rather than traced by eye: the resolved grid is keyed by the same `data_key` the `ConstraintSet` is bound to, and `Result.normalize` rescales exactly the column `Constraint.index` returns. - **`fit_type = am` now writes the marginal histograms and the credible intervals it had been silently skipping (#771).** Adaptive MCMC overrode the method that writes them with a bare `pass`, so every run of it produced neither, at any number of parameters. Nothing else in the sampler agreed with that: it called the method on the usual stride, it created the `Results/Histograms` directory the files belong in, and it wrote the standard samples file the method reads. It also accepted `credible_intervals`, `hist_bins` and `output_hist_every` without comment -- those keys are shared across the MCMC fit types, so nothing reported them as doing nothing -- and the documentation presents these files as what a Bayesian fit produces, with no exception for `am`. What a user got was a completed run, an empty histogram directory, and no credible interval for any parameter, from a sampler whose samples were fine all along. The override is gone, and the final pair of files is now written at the end of a run that reaches `max_iterations`, next to the constraint report that was already written there; a run that stops early on convergence already had that call. On a Gaussian target of standard deviation 0.2 the recovered 68% intervals are [0.276, 0.674] and [0.275, 0.706] against an exact [0.30, 0.70]. Every other sampler was unaffected, and now a test says so for all of them rather than for each in isolation. - **A fit with exactly one free parameter and the whitened DREAM proposal no longer dies halfway through burn-in (#767).** `p_dream`, and `dream` with `proposal = whitened`, ended with `ValueError: diag requires an array of at least two dimensions`. The preconditioner estimates the covariance of the pooled chain history, and `numpy.cov` of a *single* column returns that column's variance as a 0-d array rather than a 1x1 matrix, so taking its trace raised; regularizing the diagonal would have raised next. It was not a startup failure, which is what made it expensive. The covariance is first refreshed at `precondition_adapt` iterations, and that key defaults to `burn_in // 2` -- always before `burn_in`, so the crash always landed before the first sample was recorded. Half the burn-in was simulated, `Results/samples.txt` held its header line and nothing else, and the run was a total loss with nothing in the message to act on. The covariance is now shaped to 1x1 before it is used. One parameter is a shape to repair, not a case to skip: its 1x1 covariance is a perfectly good preconditioner -- a scale -- and the rest of the whitened path, the Cholesky factor and the transforms built from it, already handled it. A one-parameter fit now runs to completion with preconditioning active, and on the flat-box target #709 was measured against it samples what it should: over a flat posterior on `[0, 1]` the variance comes back 0.0843 where uniform is 0.0833, with 0.210 of the mass within 0.1 of a wall where uniform is 0.200. Only the whitened proposal was affected, and only at one parameter. Plain `dream` (`proposal = de`) never estimates a covariance, `proposal = kalman` builds its gain from explicit matrix products rather than `numpy.cov` and was checked at one parameter as well, and two or more parameters were always fine. Nothing in the suite had run the whitened proposal with one parameter; a one-parameter fit driven to completion through both of its entry points now runs in the default tier. - **Adaptive MCMC no longer under-samples along the walls of the box: a proposal that leaves it is rejected, not redrawn (#709).** `pick_new_pset` drew its Gaussian step again and again until one landed inside the box, and `got_result` accepted the survivor with the plain Metropolis ratio. That ratio is right only for a symmetric proposal, and a redrawn one is not: its density is the Gaussian renormalized over the box, `q(x -> y) = phi(y - x) / Z(x)`, where `Z(x)` is the share of the Gaussian centred on `x` that lies inside, and `Z` is smallest at a wall. Detailed balance with the plain ratio then holds for `pi(x) Z(x)`, not `pi(x)`, so the chain thinned out in a shell a few proposal widths deep along every wall -- where the edge of a credible interval sits whenever a parameter presses against its bound. Each half was tested and correct on its own, the Gaussian step and the Metropolis ratio; nothing held the two to one target. Measured through the real run loop, 100,000 iterations each way. On a flat posterior over `[0, 1]` with the default `step_size` of 0.2, the density came out at 0.71 of what it should be in the outermost tenth at either end of the box and 1.17 in the middle, and the 68% interval `[0.199, 0.805]` where `[0.16, 0.84]` is exact; it is now flat to within 3% and `[0.161, 0.843]`. That is the fixed-step phase, which every run passes through. In the adaptive phase, on a posterior whose mode sits on a wall (a Gaussian of sd 0.05 at the lower bound), the mass within 0.2 sd of the wall was 8-11% low and the median 5-10% high; both are now within 2%. The proposal is drawn once. If any component leaves its parameter's box the move is rejected then and there: the posterior is zero outside, so no simulation could change the outcome, and none is run. The iteration is spent like any other rejection -- the current point is recorded again in the history, the samples and the trajectory files -- and the chain waits for the others at the generation barrier. A generation in which every chain's proposal leaves the box has nothing to submit, which the scheduler would read as a job pool run dry and end the run on, so the barrier goes around instead, as DREAM's already did. Folding the proposal back inside, which `mh` and `pt` do, is not an alternative here. The fold acts on each coordinate separately, and the folded proposal is symmetric only if the Gaussian is unchanged by flipping the sign of one coordinate: true of their isotropic step, false of the correlated covariance this sampler adapts to. On a flat posterior over the unit square with a proposal correlation of 0.9, folding keeps both marginals flat and gets the joint badly wrong -- 0.209 of the mass in the two corner squares of side 0.2 that the correlation points at and 0.003 in the other two, where each pair should hold 0.080. Rejection gives 0.080 and 0.083. What a user will see change: - The **acceptance rate** reads lower where a wall is in reach, because a proposal that left the box used to be redrawn out of sight and is now counted as the rejection it is. At `verbosity = 2` the rate is printed with the number of proposals that left the box. - The **adapted step is shorter** there, and settles. The scale is steered to an acceptance rate of 0.234, and it never saw a wall rejection, so against a posterior that is broad for its box it grew without limit: on a correlated Gaussian cut off by the unit square the proposal sd passed 11 box-widths by iteration 30,000 and was still climbing, at a cost of 590 Gaussian draws per iteration and rising (the loop allowed 10,000 per parameter). It now settles at 0.58 box-widths, at one draw per iteration. A 150,000-iteration run of that example takes 96 seconds; before, two such runs were a third and two-thirds done after 17 minutes. - A run performs **fewer simulations than iterations**: 16% fewer on the flat example above, 60% fewer on the correlated one. `max_iterations` counts iterations, as before. - A run that never proposes outside the box is unchanged to the bit. One that does follows a different chain from the same `random_seed`, since the redraws are no longer made. Gone with the loop: with exactly one free parameter its fallback could not run (`while num != 10000 * len_params` ends at the 10,000th failure, which is where the fallback began), so a chain that failed 10,000 draws returned no proposal, was dropped from its generation without a word, and left the run waiting on it until the job pool emptied. `mh`, `pt`, `dream` and `p_dream` were held to the same flat-box target and sample it evenly; the defect was `am`'s alone. - **A Bayesian fit with one free parameter, or one constraint, now gets its histograms, credible intervals and constraint report instead of silently getting none (#769).** Both summaries are built by reading a results file back with numpy and checking the array's shape, and numpy drops any axis of length one. A single free parameter made the samples file read back as a flat list of numbers rather than a column, a single constraint did the same to the constraint samples, and a run that had recorded exactly one sample did it along the other axis; each was then taken for a file with nothing in it. The fit itself ran and wrote its samples normally, so what a user saw was a completed run whose `Results` directory was missing the credible intervals -- with only a `No samples collected` line in the log, next to a samples file that plainly had samples in it. Both readers now ask numpy for a two-dimensional array outright, so the number of rows is what decides whether there is anything to summarize. The guard still does its real job: a file holding only its header is still skipped, which matters because these run on a stride from the fit loop and fire before the first sample exists. Affects `mh`, `pt`, `dream` and `p_dream`; `am` writes its own per-chain output and was never affected. Two or more parameters, and two or more constraints, always worked. - **An adaptive-MCMC run no longer ends in `FileNotFoundError` when the combine step meets a per-chain trajectory file it never wrote (#760).** `combine_chains_traj` concatenates each key's per-chain files by loading every one of them by name, with no check that it exists. Both call sites are run terminations, so the failure landed after all the sampling work was done, on the last step before the run reported `'STOP'`. Two configurations reach it, both needing `population_size > 1` and `output_trajectory`: a name no simulation produced a column for is never written for any chain (#755), and a run that converges inside the adaptive window stops before `valid_range`, which is the first iteration at which any trajectory is written -- `check_convergence` needs only `iteration > burn_in`, while the writers need `iteration >= burn_in + adaptive`. The second requires no mistake by the user, and turns a successful early convergence into a crash; it can also leave the chains in different states, since it fires on one chain's iteration count while the others are behind it. (A run too short to reach `valid_range` is not a third way in: the constructor already refuses `max_iterations` below `burn_in + adaptive + 2`.) A source file that was never written is now skipped, and the run says once which keys and chains are missing from the combined output and what the two causes are, rather than either crashing or quietly producing less than the conf asked for. A key with no source files at all leaves no empty combined file behind either, because the combined file is opened only once there is something to put in it. Fixed in the same place, because the rewrite turned it up: `np.loadtxt` drops a one-row file to one dimension and `np.savetxt` then writes that sample down the page as one value per line. A chain with exactly one sampling iteration past `valid_range` had its trajectory transposed in the combined file, mixing a column into a table of rows. - **A recorded adaptive-MCMC sample is written whatever its values, so a trajectory that sits at zero is no longer dropped from the output file (#758).** `write_out_trajactorys` chose what to write by discarding every all-zero row of its buffer. Each buffer holds one slot per chain -- `arr_length` is the constant 1 and nothing increments `factor` -- rewritten with the chain's current state on each sampling iteration and appended to the file, so that filter ran once per sample. It was standing in for "was this slot recorded", which is not what it measures: a recorded trajectory whose observable is zero at every recorded time point looks exactly like a slot nothing was written into. An absent species under a knockout, a perturbation that removes the entity being observed, a `_Cum` counter over a window in which nothing fired -- these are ordinary, and they are often the conditions a fit exists to compare against. The consequence was not a missing file but a quietly incomplete one. `traj__chain_.txt` carries one row per sampling iteration and nothing in it marks a gap, so the row count stopped matching the sample count with no way to tell from the file which samples were absent, and the samples that went missing were exactly those at which the observable was switched off. A posterior predictive band, a mean trajectory or a credible envelope taken from that file was computed over a biased subset. Where every sampling iteration was zero no file appeared at all, which is how this was found. Each slot now records whether a result column supplied it, at the point the accumulating loop fills it, and the writers read that instead of the values. A name that matched no column is therefore still never written (#755) rather than producing a file of zeros. The same substitution of a sentinel for a fill marker is fixed in `write_out_params`, where the discarded row is a parameter vector of all zeros: out of reach for a continuous proposal over several parameters, but a box whose lower bound is zero can be hit exactly by the reflecting fold, and the flag costs one line. `write_out_scores`, which never filtered, was the precedent. - **An `output_trajectory` name that no simulation produces a column for is named once, instead of silently producing no file (#755).** Neither `output_trajectory` nor `output_noise_trajectory` (`fit_type = am`) was validated against anything. A name that matches no column is filled by nothing -- the accumulating loop runs only when the name is a column of the result -- so its buffer stays as allocated, all zeros, and the write step skips it, because it writes only the rows that are not all zero. What the user saw was a `Results/A_MCMC/Runs/` directory with fewer `traj_*.txt` files than the conf asked for, and nothing saying which name was dropped or why: a typo looked exactly like a fit that never asked for that trajectory. The first completed simulation of a run now reports every configured name it carries no column for, once, and lists the columns it does carry. The check is against a real result rather than at config load, where the rest of the codebase refuses an undeclared name (`profile_likelihood_params`, `design_target`), for two reasons: no model class exposes its observable names -- `BNGLModel` records only whether an observables block exists -- and a column need not come from the model file at all, since the measurement layer contributes its own. A completed result is authoritative and carries every model in the fit at once, so a name valid for one model is not reported against another. It is a warning and not a refusal because the fit itself is sound; only an output file is missing, and by the time this can be known the simulations are already running. #751 is why this went unnoticed for so long. The parser used to leave list separators inside the value, and the sampler stripped a comma out of each name it was handed, so `output_trajectory = A,B,C` -- one token to the parser -- reached the write step as the single name `ABC`, which no model has. Had a name that matched nothing ever said so, that would have surfaced the parsing defect years earlier; instead each defect hid the other, and the only symptom of either was output that was not there. Not addressed here: the same write step cannot tell a buffer that was never filled from one filled with zeros, so a trajectory that is legitimately zero across the recorded window is also dropped without a word. That one is data-dependent rather than a property of the conf, and unpicking it means giving the buffers a fill marker instead of reading zero as one. - **A key that takes several values accepts commas between them, and no longer takes a trailing comment as one of the values (#751).** The items of a multi-string value (`profile_likelihood_params`, `output_trajectory`, `output_noise_trajectory`, `design_target`, `design_observables`, `worker_nodes`, `postprocess`, `qualitative_scale`) were split on whitespace by a grammar whose item pattern accepts every punctuation character, so both separators travelled into the value. `profile_likelihood_params = K_par, Kf` parsed as `['K_par,', 'Kf']` and the fit was refused for naming a parameter `K_par,` -- in a message that printed the name with its comma still attached and then advised listing parameter ids, rendering the declared parameters as a comma-separated list, which is the form that had just failed. The comma form is not a guess at what a user might try. `config_keys.rst` documents `output_trajectory = ObservableA, ObservableB, FunctionA` and says in as many words that "multiple values can be defined separated by a comma"; `gradient_fitting.rst` documents `profile_likelihood_params = k1, k2`; four shipped example configs use it. Those two `output_*` keys worked only because the adaptive MCMC sampler stripped a comma out of each name as it read them, which repaired `A, B` and silently joined `A,B` -- a single token to the parser -- into one observable name `AB`, which no model has, whereupon that trajectory was skipped without a word. The hand-strip is gone: the separator is the parser's business now, for every key and every spelling, so `A, B`, `A,B`, `A , B` and `A B` are the same list. Separately, this grammar ran to the end of its value list with no comment rule, and `#` is one of the punctuation characters its item pattern matched, so `profile_likelihood_params = k1 k2 # profile these` parsed as five parameter names and reported `#`, `profile` and `these` as undeclared parameters. - **A profiled noise scale and a profiled linear coefficient are averaged over the same simulations of the best fit as the log-likelihood, and a simulation that could not be scored no longer contributes the previous one's value a second time (#743, ADR-0131 corrected).** Under `noise_profiling = 1` (ADR-0108) or `linear_profiling = 1` (ADR-0132) the value in `Results/profiled_noise.txt` / `Results/profiled_linear.txt` is the only place the estimate for that parameter appears — it is fitted but never proposed, so it is a coordinate of no parameter set and appears in no `sorted_params_*.txt` row. On a stochastic fit that value is averaged over the replicate runs of the best fit, and it was read from the objective before the guard that drops a run which could not be scored. Two consequences, both of them a wrong reported parameter value rather than a wrong caveat: - **A run the objective could not score contributed the *previous* run's values.** The objective assigns its profiled values only once a whole evaluation has succeeded, so a degenerate profile leaves the last successful run's values sitting on it; that run was then counted twice and the failed one contributed a number that was never its own. Over runs worth sigma = 4, unscoreable, and 6, the file reported 4.67 where the two runs that scored give 5. At a single-simulation fit the stale value came from some other parameter set altogether. A stale `at_bound` flag was counted the same way, in a column reported as "how many runs held this coefficient at a bound". - **A run dropped for scoring a different number of points kept contributing its profiled values.** A sum over a different `n` is a different quantity, which is why its log-likelihood is dropped; a scale profiled over a different set of scored points is a different quantity by the same argument. So the two averages were over different sets of runs. Each run's profiled values now travel beside its log-likelihood and are read only for a run that produced one, so every rule that drops a run from the reported average drops its profiled values with it. A fit that profiles nothing is untouched, and so is one whose replicate runs all scored. `_resolve_profiled_noise` still leaves its previous values in place when it refuses a degenerate group: every caller bails out on that refusal, and the guard above closes the one path that read them. Following from that, when *no* simulation of the best fit can be scored the two files are not written at all, where one of them used to be written from whatever the objective was carrying. Neither file states what its value is averaged over — there is no `replicates` or `n` row in them — so a value written when nothing scored could not be told apart from one that was: it is the previous evaluation's whenever the profile was degenerate, and the run's own only when the run scored points and merely summed to a non-finite likelihood. They are also the siblings of `information_criteria.txt`, which is not written either in that case. Because an absent file otherwise reads exactly like a fit that profiled nothing, the run says why in the log and names the profiled parameters it is not reporting. - **`information_criteria.txt` reports how many simulations of the best fit were run beside how many produced a usable log-likelihood, so an average over 3 of 10 runs no longer reads as an average over 3 of 3 (#741, ADR-0131 amended).** A simulation that fails, scores nothing, or scores a different number of points from the rest is left out of the mean, and `replicated_information_criteria` set `replicates` to the surviving count. That number was honest, but nothing in the file said a run had been lost, and the header explained `replicates` as though none could be: "the best fit is run `best_fit_replicates` times and log_likelihood is the mean over those runs". The `%d of %d simulation(s)` warning went to the log alone, and the console's only mention of the count rides on a clause that appears once two runs have produced a value — so the case where nine of ten were lost was the case it hid. This matters where the criteria are used. The writer's own docstring says they "rank this fit against competing models"; that comparison is made by a reader holding two of these files side by side. A model whose best fit lost seven of its ten runs is scored on a mean over the three that worked, optimistic in the way #720 described, and could win the AIC comparison on it. The file now carries `replicates_requested` beside `replicates`, says in words how many runs produced nothing and that the mean is optimistic by however much they would have pulled it down, and prints the same on the console — including when a single run survived. `InformationCriteria` carries the count and `replicated_information_criteria(..., requested=N)` sets it; a criteria object built without one is written exactly as before. The averaging is unchanged and no minimum-success threshold was added, unlike #720's confirmation ranking. That bar exists because candidate means were ranked against each other inside one run and a candidate measured over one run could beat one measured over ten; this file holds one parameter set the run has already settled on, so there is nothing for a survivor mean to beat unfairly, and refusing to emit criteria over a low count would remove information rather than add it. - **The best-fit confirmation stage ranks a candidate on all of its replicate runs, so one that fails most of them can no longer be pinned as the run's answer (#720, ADR-0146).** The stage exists to stop a fit reporting whichever parameter set got a lucky simulation: it runs the search's top `best_fit_candidates` parameter sets `best_fit_replicates` more times each and reports the one with the best average. But a replicate that fails, or that returns a non-finite objective value, produced no number to average, and `mean_objective` averaged only the runs that did — a survivor-only mean over a sample size that varied per candidate from 1 to `best_fit_replicates`. `failures` was carried on the row and printed in a column, but it entered no comparison: the ranking keyed on that mean alone and disqualified a candidate only when every one of its runs had failed. A parameter set that failed nine runs of ten and survived one was therefore scored on a single simulation again, the unlucky draws deleted rather than averaged in, and it beat a parameter set that returned a slightly worse value in ten runs of ten. `_emit_best_fit_confirmation` pins the winner on the trajectory by design, so the saved simulations, the best-fit model file, the information criteria, a refine's start point and a bootstrap replicate's answer all then described the failure-prone set. A candidate is now **confirmed**, and eligible to win, only when more than half of the runs that came back produced a usable value — exactly the condition for the middle of all its runs to be a real number, since a failure is known to be worse than every value. The unconfirmed ones are ranked below every confirmed one whatever their averages say and are shown in the table with a `confirmed` column, their failure counts and an `unconfirmed` line naming them; nothing is hidden, and a candidate that loses a minority of its runs still competes as it did. When no candidate clears the bar the stage pins nothing, the search's own pick stands, and the report and the log say that the simulations are failing rather than offering the best of a bad lot. It was also silent at the point of decision, which is the half a reader would have caught: the console said nothing about failures, and with a single survivor `standard_error` was `None` so the documented "raise `best_fit_replicates` if the standard errors overlap" guidance could not flag it either. The console now says how many of the winner's runs produced nothing and how many candidates were not confirmed. Its one line about a changed answer used to assert flatly that "the one the search liked best does worse when it is run again", which is false when the search's own top pick is the candidate that had the better average and lost on reliability; it now says which of the two things happened. `winner_runs` and `winner_failed` put the winner's sample size in the file beside its average. - **A `time_error` fit is refused by the PEtab export instead of being written as an exact-time one (#738).** A `time_error` clause on a `noise_model` line says the reported measurement times are not exact: the objective integrates each observation's density over a prior on its true sampling time, and loading the config swaps the whole per-point objective for a `MarginalizedTimeObjective`. The exporter had no awareness of the clause at all, so such a job exported without a warning and the emitted problem carried the measurements at their nominal times as though exact. Round-tripped, the fit came back scored by `LikelihoodObjective` — a different statistical model, silently. Its sibling in the very same `noise_model` grammar was guarded all along: `cumulative` is stored under the same kind of structural key and refused by `_reject_cumulative`, for exactly this reason. That function's twin for `time_error` was never written; it is now, beside it. PEtab cannot express the clause — a measurements row carries one exact `time` and has no field for a distribution over it — so, as with the `U` tag, the fix is a refusal naming the boundary rather than a richer export. Walking the whole noise-model field grammar, `time_error`/`sigma_t` was the last of its shape: `cumulative` and a mean-centred `location` were already refused, and the prediction formula, measurement formula and noise source all export and round-trip under test. - **A `U`-tagged free parameter is refused by the PEtab export instead of being written as a hard-bounded one (#736, ADR-0025).** `uniform_var = v 0 10 U` and its `loguniform_var` twin mean the box is enforced only during initialization: it seeds the first population and the search is then free to leave it, so the parameter's box is `(-inf, inf)`. The exporter read the two numbers and dropped the tag, and `lowerBound`/`upperBound` in PEtab are hard box constraints — so the job exported to a table byte-identical to the bounded spelling, and a re-import produced a fit constrained to `[0, 10]` where the original was not. No warning, no exception. It is the third of the same shape as #719 and #733: the export read part of a declaration and discarded the rest without saying so. This one has no mapping waiting to be written. PEtab's nearest shape, blank bounds plus an explicit uniform `priorDistribution` over the box, says something else — PEtab bounds truncate a prior rather than seed a draw — and does not survive the trip either, since the importer resolves a uniform prior to a bounded parameter whatever the bounds say. That is correct for a PEtab row, so nothing changed on the import side. The export now refuses, in code and naming the boundary, which is the contract it already holds itself to for every other construct PEtab cannot state; the message says to drop the tag to export the box as real bounds. Exporting the box with a warning was considered and rejected: the exporter emits no warnings anywhere, and a warning is the signal the start-point work found people do not act on. Nothing measurable is lost — no `U`-tagged declaration exists in the examples, tests or benchmarks. The untagged and `B`-tagged spellings export exactly as before. - **The PEtab export reads edition-2 `parameter:` records, so a free parameter written the new-era way no longer vanishes from the exported problem (#733).** The exporter picked its free parameters out of the config by matching key names against `(_var$|^var$|^logvar$)`. A `parameter:` record is stored under a `('parameter', id)` key, which matches none of the three, so it was skipped rather than refused: no row was written, nothing downstream noticed the id was missing, and the export finished cleanly having dropped the parameter. A conf whose parameters were *all* records reported "No exportable free parameters found", naming only the `*_var` keywords — pointing the user away from the syntax the edition-2 documentation teaches. Both contradict the exporter's own contract, that everything it cannot write raises `NotImplementedError` with the boundary named in code. The damage was widest on truncated priors. A record is the only grammar carrying `lower`/`upper`, so it is what the importer emits for a prior truncated to a box (ADR-0020/0047) — meaning a PEtab problem with bounded priors imported fine and then exported to a table missing exactly those parameters. On the tutorial's own PEtab priors problem (lesson 15) that was three of four: a log-normal, a gamma and a normal all disappeared, leaving one plain uniform, and the re-import produced a fit over one parameter instead of four. The exporter now reads both declaration spellings in one pass, in declaration order, and builds each record through the same mapping the fitter loads a job with, which moved to `pybnf/parameter_record.py` so the two cannot drift. A record therefore meets the same boundaries the positional line does, reached by a different spelling: a no-prior point start, a natural-log sampling scale, a three-parameter family such as student_t, and the log forms PEtab defines for no family are each refused with the keyword the record built. `initial_value:` is honoured as the start point it is, alongside a `start_point` line when the two agree and refused when they disagree, and an out-of-box one is a `PybnfError` rather than the bare `OutOfBoundsException` that reaches users as "an unknown error". - **The PEtab export writes the fit's start point, so a round trip no longer moves the fit back to a sampled draw (#719).** `nominalValue` is where a PEtab problem states the point a fit starts from, and since #583 PyBNF reads it: an imported problem's nominal becomes a `start_point` line. The export direction dropped it twice over. `write_parameter_table` built each record as four fields and had no `nominalValue` column, so a row that carried one wrote a file that could not reproduce it; and `_free_parameters_from_conf` built its free parameters from the `_var` lines alone, never reading the `start_point` lines beside them, so on the plain "publish my job as PEtab" path there was no value to drop in the first place. Both are fixed. A PEtab problem now survives import -> export -> import with its start point intact, and a native conf's `start_point` line exports as that parameter's `nominalValue`. The column is written only when some parameter declares a start, and a parameter that declares none writes an empty cell, so a job with no start point exports exactly the four-column table it did before. A `start_point` outside the parameter's own box is refused at export, since PEtab has no way to state one and the re-import rejects exactly that; so is a `start_point` naming a parameter no exported declaration claims, which was the other way the value could vanish without a word. This closes the gap between the export and the round-trip identity the PEtab documentation states. - **A proposal within a rounding of a parameter's bound no longer ends the fit (#706, #750).** `FreeParameter._reflect` folds an out-of-box proposal back between a parameter's walls, and every optimizer and sampler proposal reaches it through `set_value`. It could return a value a hair outside the box, and the constructor's bounds check then raised an `OutOfBoundsException` that nothing on the proposal path catches, so the run stopped with "Sorry, an unknown error occurred". Two roundings did it. With two finite walls the fold reconstructed the descending leg of its triangle wave from the far wall, so a proposal one ulp below the lower wall came back below it -- 0.09999999999999964 on `[0.1, 5]` -- and on `[1e-09, 1e+09]` it came back as 1e-07, inside the box but a hundredfold from where it was proposed, because the ulp being dropped belonged to the larger wall. The wave is now carried as the distance travelled from the lower wall, and the folded value is clipped into the box last, which is what the second rounding needs: on a log scale the theta-to-u round trip leaves the box on its own, since `10 ** log10(20)` is 20.000000000000004. The fold may now return a bound exactly, which the constructor has always accepted; a NaN still raises there, as before. Differential evolution met the first rounding, on `k_elim` against its lower bound of 0.1 in tutorial lesson 25 (#703), where a mid-run step left the box by an ulp. The second is worse, because on a log scale a bound only has to be *reached*, not overshot: a bounded optimizer works in sampling space over `[to_sampling_space(lb), to_sampling_space(ub)]` and projects every iterate into that box, so a bound that is active at the optimum arrives back at the parameter as the image of the box's own corner -- which lands outside the box for 91 of the 240 walls the tests now sweep, 86 of them then refused. That is how `job_type = profile_likelihood` met it (#750): the inner re-optimization of the nuisance parameters at a grid point lands on an active bound as a matter of course, so the job could not finish on any problem where some parameter pressed its box -- and reaching a bound is a documented profile outcome (a practically non-identifiable parameter), which was raising instead of being recorded. - **Tutorial lesson 25 no longer promises rates its one curve cannot settle, and its check no longer passes on the luck of one seed (#703).** The lesson fits a transit-compartment model to a single plasma curve, and its README said island DE "recovers all three" rates. The curve's objective has two minima: the true `k_transit` = 12.76 and `k_abs` = 9.11, and 10.09 and 14.65, whose absorption delay has the same mean to 0.01% and the same variance to 0.1%, so the second pair's curve never strays more than 0.0002 from the first on a peak of 1.7. The true pair scores lower (1e-9 or less, against 2.2e-8) only because the data are noise-free, and a population settles into one basin while its best objective is still ten to a thousand times that difference: the lesson's fit reached the true pair from 23 of seeds 1 to 60 and 1234. Its check ran from seed 1234 alone, which happens to be one of them, so any change to differential evolution's draws could fail it, and #698 and #700 each did. The README now teaches what one plasma curve can and cannot determine, and what data would separate the rates: with a second experiment whose dose starts in the absorption compartment the objective has a single minimum. The check accepts either minimum, `k_elim` within 3% in both, from seeds 1234, 1 and 2 (the first reaches the true pair and the other two the second), and a new default-tier test holds it to rejecting a fit between the minima, such as the early stop of #648. Tutorial checks can now list other minima and several seeds. The conf now polishes with `refine_method = trf`: DE can stall in the narrow valley the minima lie in, and from 2 of the 61 seeds the Simplex polish stopped 10% and 15% away from both; with `trf` every seed ends within 0.03% of one. No island setup that was measured made the true pair dependable. Islands of ten kept apart stall on their own under the default crossover; crossing with the target (`de_cross_with_target = 1`, #700) lets them converge, and with eight islands and a single migration late in 160 generations, four times the budget, 37 of 39 seeds reached the true pair. - **Differential evolution no longer proposes a candidate that is an exact copy of the parameter set it was built from, under `edition = 2` (#698, ADR-0143).** A candidate is its base with each parameter moved by the donors' difference with probability `mutation_rate`, and nothing guaranteed that any parameter moved: at the default rate one candidate in eight on a three-parameter model was a copy. Under the default seed policy a copy scores exactly what its base scored and takes the place of any worse member, and in `ade` the copies built more copies until the population was one parameter set and the run stopped, reporting convergence; on the stochastic recovery benchmark that stopped 24 of 30 `ade` fits early. One parameter chosen at random is now always mutated, as binomial crossover guarantees, and because PyBNF's candidates share values with the members they are built from, the choice is made again among the parameters the donors' difference actually moves when it does not move the one drawn, and other donors are drawn when it moves none. On the benchmark no fit stops early any more, and `ade` succeeds from 12 of 30 seeds against 4 (8 seeds only with the guarantee, none only without); `de`, which never collapsed, goes from 6 to 8, within noise. The learned mutation settings (`de_adapt_mutation = 1`) now use the same guarantee, which removes the early stops they still had. The legacy edition keeps the original proposal, draw for draw, and `de_force_mutation` sets the guarantee either way. Tutorial lessons 24, 25 and 44 run more iterations, since at their old budgets the new proposal left the fits of their recovery checks just short of the documented values. - **Scatter search no longer archives a lucky draw as a local minimum (#660 step 3, ADR-0136).** Scatter search decides everything by ranking, and for a stochastic model each objective value is one draw. The reference set stored that draw as fact, so a member whose value was a lucky draw could never be beaten by an honest child, its stuck counter climbed, and it was retired into the archive of local minima as one it never was, which is the list the run reports as its best results. When running a parameter set again would give a different answer, every reference member is now ranked on the mean of its draws, the draw-to-draw spread of the objective is pooled across the fit, and a decision the spread leaves in doubt (a child against its parent, or two neighbours whose rank gap sets a combination's step size) is not made: both sides are simulated again at fresh seeds, up to `ss_noise_max_draws` draws each, and the decision waits. A member counted stuck is drawn again too, so a lucky value regresses to the truth. The pooled spread and the separation test live beside the CMA-ES measurement in `pybnf.algorithms.noise_handling`. `ss_noise_handling = 0` turns it off; a deterministic fit is unchanged. - **CMA-ES no longer ranks a stochastic model's population on single noisy simulations (#661, ADR-0135).** CMA-ES reads only the ordering of its population, and for a stochastic model each objective value is one draw, so when the noise was comparable to the real differences between candidates the ordering was partly random: the distribution update was pulled in arbitrary directions and the step-size adaptation, reading noise as stagnation, shrank a step that should not have shrunk. When running a parameter set again would give a different answer, each generation now ends by simulating a few of its candidates again at fresh seeds and measuring how far they move in the ranking against what pure noise would do, the uncertainty handling of Hansen, Niederberger, Guzzella and Koumoutsakos (IEEE Transactions on Evolutionary Computation 13(1), 180-197, 2009). While the ranking is unreliable every candidate is simulated more times, up to `cmaes_noise_max_evals`, and ranked on its average, and the step size is held up; when it is reliable the extra simulations are dropped again. `cmaes_noise_handling = 0` turns it off. A deterministic fit is unchanged. The measurement lives in its own module so scatter search can reuse it (#660). - **The information criteria of a stochastic fit no longer come from a single simulation (#676).** `Results/information_criteria.txt` reports AIC, BIC and AICc from the log-likelihood of the best fit, and that log-likelihood came from simulating the best fit once more at the end of the run. For a stochastic model that one simulation is a draw, so the reported AIC was a noisy number that changed from one run to the next for the same parameter set, and under the default seed policy the draw it reported was the very trajectory the search had scored. Since the best-fit confirmation stage above landed, the averaged objective value in `Results/best_fit_confirmation.txt` and the single-draw log-likelihood in `Results/information_criteria.txt` also disagreed about the same parameter set, with nothing to say which to trust. At the end of a stochastic fit the best fit is now simulated `best_fit_replicates` times at fresh seeds, all submitted at once, and the log-likelihood behind the criteria is the mean over those runs. The file gains a `replicates` line saying how many runs that was and a `log_likelihood_standard_error` line saying how far they spread; AIC, BIC and AICc each carry twice that, so two models whose AIC values differ by less than it have not been told apart by the runs. The mean of the log-likelihoods is the same average the confirmation stage reports for the objective value, so the two files now describe the same parameter set the same way (ADR-0131). A deterministic fit, a legacy-edition fit, a fit whose seed policy pins every trajectory, a fit whose wall-time budget is spent, and the periodic checkpoint all keep their single simulation; the two new lines then say `1` and `n/a`. - **A bootstrap replicate of a multiple-shooting fit no longer reports a start belonging to the replicate before it.** A bootstrap run reuses the algorithm object across replicates, and `job_type = ms` kept adding each start's ladder result to a list that was never cleared. Everything that reports on the ladder reads that list, so the second replicate could pick a start from the first, fitted to different resampled data, as the one behind `Results/continuity_defects.txt` and the best stage trace. Found while adding the per-start summary above, which reads the same list. - **The startup parallelism report now describes how busy a fit will actually be, rather than how busy its first round is (#655).** The report added in v1.8.0 measured "how many jobs the fit runs at once" from the first batch of jobs submitted. For scatter search that first batch is the initialization round, which is `init_size` parameter sets (default ten per free parameter) and has nothing to do with the population. Every round after it runs `population_size` x (`population_size` - 1) simulations. A fit with seven free parameters and `population_size = 20` on 384 processors was told 314 processors would sit idle and that it should lower `population_size`, while it was in fact about to run 380 simulations at a time. Following that advice would have reduced the number of processors the fit could use. An algorithm can now state how many parameter sets it keeps out for evaluation once it is under way, and the report uses that number. Scatter search reports its population pairs, which is the same number it already prints as "simulations per iteration". Profile likelihood reports its directional tracks, since its first batch is a single preflight evaluation. When the first round differs from the steady state, the report says so instead of warning about it. The count is scaled by `smoothing` and `parallelize_models`, which turn one parameter set into several jobs. The advice also names the setting each fit actually reads. Powell and simplex are told about `n_starts` rather than `population_size`, which they do not use, and profile likelihood, whose concurrency follows how many parameters it profiles, is advised only about how many processors to reserve. `init_size` is now documented in the cluster guide, with a note that raising it fills a large allocation during scatter search's initialization round. The table of how many simulations each fit type runs at once listed profile likelihood as `pl`, which is not a valid `job_type`. It now reads `profile_likelihood`. ## [v1.8.1] - 2026-08-23 ### Fixed - **The CMA-ES stagnation restart now fires, because its tolerance is calibrated from the objective rather than inherited from a step length (#653, ADR-0128).** ADR-0106 gave `cmaes_tolfun` its own key precisely because it is a range in objective units while `cmaes_stop_tol` is a step length in sampling space, and then had an unset `cmaes_tolfun` fall back to it anyway. `cmaes_stop_tol` defaults to 1e-11, so the stagnation range was 1e-11 in objective units, which on an objective of any ordinary magnitude never fires. This is #648 in the mirror. There a ratio read as a range was far too loose and stopped fits early at a wrong answer. Here a step length read as a range is far too strict, and what it costs is the restart trigger the battery exists for: without it a run polishes a local basin and never yields to a restart, which is the failure the battery was built to prevent. The documentation had already told readers the default was "rarely what you want", which described the defect rather than fixing it. An unset `cmaes_tolfun` is now 1e-11 times the objective spread measured across the first generation's population. That spread is a real measure of how much the objective varies over the search box, it is in the units the tolerance needs, and it is taken before anything has converged, so it does not drift with the objective the way a fraction of the current value would. It is calibrated once and every restart reuses it, so a late restart is not held to a stricter bar than an early one. The fraction is chosen so a problem whose initial population spans one objective unit gets exactly the 1e-11 this key always defaulted to, leaving a reference-scaled problem unchanged. The run logs the value it picked. A population that cannot supply a spread keeps the old fallback rather than inventing a number, and an explicit `cmaes_tolfun` is never touched. Only affects `cmaes_restarts > 0`, which is not the default. ## [v1.8.0] - 2026-08-23 ### Added - **A cluster fit now says how many of the reserved processors it is actually using (#621).** How many simulations a fit runs at once follows from its settings, mainly `population_size`, and not from how many processors were reserved. When the two do not match, the extra processors can sit idle for the whole run and nothing said so, so a user could reserve several machines and quietly use a fraction of one. After the first set of jobs is submitted, PyBNF now compares the number of jobs running with the number of workers that connected, writes both numbers to the log so a finished run can be checked afterwards, and warns when the two differ by a large margin in either direction. The warning names both numbers and points at `population_size`. Differential evolution, CMA-ES and scatter search finish a whole generation before starting the next, so some idle time near the end of each generation is expected with them. Those three say so in the message, to save a user looking for a fault that is not there. Only cluster runs are reported, because a local run's worker count is exactly the number that was asked for. Reading the worker count never stops a fit if it fails. - **On a cluster whose machines are not all the same size, each machine now runs a worker per CPU it was granted rather than one count for all of them (#617, ADR-0124).** When PyBNF started workers across several machines it sent the same worker count to every one. On a cluster whose machines differ in size, that single count is too many workers for a small machine, where they compete for its processors, and too few for a large one, which sits partly idle. The `srun` launcher (`-t slurm-srun`) now works out each machine's count from what SLURM granted it and starts that many there. A single worker count cannot express this, because one `srun` job step binds every machine in it to the same number of CPUs, so PyBNF starts one job step per distinct machine size, each on its own machines, which SLURM runs at the same time. A run on two 40-processor machines and one 96-processor machine starts 40 workers on each of the first two and 96 on the third. The arrangement is written to the log, and each job step writes its own worker log (`dask_workers.log`, `dask_workers_2.log` and so on) so their output does not interleave. An allocation whose machines are all the same size, which is the common case, still runs as one job step exactly as before. If SLURM does not report a per-machine list PyBNF can line up with its machines, it falls back to sizing every machine the same, with a warning saying so. This applies only to the `srun` launcher. The SSH launcher (`-t slurm`) still uses one count for all machines, because `dask ssh` takes only one. Setting `parallel_count` still splits that total evenly across same-size machines on either launcher; on machines of different sizes the `srun` launcher now splits it in proportion to machine size (#643, see below). - **A multi-machine fit can start its workers without logging in anywhere, so clusters that use host-based or Kerberos SSH can run one at all (#614, ADR-0122).** PyBNF had exactly one way to run across several machines: `dask-ssh`, which logs in to every node with **paramiko** rather than with the operating system's `ssh`. paramiko offers a public key or a password and nothing else. A cluster whose nodes authenticate to each other by **host-based** SSH (no user credential exists to offer) or by **Kerberos/GSSAPI** (paramiko can, dask never enables it) therefore refused every login PyBNF attempted, on a machine where `ssh OTHER hostname` from the same shell succeeds — and the advice in `docs/cluster.rst`, to set up SSH keys, could not fix a cluster that is not asking for a key. The run stopped before a single simulation. **`cluster_type = slurm-srun`** (or `pybnf -t slurm-srun`) starts the workers with SLURM's own `srun` instead. No credential is involved because the scheduler granted the allocation before PyBNF started: PyBNF runs a dask scheduler on its own node, has `srun` place one worker process group on each node of the allocation, and connects through the scheduler file the scheduler writes — a file it now also *creates*, so under this launcher `scheduler_file` chooses where it goes (default `dask_scheduler.json` in `output_dir`) rather than naming a cluster to attach to. `-t slurm` and every other path construct exactly the calls they did before; this is a second launcher, not a change to the existing one. Both bring-up steps wait on a real signal rather than a fixed sleep, and watch the process they started while waiting: the scheduler is ready when its connection file *parses* with an address (dask writes that file in place, so a reader can catch it half-written), and the workers are ready when one has registered. That second check is the one that cannot be dropped — connecting to our own scheduler always succeeds, so a failed placement would otherwise turn into a fit that submits jobs and never gets one back, with `srun`'s explanation sitting unread. Both failures quote that output. Running outside an allocation is refused up front, since `srun` there does not place a task but submits a job and waits, which would read as PyBNF hanging. The workers are placed one task per node with the CPUs their processes need, taken from what SLURM granted (`$SLURM_CPUS_ON_NODE`): under cgroup binding a task that took the default single CPU would confine every worker it forked to that one CPU and quietly serialize the node. Two logs, `dask_scheduler.log` and `dask_workers.log`, are written to the output directory. - **A fit's start point is a supported, validated, per-parameter fact that every optimizer reads, and the resolved start is recorded beside the results (#583, #559, ADR-0117).** There was no supported way to say "start this fit at exactly this point, inside the declared box." Every route failed, and every one failed **silently** — the fit ran, converged, and reported a plausible number from a point that was not the one that was asked for. The new **`start_point = `** key says it directly: one line per parameter, at every edition, alongside the legacy `*_var` declarations, in the parameter's own units whatever its scale. The edition-2 `parameter:` record's `initial_value:` field now means the same thing and resolves through the same carrier. This closes a defect, not just a gap. `initial_value` has existed since ADR-0043 and is honored by twelve fit_types — but `StartPointOptimizer._resolve_start_pset` read `FreeParameter.value` on none of its branches, so `cmaes` / `powell` / `sim` / `gntr` / `lbfgs` / `trf` / `ms` — exactly the optimizers both issues are about — discarded it and started at the box centre. Declared at `k = 0.3, S0 = 100`, a fit started at `k = 1.505, S0 = 89.44`. ADR-0043's own text claims `initial_value` is "respected in every algorithm family"; that held only for the no-prior form, and was false the instant a bound or prior was present. A start point is **partial by design** — name only what you want pinned — and for a multi-start fit it pins start 0 while the rest stay independent draws, so the scatter survives. An out-of-box value is **refused, not moved**: PyBNF never clamped such a value to the nearest bound, it *reflected* it back inside with a periodic triangle-wave fold, so `250` into `[1, 100]` became `50` and the fit proceeded from an arbitrary interior point with a `DEBUG` line. Every refusal is a proper configuration error rather than the bare `OutOfBoundsException`, which reached the user as "an unknown error … please report this bug" — including on PEtab import, where an out-of-box `nominalValue` crashed that way on a file the user never wrote. **`Results/start_point.txt`** records where the run actually began — value, source, and the declared box, per parameter — before anything is scored. That artifact alone would have caught all four of the reported failures, and it is the only way to recover a CMA-ES start at all, since that start seeds the distribution mean and is never itself evaluated. The PEtab importer now emits `start_point` from `nominalValue`, a mapping ADR-0043's field table advertised and neither direction implemented. The two spellings are synonyms **except on `profile_likelihood`**, where a complete `initial_value:` specification keeps its established meaning — those values are θ\*, so the polish is skipped. `start_point` never carries that meaning; it names where the polish starts. The distinction is load-bearing: a `nominalValue` is not a claim of optimality, and unifying the two would have made every PEtab problem with a full `nominalValue` column profile around the nominal point rather than the optimum, putting every confidence bound in the wrong place without a word. - **A `model:` declaration carries that model's own CVODE tolerances, so a per-species absolute tolerance can finally be written by hand (#586, ADR-0116).** `sbml_atol` and `sbml_rtol` are one key each over *every* SBML/Antimony model in a fit, and that — not the grammar — is why neither could ever take a vector: a positional one is ordered against a species list a conf author cannot see, and a species-keyed one has no reading across models that do not share species names. Give the statement one model to be about and both objections disappear. Under `edition >= 2` a single-model `model:` line now takes three optional labeled fields, in any order: - **`atol:`** takes exactly what `sbml_atol` takes (a number, `auto`, or `tracking [decades]`), for that model alone. - **`rtol:`** is the per-model `sbml_rtol`. A scalar, and only ever a scalar — CVODE takes one relative tolerance and there is no per-species one to route. - **`species_atol:`** is the hand-written vector, written as ` ` pairs: `model: weber.xml, atol: auto, species_atol: PKD 1e-3, CERT 1e-2, rtol: 1e-9`. **The map is a set of exceptions, not a replacement.** Every species it names takes the number stated, verbatim — neither ADR-0103's `1e-8` ceiling nor ADR-0105's model-scalar floor binds a number a person wrote, for the same reason a pinned `sbml_atol` number is not clamped either. Every species it does *not* name keeps whatever it would have had: the derived vector, a pinned scalar broadcast, or the backend default. The steady-state cutoff stays the model-wide number, and under `atol: tracking` the map becomes the ceiling the trajectory-following weights sit beneath. **The species names are bngsim's, and an unknown one is an error at config load** that lists the names the model does have. That matters twice: bngsim renames a species id colliding with an Antimony reserved word (`NULL` integrates as `_ant_NULL`, and the error says so), and a parameter driven by a rate rule becomes a state the `` never mentions — nameable here and reachable by nothing else. On `Smith_BMCSystBiol2013`, the one subset-I slug where the two name lists disagree, the derived vector declines outright and a hand-written map is the only route to `CVodeSVtolerances` there. **Nothing changes without one of those fields.** `sbml_atol`/`sbml_rtol` keep every meaning and remain the fit-wide default; a job that states nothing gets the same `Simulator.run` call it got before, argument for argument. A `model:` line carrying a tolerance field must declare exactly one model, and that model must be an `.xml`/`.ant` one on `sbml_backend = bngsim` integrated by CVODE. A BNGL model states `atol`/`rtol` in its own `begin actions` block, the RoadRunner backend has its own integrator settings, and `sbml_integrator = gillespie` runs every action stochastically — all three are refused rather than accepted and ignored, as is a `species_atol` on a bngsim without lanl/bngsim#196. Measured on `Weber_BMC2015`: writing bngsim's own `rtol * y_i` rule by hand — the rule ADR-0105 measured and declined to apply automatically — costs 192 integrator steps against `auto`'s 184 and the clamped default's 210, agreeing with all of them to ~3e-08 at the final state. The derivation is still the better default; what it could not do was be overridden for one species of one model. - **The information criteria are checkpointed alongside the parameter sets, so a run is scoreable before it ends (#560).** `sorted_params_backup.txt` has been written throughout a run since forever; `information_criteria.txt` was written **only** on the terminal path. Every downstream consumer of a fit's result needs both, so a run was un-scoreable at every moment except its last — even though the best parameter vector had been on disk the whole time and the number was fully determined long before the process exited. That gap is not cosmetic: `log_likelihood` in `information_criteria.txt` is the only place PyBNF reports the **full normalized** log-likelihood (every dropped per-point constant restored), while the minimized `Obj` column of the parameter table is the *reduced* objective, so an absolute AIC/BIC — or any benchmark score built on one — could not be computed from the checkpoint alone. Each checkpoint now also writes `Results/information_criteria_backup.txt` (and `information_criteria_refine_backup.txt` during a refine, mirroring the parameter file beside it). Same format as the final artifact, differing only in `#` comments that mark it a snapshot and name the parameter set it describes, so one parser reads either file. The final `information_criteria.txt` is unchanged, in name, content, and meaning. Cost is one extra simulation per checkpoint, and only when the checkpoint has something new to say. Nothing is spent while the best fit is unchanged — the file on disk already describes it, which is exactly the long converged tail that motivated the issue: a `gntr` fit of `Brannmark_JBC2010` (100 starts × 1000 iterations) reached its final objective with ~40 minutes left, all of it spent waiting for a file rather than for an answer. Nothing is spent at all unless the objective is a proper likelihood, since no information criterion is defined for `sos` / `sod` / `norm_sos` / `kl` / `wasserstein` / `direct_pass`. Otherwise it is one simulation per `backup_every * population_size * smoothing` returns — 1 in 1000 at the common `backup_every = 10, population_size = 100`. The new key **`backup_information_criteria`** (default 1) turns it off for a model where even that is too expensive. A killed or crashed long run is worth what it should be, too: the parameters survived it already, and now the score does. - **`sbml_atol` takes `auto` and `tracking`, so a model whose own scale asks for a *looser* absolute tolerance can have one (#557, ADR-0114).** ADR-0103's derivation is allowed only to tighten. That is a no-regression rule and it is why the derivation could be applied to every model without a flag — but a model whose species all sit far above one has a real, computable tolerance need, and the derivation computes it and then discards it. `Weber_BMC2015`'s seven species live at `1.24e+02 .. 4.21e+07`, so it asks for `4.665e-03` and is handed `1e-08`, 5.7 decades tighter. ADR-0105's per-species vector cannot rescue that either, and the way it fails is the point: each entry clamps into `[scalar_atol, default_atol]`, the scalar has itself been clamped to `default_atol`, that interval collapses to a point, and the vector correctly declines — so the per-species mechanism silently declines exactly the models that span the most decades. Measured over the subset-I corpus, the clamp binds on **10 of the 22 slugs with a readable nominal state**. - **`sbml_atol = auto`** lifts the ceiling on both derivations and nothing else. ADR-0105's *floor* — no species resolved below the model's own scalar — stays exactly where its measurement put it (releasing it killed 91 of 100 `Brannmark_JBC2010` box points against 39), and the `1e-16` floor stays too. What is left, `rtol * max(y_i, median)`, is bngsim's own `derive_atol` with the model's median as its `floor`, and is asserted against the library rather than against a second copy of the arithmetic. - **`sbml_atol = tracking [decades]`** wires lanl/bngsim#213's `CVodeWFtolerances` — an absolute tolerance re-evaluated against the state being integrated, so a species that starts at order one and decays to nothing keeps a tolerance that means something for it. This is the half ADR-0105 named as out of reach. The ceiling is `auto`'s vector, held from the **nominal** state, so `tracking 0` is `auto` exactly and a fit that moves initial conditions does not move its own tolerance; bngsim's bare `"auto"` ceiling, which re-derives from the live state at every run, is deliberately not used. An unstated depth is left to bngsim rather than copied. **Nothing changes without one of those two words.** Unset is byte-identical — same clamps, same vector, same steady-state pairing — and a number remains the documented off-switch that pins `CVodeSStolerances` ulp for ulp. `tracking` on a bngsim without the capability is *refused*, at config load and again in the model constructor, rather than silently integrating at something else. **One thing #557 claims does not reproduce, and it is recorded rather than repeated.** The issue's headline is that Weber integrates 6 of 30 sensitivity-applied box points at the clamped tolerance, against 22 at `1e-04`. On the current stack all six arms — unset, `1e-08`, `1e-04`, `auto`, `tracking`, `tracking 6` — integrate **30 of 30**, within 7.2 s of each other. bngsim moved from 0.12.2 to 0.13.0 in between, and the documented Weber-specific change is lanl/bngsim#305, whose own entry measures that slug's `t = 24` crossing and reports the step count roughly halving. So this ships as a capability with a measured cost rather than as a rescue, and the instrument that discriminates is **integrator steps** rather than pass/fail. Over 20 box points per slug with the gradient sensitivity request applied, `auto` costs 0.38x–0.67x the CVODE steps on the six slugs the clamp binds hardest (`Perelson` 8 440 → 3 216; `Weber` 78 641 → 36 154; `Laske` 236 629 → 158 543), is bit-identical on `Giordano_Nature2020` — the control, whose derivation tightens and so cannot see the ceiling — and moves `J_paper` at the PEtab nominal point in its sixth decimal. It is not free in the other direction either: error-test failures rise as steps fall, and `Laske` lost one box point of twenty at one seed (both arms lose one at a second seed). ADR-0114 has every arm and says what it has not bisected. - **Measurement-time uncertainty, phase 2: a marginalized-time likelihood can now be fit by `job_type = lbfgs` and `job_type = gntr` (#588, ADR-0113).** Phase 1, below, refuses every gradient job type, because the gradient of the marginal `-log z_k` did not exist yet. Phase 2 supplies it and lifts the refusal for the two job types that can consume it. Nothing is added to the model file and the differential equations are not augmented. The published method augments the model with one extra state per observation, because the solver it uses returns sensitivities of model states only. PyBNF does not need that. Its forward-sensitivity engine already stores `∂y(τ)/∂θ` at every node of the same trajectory phase 1 integrates over, so `∂z_k/∂θ` is a second quadrature over a trajectory PyBNF already holds, reusing each noise family's existing derivative and the time prior's own `∂p/∂σ_t`. The trapezoid rule is a linear functional that does not depend on the parameters, so what comes out is the exact derivative of the number phase 1 reports rather than an approximation of it. A finite-difference check agrees to about 1e-9 in every column. `-log z_k` is the log of an integral rather than a sum of squares, so it takes the scalar gradient path. `lbfgs` consumes the scalar gradient and `gntr` builds its Gauss-Newton Fisher matrix from the per-observation scores. Three job types stay refused, each with its own reason: `trf` needs an exact residual vector and the log of an integral is not one, `hmc` runs on the analytical model rather than the simulator, and `ms` runs through the shooting layer. Integration error controlled by the solver, which is the other thing augmenting the equations would buy, is a separate question, and phase 2 keeps phase 1's fixed grid. Measured end to end on tutorial lesson 49 through the real simulation backend, `lbfgs` and `gntr` both recover `k = 1.02` against a truth of 1.0, with every multi-start reaching the same optimum. The lesson carries that arm as `marginal_gradient.conf`. - **Measurement-time uncertainty via posterior marginalization, phase 1 (#587, ADR-0112).** A new `time_error` clause on the `noise_model` line treats the latent sampling time as a random variable and *integrates it out* of the likelihood, instead of assuming each datum was collected at exactly its reported time — an assumption that biases estimates and makes posteriors overconfident when sampling times actually drift (handling delays, imperfect synchronization, reporting error). Written whole-fit as `noise_model = , = , time_error = truncated_normal, sigma_t = fit st__FREE` (or `uniform`; `sigma_t = fix_at `), it replaces the per-point likelihood with a `MarginalizedTimeObjective` whose per-observation contribution is `−log ∫ p(ȳ_k | y(τ)) p(τ | t_k) dτ` — the `n_t`-dimensional marginal factorizes into one-dimensional integrals (the method of Vanhoefer, Nakonecnij, Binder & Hasenauer, bioRxiv 2026.05.09.724053; the temporal analogue of Raimúndez et al. 2023 nuisance marginalization). The search stays `n_θ` (+ one `σ_t`), not `n_θ + n_t`. Phase 1 evaluates each integral by log-space quadrature over the stored trajectory, reusing every noise family's normalized `log_density` (ADR-0056) as the integrand and the gradient-free optimizers/samplers (`de`/`pso`/`ss`/`mh`/ `dream`/…) unchanged — nothing is added to the model file. Edition-2 only. The `σ_t → 0` limit is the standard likelihood (a `fix_at 0` clause short-circuits to it). LOO/WAIC and `information_criteria.txt` work out of the box: the marginal per-observation `log z_k` **is** a normalized per-observation log-likelihood, so the objective reports it through the same `evaluate_pointwise` hook the per-point families use (`Σ_k log z_k = −score`), and an estimated `σ_t` is already counted in `k`. A marginalized time course is simulated on a **dense uniform grid** over the support (`t_end:` required, `t_start:`/`n_steps:` optional on the experiment line — decoupled from the sparse reported times, which only centre each timing prior), and `sigma_t = fit …` estimates the timing scale jointly (recognized as a declared nuisance). Worked end to end in **tutorial lesson 49** (`examples/tutorial/49_measurement_time_uncertainty/`): ignoring the timing spread biases the decay rate to `k ≈ 1.36` (truth 1), marginalizing recovers `k ≈ 1.06`, and estimating `σ_t` recovers `k` while re-discovering a non-zero timing error. Deferred and refused at build with a reason: a per-observable time prior, a prediction-dependent `σ`, the count family, and, in phase 1, every gradient `job_type`. Phase 2 above lifts that refusal for `lbfgs` and `gntr` in this same release, and does it by chaining the stored forward sensitivities rather than by the augmented equations phase 1 expected to need. `trf`, `hmc` and `ms` stay refused. `noise_profiling`, which *maximizes* a scale out, is refused as ill-defined alongside marginalization, which *integrates* the time out. - **Multiple shooting, `job_type = ms` (#563, ADR-0110).** The consumer of the constrained-transcription layer below, and the thing #563 was actually asking for. Each scored experiment's time course is cut at knots; segment *j* is integrated from its own start state — segment 0's is the model's own initial conditions, and each interior knot carries an auxiliary state that is searched, bounded and differentiated but is **never** a reported fit result — and continuity `Phi_j(z_j, theta) - z_{j+1} = 0` is enforced by an augmented Lagrangian whose subproblem is solved by `gntr`'s own Gauss-Newton trust-region step machine. Every reported score comes from discarding the auxiliary states, re-simulating theta with ordinary single shooting, and scoring *that*, so a run that leaves continuity unconverged scores as what it actually is — and every certified iterate lands in the ordinary trajectory at that score, so `sorted_params`, the best-fit simulations, the information criteria and the inference-data sidecar are produced by the same code every other `job_type` uses. Why it is worth the machinery: on `Borghans_BiophysChem1997` a correctly-shaped oscillator whose period is wrong by more than about 3 % scores *worse than fitting no dynamics at all*, so under single shooting the flat line is the ceiling on essentially the whole box and fifteen independent global searches terminate at it. Over one short segment a period error cannot accumulate: the information moves out of a residual term that saturates and into continuity defects, which carry a direction. The prototype's structural finding is what keeps the implementation small — a segment-start state is an `IC` route with chain-rule factor 1, so the existing gradient/Fisher assembly builds its column with no new residual math, and each segment is presented to it as an ordinary *experiment*. The continuity block is the only new assembly surface. Knots are named by their exact fraction of the horizon, so a coarser rung recognises a finer one's knots and the `4-2-1` ladder *continues* rather than reseeds. New keys `ms_segments`, `ms_coarsening`, `ms_penalty`, `ms_penalty_growth`, `ms_max_penalty`, `ms_feasibility_tol`, `ms_optimality_tol`, `ms_inner_iterations`, `ms_aux_decades`, `ms_max_iterations`, defaulted from ADR-0109's measurements rather than from taste. Requires the bngsim backend — a knot carries the model's *state*, so both a generated network (`.net`) and an SBML/Antimony model are supported, through two backends that differ only in what a simulation returns: on the SBML path the columns an experiment scores and the columns a continuity row differences are the same columns, and on the `.net` path they are not, so that backend asks for the observable and species selector families together and one integration still serves both (#577). A network-free (NFsim) model enumerates no state and is refused. It also refuses, by name, a fit whose scored quantity is a function of a whole series — an analytic per-series scale, a data normalization, a cumulative-to-incident difference — since cutting the series would change it. An analytically profiled noise scale (ADR-0108) is deliberately fine: it is profiled over pooled residuals, so the segments pool the same ones, and the constraint terms never enter the likelihood. What is *not* claimed: that multiple shooting improves the typical fit (48 paired starts: 24–24, medians tied at every radius, at 2–7x the simulations), or that it solves Borghans from an uninformed start (0/24). The measured case is the tail and the robustness. Segment simulations run serially on the master in this cut; parallelising them, and the acceptance benchmark, are the follow-on work. - **A constrained-transcription layer, `pybnf.transcription` (#563, ADR-0109).** Infrastructure for restating a fit as a larger, better-conditioned problem with internal auxiliary variables and equality constraints that tie them back together — the reusable half of multiple shooting, which will be its first consumer. Four pieces: an **augmented variable layout** that carries the fit's reported free parameters and the transcription's internal blocks in one vector while keeping them rigorously apart (an auxiliary state is searched, bounded, and differentiated; it is never a reported fit result); an **equality residual/Jacobian interface** whose Jacobian is block-sparse with a condensing seam left open, and whose defects are scaled so one penalty means one thing across states of different magnitude; the augmented Lagrangian offered in all three forms PyBNF's optimizers consume (scalar for `lbfgs`, an *exact* stacked least-squares residual for `trf`, Gauss-Newton for `gntr`); and an **optimizer-agnostic augmented-Lagrangian outer loop** with a transcription homotopy and best-iterate certification through the ordinary single-shoot path. No behaviour change to any existing fit: nothing imports it yet, it defines no configuration key and no `job_type`, and it makes no simulator call — which is what lets the whole layer be verified offline (93 tests, ~1.4 s) against an equality-constrained quadratic whose multiplier is known analytically and a closed-form linear-ODE shooting problem measured against an independently computed single-shoot optimum. Three measurements from the #563 prototype are baked into the defaults rather than left as tuning advice, each contradicting the plan that preceded it: the penalty schedule starts **tight** (`rho0 = 10`, `gamma = 5` beat `0.1`/`3` on quality *and* halved the cost), the segment ladder is the **mechanism** rather than a later refinement and starts in the middle (`4-2-1`, not `8-4-2-1` — many short segments certified worse than their own start under partial observability), and a run reports its **best certified** iterate rather than its last (on one start the final stage held `-147.0` while an earlier iterate certified at `-196.3`). - **`noise_profiling = 1`: profile an estimated noise scale out of the search analytically (#562, ADR-0108).** ADR-0066 already profiles a declared column's optimal multiplicative **scale** out of the fit; this is the other half of the same classical trick. Every noise parameter declared `= fit ` is removed from the search and replaced, at each evaluation, by its closed-form maximum-likelihood value over the scored points that share it — the weighted residual RMS `sqrt(sum w r**2 / sum w)` for the Gaussian families (`normal` / `lognormal` / `lnnormal`), the weighted mean absolute residual `sum w |r| / sum w` for `laplace`. Opt-in; `0` (the default) is an exact no-op. Why it matters beyond the dimension count: at a random point in the box the sampled scale is nowhere near its optimum, so the `log sigma` term dominates and a global sampler ranks candidates mostly by *how wrong their sigma happens to be* rather than by how well their dynamics fit. On `Borghans_BiophysChem1997` every optimizer that can run it converges to the same attractor, and that attractor is exactly the **no-dynamics solution** (a flat line at the best constant with sigma at the residual RMS, `-51.204092` analytically and `-51.204092` reported). Across the Grein et al. 2026 subset-I corpus a plain free-parameter sigma accounts for **32 parameters in 13 of 23 slugs** — 4 % to 33 % of the search. A profiled scale also has no box to run into, so a fit can no longer optimize its sigma into an upper bound and absorb model misfit as "measurement noise" (`Schwen_PONE2015`'s `IR_obs_std`). The switch is all-or-nothing within a fit and is refused *before the run starts*, naming the reason, for anything without a closed form: a `formula` / `prediction_formula` / per-measurement sigma, a `student_t` `df`, the `neg_bin` dispersion, a `location = mean` prediction on a log scale, one free parameter serving as the scale of two different families, or a fit with nothing to profile. A **fixed** scale (a data column, `fix_at`, `relative`) is not searched, so it is simply left alone. Refused for the Bayesian samplers too: profiling *maximizes* the nuisance out where a posterior *integrates* it out, so the draws would not be posterior draws. Profiled parameters stay declared (the same `.conf` runs with and without the key; their bounds and prior become inert) and stay **estimated**, so they keep counting in `k` in `information_criteria.txt` — otherwise every AIC/BIC would shift between the two runs. Their fitted values are written to the new **`Results/profiled_noise.txt`** and echoed on the console, since a value the fit solves for rather than proposes is not a coordinate of the best parameter set and appears in no `sorted_params_*.txt` row. Gradient support comes free by the envelope theorem — `job_type = lbfgs` and `gntr` consume the exact scalar gradient with the sigma columns dropped, and no new forward sensitivity is needed. `job_type = trf` refuses a profiled fit (as it already refused a searched free scale): under profiling the least-squares residual norm is identically constant, so a trust-region residual model carries no information about the parameters. - **`Results/method_chain.json`: which methods a run actually executed (#564, ADR-0107).** Written by every run — budget or no budget — it carries the chain the conf requested (`requested_methods`), the chain that ran (`executed_methods`), and one entry per phase (the fit, the refine, each bootstrap replicate) with its status (`completed` / `wall_time_expired` / `skipped`), its stop reason, its elapsed seconds, its completed simulations, and the best objective it reached. `requested_methods` longer than `executed_methods` is the machine-readable form of a downgrade, so a scoring harness can assert on the method it measured in one line instead of parsing stdout. A `bootstrap` phase records `replicates_requested` / `replicates_completed`, because `bootstrap = 30` in a conf is worth nothing if the budget stopped the run at 11. The file is rewritten after every phase (so a run whose refine raises still leaves the record of the fit that happened), is strictly valid JSON (a non-finite objective is recorded as `null`, never `Infinity`), and — like `stop_reason.txt` and `information_criteria.txt` — a failure to write it is logged and swallowed rather than taking a finished run down. ### Fixed - **Setting `parallel_count` on a cluster whose machines are not all the same size no longer stops the `srun` launcher from starting any worker (#643, ADR-0126).** `parallel_count` gives a total number of workers over all machines. The `srun` launcher (`-t slurm-srun`) split that total evenly and asked SLURM for the per-machine share on every machine in one job step. On a mixed allocation the even share can be more than a smaller machine was granted, so SLURM refused the step with `srun: error: Unable to create step for job NNNNN: More processors requested than permitted`, and no worker started. For example, on a 96-CPU and a 40-CPU machine with `parallel_count = 136` the even share is 68 per machine, and 68 is more than the 40-CPU machine holds. The total is now split in proportion to each machine's granted CPUs, one job step per distinct size, the same way the automatic (unset `parallel_count`) sizing already works, so each step asks only for what its machines hold and none is refused. On the 96-CPU and 40-CPU example that is 96 workers on the larger machine and 40 on the smaller. An allocation whose machines are all the same size is unchanged: `parallel_count` is still split evenly in one step. The SSH launcher (`-t slurm`) is unchanged, because `dask ssh` takes one worker count for all machines. - **A multi-machine fit started from the shell `salloc` opens no longer asks SLURM for more processors than it granted (#642, ADR-0125).** On many clusters `salloc` returns a shell on the **login node** while the allocation is held on a compute node, and that shell is exactly where PyBNF is meant to run: it holds the allocation. Started there, the `srun` launcher (`-t slurm-srun`) stopped at once with `srun: error: Unable to create step for job NNNNN: More processors requested than permitted`, and no worker ever started. PyBNF was sizing the run by `$SLURM_CPUS_ON_NODE`, which SLURM sets only inside a job step running on an allocated node, so on the login node it is absent — and the two remaining numbers describe the machine asking, which there is the login node, not in the allocation and usually several times larger. A job granted 20 CPUs was therefore sized as though it held 128, and 128 is a request SLURM refuses. The count now falls back to `$SLURM_JOB_CPUS_PER_NODE`, the per-node list SLURM publishes for the **job**, which is set correctly in that shell (its smallest entry, since one number has to be acceptable on every machine in the step). Inside the allocation nothing changes: `$SLURM_CPUS_ON_NODE` is still preferred where SLURM sets it, and the same `srun` command is built. The default path also hands the per-machine counts it already read straight to the command rather than having it read the environment a second time, so a stale `$SLURM_CPUS_ON_NODE` — including one exported by hand as the workaround for this bug — no longer sizes a later run. The SSH launcher (`-t slurm`) reads the same count, so it too stops starting a login node's worth of worker processes on each machine when it is launched from the login node. - **When the workers cannot be started, the message says what went wrong and what to try instead (#618).** A multi-machine run whose workers failed to start stopped with `Failed to start the dask-ssh cluster (dask-ssh exited with code 1)` and, on the cluster this was reported from, nothing else. The real cause was that the login to the other machines had failed, and no part of the message said so, named a cause, or named a way of running that needs no login. Stopping was right — carrying on with fewer machines than were asked for wastes the whole run — but a hard stop makes that message the entirety of what the user gets. It was empty because the half of dask's output that explains the failure was discarded twice over. `dask ssh` prints its own account of a refused login — the node it was connecting to, and the exception paramiko raised — to **stdout**, and lets only the traceback fall to stderr; PyBNF captured stderr and sent stdout to `DEVNULL`. And dask ends a failed bring-up with `os._exit(1)`, which does not flush Python's buffers, so its few hundred bytes of stdout never reached the 8 KB that would have forced a write to the file. Measured against dask 2026.7.1 on a login that fails: **0** of dask's own lines survived; **15** survive now that PyBNF captures both streams into one file and runs dask unbuffered. The message now quotes what dask said, and says so plainly when there was nothing to quote rather than falling back to "Check the cluster log directory" without naming a directory. The traceback frames are folded out of it — a failed login writes one traceback per node per retry, three retries each, and the sentences that say what happened are buried in dask's and paramiko's own source: **137** captured lines became **32**, losing none of those sentences. The log still keeps every line. When the output reads as a refused credential — "Authentication failed", "No authentication methods available", an encrypted key, a host key that did not match — the message says the login is the likely cause and says what PyBNF logs in with: paramiko, which can offer a public key or a typed password and nothing else, so a cluster that authenticates its nodes to each other by host-based or Kerberos (GSSAPI) SSH refuses it however it is configured, `ssh` from the same shell succeeds anyway, and `ssh-keygen` cannot help. A machine that could not be reached at all is deliberately *not* answered that way. Whatever the cause, the message names both ways of running on several machines that need no login: `cluster_type = slurm-srun` (#614), which starts the workers inside the allocation SLURM already granted, and a `scheduler_file` naming a cluster that is already up. `docs/cluster.rst` and `docs/troubleshooting.rst` now say the same. - **The SSH cluster launcher waits for its workers to register instead of sleeping for ten seconds (#398).** `cluster_type = slurm` slept ten seconds after starting `dask ssh`, on the assumption that the workers were up by then, and ten more after asking the cluster to stop, on the assumption that it had stopped by then. Both assumptions were wrong in both directions. Measured on a real SLURM cluster the workers took 26 to 59 seconds to register, so the fit began before its cluster was ready, and on a fast cluster ten seconds is longer than needed and every run paid it twice. Startup now polls the scheduler until every expected worker has registered, up to a time limit, and watches the `dask ssh` process on every pass, so a failed login is reported the moment `dask ssh` exits, quoting what it said, rather than after a fixed wait. Requiring the full worker count also turns a quietly undersized cluster into a clear error, which is what #200 asked for. Teardown asks each process it started to stop, waits until it has actually exited, and kills one that will not stop within a bounded time, so teardown returns as soon as the processes are really gone. On a real two-node run that took about one second rather than ten. The two limits are named constants in `pybnf/cluster.py`, `SSH_WORKER_TIMEOUT` and `TEARDOWN_TIMEOUT`, matching how the `srun` launcher sets its own. - **The cluster tests now notice when an outside program is renamed (#619).** The tests for starting a cluster checked that PyBNF built a particular command, against a copy of that command written into the test file. Nothing checked that the command could be run. When distributed stopped installing `dask-ssh` (#615), every test kept passing while every real multi-machine run died on `FileNotFoundError` — and because the outdated name was written in as the expected answer, *correcting* PyBNF would have read as a test failure. A handful of checks now ask the installed programs themselves. They take each command from the code that builds it for a real fit — no argument list is written down a second time — and confirm that the dask command line interface PyBNF invokes runs, that it still has the `ssh`, `scheduler` and `worker` subcommands, and that its `--help` still declares every option PyBNF passes it, so a renamed *option* fails as loudly as a renamed command. The `ssh` command is checked the stricter way, by handing dask the whole command PyBNF builds and letting dask parse it. The same questions are asked of `srun` and `scontrol`, skipped wherever SLURM is not installed — which is every developer machine, but not the clusters where PyBNF's tests are also run, and where a renamed `srun` option is worth catching before a fit walks into it. Nothing was taken away: the existing tests go on pinning the argument lists PyBNF is supposed to build, and the one copy of the dask invocation they keep is now compared against `cluster.DASK_CLI`, so it cannot drift away from the original unnoticed either. - **A multi-machine fit sizes its worker pool by what the job was granted, not by how big the machine is (#616).** PyBNF decided how many worker processes to start on each node by calling `multiprocessing.cpu_count()`, which reports every processor the machine has whatever the job scheduler granted. On a cluster where this was measured, a job that asked for **4** CPUs was told the node had **128**: PyBNF would have started one worker per processor and overshot the job's real capacity **32-fold**. Every worker is a separate process, so that multiplies memory use and leaves the workers competing for the same four CPUs — a fit that runs slower than it would have on the share it was given, or that runs out of memory. The defect could hide because a job that asks for *whole* nodes gets the right answer by coincidence: there the two numbers are equal. Both launchers now take the count from `Cluster.cpus_per_node`, the one place that decides it, which prefers **`$SLURM_CPUS_ON_NODE`** — what the allocation granted, and the only one of the three numbers that describes the *allocation* rather than the process asking, so it is still right for a worker started on another machine — then **`dask.system.CPU_COUNT`**, the machine's processors narrowed by CPU affinity and by any cgroup quota, which is what a single-machine run already sizes itself by, and only then the machine's whole processor count, which is correct only when nothing is limiting the job. The count and **which of the three it came from** are written to the log at the start of the run, so an unexpected number of workers can be traced to the number PyBNF believed; setting `parallel_count` still overrides all of it, and the log then names that key as the source. `-t slurm-srun`, which already read what SLURM granted, is unchanged apart from logging the source. - **Multi-machine fits run the command dask actually installs, so a cluster run gets past its first second (#615).** PyBNF started remote workers by running **`dask-ssh`** — one of three standalone scripts (with `dask-scheduler` and `dask-worker`) that distributed stopped installing in **2026.6.0**. On any current install, every run that used more than one machine died immediately with `FileNotFoundError: dask-ssh`, before a single simulation, on both routes through that code: `-t slurm` and a hand-set `scheduler_node` / `worker_nodes`. The same feature is now a subcommand of the unified `dask` program, and PyBNF runs it as **`dask ssh`**. The subcommand form is not a version trade: `ssh`, `scheduler` and `worker` are registered in the `dask_cli` entry point group in **2024.1.0**, the oldest dask/distributed `pyproject.toml` allows, so the fix works across the whole supported range and no floor had to move. PyBNF invokes it as ` -m dask ssh` rather than as a bare `dask` from `PATH`, because `dask ssh` passes its own `sys.executable` on to the workers it starts remotely — a `dask` picked up from `PATH` could therefore run the fit's workers under a different Python than the fit. A missing command is now **refused before anything is launched**, naming the subcommands the installation does offer and the package that provides them, instead of surfacing as `FileNotFoundError` inside "an unknown error … please report this bug". The check asks the same `dask_cli` entry point group dask's own CLI builds its command set from, so it cannot disagree with what `dask` will actually run — and the test that exercises it against the real environment is the one that would have gone red for this issue, where every mocked assertion stayed green against a command that no longer existed. `docs/cluster.rst` and `tests/full_tests/cluster_manual.sh` now give `dask scheduler` and `dask worker` for manual setups, for the same reason. - **The temporary directory dask leaves behind is now removed under the name dask actually creates (#620).** `_cleanup_dask_workspace` deleted `dask-worker-space` — the name dask used before it renamed the directory to `dask-scratch-space`. Every dask version the project supports (`>=2024.1.0`) creates the new name, so the cleanup matched nothing and the scratch directories accumulated in the working directory and in the home directory, one per run. The failure was invisible by construction: the cleanup ran, raised nothing, reported nothing, and looked exactly like a cleanup that had worked. On a cluster, where home is usually shared storage under a quota, that ends with a user unable to write at all — with no indication that a fit was the thing filling the quota. Both names are now removed, from both locations, so directories left by earlier runs are cleared alongside new ones. There is no public API that reports the name: `distributed` hardcodes the `dask-scratch-space` literal itself, so the cleanup carries both names deliberately rather than deriving one. - **A multiple-shooting segment that comes back short of its end knot is refused instead of being read at the wrong time (#584).** `job_type = ms` builds every continuity row from the **last row** of a segment's trajectory, which is the end knot only if the integration reached it. An integrator that stops early and returns a partial result rather than raising breaks that quietly: the trajectory is finite, in the right columns, and its final state belongs to an earlier instant, so the defect `Phi_j(z_j, θ) − z_{j+1}` becomes a difference of states taken at two different *times* — a different constraint, satisfied by a different trajectory. Nothing downstream could see it, because the symptom is a nonzero defect, which is what an honest stage shows after θ has moved. Measured on the offline fixture: a stage seeded to be feasible at iteration zero, whose defect norm is exactly `0`, reported `0.0399` and would have optimized against it. The segment seam now checks that a span reached the knot it was asked for, and treats one that did not the way it treats any point that did not integrate — the local model goes non-finite and the search backs off, rather than the run dying on a point-specific failure. - **A segment simulation that carries no forward sensitivities now names the segment seam (#584).** It stopped on the gradient assembly's own message — *enable the gradient path (apply_routing) on every scored model* — which tells a `job_type = ms` user to enable a path they did enable. A sensitivity request is applied **per action** (#475/#482), so a segment run under a suffix the request never reached comes back with a perfectly good trajectory and no tensor at all; that is an internal wiring error and every point in the fit hits it, so it is now refused where it happens, saying which segment of which experiment returned what. - **The shooting suite now tests behaviour at pathological parameter points, not only the arithmetic at well-behaved ones (#584).** The two defects above were found by these tests. The two before them (#578, #581) were not: both reached `main` through a suite of 51 passing offline tests and were caught by pointing `job_type = ms` at a real model. The offline fixtures are a closed-form flow and an exponential decay, deliberately well-behaved so that every derivative is checkable exactly — which is the right design for what they verify, and leaves nothing asking what the method does when a *point* misbehaves. The offline backend's failures are now switchable — a region that does not integrate, a region whose trajectory is finite while its forward sensitivities overflow (#581's exact shape), a span longer than the "model" can carry, a one-off refusal on the n-th call, a missing sensitivity tensor, and a span that stops short of its end knot — and the new tests assert the three properties ADR-0110 states as design and nothing checked: an unusable local model **backs the search off** rather than ending the run; a run that stops early still **reports what it has already earned**; and an unusable point **never becomes a reported fit**, including the corner the method's own advantage creates, where every segment integrates, the whole horizon does not, and the run therefore reports no fit rather than a segmented score no ordinary run could reproduce. Each runs in milliseconds with no simulator. - **A gradient fit says, at job start, when bngsim declined the analytic `∂f/∂θ` for one of its models and CVODES' difference quotient is carrying every sensitivity column instead (#606, ADR-0121).** `CVodeSensInit1` takes one sensitivity-RHS callback for every column, so a single rate law bngsim cannot differentiate — an `abs()`/`floor()`/`erf()` term, a comparison it cannot solve, or a derivation that ran out of its build-time budget — declines the analytic path for the **whole model**. The substitution is correct and costs an extra right-hand-side evaluation per column per step, so an N-parameter fit pays roughly N× the sensitivity cost. On a fit measured in hours that is not a slower answer but no answer: on `Smith_BMCSystBiol2013` all 25 columns fell back, every start timed out to `inf`, and thirteen hours produced nothing. PyBNF surfaced none of it — no console line, no refusal, nothing distinguishing a gradient fit on the analytic path from one on the fallback. The decline did reach `.log`, because bngsim's logger propagates to root, but as one line per model written mid-run from N worker processes into a shared, noisy file: discoverable by someone who already suspects the problem, which is the wrong order. PyBNF now checks **once per model at gradient-path setup**, on the head node, before the fit has evaluated anything, and names the model, its column count, the expected cost multiplier and bngsim's own reason. The check is not the log line: the verdict is read off the compiled codegen artifact — whether it exports the analytic sensitivity RHS symbol, the exact symbol bngsim's C++ resolves to choose between the two paths. That matters because bngsim reports a decline while *generating* codegen source, and since lanl/bngsim#174 a warm structural cache skips source generation entirely: measured on 0.14.0 and 0.13.0 alike, the same declining model reports its decline on the first construction and says nothing on the second, while running on the same fallback both times. Since the cache persists on disk, the run that hears nothing is typically the second run of a fit — the one made after the first came back empty. bngsim's reason is still captured and reported when it is heard, but only ever as prose; nothing keys off its absence. The new `sensitivity_fallback` key chooses what happens next: `warn` (the default — today's behaviour plus the sentence), `error` (refuse the fit, for a long unattended run), or `ignore` (skip the check, including the one simulator construction per model it costs). The policy keys off the verdict rather than the reason, so a fit refuses or does not refuse reproducibly whatever state the codegen cache is in. A model PyBNF cannot read an answer for — a `codegen=False` run has no artifact — reports **no opinion**: logged, never warned about, and never refused, because guessing is wrong in both directions. One decline is worse than slow: a model that also branches at a crossing whose time moves gets a difference quotient that integrates straight through it, so every column is wrong at and after the crossing. bngsim ≥ 0.14.0 refuses such a run rather than return a gradient it has flagged as wrong; on 0.13.0, which PyBNF's floor still admits, the same model only warns and returns it, and PyBNF now says so at verbosity 0 whenever bngsim's reason reaches it. - **A comment is never part of a filename, so a file list stops at the `#` (#599, ADR-0120).** The three shared file tokens in `parse.py` were unanchored and lazy. Both properties are deliberate and both are kept. Unanchored is what lets a path begin with anything a filesystem allows, and lazy is what makes a comma list stop at the first extension rather than swallowing the line. The defect was that nothing stopped either property from operating across a comment, so a stray trailing comma followed by a comment that happened to mention a file declared that comment as a second file: model: a.xml, # note about b.xml -> models = {'a.xml', '# note about b.xml'} model = a.xml : d.exp, # note about e.exp -> exp_data = {'d.exp', '# note about e.exp'} Measuring the reach found three more declarations with the identical defect that the issue had not named, `mutant`, the edition-2 `experiment:` record's `data:` field, and the `data =` key, which is the argument for fixing the shared token rather than each declaration. The bogus entry always died somewhere rather than being fitted, but which error the user got depended on which extension their comment happened to contain, and none of them pointed at the comma. A bogus `.bngl` reported a missing model file, a bogus `.xml` reported a parse error in SBML the user never wrote, a bogus `.ant` demanded an optional dependency they never asked for, and a bogus `.exp` sent them into their BNGL to add a `suffix=>` action for a measurement that was a sentence in English. - **A malformed `data` line now says what the key takes instead of denying that the key exists (#609).** The per-key format-hint chain had no branch for `data`, so any malformed `data = ...` line fell through to the generic fallback and was told `data is not a valid configuration key`, which is untrue. `data` is a real key: a comma list of `.exp` files bound to a bring-your-own callable objective (ADR-0050), valid alongside `objective = callable`. The old message sent the reader hunting for a typo in the key name rather than in the file list, which is the one place the error actually was. It now reports its real format, in the style of the neighbouring branches. Found while fixing #599. - **The discrete-event gradient gate reads a bngsim capability instead of a version floor, so a from-source build that merely *declares* a new enough version no longer passes it (#558, ADR-0119).** `BNGSIM_HAS_EVENT_SENS` gates forward sensitivities that survive a discrete event, and it does not gate a missing feature — it gates **silent wrongness**: a build below the line does not refuse an event it cannot differentiate, it returns a finite tensor with the event's contribution missing. It decided by comparing `bngsim.__version__` against exactly `0.12.2`. bngsim bumps its version at the *start* of a release cycle, so every from-source build made between that bump and the fixes that set the line (lanl/bngsim#144, #146) declares the same string as the release that carries them, clears the floor, and is reported as carrying them. The two failure directions are not symmetric — a false *absent* is a refusal and a metaheuristic fit; a false *present* is a gradient fit that runs to completion and reports a converged wrong number — and a version compare could only ever be wrong the second way. The neighbouring `BNGSIM_HAS_PER_SPECIES_ATOL` already carried the argument, for the same version string, in a comment on the wrong flag. The gate now resolves through published capabilities: `features['event_sensitivities']` if bngsim ever publishes a dedicated key (both directions, so the flag starts reading the real answer on the first build that grows one, with no PyBNF release), otherwise `features['effective_ic_sensitivity']` — a **witness**, usable because lanl/bngsim#155 added it three commits after #146 inside the same release window, so a build that publishes it necessarily carries the fixes. The version survives only as a veto: it can no longer report the capability present on its own, because the witness shipped *in* 0.12.2 and a build claiming 0.12.2-or-newer without it is exactly the pre-release build at issue. **No install that works today is refused** — every released bngsim at or above the floor publishes the witness — and a refusal now names the route that decided (`event_sens_probe()`) rather than telling a reader who already has 0.12.2 to upgrade to 0.12.2. - **A fit whose bngsim loaded a compiled core older than its own C++ says so at job start, not in import noise (#558, ADR-0119).** An editable bngsim serves live Python from the source tree while loading `_bngsim_core*.so` from a separately built artifact with auto-rebuild off, so the two halves drift — one install reporting `0.12.2` was found with a core binary three days older than the `.cpp` beside it. Every version, metadata and feature-key check passes there, because nothing in the Python layer moved. bngsim detects it by mtime and warns, but it warns at *import*, which for PyBNF is while the `pybnf` package loads: before `init_logging`, before the config is read, and before the user has committed to anything. PyBNF now repeats it at job start — the core's identity line (path, build commit, mtime) to the log unconditionally, the staleness report promoted to a console warning at verbosity 0 — where a reader can still stop a run that would otherwise spend hours producing statements about code that is no longer in the tree. `bngsim_build_id()` exposes the commit the core was built from, which is the only thing on hand that tells two installs declaring one version apart. Every read is guarded and memoized; an install that cannot answer reports no opinion rather than taking the fit down. - **A `parameter:` record is now held to the same declaration rules as the equivalent `*_var` line (#603, ADR-0118).** `_check_variable_keyword_combination` refuses an incoherent pairing of free-parameter declarations and `job_type` — an unbounded prior handed to a box-mode optimizer, `var`/`logvar` handed to a method that draws a population, a mix of point starts and boxes. It decided what it was looking at by pattern-matching **config key names** with `re.search('var$', k[0])`, which a `('parameter', )` key never matches, so the whole rule was silently bypassed by the edition-2 syntax: legacy normal_var = p1 0 1 -> refused: Box-mode optimizer requires a bounded prior record parameter: p1, prior: normal -> ACCEPTED That matters more than an ordinary validation gap, because the record syntax is the *only* one that can express `initial_value` — so the surface most likely to be used for careful, seeded work was the one with no coherence checking at all. The rule now keys on the **built `FreeParameter`** rather than on the key that declared it, and runs after the variables exist rather than before. Both syntaxes produce the same `FreeParameter`, so both now get the same answer. The obvious repair — re-deriving the keyword set from `{v.type}` — was tried and rejected: a *truncated* prior carries a real finite box while its family does not, so it would have falsely refused `prior: normal, ..., lower: X, upper: Y` on `job_type = gntr`, which is an entire benchmark corpus. Verified equivalent to the old rule for every untruncated declaration (family-level and parameter-level bounded-support agree across all 48 registered prior keywords), and run against **1049 real `.conf` files** with zero refusals. The dead branch for ADR-0015's third fit_type category is deleted: every registered refiner now also carries `start_from_box`, so the "point-only start optimizer" category is empty and its code was unreachable. Error messages now name the offending **parameters** rather than a keyword the record user never wrote, and point at `start_point` for the case they usually mean — a bounded box searched from a chosen point. - **`starting_params` is now a configuration error on a `job_type` that has never read it (#559, ADR-0117).** It has exactly one read site — the Bayesian sampler base — so on the other fourteen `job_type`s it was accepted, validated against nothing, and then discarded without a word: a `gntr` job seeded with it produced **bit-identical** output to the same job with the line deleted. The error names `start_point` as the replacement, which every `job_type` reads and which is matched by **name** rather than by position (`starting_params` is positional against declaration order, while every result file PyBNF writes is alphabetical, so round-tripping a result row back into it silently permutes the values). Unchanged for the six samplers that do read it, which is every shipped conf that sets it. - **A mixed bounded/unbounded parameter set no longer starts every parameter at the wrong place (#583, ADR-0117).** Start resolution was all-or-nothing, so one unbounded parameter sent *every* parameter down the point-start branch — where a bounded parameter's `p1` is its **lower bound**, read as if it were a sampling-space start value. A `loguniform_var` over `[1e-3, 1e3]` started at `10**1e-3 = 1.0023`, its lower corner, with nothing logged at any level. Resolution is now per parameter. - **CMA-ES no longer freezes a coordinate whose prior is truncated (#583, ADR-0117).** The per-coordinate box width was `p2 - p1`, which is the box only for a `uniform`/`loguniform` declaration; for a **truncated** prior those are the family's location and scale, so the width came out as the scale and, for the entirely ordinary `sd == mean`, as exactly `0.0`. CMA-ES squares these into its initial covariance diagonal, so that coordinate got a singular covariance and could never move for the whole run. Widths now come from the prior's own support — bit-identical for every `uniform`/`loguniform` box. - **`FreeParameter` no longer skips its bounds check for a value of exactly `0`, and no longer mutates the shared template on a fold (#583, ADR-0117).** The check was guarded by truthiness rather than `is not None`, so `initial_value: 0` — a legitimate value for a linear parameter — was stored unvalidated. Separately, folding an out-of-box value wrote `self.value = self.lower_bound` onto the template `FreeParameter` living in `Configuration.variables`, which every Algorithm aliases and which rides the algorithm's pickle, so the contamination survived a checkpoint and a `--resume`. Nothing read it. - **The concurrent multi-start scatter honors `initialization`, and says so when it cannot scatter (#583, ADR-0117).** It called the Latin-hypercube sampler unconditionally, so `initialization = rand` was a silent no-op for the whole gradient/CMA-ES multi-start family; and a `population_size > 1` on a fit with no box to scatter across was silently reduced to a single start — the same "accepted, does nothing, says nothing" shape as `starting_params`. - **`job_type = profile_likelihood` no longer reports that no start point was supplied when some were (#583, ADR-0117).** A partial specification cannot be θ\*, so it correctly falls through to the polish — but it said "No initial_value supplied" while doing so. It now names the parameters that were left undeclared. A complete specification keeps its established meaning there: those values are the optimum, and the polish is skipped. - **`sbml_rtol` is checked for finiteness, not just for sign, so `sbml_rtol = inf` no longer reaches CVODE (#586).** The conf grammar's number token also matches `inf` (ADR-0047's open truncation side), so `sbml_rtol = inf` parsed to a float that cleared the config check's bare `tol <= 0.` and was installed as the relative tolerance — which turns relative error control off rather than erroring. `sbml_atol` never had the hole, because `parse_atol_setting` has always demanded finiteness. Found while adding the per-model `rtol:` field of #586, which would otherwise have inherited the same check. - **A differential evolution fit whose objective is small no longer stops before its parameters have separated (#648, ADR-0127).** The convergence fix below replaced a dimensionless ratio with an absolute objective range, and let an unset `de_tolfun` keep using `stop_tolerance`'s number so existing configurations kept the threshold magnitude they had. ADR-0115 justified that by arguing an absolute range at the same magnitude is a stricter stop. It is stricter only above an objective of 1. Below that it is looser, and without limit. At the 2e-05 a well scaled sum of squares fit reaches, a range of 0.002 stops the run at a spread fifty thousand times wider than the ratio of 0.002 it replaced, so the population satisfies it almost as soon as it is scored. What that costs is a wrong answer that looks right. Tutorial lesson 25 fits three pharmacokinetic rates from one observed curve, and stopped early reporting `k_transit` as 11.18 against a true 12.76 and `k_abs` as 11.50 against a true 9.11, at an objective of 2.19e-05 and with the third rate correct to three digits. The two affected rates trade against each other, so the early stop fits the data well and nothing in the output looks wrong. Whether `de_tolfun` was set now decides what the threshold means. An explicit `de_tolfun` is an absolute range in the objective's own units and is honoured as written, at any sign and any scale, exactly as it was. An unset one falls back to `stop_tolerance`, which keeps the meaning that key has always had: a dimensionless ratio where the objective is positive, applied as `max - min <= tol * min` so an all-zero population never divides, and an absolute range only where the objective is not positive and a ratio would mean nothing. So a positive-objective fit converges as it did in 1.7.0, a likelihood fit keeps the fix below, and one value of `stop_tolerance` means the same thing on fits whose objectives are decades apart. Ignoring failed simulations, and the island guard, are unchanged. - **`job_type = de` and `job_type = ade` no longer stop after generation 0 on a negative objective (#561, ADR-0115).** The Differential Evolution family tested convergence with a *ratio* of objectives — `max(fit) / min(fit) < 1 + stop_tolerance` — which reads as convergence only on a positive objective bounded below by 0 (a χ², an SSE). On a likelihood objective (a negative log-likelihood, unbounded below) it fired at generation 0: an all-negative population lands the ratio in `(0, 1]`, and a single `inf`-scored failed simulation makes it `-inf`, below *every* threshold — so no value of `stop_tolerance` disabled it, and both members of the family were unrunnable on any estimated-σ likelihood fit (the whole Grein et al. 2026 benchmark subset-I corpus, 23/23 slugs). On `Borghans_BiophysChem1997` (`islands = 4`, `population_size = 400`, `max_iterations = 600`) the run terminated inside the first generation, spending a 240,000-evaluation budget on 0 generations of search — and `stop_tolerance = -1e9` still fired, because the failed-sim `-inf` is below that too. The convergence test is now an **absolute range in objective units**, `max - min <= de_tolfun`, assessed over the **finite** fitnesses only — sign-agnostic, and with failed simulations (`inf`) ignored so one dead candidate can neither trigger nor block the stop. It gets its own key, `de_tolfun` (a range in objective units, where `stop_tolerance` was a dimensionless ratio), which falls back to `stop_tolerance` when unset so an existing config keeps its threshold magnitude — mirroring how `cmaes_tolfun` splits off `cmaes_stop_tol` (ADR-0106, the CMA-ES sibling of exactly this defect). The shared check lives in one `DifferentialEvolutionBase` helper, so `de` and `ade` cannot drift apart again (the missing `!= 0` guard that let `ade` divide `0/0` on an all-zero population was such a drift; the range form has no division). In an island run, convergence is assessed only once every island has completed an iteration, so ignoring `inf` cannot let one finished island stop the whole search before the others have run. Regression tests pin each failure mode at its decision point (an all-negative spread and an `inf`-defeated threshold are *not* converged; a collapsed finite population *is*; `de_tolfun` is its own knob; the all-zero `ade` population no longer divides) plus an end-to-end guard that both optimizers advance past generation 0 on an objective that is negative everywhere the population lands. - **`job_type = ms` no longer dies on a fit with a measurement-model formula observable (#578).** Multiple shooting was unusable on essentially the whole PEtab-imported corpus — including its own motivating problem, `Borghans_BiophysChem1997` — failing on the first outer iteration with `Measurement model 'Ca' would shadow an existing simulation-output column`. The setup was all correct (noise profiling, the `4-2-1` ladder, knot placement); it died the moment the loop asked for a second evaluation. The measurement layer materializes each `observable: , formula: ...` column *into the trajectory in place* (ADR-0036) and deliberately refuses a column that already exists. Every ordinary fit satisfies that for free, because the propose/score loop scores a freshly simulated `Data` every time. Multiple shooting caches its segment trajectories per point — so that one augmented-model evaluation costs one pass of segment simulations rather than two — and the outer loop then re-evaluates at the point the inner solver finished at, which is a cache hit on those very objects. Fixed at the cause: the assembled objective is now memoized on the point, so each simulated trajectory is scored exactly once. That also removes a redundant gradient/Fisher assembly per outer iteration, the larger of the two costs. Sound because the objective at a point does not depend on the multipliers — only the augmented model combines them. The shooting suite structurally could not see this: its fixtures score native columns (a species, an observable), so the measurement layer never ran. The regression tests use a formula-observable fixture, and both fail with the exact production error when the memo is removed. - **`job_type = check` runs again (#569).** #564's method-chain record was built from `alg.res_dir` on a line that ran before the `job_type != 'check'` branch nine lines below it, so every check run died in setup with `AttributeError: 'ModelCheck' object has no attribute 'res_dir'` — no objective value at all, not merely a noisy tail. `ModelCheck` deliberately does not subclass `Algorithm`, so it has no `res_dir`, and none of `stop_reason`, `completed_simulations` or `trajectory` either; the fit-phase recording that reads all three was above the branch too. Both now sit inside it, with the boundary stated in a comment so the next addition to `main()` lands on the correct side of it. A check run still writes no `method_chain.json`: it is one evaluation of the parameters as given, not a chain of search phases. The regression test that was missing — a `job_type = check` job driven end to end through `main()` — now guards the path; every existing check test called `run_check()` directly, which is why two commits landed on top of the break. - **`wall_time_fit` no longer silently downgrades `refine = 1` to no refine at all (#564, ADR-0107).** `refine = 1` requests a *method* — search globally, then polish the result with a local optimizer — but a wall-clock-budgeted search runs until the clock stops, so it has no reason to leave anything behind, and the polish (new work, forbidden once the budget is spent, ADR-0093) never started. Not occasionally: **15 of 15 runs** in a benchmark campaign (`Borghans_BiophysChem1997`) configured as `cmaes` + `refine = 1, refine_method = gntr` actually ran plain `cmaes`. And the downgrade was invisible — ADR-0093's whole promise is that a budgeted run "writes exactly what a converged one writes", so `sorted_params_final.txt` and `information_criteria.txt` looked identical either way and the only trace was one line on stdout. A harness that scores a directory could not tell which method it had measured. A new global key **`wall_time_refine_frac`** (default `0.1`) holds that share of `wall_time_fit` back from the search, so the refine runs on a slice the search was never allowed to spend. The run's total is unchanged — one deadline still bounds the whole run, it is just partitioned rather than first-come-first-served — and the split is stated on the console before the search starts. The reserve is a floor, not a cap: a search that converges early hands everything it did not spend to the polish. No reserve is taken when there is no refine to protect (no `refine`, no budget, or a `refine_method` naming the algorithm the fit itself ran), so a run that asks for no polish is byte-identical to before. `wall_time_refine_frac = 0` restores the old split, and the resulting skip is now a `print0` warning that names the method that did *not* run, the method that ran alone, and the key that would have made room. - **A refined run's `sorted_params_final.txt` describes the refined point (#564).** The refiner wrote its result only to `sorted_params_refine_final.txt` while rewriting `information_criteria.txt` from the same end-of-run tail, so two files in one `Results/` disagreed about which parameter set they described — and the *conventional* name carried the point the requested method chain did not end on. A refine's end-of-run output is the run's end-of-run output, and is now written under both names; the `refine_`-prefixed file is unchanged. - **A bootstrap replicate's refine no longer writes into the main run's `Results/` (#564).** `_refine_best_fit` redirected a replicate's `sim_dir` and `failed_logs_dir` to the `Results-boot{N}` peers but not its `res_dir`, so every replicate's polish overwrote the *main* fit's `sorted_params_refine_final.txt` and `stop_reason.txt`. The refiner now writes where the fit it is polishing wrote. - **A refine's wall-time stop reason is appended to `Results/stop_reason.txt`, not written over the fit's (#564).** Both phases share one Results directory; a run where the search hit the deadline *and* the polish did has two facts to report, not one that replaces the other. ## [v1.7.0] - 2026-08-12 ### Changed - **An SBML model on `sbml_backend = bngsim` no longer charges every species for the tolerance its smallest one needs (#549, ADR-0105, supersedes ADR-0103).** ADR-0103 derived a single `atol` from the median species value because bngsim's `Simulator.run` took only a scalar, and it wrote down what that cost: `Brannmark_JBC2010` reads `3.3e-10`, which holds its `IR`/`IRS`/`X` species at ~10 to `3.3e-11` *relative* — three decades tighter than the `rtol` that governs them — to buy a resolution for a `1.76e-9` transient that the same ADR had already decided not to chase. Now that lanl/bngsim#196 routes a vector to `CVodeSVtolerances`, that over-tightening is given back per species: `atol_i = clamp(sbml_rtol * y_i, the model's scalar, 1e-8)`. Measured on 100 points sampled from Brannmark's own fit box with the fit's sensitivity request applied, 39 dead simulations become **33**, in 428 s rather than 576 s. **The lower clamp is the change, and it is there because the obvious rule loses.** "Resolve each species to `rtol` of its own magnitude", full stop — which is what #549 proposes — puts that transient at `1.76e-17`, and on the same 100 points killed **91 of 100**: ADR-0103's withdrawn *minimum* rule reappearing one species at a time. The `1e-16` floor #549 asks about rescued nothing (91 either way), because the damage is done well above it. A tolerance below `rtol*|y|` is inert until a species has decayed far below its nominal value, and what it then demands is that a species which has decayed into nothing be resolved as if it had not; telling that apart from a genuinely tiny species needs the trajectory, not the initial values. So every entry now lies in `[the model's scalar, 1e-8]`: no species is integrated more tightly than PyBNF integrates it today, and no model that runs today can start failing. A control confirms the mechanism is the values and not the plumbing — the same scalar sent as a uniform vector reproduces the baseline exactly, 39 and 39. Over the 23-slug subset-I corpus 19 models take the same scalar call as today and the 4 that #546 tightened take vectors, which is the shape of a refund: only a model that was charged can receive one. A species declared at zero has no magnitude of its own and falls out of the same expression at the model's scalar, leaving it where ADR-0103 put it. The derivation stays a property of the **model file** — read off the SBML document at load, held for the whole fit, never re-derived from the fit point, which is why bngsim's state-reading `AUTO` token is not used: a tolerance that moved with a fitted initial condition would put a step in the objective wherever the derivation crossed a rounding boundary, and it would be invisible, since the objective still looks correct and only the search behaves oddly. ADR-0103's median-derived scalar does not retire either — it becomes the steady-state convergence cutoff, passed explicitly whenever the vector is in force, because bngsim's own fallback for that cutoff is the Simulator's `1e-8` rather than anything derived from the vector, and taking it would silently return every `time = inf` measurement and every pre-equilibration phase to "equilibrium at t = 0" on a small-scale model. `sbml_atol` remains a single number and remains the off-switch: stating it integrates every species at that value and pins the pre-#196 code path bit-for-bit. A bngsim without the capability keeps ADR-0103's scalar unchanged, detected by name (`bngsim.AUTO`) rather than by version, because the build that first carried #196 declares the same version string as the wheel that predates it. ### Fixed - **The CMA-ES restart battery's TolFun trigger no longer gets more eager as the fit gets better (#550, ADR-0106, amends ADR-0082).** ADR-0082 made TolFun's stagnation threshold *relative* to the current objective, `frange <= cmaes_stop_tol * max(1, |f|)`. PyBNF minimizes a negative log-likelihood, which is unbounded below, so `|f|` **grows** as the fit improves and that threshold rises as CMA-ES approaches the optimum — while Hansen's window `10 + ceil(30N/lambda)` shrinks as IPOP grows the population (30 generations at `lambda = 32`, 11 at `lambda ≈ 1900`). The two move in opposite directions and compound, so a late restart must improve by *more* within *fewer* generations than an early one, and IPOP's large-population restarts — the ones grown to do the heavy lifting — were the ones cut off mid-descent. On `Elowitz_Nature2000` (Grein subset-I, k=21) restart 3 was descending `OG` 53.0 → 26.4 → 5.105 and was killed at the bottom of that descent because `0.001/generation × 11 generations = 0.0105` fell just under `1e-4 × 121.06 = 0.0121`; across two fits **all 14 restarts** fired on TolFun and not one run ever converged. TolFun now compares an **absolute** range in objective units, as Hansen's `tolfun` does (pycma's relative variant `tolfunrel` normalizes by the run's *initial* median, a scale that does not drift with fit quality either; neither reference form uses the current `|f|`). On that fit the threshold becomes the configured `1e-4` and the descent clears it by a factor of 100. **TolFun also gains its own key, `cmaes_tolfun`.** `cmaes_stop_tol` is a step length in the parameter sampling space and TolFun is a range in objective units; no single value is right for both, and reaching a TolFun that fired at all on the fit above meant declaring the search distribution converged at `1e-4` in `u`, seven orders looser than the default. Unset, `cmaes_tolfun` follows `cmaes_stop_tol`, so an existing config keeps the threshold magnitude it had — the only change is dropping the `|f|` factor, which can only make TolFun fire less, and a fit whose objective satisfies `|f| <= 1` is unchanged outright. The restart reason now also reports the tolerance it used (`range 0.0105 over the last 11 generations, tolerance 0.0001`), so a restart's arithmetic is checkable from the log. `cmaes_restarts = 0` (the default) is untouched: the battery is still restart-gated (ADR-0070). - **A PEtab v1 problem whose parameter table merely *has* a prior column no longer loses its log estimation scale (#548).** `petab1to2_preserve_scale` re-injects the `parameterScale` that `petab.v2.petab1to2` drops, skipping any row that already carries a prior so a scale petab1to2 already folded into one is not clobbered. That guard is right in intent and impossible to implement in v2 alone: petab1to2 **materializes** v2's implicit default — `priorDistribution = uniform` over the bounds — into the converted table whenever the v1 table has a prior column at all, even an entirely empty one, and after conversion a materialized default and a declared `uniform` are the same cell. So the decider was whether the upstream TSV happened to carry a prior column, a cosmetic property of the file: `Zhao_QuantBiol2020` (four prior columns, 100% empty) lost **all 28** of its log10 parameters and `Schwen_PONE2014` (six real `parameterScaleNormal` priors, the rest blank) lost **24 of 25**, while `Giordano_Nature2020`, whose v1 table has no prior column, converted correctly. The conversion now reads which rows carried a prior from **v1**, where a blank is still a blank, and skips only those; `inject_log_uniform_priors` gains an optional `declared_prior_ids` and keeps the conservative v2-only reading when it is omitted. This was silent by construction: the re-injected prior sets only the search scale and initial sampling, and PyBNF's optimiser objective excludes the prior, so the objective, the `simulatedData` oracle check and the finite-difference gradient check all still passed — `Zhao`'s nominal `J_paper` is unchanged to 13 significant digits. Only the search was wrong, and a multi-decade parameter sampled linear-uniform presents as a fit needing more starts: `Zhao`'s `gamma_*` sit on `[1e-08, 1]` with an optimum near 0.05–0.39 and its `sd_*` on `[0.001, 1e5]` with MLEs of 186–5013, so across 28 parameters effectively no box-sampled start lands near the basin. A 100 × 1000 multi-start stalled at ~718 and was decelerating; on the corrected scale it beat that in under 90 seconds. `Schwen` is the discriminating case — its six declared priors survive as `log-normal`, its five genuinely `lin` parameters stay `uniform`. - **An SBML model whose species are far below 1 no longer integrates — and differentiates — at a tolerance larger than its own state (#546, ADR-0103).** `Giordano_Nature2020`'s assembled gradient disagreed with central differences on 41 of its 50 fitted parameters, by up to 26%, identically at every finite-difference step size, with no refusal and no warning. The model is piecewise-in-time — 110 `piecewise` expressions across 14 assignment rules, all gated on the COVID NPI stage boundaries — and the error partitioned along whether a parameter sat behind a time gate, so it read as unhandled switching. It is not: bngsim's SBML loader already registers every `time` inequality as a CVODE root (13 for this model), and the crossings are landed on exactly. The defect is the **absolute** tolerance. CVODE weights each state by `rtol*|y| + atol`, so a constant `atol` declares values beneath it to be noise — a statement about the model's units, and bngsim's `1e-8` is BNG2.pl's, right for a model in molecule counts. Giordano is a population-*fraction* model whose species sit at `1.7e-8..1`, median `3.7e-7`: its early trajectory carries no significant digits, and the forward-sensitivity solve carries fewer still, since CVODES scales the state tolerances by the parameter magnitude for the sensitivity vectors. The gate correlation is real but incidental — a gated parameter acts only inside its own stage window, and the earliest windows are where the states are smallest. Tightening `rtol` by four decades changes nothing; tightening `atol` fixes it. The bngsim SBML/Antimony path now derives `atol` as `rtol` times the model's median strictly-positive species initial, clamped to at most the backend default and at least `1e-16`, so it can only ever tighten: across the 23-model subset-I corpus 19 are untouched and 4 tighten. Giordano's worst column goes **7.7e-02 → 4.5e-04** for ~14% more wall clock; Brannmark 5.0e-05 → 3.6e-05 and Bertozzi 2.7e-05 → 2.4e-05 at no measurable cost. The median rather than the minimum, because `Brannmark_JBC2010` seeds one transient intermediate at `1.8e-9` against principal species at `0.1..10`, and resolving *that* asks for `1e-17`, which makes the model fail outright on `mxstep` at interior fit points. Unchanged: every BNGL/net model (its tolerances come from the actions block, and BNG2.pl parity is what that backend is measured against), every stochastic run, the RoadRunner backend, and every SBML model of order-one scale. - **A scalar fit no longer re-derives the analytical Jacobian on every action (#544).** #543 warmed the cached engine template so clones inherit the compiled sensitivity RHS, but it warmed only when a sensitivity request was active — and the same never-warmed-parent shape costs the **scalar** path too, one artifact over. bngsim's `Simulator.__init__` calls `model.prepare_analytical_jacobian()`; PyBNF builds that `Simulator` on the per-action *clone*, and the clone is discarded. `clone()` carries the `_jac_attempted` sentinel parent → child precisely so a derived parent yields cheap clones, but nothing ever derived it on the parent, so every scalar action re-ran the SymPy derivation from scratch. PyBNF now warms unconditionally and lets the *shape* of the warm depend on the request rather than gating the warm itself on there being one. Measured through PyBNF's own action path on the 44-species `yeast_cell_cycle` model: **0.1542 s → 0.0057 s per action** (27x), derivations 10 of 10 → 0 of 10. As reported on #544, `Smith_BMCSystBiol2013` goes 0.0401 s → 0.0224 s, 20 of 20 → 0 of 20, taking that job's shipped 64,000-evaluation `cmaes` budget from ~36 core-hours to ~14. Trajectories are bit-for-bit unchanged. Two guards this needed: a scalar warm must **not** satisfy a later gradient warm (bngsim clears a plain-RHS artifact and regenerates it at the first sensitivity request, so a scalar-warmed template would be correct and save the gradient path nothing) — the warm-state predicate and the per-shape attempt memo both answer "not yet" for the sensitivity shape; and a model with **no ODE action** warms nothing, since bngsim derives the Jacobian under ODE dispatch and nowhere else, having deliberately moved it off the load path so a stochastic run never pays SymPy. Applies to every SBML/Antimony fit through `sbml_backend = bngsim`, on every `job_type`; the larger the model the bigger the effect, since the derivation scales with the network and the solve does not. The `.net` backend clones from a held `_engine_model` in the same never-warmed shape and does re-attempt the derivation per evaluation (measured 4 of 4), but it costs it essentially nothing: a BNGL network is all-Elementary, so bngsim takes its closed-form C++ Jacobian rather than SymPy — 0.1511 s → 0.1472 s per evaluation on `egfr_ground.net` (356 species), within noise. Left alone rather than warmed on a measurement that does not justify it. - **A gradient fit no longer rebuilds the analytical sensitivity RHS on every action (#543).** `_get_engine_template` caches one loaded bngsim model per SBML text per worker process and #415 clones it per action, precisely so the parse and the derived Jacobian are paid once. bngsim's `clone()` cooperates, carrying the compiled sensitivity artifact parent → child — but the `Simulator` was built on the *clone*, so the clone was what discovered and recorded that artifact, and the clone is discarded at the end of the action. Discovery flowed child-ward only: the template's `_codegen_so_path` stayed empty forever, and every action regenerated the C source, and every symbolic derivative behind it, because bngsim keys its compiled `.so` on a hash of that source. PyBNF now warms the template once per process, with the sensitivity request the actions will use (a scalar-shaped warm would be correct and save nothing — bngsim regenerates it at the first sensitivity request), so `clone()` propagates it from then on. Measured on `Smith_BMCSystBiol2013` (133 species, 16 sensitivity columns) through PyBNF's own action path: **2.015 s → 0.537 s per action**, source generated 4 of 4 times before and 0 of 4 after; the tensor is bit-for-bit unchanged. Inside a real `gntr` run of that job — dask worker, condition applied, the experiment's own sample times — its first action goes from 3.731 s to 0.327 s. Invisible in a profile that looks at simulation — the integration is untouched and all of the difference is `Simulator(...)` construction. A **scalar** (metaheuristic) fit never generates the source at all and was left untouched here (0.194 s per action either way), which bounded who *this* entry helps; #544 above warms it for a different artifact. The `.net` backend clones from a held `_engine_model` in the same never-warmed shape but is **not** affected here: its `.net` codegen memo is keyed on the file path rather than on generated source, and it regenerates nothing (measured 0 of 6). - **A pre-equilibration condition that doses a species from a fitted parameter no longer reports a zero gradient column for it (#538, ADR-0101).** `preequilibrate:` applies its condition inline, so a species target becomes a `setConcentration` written *before* the first phase — and `Model.set_concentration` reads an assigned amount as a literal initial condition (`∂x_k(0)/∂θ = 0`), retiring whatever seeding the species' `.net` expression carried (lanl/bngsim#113). ADR-0098 supplies that row for a write between two phases; with nothing pending it had nothing to rebuild and left the write to the backend, so an amount like `"A()" = 2*k_deg` contributed **exactly zero** to `k_deg`'s derivative. Nothing failed and no refusal fired: the fit simply walked a wrong steepest direction to a plausible answer. PyBNF now *declares* the assignment's own `∂x_k(0)/∂θ` (`Model.declare_ic_sensitivity`, the API bngsim documents for a hand-assigned θ-dependent initial condition), so bngsim's own seeding starts from it — narrowly, only when a fitted parameter reaches the amount, so every protocol whose gradient was already right reaches the backend through the same calls as before. An `addConcentration` re-declares the row its constant shift left alone. Visible only with a **fixed-duration** equilibration (`equil_t_end:`, what the preincubate → wash → dose-scan protocols use); a steady-state equilibration relaxes the dose away, so the derivative is genuinely zero there. Also new: an intervention amount that reads a fitted parameter **no** requested forward-sensitivity column carries is now refused by name, on both this path and the mid-protocol one, rather than silently contributing a zero row. - **A pre-equilibrated dose-response experiment can now be fit by a gradient method (#532, ADR-0098).** The preincubate → wash → dose-scan protocol (`preequilibrate:` + `condition:` + `type: parameter_scan`) refused every scored gradient evaluation, and the refusal landed at *scoring* — so a `trf` fit of Erickson 2019's `igf1r` job "finished" with `inf` at all ten starts and said only `Unknown error during job bestfit_infocrit`. Two things were wrong. The guard itself was **stale**: bngsim 0.12.0 (lanl/bngsim#81, #111) carries the state each dose restores *together with* its `dx/dθ`, which is exactly the capability the guard said did not exist; it is now a capability gate (`bngsim >= 0.12.0`, the new `pyproject` floor), and each dose's tensor stacks down the dose axis like any other scan's. Underneath it, the protocol's **wash** was silently discarding the equilibration's derivative — `Model.set_concentration` drops the pending `dx/dθ` rather than guess an externally supplied amount's, so *no* pre-equilibrated experiment with a species intervention could be gradient-fit, a measured time course failing outright with `carry_sensitivities=True, but no matching forward-sensitivity seed from a prior phase is available`. PyBNF now supplies the row it knows: the intervention's own `∂x_k(0)/∂θ` — `0` for a literal amount, the exact derivative for one written over model parameters (differentiated through the `.net`'s derived ids), the carried row for an `addConcentration` — with the rest of the matrix preserved, and an honest refusal naming the assignment when it lies outside the arithmetic grammar. A `resetConcentrations()` that follows a `saveConcentrations()` is likewise recognised as returning to a *carried* state, which is what made a model's **second** pre-equilibration experiment refuse. All seven `igf1r` rate constants now agree with central differences to ≤ 2.3e-04 on all three experiments, and the fit reaches a finite objective. - **A refusal raised while simulating now stops the fit and states its reason, instead of returning `inf` at every start (#532).** `Job.run_simulation` swallowed a user-targeted `PybnfError` into its generic "unknown error" arm, so a property of the *setup* — a model construct this `job_type` cannot handle, a missing backend capability — was reported once per evaluation as a failed simulation and the run continued to a meaningless finish. The documented fail-fast policy (re-raise; it would fail every job) existed one layer up and was never reached. Scoring failures are unchanged: a per-point objective failure still penalizes that point (#388). - **A gradient start that reaches a point it cannot differentiate no longer takes the whole multi-start fit down with it (#528, ADR-0092).** A stiff parameter point can score finitely while its forward sensitivities diverge, leaving a finite objective with a non-finite gradient. That model went straight into the trust-region factorization, where LAPACK refuses it (`Sorry, an unknown error occurred: numpy.linalg.LinAlgError: SVD did not converge`) and the exception unwound out of the run loop — 19 healthy starts of a 20-start `gntr` fit discarded because of one, which inverts the reason multi-start exists. (`lbfgs` aborted the same fit by a different route: its NaN direction proposed a NaN point, rejected as `OutOfBoundsException: Free parameter k cannot be assigned the value nan`.) An unusable local model is now treated exactly as a failed simulation already was: **mid-search the trial is rejected** — the trust region shrinks, or the line search backtracks — and that start carries on from its current iterate; **at the start point that one start stops**, saying which model was unusable (`the Fisher model (gradient + EFIM Hessian) at the start point is not finite (the point scored, but its derivatives did not)`), while every other start keeps running and the global best is taken across the survivors. The two LAPACK calls in the step math (`svd`, `eigh`) are wrapped as well, so a factorization that fails on finite-but-pathological input routes the same way. `profile_likelihood`, which drives the same runners, was a third casualty of the same missing guard: a slice whose derivatives diverge now ends that one direction at a wall it names, with the un-optimizable grid point contributing no profile value — rather than entering an un-minimized upper bound as if it were the profile, which would inflate that point's Δχ² and could close the confidence interval too narrowly. A fit whose models are all finite is unchanged. - **A negative count is no longer scored as a perfect fit, nor counted in `n` (#523, ADR-0090).** The `neg_bin` family has a negative observation contribute nothing to the objective — right for the fit, since a negative count has no negative-binomial probability, and real surveillance data contains them (a downward revision of a cumulative total makes a negative daily increment). But a negative-binomial PMF is self-normalizing, so that zero cost became `log p = 0` — probability **one** — in the pointwise log-density, a better per-point density than the family assigns any real count, including one the model predicts exactly. Those points were also counted as scored points, entering `n` for AIC/BIC and the LOO/WAIC observation axis. An observation outside its noise family's **observation domain** is now excluded exactly as a NaN observation already was: off the observation axis, out of `n`, and reported once per observable with a count (`excluded 4 measurement(s) of 'cases' in ...: this observable's noise model scores only a non-negative count`). Scoring data containing negative counts now matches scoring the same data with those rows deleted. The **cost** path is deliberately unchanged — such a point still contributes nothing to the objective and to the gradient — and every family whose support is the whole real line (`normal`, `lognormal`, `lnnormal`, `laplace`, `student_t`) is byte-identical. - **Steady-state (`time = inf`) measurements now load and fit (#521, ADR-0086).** A PEtab problem measured only at equilibrium imported fine but crashed at configuration load (`OverflowError: cannot convert float infinity to integer`): the experiment was materialized as an ordinary time course, which derives its step count from a (here infinite) endpoint. PyBNF's only steady-state route was the dose-response `parameter_scan` of ADR-0046, which needs a swept axis a plain equilibrium observation does not have. An `.exp` whose `time` column is all `inf` is now recognized as a **steady-state experiment**: it emits `simulate({...,steady_state=>1,n_steps=>1})` — the relaxation-with-early-stop primitive pre-equilibration already used — and the objective scores the datum against the run's final (equilibrium) row. `t_end:` bounds the relaxation (default `1e6`) instead of timing a readout, and `type: steady_state` may state explicitly what the data implies. Supported on BNGL (BNG2.pl/bngsim), bngsim SBML/Antimony, and RoadRunner (which uses its own steady-state solver, falling back to the bounded integration); forward sensitivities are carried at the equilibrium, so `trf`/`lbfgs`/`gntr` fit these problems. NFsim (`method: nf`) has no steady-state solve and is refused, as is an experiment mixing `inf` with finite times. This unblocks `Blasi_CellSystems2016`, the last unimported subset-I problem of the Grein et al. 2026 benchmark collection. - **PEtab conditions measured only at `t = 0` now load and evaluate as initial-state observations (#510).** A data-derived ``TimeCourse`` previously required at least one positive output time, so one legitimate initial-state condition rejected the entire imported problem (including every ordinary time course); this blocked ``Schwen_PONE2014``. SBML/RoadRunner and SBML/bngsim now return the initialized model as a one-row ``t = 0`` trajectory without invoking an integrator. The bngsim gradient path also supplies the initial-condition identity derivative and differentiates parameter-driven SBML ``initialAssignment`` expressions, so ``trf`` / ``lbfgs`` / ``gntr`` retain correct forward sensitivities for an initial-only experiment. - **PEtab natural-log Gaussian observables now import exactly (#509, ADR-0084).** PEtab v1 ``observableTransformation = log`` and v2 ``noiseDistribution = log-normal`` previously reached PyBNF's internal ``Gaussian(LN)`` kernel but could not be serialized into the generated ``.conf``; ``import_job`` raised ``NotImplementedError``. The new explicit ``lnnormal`` noise family is ``Gaussian(additive_on=LN, location=MEDIAN)`` and is kept distinct from PyBNF's existing ``lognormal`` (log10) family. Imports now preserve the natural-log residual and sigma units, and normalized pointwise log-likelihoods use the natural-log Jacobian ``-log(y)`` (so ``information_criteria.txt`` is on the correct absolute scale). This unblocks ``Blasi_CellSystems2016`` and ``Laske_PLOSComputBiol2019``. The exact reverse mapping also exports ``lnnormal`` as PEtab v2 ``log-normal``; log-scale Laplace remains unsupported by the native configuration surface. - **PEtab import now preserves replicate-specific `observableParameters` / `noiseParameters` bindings (#508, ADR-0083).** The per-measurement sidecar was keyed only by column, time, and placeholder, so repeated PEtab cells from different replicates collided and the last replicate's token silently replaced the others. This blocked `Fiedler_BMCSystBiol2016` by orphaning the first gel's scale parameters and could silently fit other problems with the wrong per-row scaling/noise binding. Replicate-aware sidecars now add a 1-based `replicate` column, using the same row-dealing partition that creates each `_repN.exp`; configuration loading and PEtab re-export select tokens by `(replicate, time)`. Legacy four-column sidecars remain valid and retain their shared-across-replicates meaning. - **PEtab import: `observableParameters`/`noiseParameters` placeholders with a fixed noise parameter, multiple noise tokens, or an affine/prediction-scaling noiseFormula now import (#495, ADR-0075).** Three related gaps in the placeholder→parameter mapping left benchmark-collection problems unimportable. (a) `Oliveira_NatCommun2021` — a `noiseParameters` id (`sd_cumulative_*`) that is **fixed** (`estimate=0`) was emitted as a `fit` free sigma the `.conf` never declared, so the job failed to load; it now inlines as a constant sigma. (b) `Fiedler_BMCSystBiol2016` — a **multi-token** `noiseParameters` cell (`s_gel;sigma`, a `noiseParameter1 * noiseParameter2` formula) was never split (only `observableParameters` was), so the whole cell was mis-read as one id; it now splits and binds each `noiseParameter${n}` per data point when row-varying. (c) `Raia_CancerResearch2011` — an affine `noiseParameter1 + noiseParameter2 * (species…)` produced a `noise_model` line that referenced the simulated trajectory a `formula` sigma cannot see, failing to parse; it now imports as the new `prediction_formula` source (see Added). All three land on crafted simulator-free fixtures scored against a hand-derived NLL. - **PEtab SBML import: an `observableFormula` referencing an SBML `assignmentRule` variable now imports instead of being rejected as "not a model entity" (#493).** An `assignmentRule`-defined variable (a ``/`` with `constant="false"` whose value is set by an ``) is a *derived* model output — the backend recomputes it at every step — so it is exactly the kind of quantity an observable is built from (the SBML analogue of a BNGL global function, which PyBNF already accepts). `import_job`'s measurement-model importer recognized only species / parameters / observables / functions, so six PEtab benchmark-collection problems (`Giordano_Nature2020`, `Laske_PLOSComputBiol2019`, `Rahman_MBS2016`, `SalazarCavazos_MBoC2020`, `Smith_BMCSystBiol2013`, `Zhao_QuantBiol2020`) could not be imported at all — a bare-name formula raised *"has a bare observableFormula … which is not a model entity"* and an expression formula raised *"references … which is not a known model entity"*. The importer now **inlines** any referenced assignment-rule variable down to the species/parameters its rule is defined over (recursively), exactly the resolution the config-load measurement layer already performs (#465, ADR-0036); the model file is carried byte-verbatim and a formula naming no rule variable is returned unchanged (the bare-name common case stays dependency-free). - **PEtab import: `observableTransformation = log10` is no longer dropped in the v1→v2 conversion (#499).** A PEtab **v1** observable with `observableTransformation = log10` (Perelson_Science1996, Borghans_BiophysChem1997, Elowitz_Nature2000, and other multi-decade-signal benchmark problems) imported as a **linear** `gaussian` noise model, so the fit optimized the *wrong* objective — a linear residual with no change-of-variables Jacobian instead of the `log10` residual the problem (and the paper) specify. The `log10` transformation was silently dropped by `petab.v2.petab1to2` (PEtab v2 removed the `observableTransformation` column and has **no** `log10` `noiseDistribution`; it downgrades `log10-normal` to a blank distribution), and `import_job` read only `noiseDistribution`, so the observable resolved to `Gaussian(LINEAR)`. Directly parallel to the `parameterScale` drop `petab1to2_preserve_scale` (#491) already fixes: - **`pybnf.petab.petab1to2_preserve_scale` now also re-injects `observableTransformation`** as a preserved extra column on the converted v2 observables table (v2 lint-clean; other tools ignore it). Since v2 has no faithful `log10` `noiseDistribution`, this extra column is the only channel for a `log10` residual — the observable twin of the `log-uniform` parameter-scale re-injection. - **The importer selects the noise family's *additive scale* from `observableTransformation`, not just the family from `noiseDistribution`.** `log10 + normal` → the native `lognormal` family (`Gaussian(LOG10, MEDIAN)`, which already carries the correct log10-space residual **and** the `Σ log(y·ln10)` Jacobian), emitted as `objective = lognormal` (or a `noise_model = lognormal, …` line); `log + normal` maps to the natural-log `Gaussian(LN)` (v2's `log-normal`). The `pybnf.petab.observables` adapter reads it the same way (`log10` → LOG10, `log` → LN, `lin` → unchanged), with a guard against a transformation that contradicts a log `noiseDistribution`. Natural-log Gaussian now routes to ``lnnormal`` (#509, ADR-0084); log ``laplace`` families still raise ``NotImplementedError`` — the boundary stays in code, never a silent mis-recovery. A linear problem is byte-for-byte unchanged. ### Fixed - **The router now reads what the backend seeded instead of inferring it, so an IC-seeding parameter is routed exactly once (#537, ADR-0100).** bngsim 0.12.2 (lanl/bngsim#157, answering our lanl/bngsim#155) exposes `Model.effective_ic_sensitivity()`, the `{species: {param: ∂x(0)/∂θ}}` the solver will actually be seeded with, from model structure alone and with a present-but-zero entry distinguished from an absent one. The rule is now stated rather than guarded around: **route a bound id's own parameter axis, and add an initial-condition term only for a `(species, param)` pair the backend reports absent.** That subsumes both of the defects either side of it — #535 was an axis dropped that was needed, #537 an axis kept that duplicated, and both came from inferring what the parameter axis contained. `ode_rhs_symbols` is demoted to an optimization (dropping a provably identically-zero axis), so a model that cannot answer keeps its axis instead of facing two silent wrongs; the refusal added for that case is gone, as is the `lowered_ic_species` build discriminator that existed only to survive the interval. `Fiedler_BMCSystBiol2016` returns to 3.68e-06 on 0.12.2 — its value before lanl/bngsim#147 widened the seeded class — and seven slugs whose parameters seed an initial condition switch from the ic axis to the parameter axis, unchanged to the digit. **The minimum bngsim is now 0.12.2.** One limit is deliberate: the two axes are only interchangeable in a common unit convention, and the ic axis is rescaled by each species' PyBNF-value-to- concentration factor while the parameter axis is not, so a species whose factor is not 1 keeps the ic route (substituting overstates its column by `1/factor`, measured across three compartment sizes). Exactly one benchmark model has non-unit factors and it seeds nothing. - **A gradient fit is refused, rather than silently doubled, on a bngsim build that seeds a lowered `initialAssignment` into the parameter axis (#537, ADR-0100).** The backend authors confirmed (lanl/bngsim#155) that `output_sensitivities(axis='parameter')` is the **total** derivative — `d_param[θ] = (RHS path) + Σ_j (∂x_j(0)/∂θ)·d_ic[x_j]` — so a route holding both a parameter's own axis and an IC term is correct only while the backend seeds nothing for that species. lanl/bngsim#147 changes that for compound parameter-only assignments, which it lowers to a synthetic `_ic_` derived parameter. `route_for_model` now detects exactly that build-and-model combination and refuses by name, so `Fiedler_BMCSystBiol2016` — the one slug in 23 whose free parameters legitimately route both axes, correct on every build through 0.12.1 at ≤3.7e-06 — fails loudly on the first build carrying #147 instead of reporting seven columns at the wrong value. The assembly's numeric guard was also generalized: it now compares the parameter slice against the **weighted sum** of the IC terms rather than a single slice, which catches a non-unit seed (`X(0) = a*X0` gives `d_param[X0] = 3.0·d_ic[X]`, agreeing to roundoff rather than bit-for-bit — the counterexample the first cut missed). The net backend is pinned to the same contract by test, having been verified bit-identical on `e2e_ode_decay.net`. - **A parameter's own sensitivity axis is not "identically zero" when it seeds an initial condition — it is the whole derivative, and adding the IC axis to it doubles the column (#537, ADR-0100).** ADR-0097 drops the `sensitivity_params` axis of a parameter that only seeds species initial values, documented as sparing a wasted vector on an axis that "would be identically zero". Measured on `Raia_CancerResearch2011`, that axis is not zero: bngsim seeds `∂x(0)/∂p` into it as well (lanl/bngsim#43, widened to compound `initialAssignment` expressions by lanl/bngsim#147), so `d_param[init_Rec_i]` and `d_ic[Rec_i]` come back **bit-for-bit identical** across every output. The drop is therefore load-bearing for correctness, and forcing both axes into the route reproduces #537's signature exactly — `init_Rec_i` at 2.00000× its central difference, every other column untouched. Three changes. The claim is corrected wherever it appears. The **unanswerable-model fallback is now a refusal**: ADR-0097 kept the axis when `ode_rhs_symbols()` could not say, on the strength of that false premise, which means one error deletes half a derivative (#535) and the other doubles it — with no safe default left, the router names the parameter and points at a gradient-free `job_type`. Both shipped backends always answer, so no fit changes. And the assembly now **checks numerically**, per experiment: a route holding both its own parameter axis and an IC axis whose tensor slices are bit-identical is refused, since two independent derivatives of a live model do not coincide to the last bit. `Fiedler_BMCSystBiol2016`, which legitimately routes both axes for seven parameters (their columns genuinely differ — RHS path versus seeding), is unchanged at ≤3.7e-06. ### Changed - **Tutorial Lesson 6 no longer teaches a refusal that stopped being one (#536).** Its whole premise was that a hard `if(t < tau, …)` step in a rate law is not differentiable, so `trf` refuses it and the fix is to smooth the step into a sigmoid. bngsim 0.12.2 differentiates it, and the lesson's three claims all failed when measured: the fit is not refused; its gradient is *correct* (agrees with a central difference to 3e-06 at a well-conditioned step — the apparent 3e-03 disagreement at `h = 1e-6` grows as `h` shrinks, which is roundoff, not a defect); and `tau` — which the lesson said "the hard-`if()` model could never expose to a gradient" — is recovered to 1.2e-08, because the solver contributes the crossing term where the switch fires. `step_input_trf_refused.conf` is renamed `step_input_trf.conf`, now fits `tau` alongside the rest, and is checked as a *recovery* rather than a refusal. The refusal half of the lesson is kept rather than deleted, because Lesson 3 sends readers here for it and because "what does a gradient fit refuse" is still worth teaching — it moves to a new `step_input_ssa_refused.conf`, a scored `method: ssa` experiment. That reason is durable in a way the old one was not: an SSA trajectory is a random walk, so there is no derivative with respect to the rate parameters to carry, and forward sensitivities exist only for the ODE backend. The smooth-sigmoid variant stays, re-pointed from "the differentiable one" to a **conditioning** choice — both fit the same parameters to the same accuracy and differ in what they cost the integrator. Without this the nightly `recovery` job would have gone red on its own the moment bngsim 0.12.2 reached PyPI, since CI installs bngsim unpinned. ### Added - **`sbml_rtol` and `sbml_atol` config keys (#546, ADR-0103).** State the CVODE tolerances for every deterministic run of every SBML/Antimony model in a fit, on the `sbml_backend = bngsim` path — which previously had no way to ask for either, since `atol`/`rtol` are read only out of a BNGL actions block and an SBML model has none. Unset (the default) leaves `rtol` at the backend default and derives `atol` from the model's own species magnitude, as described under Fixed above. Setting either under another `sbml_backend` is a config error rather than a silent no-op, and a model whose derived tolerance hits the `1e-16` floor is told once, by name, that `sbml_atol` is where to say so. - **A chain of two or more data-level normalizations on one column is now differentiable (#539, ADR-0102).** `normalization` is an ordered chain (ADR-0066), and five of its transforms — `peak`, `init`, `zero`, `unit`, `floor` — rewrite the column in place. Stacking two of them (`normalization pStat = floor 0.03, peak`) was refused on every gradient `job_type`, because the sidecar recording *how* a column was normalized kept one record per column: the second transform overwrote the first one's facts, and its intermediate values were gone with them. The sidecar is now the **list** of a column's transforms in chain order, and the gradient folds it — each stage's chain rule is the same closed form it always was, read in that stage's own inputs (the previous stage's per-row sensitivities) rather than in the raw ones, so the fold wraps the sensitivity accessor once per stage. The values a stage produced, which its own rule reads, ride along on its record when the next transform overwrites them; a column normalized once — every fit in the wild — retains nothing and costs what it did before, and its gradient is the same arithmetic in the same order. Every normalized *value* is unchanged: this is about what is recorded alongside them, not what is computed. `floor 0.03, scale` is unaffected either way, and a chain of any length still composes with the analytic `scale` (ADR-0099), the cumulative and per-measurement transforms, and the EFIM path. - **A gradient routing that would read one sensitivity column twice is now refused, not assembled (#537).** A route's derivative is the sum over its contributions, so two of them naming the same native `(axis, key)` column add that column twice — and the result is a clean integer multiple of the true derivative, which nothing downstream can detect: the objective stays finite, smooth and plausible, and the fit simply walks a scaled surface to a wrong answer. That is the shape a `Raia_CancerResearch2011` column came back in, once, during the #535 finite-difference sweep — exactly 2× its central difference on `init_Rec_i`, the fit's only initial-condition-axis parameter — and it has not reproduced in six runs at the same point, so there is no confirmed defect to fix. What there is now is a standing check on the narrow invariant it violated. Two halves: `route_experiment` **folds** terms that meet on one column into a single contribution carrying their summed derivative tree (the same sum, one tensor read, and `at_point` still refreshes a point-dependent path), which makes one-contribution-per-column structural rather than incidental; and the gradient, EFIM and constraint assemblers each `check_column_multiplicity()` on every routing they consume, once per experiment, raising a `PybnfError` that names the free parameter, the axis and key of the repeated column, and each duplicate's factor. The fold would otherwise blind that check — two same-column terms become one term of doubled factor, indistinguishable after the fact from a legitimate single term — so each contribution records the chain-rule path(s) it came from (`origins`: `bind`, or `ref:` per condition parameter-reference) and the check reads those. One path reaching one column twice is the defect; two paths meeting on it is arithmetic, and only the labels separate them once the factors are summed. Provenance is metadata: excluded from equality, hash and repr, so a routing still compares by what it computes. Note the scope: this covers the routing and assembly half of #537's hypothesis. A doubled column arriving from the backend's own sensitivity tensor would still pass, and remains open on the issue. - **A floored or analytically scaled observable is now a gradient target (#533, ADR-0099).** The two ADR-0066 normalization primitives shipped with deliberately deferred gradients, so a fit whose experiment declared either — `normalization = floor 0.03, scale`, the chain arbitrary-unit fluorescence / blot data is fit with — was unavailable to **every** gradient `job_type`, refusing with *"Analytic per-series scaling ('scale', #479) on column '…' is not differentiable on the gradient path"*. Both are now threaded. The **floor** (`x' = x + ρ·max x`) is additive, so every row picks up the same `ρ·s_argmax` term. The **analytic scale** is the one transform that is not per-point — its `c*` is profiled out of the whole matched series — so the scored value `c*(θ)·ŷ_i(θ)` differentiates by the product rule, with `∂c*/∂θ` the closed-form derivative of the profiling condition itself (the geometric-mean ratio for a log family, the least-squares optimum for a linear one), computed once per experiment and shared by every point of the column. That term does **not** drop out by the envelope theorem: the profiling is σ-unweighted, so `c*` is not in general the objective's own minimizer over the scale. A scaled fixed-σ Gaussian fit stays an exact least-squares model, so `trf`, `lbfgs`, and `gntr` all consume it; the profiled scale is resolved per experiment, so a column scaled in one experiment stays an ordinary column in another. One boundary was stated here rather than silently mis-differentiated — a chain of two or more *data-level* transforms on one column (`floor 0.03, peak`), which kept only the last one's facts — and is closed by #539 above. - **A total wall-clock budget for a fit, `wall_time_fit`, that finalizes on expiry (#529, ADR-0093).** PyBNF's time limits were all per unit of work — `wall_time_sim` bounds one simulation, `wall_time_gen` one network generation — and nothing bounded a *run*: the only native budget was `max_iterations` × `population_size`, which is not convertible to wall time without knowing per-iteration cost in advance. The new global key sets the seconds a whole fit may run (`wall_time_fit = 10800`; `0`, the default, is unbounded). On expiry the run stops launching work, abandons what is in flight, and runs the **normal** end-of-fit path against the best point found so far — `sorted_params_final.txt`, the best-fit simulations, `information_criteria.txt`, the ArviZ sidecar, the backup rename — so a budgeted result is scoreable exactly like a converged one. Only the stop *reason* differs, and it is logged, printed, and written to `Results/stop_reason.txt` (whose presence is the signal; no existing file's format changes). The clock starts when PyBNF starts, so configuration loading and network generation count against it, and one budget bounds the whole run: no `refine` and no further bootstrap replicate begins once it is spent. This makes PyBNF runnable under wall-time-budgeted optimizer benchmarks (Grein et al. 2026), where the previous alternative — killing the process — lost the artifacts scoring needs. Two overruns are deliberate and documented: one in-flight simulation may run up to `wall_time_sim` past the deadline before it is abandoned, and finalizing re-simulates the best fit once. Refused (rather than silently ignored) for `job_type = hmc`, which runs its own in-process sampling loop. - **A model with discrete events fits on the gradient path (#536).** An SBML `event` — and so an Antimony `at (…): …`, the usual way a dosing or stimulation schedule is written — used to refuse `trf` / `lbfgs` / `gntr` at construction. That refusal (#461) was right when it was written: a forward sensitivity carried across a state jump is correct only if the solver applies the event's own jump `s⁺ = ∂h/∂x·(s⁻ + f⁻·∂t*/∂p) + ∂h/∂p − f⁺·∂t*/∂p` at each fire, and bngsim did not, so it refused sensitivities on any event-bearing model and PyBNF hoisted that refusal to a clean pre-flight gate. bngsim applies the jump now, across a fixed trigger time, a trigger whose threshold is a fitted constant (lanl/bngsim#49), and a state-dependent trigger whose crossing it differentiates in flight (lanl/bngsim#144). The gate is therefore a **capability check** rather than a blanket structural refusal: an event-bearing model is allowed through, and the subclasses the build genuinely cannot cross (an execution delay; a trigger that is not a single relational comparison) keep a per-simulation refusal naming the reason. **The floor is set by silent wrongness, not by a missing feature.** A build that refuses is safe; a build that answers an event it cannot actually differentiate, without saying so, is the thing this gate exists to prevent. Three such answers had to go first, which puts the floor *newer than bngsim 0.12.1*: a trigger reading the state came back as a finite tensor with the event's contribution missing rather than being refused (lanl/bngsim#52, through 0.11.x); an assignment reading the state — `A := A + dose`, the repeat-dosing idiom — dropped the carried `∂h/∂x·s⁻` and restarted the assigned row from zero, measured on 0.12.1 at `-10.96` against the model's own central difference of `-311.20` while the identical model built through `ModelBuilder.add_event` was right to `2e-6` (fixed after 0.12.1 by lanl/bngsim#144's jump-handler rework); and a solver root that fires nothing rewound the state but not the sensitivity history (lanl/bngsim#146, also after 0.12.1). On 0.12.1 or older the refusal stays, and stays blanket, with the message naming the upgrade. The floor lives in `pybnf._bngsim_caps` (`BNGSIM_HAS_EVENT_SENS`) and is not a dependency bump: the install floor stays `bngsim>=0.11.35`, and every scalar (metaheuristic) fit is unaffected either way. The SBML backend also gained a narrowed form of the net backend's #525 wrapper, which the lifted gate makes reachable: bngsim declining to differentiate a model's events is a *structural* verdict, identical at every parameter set, so it now surfaces as an actionable `PybnfError` instead of a `FailedSimulationError` the optimizer scores `inf` and steps around — which would have reported an unsupported event as a failed search. Every other backend failure keeps that back-off, which is the right answer for a candidate point the integrator cannot get through (#492). New coverage is a finite-difference oracle on a two-species event fixture at both levels, the backend tensor and the assembled objective gradient (`tests/test_gradient_events.py`): every sensitivity column is scored against a central difference of PyBNF's own trajectory / own loss, including the columns that cross the jump, which is the only instrument that catches a jump term that is missing rather than merely inaccurate. The backend-level oracles run on every build, since what bngsim computes is worth asserting whether or not PyBNF admits the model; the fits, and the bolus-assignment case that set the floor, are gated on `BNGSIM_HAS_EVENT_SENS`. - **Published-source organization and two curated real-world jobs.** The real-world gallery now uses `Author-Year/job_slug` paths aligned with the BNGL-Models job corpus. The former flat Kozer EGFR, Monine TLBR, Gupta FcεRI, and Mitra receptor jobs retain their tested edition-2 configurations under source-oriented collections; curated provenance, validation notes, and reproduction assets accompany the Kozer and Monine jobs. The reduced F5B-only IGF1R teaching fit is replaced here by `Erickson-2019/igf1r`, the authors' published seven-rate, three-dataset preincubate→wash→dose-scan fit (the reduced job remains under `examples/igf1r/`). Two additional workstation examples broaden the executable corpus: `Salazar-Cavazos-2019/egfr_simpull` adds an authors' multisite-EGFR ODE fit, and `Kirsch-2020/phosphoswitch_bpsl` adds a four-model, constraint-only BPSL fit. The default corpus test now distinguishes quantitative `.exp` jobs from qualitative `.prop` jobs. - **Workstation-scale exact-SSA real-world examples (#472).** The new `examples/real-world/Rijal-2025/` collection fits lacUV5/lacUD5 and 5DL1 promoter-noise data from Jones et al. (2014) with the two-state model studied by Rijal and Mehta (2025). Each edition-2 SSA job uses `method: ssa`, 200-trajectory smoothing, and measurement formulas for ensemble mean and standard deviation; paired exact moment-ODE jobs provide deterministic references, with source tables, regeneration scripts, validation notes, and reproduction figures preserved alongside them. All four configurations receive backend-free corpus checks, and the bounded SSA fits enter the opt-in bngsim recovery tier, closing the gap left by the cluster-scale FcERI SSA reference. - **CMA-ES gains an optional bounded per-run generation budget (`cmaes_run_maxgen`; #507, ADR-0085).** The global `max_iterations` budget previously left the initial run and every IPOP / BIPOP large run unbounded, so one run making slow progress in an ill-conditioned local basin could consume nearly the entire fit before later restarts launched. Setting the new positive-integer cap applies it to every run and turns reaching it into the existing per-run restart trigger; the final run then stops at the same cap. BIPOP small runs use the smaller of this user cap and their existing automatic evaluation-balancing cap. The default is unset, preserving the prior schedule and results. - **Prediction-dependent noise: a `noise_model … = , sigma = prediction_formula ` source whose σ scales with the simulated output (#495, ADR-0075).** The honest combined additive+proportional error model `σ = σ_abs + σ_rel · y` — where `y` is the observable's *predicted* value — previously had no native form: `formula` (ADR-0044) evaluates only over free parameters, and `relative` / `column_mean` (ADR-0031) read the *data*, not the simulation. The new `prediction_formula` verb builds a `PredictionFormulaSigma` whose symbols resolve either from the PSet (the estimated coefficients) or from the current simulation column of that name (a model species / observable / function), evaluated per scored point. Any new-era job may author it; PEtab import is one way to reach it. (Gradient-free score path only — a prediction-dependent σ raises `GradientNotSupported` on the #385 gradient/EFIM path, a later sub-layer.) - **Scale-preserving PEtab v1→v2 conversion: `pybnf.petab.petab1to2_preserve_scale`.** The official `petab.v2.petab1to2` **drops** the v1 `parameterScale` column (PEtab v2 removed it) and only *warns* — so a `parameterScale = log10` estimated parameter carrying no objective prior (the common case for a multi-decade kinetic parameter) converts to a *linear* `uniform_var` over the raw bounds: the same argmin, but a far harder, worse-conditioned optimization than the log10 search the modeler specified. This wrapper runs the standard converter and re-injects the dropped estimation scale in the **v2-native** form — `priorDistribution = log-uniform` over each such parameter's bounds — which PyBNF imports as a `loguniform_var` on the Log10 scale. Because the optimizer objective excludes the prior, this sets only the search scale and initial sampling, not the objective, so the fit stays the pure-MLE problem v1 specified; parameters petab1to2 already folded into a prior (`parameterScale*Normal` → `log-normal`, …) are left untouched, as are linear ones. `import_job` stays a pure v2 importer — the conversion is an explicit, named step, not a reach-back to v1 in the read path. Intended as the opt-in migration `petab1to2` itself should offer. - **General-objective trust-region optimizer: `fit_type = gntr` (ADR-0068, #481).** Fills the last empty cell of the gradient-fitting (objective × curvature-model) matrix. `trf` gives a trust-region step with a `JᵀJ` (Gauss-Newton / empirical-Fisher) Hessian but only for an exact least-squares objective; the moment the objective stops being a pure sum of squares — an estimated noise scale, a Laplace / count likelihood, or an active constraint — the gradient path dropped to `lbfgs` (limited-memory quasi-Newton). `gntr` extends `trf`'s well-conditioned trust-region step to those general-NLL objectives: its Hessian is the **expected-Fisher / Gauss-Newton information** `H = Σ κᵢ sᵢsᵢᵀ` (+ estimated-noise and constraint blocks), built from the same #385 forward sensitivities `sᵢ = ∂predᵢ/∂θ` plus small analytic per-family curvature factors (new `NoiseModel.location_fisher` / `noise_param_fisher` and `Constraint.penalty_curvature` seams) — no second-order sensitivities. It consumes the *same scalar gradient* as `lbfgs`; only the curvature differs. Internally it reuses `trf`'s Coleman–Li reflective machinery unchanged by feeding `(g, H)` through a ridge-regularised pseudo-Jacobian, so on a Gaussian least-squares fit it reduces to `trf`'s step exactly. It runs natively in the distributed propose/score loop (picklable, no `run()` override, concurrent `N`-start multi-start, registered as a box-start refiner) like every other `fit_type`. New config keys `gntr_grad_tol` (1e-8), `gntr_step_tol` (1e-8), `gntr_ridge` (1e-10), and the runtime-guarded `gntr_max_iterations`. This cut supports an estimated-σ Gaussian (`chi_sq_dynamic`), a fixed-scale Laplace, a fixed-dispersion mean-centered negative-binomial, and a Gaussian fit with static-hinge constraints; the coupled corners it cannot yet build the Fisher Hessian for (a mean-on-log estimated scale, a free-dispersion / median count family, an estimated Student-t df, or an estimated constraint scale) refuse with a pointer to `lbfgs`, which fits them. `trf` / `lbfgs` are byte-identical (they never form the Hessian). - **Kalman-inspired DREAM proposal: `proposal = kalman` (ADR-0067, Stage 3; DREAM(KZS), #358).** The third proposal operator on the unified DREAM engine (Zhang, Vrugt et al. 2020). During a burn-in window each proposal is steered toward the data by a Kalman gain `K = C_ZY (C_YY + R)⁻¹` built from the archive's parameter↔model-output cross-covariance, with the innovation `d - f(xᵢ) + ε` taken at the chain's current state (`ε ~ N(0, R)`), which accelerates burn-in on informative, mildly non-linear problems; after the window the chain reverts to `de` for a reversible sampling phase (the Kalman jump breaks detailed balance by design, so its samples are burn-in and discarded). The gain reads each archive entry's *model output vector* `f(Z)` — surfaced by the new `LikelihoodObjective.aligned_prediction_data` seam and carried in an output-augmented archive that turns on only for this proposal (the "implied axis 2b"; dormant and byte-identical for `de` / `whitened`). `kalman` requires a linear-scale Gaussian likelihood (`chi_sq` / `chi_sq_dynamic`, the source of `R = diag(σ²)`) and `n_try = 1`, and refuses any other objective or `n_try > 1` *before the run starts*. The internal ensemble size is fixed (`M = 20`, clamped to the available archive, falling back to `de` before enough outputs accrue — no new user key); one new proposal-scoped key `kalman_burnin_frac` (default `0.3`) sets the window as a fraction of `burn_in`. Validated end-to-end against a closed-form linear-Gaussian posterior (`f(x) = A x` scored by real `chi_sq`), plus pinned gain-math and burn-in-switch unit tests. - **Multi-Try DREAM: the `n_try` count (ADR-0067, Stage 2; MT-DREAM(ZS), #357).** A new integer `n_try` config key turns each chain-generation into a multiple-try step (Liu, Liang & Wong 2000; Laloy & Vrugt 2012): with `n_try = k > 1` a chain draws `k` candidate proposals, selects one in proportion to its posterior importance weight, and accepts it over the current state with a multiple-try Metropolis ratio evaluated against a `k - 1`-point reference set drawn from the winner plus the current state (`2k - 1` evaluations per chain per generation). Multiple tries per generation raise the per-generation acceptance rate and help parameter-rich / strongly correlated posteriors mix. It is the second orthogonal axis of ADR-0067 and **composes with every `proposal` value** (`de`, `whitened`) and with the snooker update — MT-DREAM(ZS) is literally multi-try parallel-DE. `n_try = 1` (the default) is the classic single-try engine and is **byte-identical** to before (verified against the DREAM/P-DREAM oracle suites and the effective-config goldens; the only change is the additive `n_try` key). The snooker proposal is non-symmetric, so under multi-try its candidate and reference weights carry the ter Braak & Vrugt (2008) Jacobian `||p - z||^(d-1)`; the current-state reference slot uses the current state's distance to the **selected candidate's** anchor — the unique choice that reduces to the published single-try snooker ratio at `k = 1` (derived from first principles and confirmed by a stationary-distribution test; both the DREAM-Suite and PyDREAM reference implementations differ on this term). The multiple-try acceptance is validated to preserve a known Gaussian target with the snooker update active. - **DREAM `proposal` operator key; P-DREAM folded into one DREAM engine (ADR-0067, Stage 1).** DREAM(ZS) and Preconditioned DREAM are now one `DreamAlgorithm` engine selected by a new `proposal` config key: `proposal = de` (default) is the classic parallel-direction proposal, and `proposal = whitened` is the covariance-preconditioned proposal that used to be a separate algorithm. The `p_dream` job type is unchanged for users — it is simply `dream` with `proposal = whitened` pinned — and `whitened` can now also be requested explicitly on a `dream` run. This is a pure refactor: `dream` at defaults and `p_dream` are **byte-identical** to before (verified against the existing DREAM/P-DREAM oracle suites and the effective-config goldens). It is the first step of ADR-0067's unification of the DREAM family into two orthogonal axes (`proposal` × `n_try`), which will absorb the requested MT-DREAM (#357) and DREAM-KZS (#358) without new sampler subclasses. - **Composable floor normalization + analytic per-series scaling for relative / arbitrary-unit data (#479).** Two composable, per-series normalization primitives so a log/relative objective on arbitrary-unit data (fluorescence, blots) can be spelled with standard tokens instead of a bespoke objective class. (1) **`floor `** — an additive measurement-noise floor `x' = x + rho*max(x)` (default `rho = 0.03`) applied **identically to the simulated and the experimental** column, so a log objective stays finite where a series legitimately touches zero. (2) **`scale`** — analytic per-series **optimal** multiplicative scaling profiled out at scoring time (hierarchical / profiled scaling; Weber et al. 2011, Loos et al. 2018), family-appropriate: the geometric-mean ratio for a log family (`lognormal`) and the least-squares optimum `c* = Σ w s d / Σ w s²` for a linear one — so an overall model-vs-data scale difference is not penalized (no per-series `scale` parameter needed). They compose as an ordered chain (`normalization = floor 0.03, scale`, per-observable / `.` / whole-fit), and together with `objective = lognormal` spell the exact sum-of-squared-log-differences-of- geometric-mean-normalized-trajectories objective of Jaruszewicz-Błońska et al. (*PLoS ONE* 2023; 18(6):e0286416). Legacy `normalization = peak` / `normalization x = peak` round-trip byte-identically; `peak`/`init`/`zero`/`unit` stay sim-only. The `peak`, `unit`, and `floor` column reductions are **NaN-aware** (`np.nanmax`/`nanargmax`, etc.), so a sparse multi-observable target — NaN in the rows where a given observable is unmeasured — is reduced over its measured points only rather than collapsing the whole column to NaN (which had silently zeroed the objective); a dense column is byte-identical. Both new primitives have a **deferred gradient** (they raise `GradientNotSupported`, so a gradient fit falls back to a gradient-free step; the motivating fits are evolutionary), and both are **refused on PEtab export** (a whole-trajectory reduction has no pointwise PEtab v2 operator; `scale`'s `observableParameters` mapping is a future direction). See ADR-0066. - **Gradient-based fitting extends to `parameter_scan` (dose-response) objectives (#476).** A gradient fit (`fit_type = trf`/`lbfgs`) can now target a dose-response objective, not just a time course. The default dose-response path already computed the per-dose forward sensitivities `∂obs(dose)/∂θ` — one sensitivity-configured ODE `run()` per swept dose — and then discarded them at row assembly; PyBNF now stacks those per-point final-row sensitivities down the dose axis into the scan `Data`, so the existing gradient assembly produces `d(objective)/dθ` for dose-response fits. The swept dose is the data's independent variable (not a fitted parameter), so the per-dose sensitivity is well-posed and consumed exactly as a time-course row is. Supported for the **reset-to-seed** strategies — the parity / integrate-to-steady-state default and the independent fixed-time scan — on both the native BNGL and SBML/Antimony backends. Newton/KINSOL (`ss_method=>"newton"`, now supported — see #478 below), continuation/bifurcate (`reset_conc=>0`), `method=>"protocol"`, and carried-state (pre-equilibration, ADR-0062) scans refuse cleanly on the gradient path with an actionable message (an *incidental*, unscored scan of the same shape still runs sensitivity-free, #475). The scalar (metaheuristic) path is byte-identical. See ADR-0064. - **Scored Newton/KINSOL (`ss_method=>"newton"`) steady-state dose-response scans are now differentiable (#478).** The KINSOL accelerator solves each dose point's steady state as an algebraic `f(x)=0` (no forward-sensitivity *integration*), so #476 (ADR-0064) kept a *scored* Newton scan gradient-free and pointed at the parity default. It is now a real speed win under a gradient fit: the KINSOL solve returns `dY_ss/dp` **exactly** (the implicit-function-theorem derivative on the analytical Jacobian, not a finite difference), and bngsim ≥ 0.11.35 (lanl/bngsim#12) maps it through the observable/function Jacobian `∂g/∂x` and exposes it as `SteadyStateResult.output_sensitivities`, mirroring the CVODE `Result`. PyBNF stacks those per-dose slices down the dose axis exactly as the parity path does — no gradient-assembly change. On the gradient path the scan runs sequentially (the KINSOL sensitivity solve is kept off the thread pool) and the KINSOL→CVODE non-convergence fallback is itself differentiable and consistent with the converged path. **Requires bngsim ≥ 0.11.35**; a build lacking the accessor refuses a scored Newton scan cleanly with an upgrade hint (a scalar Newton scan is unaffected). Continuation/bifurcate, `method=>"protocol"`, and carried-state scans still refuse on the gradient path. See ADR-0065. - **Edition-2 preincubate → wash → dose-response scan protocol (#474).** The new-era `experiment:`/`condition:` surface now expresses the full **equilibrate → intervene → measure a dose-response** protocol, so a published fit that needs it (the Erickson-2019 IGF1R competition/dissociation fit — 7 rate constants to 3 datasets, two of them a 2 h-preincubate → wash → cold-competition scan) runs in `edition = 2` with **no in-model actions block**. Two capabilities: (A) a `parameter_scan` may be the measured phase of a `preequilibrate:` experiment — the synthesizer emits `saveConcentrations()` + `parameter_scan(… reset_conc=>1)` so each dose resets to the carried post-intervention state; (B) a `condition:`'s `perturbations:` accepts a **quoted BNGL species pattern** target with a number *or a parameter-expression* value — a species `setConcentration` (a wash `"IGF1(ds,hs,label~hot)" = 0`, or a dose-tracking bolus `"IGF1(ds,hs,label~cold)" = IGF1_cold_conc*(NA*Vecf)`), vs. a parameter `setParameter`. The bngsim backend routes such a carried-state scan to its native reset-conc-to-snapshot `parameter_scan`/`bifurcate` (**requires bngsim ≥ 0.11.34**, lanl/bngsim#11), reproducing BNG2.pl exactly; the fresh-from-seed dose-response paths (ADR-0046) are unchanged. See ADR-0062. - **PEtab v2 export/import of the preincubate → wash → dose-scan protocol (#477).** The two shapes ADR-0062 added to the edition-2 fitter now export to PEtab v2, import back, and round-trip byte-for-byte (validated by petab's full `default_validation_tasks`): (1) a **species `setConcentration`** condition target — a BNGL species pattern is not a valid PEtab id, so it is aliased through the **mapping table** (`petabEntityId` → the pattern) and the condition targets the synthesized `species_<…>` id with a number or a parameter-expression value; (2) a **pre-equilibrated dose-response** — each dose becomes a two-period Experiment (a `time = -inf` pre-equilibration period + a measurement period applying both the shared wash condition and a per-dose swept-parameter condition), the combination of ADR-0052 and ADR-0046. The exporter's previous "deferred" refusals are lifted. The surrogate split × a pre-equilibrated scan (an empty surrogate set M is required) and a whole-fit `normalization` transform (the real Erickson-2019 IGF1R job) stay out of scope, raised in code. See ADR-0063. - **`examples/real-world/` — the 2019 PyBNF-paper case studies on the edition-2 surface.** The biological models from Mitra et al. (iScience 2019) — Kozer's EGFR (ODE and network-free), the ligand/receptor model (ODE and NFsim), IGF1R competition binding, the FcεRI γ-chain SSA network, and the trivalent-ligand aggregation model — re-expressed on the new-era `experiment:`/`condition:`/`data:` config surface, spanning the three simulator paths (deterministic ODE, Gillespie SSA, network-free NFsim). These validate PyBNF's bngsim-backed default path on representative, paper-scale models (issue #380), with `tests/test_real_world_examples.py` running a backend-free well-formedness tier in default CI and a real-bngsim end-to-end tier under `-m recovery`. - **Network-free (NFsim) options on the edition-2 experiment surface.** A `method: nf` experiment now accepts `gml:` (global molecule limit) and `complex:` (track molecular complexes) — the network-free counterparts of `atol`/`rtol` — carried into the synthesized NFsim `simulate`/`parameter_scan` so a large aggregating model (e.g. the EGFR clustering fit) can raise its molecule limit and track complexes as its classic hand-written action did. - **Fixed-time NF pre-equilibration (`equil_t_end:`).** Edition-2 pre-equilibration (ADR-0052) equilibrates *to steady state*, but NFsim has no steady-state solve. A `method: nf` pre-equilibration now takes `equil_t_end: