Best-fit confirmation ranks candidates on surviving replicates only, so a parameter set that fails 9 of 10 runs is pinned as the run's best fit over one that succeeds every time

[wshlavacek]

What happens

mean_objective (lines 50-52) is the aggregate the whole stage ranks on:

if not candidate.scores:
    return math.inf
return fmean(candidate.scores)

candidate.scores contains only the replicates that produced a usable value — Algorithm._run_confirmation_replicates in pybnf/algorithms/base.py increments failures[index] and continues for a FailedSimulation, and for score is None or not np.isfinite(score). failures is stored on the Candidate and printed as a report column, but it appears in no comparison: ranked (line 83) keys on mean_objective alone, and winner (86-95) disqualifies a candidate only when it has zero usable scores. There is no minimum-success threshold and no penalty, so the aggregate is a survivor-only mean with a sample size that varies per candidate from 1 to best_fit_replicates.

The consequence is what the stage was built to prevent. Its own docstring (lines 6-14) says the search goes wrong because "every one of those values came from a single simulation", so "the winner of that comparison is very often the parameter set that happened to get a lucky simulation". A candidate that fails nine replicates and survives one is scored on exactly one simulation — the same max-of-lucky-draws estimate, now with the unlucky draws deleted rather than averaged in — and it beats a candidate measured honestly over ten runs. Algorithm._emit_best_fit_confirmation then calls trajectory.pin_best(best.pset, mean_objective(best), best.name), which is by design how "everything downstream ... reports the same parameter set", so the saved simulations, the best-fit model file, the information criteria, a refine's start point and a bootstrap replicate's answer all describe the failure-prone set.

It is also silent at the point of decision. console_lines (169-195) never mentions failures, and for a one-survivor winner standard_error is None, so the documented "raise best_fit_replicates if the standard errors overlap" guidance cannot flag it; the console instead asserts "the one the search liked best does worse when it is run again", which in this scenario is false. _emit_best_fit_confirmation logs a warning only when no candidate has a usable score.

For contrast, PyBNF's other replicate-averaging path takes the opposite line: JobGroup.job_finished/average_results (pybnf/algorithms/core.py:442-465) turns a whole smoothing group into a FailedSimulation the moment one sub-run fails, rather than averaging the survivors.

Reproduction

Read-only repro, run from /Users/l119605/Code/PyBNF with
uv run --extra tests --extra petab python <file>
(script saved at /private/tmp/claude-503/-Users-l119605-Code-PyBNF/eb55c5a3-1006-4c5d-9d1e-9d14fe8045f0/scratchpad/exp.py):

import sys, os, tempfile
sys.path.insert(0, '/Users/l119605/Code/PyBNF')
from tests.test_best_fit_confirmation import (_algo, _ps, _FakeClient, _FakeAsCompleted, NOISE)
from tests.context import algorithms
algorithms.core.as_completed = _FakeAsCompleted   # synchronous dask

NOISE.clear()
for r in range(1, 11):
    NOISE[(1.0, r)] = float('inf')   # candidate A: 9 of 10 replicates unusable
NOISE[(1.0, 7)] = 5.0                #             one lucky survivor at 5.0
for r in range(1, 11):
    NOISE[(2.0, r)] = 5.1            # candidate B: 10 of 10 succeed at 5.1

algo = _algo(tempfile.mkdtemp(), candidates=2, replicates=10)
algo.trajectory.add(_ps(2.0), 5.05, 'reliable')
algo.trajectory.add(_ps(1.0), 5.20, 'flaky')
algo._confirm_best_fit(_FakeClient())
print(algo.trajectory.best_fit_name(), algo.trajectory.best_score(),
      algo.trajectory.best_fit()['v1__FREE'])
print(open(os.path.join(algo.res_dir, 'best_fit_confirmation.txt')).read())

Observed:
flaky 5.0 1.0
winner flaky
winner_mean_objective 5
winner_standard_error n/a
# rank name mean_objective standard_error std_deviation runs failed search_objective
1 flaky 5 n/a n/a 1 9 5.2
2 reliable 5.1 0 0 10 0 5.05
and on the console, with no mention of the 9 failures:
Best average objective 5 (flaky).
This is not the parameter set the search would have reported. It was number 2 in the
search ranking, and the one the search liked best does worse when it is run again.

Expected: the parameter set that produced a usable objective in 1 of 10 runs must not be
pinned as the run's best fit over one that produced 5.1 in 10 of 10 — or, at the very least,
the stage must not report it as the confirmed answer without saying that 9 replicates
produced nothing. The reliable candidate's expected objective is 5.1; the flaky one's is
effectively infinite.

Note: the real-run trigger for the dropped values is either a FailedSimulation (crash or
wall_time_sim timeout, both genuinely intermittent for SSA/NFsim) or a non-finite objective
(base.py if score is None or not np.isfinite(score)), which is easy to hit for a stochastic
model whose trajectory occasionally hits zero counts under a log-based objective. The inf
values above stand in for that second case exactly as the production path treats it.

Verification notes

CONFIRMED by running the code path end to end, not by reading alone.

The aggregation is mean_objective (best_fit_confirmation.py:50-52): if not candidate.scores: return math.inf / return fmean(candidate.scores). candidate.scores holds only the replicates that produced a usable value — Algorithm._run_confirmation_replicates (base.py ~1833-1843) does if isinstance(res, FailedSimulation): failures[index] += 1; continue and if score is None or not np.isfinite(score): failures[index] += 1 else: scores[index].append(...). ranked (line 83) sorts on that mean; winner (86-95) rejects a candidate only when scores is entirely empty. failures is carried on the namedtuple and printed in a report column, but it enters no comparison anywhere. So the ranking that decides the run's answer is over survivors only, with no minimum-success guard.

Experiment (ran it, /Users/l119605/Code/PyBNF harness): two candidates, 10 replicates each. "flaky" (v1=1.0) returns a non-finite objective on 9 of 10 replicates and 5.0 on one; "reliable" (v1=2.0) returns 5.1 on all 10. Result: winner flaky, table row 1 flaky 5 n/a n/a 1 9 5.2 above 2 reliable 5.1 0 0 10 0 5.05, and _emit_best_fit_confirmation calls trajectory.pin_best(...) so best_fit_name() == 'flaky', best_score() == 5.0, best_fit()['v1__FREE'] == 1.0. Everything downstream of _confirm_best_fit in run() — saved sims, best-fit BNGL, information criteria, refine start, bootstrap answer — then describes that parameter set.

It is silent, not merely a disclosed tradeoff. console_lines never mentions failures; with one survivor standard_error is None so it printed only "Best average objective 5 (flaky)." — no hint that 9 of 10 runs produced nothing. Worse, it printed "This is not the parameter set the search would have reported ... the one the search liked best does worse when it is run again", which is the opposite of the truth here. _emit_best_fit_confirmation warns only when no candidate has any score. And _compute_information_criteria (base.py:2060, "A replicate that fails or scores nothing is left out") applies the same rule, so the wrong winner is not caught downstream either.

Evidence considered on the other side, and why it does not settle it as intended behavior:

• The module docstring (lines 38-39) and mean_objective's docstring do say failures are left out, and tests/test_best_fit_confirmation.py:479 (test_a_failed_replicate_is_counted_and_the_rest_still_decide) deliberately pins 1-of-3 failing -> the surviving 2 decide. But that test has a SINGLE candidate. Nothing in tests/ or docs pins or even mentions cross-candidate ranking under unequal failure counts; test_a_candidate_with_no_usable_value_sorts_last_and_never_wins covers only the all-fail extreme. "Don't let one bad draw kill a good candidate" is a sound intent that the code overshoots into "let a 1-in-10 candidate beat a 10-in-10 one".
• The codebase's own convention for averaging replicates of one parameter set is the opposite: JobGroup.job_finished/average_results (pybnf/algorithms/core.py:442-465) makes a whole smoothing group a FailedSimulation as soon as any sub-run fails. So survivor-only averaging is not a house-wide rule this stage is following.
• Dropped values are not neutral noise: a non-finite objective or a crashed/timed-out stochastic run is an informative, bad outcome, and dropping it biases the estimate downward. Selection on max-of-draws is exactly what the stage exists to undo (module docstring lines 6-14), and this reintroduces it through the failure channel for the one population it operates on — the search's top-N, which over-selects lucky, boundary-adjacent parameter sets in the first place.

Confidence is "likely" rather than "certain" only because the fix involves a policy choice (require a minimum number of successes, rank by (failures, mean), or count a failure as inf); a maintainer could call the current rule deliberate. The mechanism itself is not in doubt — it is demonstrated above.

Where

pybnf/algorithms/best_fit_confirmation.py:50 — severity medium, confidence likely. Repro executed: True.