diff --git a/CHANGELOG.md b/CHANGELOG.md index a7a17329..56f2b032 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,26 @@ All notable changes to PyBNF are documented below. This project adheres to ## [Unreleased] ### Added +- **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 @@ -37,6 +57,13 @@ All notable changes to PyBNF are documented below. This project adheres to by default. Both surfaces are documented under gradient-based fitting. ### Fixed +- **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 diff --git a/docs/algorithms.rst b/docs/algorithms.rst index 5dd910b1..7f3ea63b 100644 --- a/docs/algorithms.rst +++ b/docs/algorithms.rst @@ -171,6 +171,12 @@ asynchronous variant ``ade``, ``ss``, and ``pso``). ``n_starts = 1`` (the defaul single run. For ``cmaes`` the equivalent knob is its own ``cmaes_restarts`` (IPOP/BIPOP restart). +A fit with more than one start writes ``Results/multistart_summary.txt``, one row per start +sorted by objective value from best to worst, and prints a short version of it at the end of +the run. That is how you tell a search whose starts all agreed from one whose starts all +landed somewhere different, which is the difference between an answer worth believing and +the least bad of several poor ones. See :ref:`the per-start summary `. + .. _alg-ss: Scatter Search diff --git a/docs/config_keys.rst b/docs/config_keys.rst index 5a25ad64..ea0ab02b 100644 --- a/docs/config_keys.rst +++ b/docs/config_keys.rst @@ -650,6 +650,10 @@ Required Keys **n_starts** Number of independent multi-start runs for the metaheuristic optimizers (``de``, ``ade``, ``ss``, ``pso``). A single run collapses its population into one basin, so on a multimodal objective it returns only a local minimum. With ``n_starts > 1``, that many independent searches are run one after another -- each a fresh random / Latin-hypercube population, each up to ``max_iterations`` iterations or until it converges -- and the best fit over all of them is kept. ``1`` (the default) is a single run, identical to the historical behavior. (``cmaes`` has its own multimodal restart, ``cmaes_restarts``; the gradient optimizers use ``population_size`` as their start count.) + ``n_starts`` also selects the number of starts for the concurrent local optimizers ``powell`` and ``sim``, which run their starts at the same time rather than one after another. + + A fit with more than one start writes ``Results/multistart_summary.txt``, one row per start sorted by objective value from best to worst, and prints a short version of it at the end of the run. See :ref:`the per-start summary `. + Example: * ``n_starts = 10`` diff --git a/docs/gradient_fitting.rst b/docs/gradient_fitting.rst index 6625084f..be24bbe8 100644 --- a/docs/gradient_fitting.rst +++ b/docs/gradient_fitting.rst @@ -124,6 +124,20 @@ to a standalone box-start fit: when the optimizer runs as a **refiner** (an expl injected) it always runs a single start, since the job there is to polish the one best fit, not to re-scatter. ``max_iterations`` is the per-start iteration budget. +**How the individual starts did.** The reported best fit is one number, and on its own it does not +say whether the search can be trusted: twenty starts that all reached about the same objective value +and twenty starts that all landed somewhere different look identical in it. So the run writes +``Results/multistart_summary.txt``, one row per start sorted by final objective value from best to +worst, and prints a short version at the end:: + + Multi-start summary: 20 starts. Best objective 143.21, median 143.22, worst 891.4. + 17 of the 20 came within 0.1% of the best. The more starts that reach the same low value, + the more likely the fit has found the best answer available. + Per-start table, best first: output/Results/multistart_summary.txt + +See :ref:`Should I believe my fit, or run it again with more starts? ` for how +to read the table. + **A bad start is survivable.** Stiff corners of a parameter box can defeat the solver, and a fit is expected to walk into them: a point may fail to integrate at all, or integrate while its forward sensitivities diverge, leaving a finite objective with a non-finite gradient. Neither ends the fit. diff --git a/docs/troubleshooting.rst b/docs/troubleshooting.rst index 0966ab8c..f5b4742b 100644 --- a/docs/troubleshooting.rst +++ b/docs/troubleshooting.rst @@ -140,6 +140,49 @@ it leaves the original run's record intact, since a resume has no start of its o stopped at the next polling boundary. +.. _multistart_summary: + +Should I believe my fit, or run it again with more starts? +---------------------------------------------------------- +A fit that runs several searches from different starting points reports the best of them, and +that one number cannot be checked on its own. Two very different runs look identical in it: +one where every start reached about the same objective value, and one where every start landed +somewhere different and the reported answer is only the least bad of them. + +So a fit with more than one start writes ``Results/multistart_summary.txt``, one row per start +sorted by objective value from best to worst:: + + starts 20 + best_objective 143.2085 + worst_objective 891.4471 + median_objective 143.2166 + starts_near_best 17 + # rank start objective iterations evaluations reason + 1 7 143.2085 38 402 gradient is flat (‖v·Jᵀr‖∞ ≤ 1e-08) + 2 2 143.2091 41 437 gradient is flat (‖v·Jᵀr‖∞ ≤ 1e-08) + ... + 20 13 891.4471 60 631 reached max_iterations (60) + +Read the objective column downward. A long flat run of starts at the same low value means the +search found a consistent answer and more starts would very likely not change it. A staircase +with no flat section at the top means it did not, so run more starts, or take another look at +the model and at the parameter bounds. ``starts_near_best`` counts the flat run for you, and +a short version of all of this is printed at the end of the run. + +Two rows deserve a second look wherever they appear. ``inf`` in the objective column is a start +that produced no usable fit at all, usually because its start point failed to simulate; if every +row says that, the parameter box is somewhere the model cannot be integrated. A reason of +"did not finish before the fit ended" or "did not start before the fit ended" means the run +stopped before that start was done with, which is what ``wall_time_fit`` expiring looks like +here. + +This applies to ``trf``, ``lbfgs``, ``gntr``, ``powell``, ``sim``, ``ms``, the polishing phase +of ``profile_likelihood``, and the metaheuristics ``de``, ``ade``, ``ss`` and ``pso``. The +number of starts comes from ``population_size`` for the gradient methods and from ``n_starts`` +for the rest. A refine writes ``Results/multistart_summary_refine.txt``; a fit with a single +start writes nothing, since there would be nothing to compare it against. + + Unexpected behavior when generating SBML files in COPASI -------------------------------------------------------- While COPASI is a useful tool for generating SBML files, it is important to note that some settings in COPASI do not get converted into SBML. This can lead to unexpected model behavior in PyBNF. diff --git a/pybnf/algorithms/base.py b/pybnf/algorithms/base.py index 90003e5c..abdb7a4b 100644 --- a/pybnf/algorithms/base.py +++ b/pybnf/algorithms/base.py @@ -8,6 +8,7 @@ from . import core +from . import multistart_report from .core import ( FailedSimulation, JobGroup, @@ -1444,6 +1445,10 @@ def _finalize_run(self): :attr:`trajectory` and :attr:`stop_reason` and nothing about how they were filled. """ self._announce_stop_reason() + # How every start of a multi-start fit did, next to why the fit stopped -- both are + # statements about the run rather than about its best parameter set, and a fit whose + # starts all failed still has this to report even though it has no best fit below. + self._emit_multistart_summary() # Write the final parameter sets, then copy the best simulations into the results # folder. A run that stopped before any result came back (an expired budget, an @@ -2005,6 +2010,56 @@ def _emit_information_criteria(self, ic, name='', preamble=()): % (ic.aic, ic.bic, aicc_str, ic.k, ic.n, ic.log_likelihood)) return True + def multistart_records(self): + """One :class:`~pybnf.algorithms.multistart_report.StartRecord` per start of this + fit, in any order -- :meth:`_emit_multistart_summary` sorts them (#658). + + Empty here, which is right for the great majority of fit types: they run one + search and have nothing to compare it against. The three families that run several + starts override this, each reading the numbers off wherever it happens to keep + them (see :mod:`pybnf.algorithms.multistart_report`). + """ + return () + + def _emit_multistart_summary(self): + """Write ``Results/multistart_summary.txt`` and print a short version of it (#658). + + A fit that runs several searches from different starting points reports the best + of them. On its own that number says nothing about whether the search can be + trusted, because a run whose starts all agreed and a run whose starts all + disagreed print the same thing. This is the table that separates them: one row per + start, best objective value first. + + Nothing is written for a fit with fewer than two starts, which includes every fit + type that does not run a multi-start at all and every refine (a refine polishes + the one point it was handed). A refine that does run several starts writes to the + ``_refine`` name, as the other artifacts of a second phase do, rather than + overwriting the searching fit's own table. + + Every failure is logged and swallowed. The fit has finished by this point, and a + report must never be the reason a finished run dies. + """ + try: + records = list(self.multistart_records() or ()) + except Exception: + logger.exception('Failed to collect the per-start results for the multi-start summary') + return + if len(records) < 2: + return + name = 'multistart_summary_refine.txt' if self.refine else 'multistart_summary.txt' + path = str(Path(self.res_dir) / name) + lines = multistart_report.summary_lines( + records, job_type=self.config.config.get('fit_type')) + try: + with open(path, 'w') as f: + f.write('\n'.join(lines) + '\n') + except Exception: + logger.exception('Failed to write %s' % name) + return + logger.info('Wrote the multi-start summary %s' % path) + for line in multistart_report.console_lines(records, path): + print1(line) + def _emit_profiled_noise(self): """Write ``Results/profiled_noise.txt`` for a fit that profiles a noise scale out of the search (``noise_profiling = 1``, ADR-0108, #562). diff --git a/pybnf/algorithms/multistart_report.py b/pybnf/algorithms/multistart_report.py new file mode 100644 index 00000000..77bbc071 --- /dev/null +++ b/pybnf/algorithms/multistart_report.py @@ -0,0 +1,185 @@ +"""The per-start table a multi-start fit writes when it finishes (#658). + +Several fit types run more than one search, each from a different starting point, and +report the best result. That single number cannot be checked. Twenty starts that all +reach about the same objective value say the fit has very likely found the best answer +available, and that running more starts would not help. Twenty starts that all land +somewhere different say 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. +Before this, both cases printed one number and nothing else, so a user had no basis for +telling them apart. + +Every multi-start fit now writes ``Results/multistart_summary.txt``: one row per start, +sorted by final objective value from best to worst, plus a short version on the screen. +Sorting the values that way and looking at the shape of the resulting curve is the usual +way to answer the question in the parameter fitting literature, where the picture is +called a waterfall plot. + +Three families of fit type run more than one start, and each keeps the numbers somewhere +different: the concurrent local optimizers (``trf`` / ``lbfgs`` / ``gntr`` / ``powell`` / +``sim``, and the polish phase of ``profile_likelihood``) hold them on their per-start +runner objects, multiple shooting (``ms``) on its per-start homotopy results, and the +metaheuristics (``de`` / ``ade`` / ``ss`` / ``pso``) in the mixin that drives their +restarts. So the row type and the formatting live here, apart from all three. Each family +hands rows up through :meth:`pybnf.algorithms.base.Algorithm.multistart_records`, and +:meth:`pybnf.algorithms.base.Algorithm._emit_multistart_summary` writes whatever it is +given. This module holds no state and touches no files, so it is testable on its own. +""" + +import math +from collections import namedtuple +from statistics import median + +#: One start's outcome, the row type of the summary. +#: +#: ``start`` numbers the starts from 1, matching how the run log names them. +#: ``objective`` is the best objective value that start reached; ``inf`` where the start +#: produced no usable fit (its start point failed to simulate, say) and ``None`` where it +#: never ran at all. ``iterations`` counts the steps the start took, in whatever unit the +#: method counts, and ``evaluations`` counts the simulations charged to it; either is +#: ``None`` for a method that does not keep that count. ``stop_reason`` is the method's +#: own explanation, or ``None`` for a start the fit ended before. +StartRecord = namedtuple('StartRecord', + 'start objective iterations evaluations stop_reason', + defaults=(None, None, None, None)) + +#: What the reason column says for a start that had not stopped on its own when the fit +#: ended. A wall-time budget is the usual cause. +UNFINISHED = 'did not finish before the fit ended' + +#: The reason for a start that never began, which happens when a fit whose starts run one +#: after another ends before it reaches them. +NOT_STARTED = 'did not start before the fit ended' + +#: How close to the best objective a start has to land to count as agreeing with it, as a +#: fraction of the best value (so 0.001 is a tenth of a percent). Relative rather than +#: absolute because objective values here span many orders of magnitude, with a floor of +#: 1.0 on the scale so a best value at or near zero compares sensibly. +PLATEAU_REL_TOL = 1e-3 + + +def _has_value(record): + """Whether this start reported a real, usable objective value.""" + obj = record.objective + return obj is not None and not math.isnan(obj) + + +def sorted_records(records): + """The records ordered by objective value, best first. + + A start with no value at all sorts last, and ties keep start order, so the table is + the same on every run of the same fit. + """ + return sorted(records, + key=lambda r: (0, float(r.objective), r.start) if _has_value(r) + else (1, 0.0, r.start)) + + +def _finite_values(records): + return [float(r.objective) for r in records + if _has_value(r) and math.isfinite(r.objective)] + + +def plateau_count(records, rel_tol=PLATEAU_REL_TOL): + """How many starts finished within ``rel_tol`` of the best objective value. + + This is the number that answers the question the table exists for. A count close to + the number of starts means the starts agree on an answer; a count of one means the + best result stands alone and the search has not shown that it found anything. + Returns 0 when no start produced a finite value. + """ + values = _finite_values(records) + if not values: + return 0 + best = min(values) + tolerance = rel_tol * max(abs(best), 1.0) + return sum(1 for v in values if v - best <= tolerance) + + +def _objective_text(value): + if value is None or (isinstance(value, float) and math.isnan(value)): + return 'none' + return '%.10g' % value + + +def _count_text(value): + return 'n/a' if value is None else '%d' % value + + +def _reason_text(reason): + if not reason: + return UNFINISHED + # The table is tab separated, so a reason carrying a tab or a newline of its own would + # break the row it sits in. + return ' '.join(str(reason).split()) + + +def _percent_text(rel_tol): + return '%g%%' % (rel_tol * 100.0) + + +def summary_lines(records, job_type=None, rel_tol=PLATEAU_REL_TOL): + """The whole text of ``Results/multistart_summary.txt``, as a list of lines. + + ``job_type`` names the fit type in the header when it is known, so the file can be + read without the configuration file beside it. + """ + rows = sorted_records(records) + values = _finite_values(rows) + lines = [ + '# One row per start of this fit, sorted by final objective value, best first.', + '# This fit ran several searches from different starting points and reported the', + '# best of them. This table says how the others did, which is what tells you', + '# whether to believe the reported best.', + '# If many starts sit at about the same low objective value, the fit has very', + '# likely found the best answer available and more starts would not help. If the', + '# values keep climbing with no group of similar ones at the top of the table,', + '# the reported best is only the least bad of these starts, so run more starts or', + '# take another look at the model and the parameter bounds. Plotting this column', + '# in this order is often called a waterfall plot.', + '#', + '# objective: the best objective value that start reached. "inf" means the start', + '# produced no usable fit, and "none" means it never ran.', + '# iterations: steps that start took, in whatever unit the method counts.', + '# evaluations: simulations charged to that start.', + '# Either count is "n/a" for a method that does not keep it.', + '#', + ] + if job_type: + lines.append('job_type\t%s' % job_type) + lines.append('starts\t%d' % len(rows)) + if values: + lines.append('best_objective\t%.10g' % min(values)) + lines.append('worst_objective\t%.10g' % max(values)) + lines.append('median_objective\t%.10g' % median(values)) + lines.append('# starts_near_best counts the starts that came within %s of the best ' + 'objective value.' % _percent_text(rel_tol)) + lines.append('starts_near_best\t%d' % plateau_count(rows, rel_tol)) + lines.append('#') + lines.append('# rank\tstart\tobjective\titerations\tevaluations\treason') + for rank, r in enumerate(rows, start=1): + lines.append('%d\t%d\t%s\t%s\t%s\t%s' + % (rank, r.start, _objective_text(r.objective), + _count_text(r.iterations), _count_text(r.evaluations), + _reason_text(r.stop_reason))) + return lines + + +def console_lines(records, path, rel_tol=PLATEAU_REL_TOL): + """The short version printed at the end of the run, as a list of lines.""" + rows = sorted_records(records) + values = _finite_values(rows) + if not values: + return ['Multi-start summary: none of the %d starts produced a usable fit. ' + 'Per-start table: %s' % (len(rows), path)] + lines = ['Multi-start summary: %d starts. Best objective %.6g, median %.6g, ' + 'worst %.6g.' % (len(rows), min(values), median(values), max(values))] + unusable = len(rows) - len(values) + if unusable: + lines.append(' %d of the %d produced no usable fit.' % (unusable, len(rows))) + lines.append(' %d of the %d came within %s of the best. The more starts that reach ' + 'the same low value, the more likely the fit has found the best answer ' + 'available.' + % (plateau_count(rows, rel_tol), len(rows), _percent_text(rel_tol))) + lines.append(' Per-start table, best first: %s' % path) + return lines diff --git a/pybnf/algorithms/optimizers/concurrent_multistart.py b/pybnf/algorithms/optimizers/concurrent_multistart.py index fc7ad7f9..9e648527 100644 --- a/pybnf/algorithms/optimizers/concurrent_multistart.py +++ b/pybnf/algorithms/optimizers/concurrent_multistart.py @@ -33,6 +33,8 @@ count -- all plain list/dict/int, so the optimizer pickles for backup/resume (ADR-0007). * **The run loop** (:meth:`start_run` / :meth:`got_result` / :meth:`_route`): seeding, routing, per-iteration reporting, and the ``STOP``-on-last-start coordination. +* **The end-of-run per-start table** (:meth:`multistart_records`, #658): what every start + reached, so a reader can tell a fit whose starts agreed from one whose starts did not. * **Resume plumbing** (:meth:`add_iterations`) and the ``__init__`` / :meth:`reset` skeleton. @@ -76,6 +78,7 @@ import logging from .local_base import StartPointOptimizer +from ..multistart_report import StartRecord from ...printing import print1, print2 logger = logging.getLogger('pybnf.algorithms') @@ -227,6 +230,7 @@ def _init_orchestration(self): self.pending = {} # dispatched pset name -> owning runner index (routing map) self.probe_counter = 0 # global submission counter -> unique pset names self.active = 0 # starts not yet terminated + self.start_evals = {} # runner index -> completed evaluations charged to it (#658) def add_iterations(self, n): """Extend every start's per-start iteration budget by ``n`` (the ``-r`` resume path). @@ -250,6 +254,7 @@ def start_run(self): self.active = len(self.runners) self.probe_counter = 0 self.pending = {} + self.start_evals = {} out = [] for idx, runner in enumerate(self.runners): out.extend(self._seed(idx, runner)) @@ -261,6 +266,10 @@ def got_result(self, res): starts keep going), or ``'STOP'`` only when the last live start finishes.""" idx = self.pending.pop(res.name) runner = self.runners[idx] + # Charge the completed simulation to its start, for the end-of-run summary (#658). + # Counted here rather than at submission so a job abandoned when the fit stopped is + # not billed to a start that never saw its result. + self.start_evals[idx] = self.start_evals.get(idx, 0) + 1 prev_iter = runner.iteration out = self._advance(idx, runner, res) if runner.iteration > prev_iter: @@ -274,6 +283,26 @@ def got_result(self, res): return [] return out + def multistart_records(self): + """One row per start for ``Results/multistart_summary.txt`` (#658). + + Everything the table needs is already on the per-start runners: the objective each + one reached (``fval``), the steps it took (``iteration``), and why it stopped + (``stop_reason``). Read at the end of the run rather than collected as each start + finishes, so a start that was still going when the fit stopped -- a wall-time + budget is the usual reason -- appears in the table too, with whatever it had + reached and no stop reason. Leaving it out would make the table look like a + complete set of starts when it was not. + + Empty before :meth:`start_run` builds the runners. + """ + return [StartRecord(start=i + 1, + objective=runner.fval, + iterations=runner.iteration, + evaluations=self.start_evals.get(i, 0), + stop_reason=runner.stop_reason) + for i, runner in enumerate(self.runners)] + def _route(self, idx, pset): """Record ``pset`` (already uniquely named by the leaf) as owned by start ``idx`` and return it for submission. The name is the routing key ``got_result`` reads back.""" diff --git a/pybnf/algorithms/optimizers/multiple_shooting.py b/pybnf/algorithms/optimizers/multiple_shooting.py index 4c73a97c..22370e0c 100644 --- a/pybnf/algorithms/optimizers/multiple_shooting.py +++ b/pybnf/algorithms/optimizers/multiple_shooting.py @@ -62,6 +62,7 @@ from pydantic import Field from .gradient_base import GradientOptimizer +from ..multistart_report import NOT_STARTED, StartRecord from ...config_schema import PyBNFConfigModel from ...printing import PybnfError, print0, print1, print2 from ...registry import register_fit_type @@ -253,6 +254,18 @@ def __init__(self, config, refine=False): feasibility_tol=config.config['ms_feasibility_tol']) self.homotopies = [] + def reset(self, bootstrap=None): + """Clear the per-start ladder results along with the rest of the run's state. + + A bootstrap replicate reuses the algorithm object, and everything that reports on + the ladder reads this list: the reported best start (:meth:`_best_homotopy`, behind + ``continuity_defects.txt`` and the best stage trace) and the per-start summary. Left + uncleared, a replicate could report a start belonging to the replicate before it, + fitted to different resampled data. + """ + super().reset(bootstrap) + self.homotopies = [] + def _start_banner(self): return ('Running multiple shooting: coarsening ladder from %i segment(s) placed by ' '%s, up to %i outer iteration(s) per rung, from %i start point(s)' @@ -361,6 +374,30 @@ def _run_starts(self, specs, rungs): 'ordinary single-shoot path and are not comparable with an ' 'ordinary fit\'s.') + def multistart_records(self): + """One row per start for ``Results/multistart_summary.txt`` (#658). + + This fit type drives its own search rather than the per-start step machines the + base's summary reads, so the numbers come from the homotopy result each start + produced: its best certified objective, the outer iterations it took over the whole + ladder, and why the ladder stopped. A start the run never reached -- the loop + breaks when the wall-time budget goes -- is listed as one that never ran, so the + table does not read as a complete set of starts when it is not. + """ + rows = [] + for i in range(max(len(self.start_psets), len(self.homotopies))): + if i >= len(self.homotopies): + rows.append(StartRecord(start=i + 1, stop_reason=NOT_STARTED)) + continue + result = self.homotopies[i] + rows.append(StartRecord( + start=i + 1, + objective=result.best_score, + iterations=sum(len(stage.outer.iterates) for stage in result.stages), + evaluations=result.n_evaluations, + stop_reason=result.stop_reason)) + return rows + def _record_iterate(self, record): """Enter one certified outer iterate in the ordinary trajectory. diff --git a/pybnf/algorithms/optimizers/multistart.py b/pybnf/algorithms/optimizers/multistart.py index 7c750be7..396a2661 100644 --- a/pybnf/algorithms/optimizers/multistart.py +++ b/pybnf/algorithms/optimizers/multistart.py @@ -56,10 +56,22 @@ set empties. For a synchronized method the set is already empty, so the next start begins immediately. -All added state is plain ``int`` / ``set`` (the start index, the in-flight names, the -draining flag), so the optimizer pickles for backup/resume exactly as before (ADR-0007). +All added state is plain ``int`` / ``set`` / ``list`` (the start index, the in-flight +names, the draining flag, the per-start tallies), so the optimizer pickles for +backup/resume exactly as before (ADR-0007). + +Reporting the starts (#658) +--------------------------- +Keeping the global best is not the whole job: a user also needs to know whether the +starts agreed, since a run whose starts all reached the same value and a run whose starts +all landed somewhere different otherwise print the same single number. The mixin is the +only place that knows where one start ends and the next begins, so it tallies each start's +best objective value and evaluation count as the results come back and hands them to +``Results/multistart_summary.txt`` through :meth:`~MultiStartOptimizer.multistart_records` +(see :mod:`pybnf.algorithms.multistart_report`). """ +from ..multistart_report import NOT_STARTED, StartRecord from ...config_schema import PyBNFConfigModel from ...printing import print2 @@ -67,6 +79,11 @@ logger = logging.getLogger('pybnf.algorithms') +#: Why a start of a metaheuristic ended, for the end-of-run per-start table (#658). These +#: methods return a bare stop signal without saying which of their two stopping conditions +#: fired, so the row states both rather than guessing at one. +_SEARCH_ENDED = 'the search ended on its own (converged or reached max_iterations)' + class MultiStartConfig(PyBNFConfigModel): """The shared ``n_starts`` config field (#498), mixed into the schema of every @@ -100,6 +117,7 @@ def __init__(self, *args, **kwargs): self._start_index = 0 self._inflight = set() self._draining = False + self._start_stats = [] self.n_starts = self._resolve_n_starts() def _resolve_n_starts(self): @@ -116,6 +134,7 @@ def start_run(self): self._start_index = 0 self._inflight = set() self._draining = False + self._start_stats = [] if self.n_starts > 1: print2('Multi-start: up to %i independent starts, keeping the global best' % self.n_starts) @@ -126,6 +145,7 @@ def got_result(self, res): # on emission, and add_to_trajectory already recorded it); untrack it before any # stripping. self._inflight.discard(res.pset.name) + self._record_for_summary(res) if self._draining: # A straggler from the just-finished start: its score is already in the # trajectory (recorded before got_result), and the inner search has been @@ -136,6 +156,7 @@ def got_result(self, res): self._strip_prefix(res) response = self._search_got_result(res) if response == 'STOP': + self._current_start_stats()['stop_reason'] = _SEARCH_ENDED if self._start_index + 1 < self.n_starts: logger.info('Multi-start: start %d/%d finished; %d job(s) to drain ' 'before the next start', self._start_index + 1, self.n_starts, @@ -159,6 +180,62 @@ def reset(self, bootstrap=None): self._start_index = 0 self._inflight = set() self._draining = False + self._start_stats = [] + + # --- the end-of-run per-start table (#658) ------------------------------- # + def _current_start_stats(self): + """The running record of how the start now in progress is doing, created on + first use. A plain list of dicts, so it rides the backup pickle like the rest of + the mixin's state.""" + while len(self._start_stats) <= self._start_index: + self._start_stats.append({'objective': None, 'evaluations': 0, + 'stop_reason': None}) + return self._start_stats[self._start_index] + + def _record_for_summary(self, res): + """Charge one completed evaluation to the start it belongs to, and keep the best + objective value that start has reached (#658). + + A metaheuristic keeps no single "final objective" of its own -- its answer is the + best member it ever produced -- so the mixin takes the best objective value seen + while that start was running. Every result passes through here, including the + stragglers of a start that has already finished, which belong to that start and + not to the next one. + """ + stats = self._current_start_stats() + stats['evaluations'] += 1 + try: + # Defensive: this runs on every completed result, in the run loop, purely to + # fill in a report. A result that somehow carries no usable objective costs + # the fit a row of a table, not the fit. + score = float(getattr(res, 'score', None)) + except (TypeError, ValueError): + return + if stats['objective'] is None or score < stats['objective']: + stats['objective'] = score + + def multistart_records(self): + """One row per start for ``Results/multistart_summary.txt`` (#658). + + These starts run one after another, so a fit that stops early (a wall-time budget, + say) may never reach some of them. Those are listed too, as starts that never ran, + because a table that quietly showed six rows for a twenty-start fit would read as + a twenty-start fit that agreed with itself. + + The iteration count is left out. Each of these methods counts its own progress + differently -- generations per island, unproductive iterations, and so on -- and + there is no shared number to put in the column. + """ + rows = [] + for i in range(max(self.n_starts, len(self._start_stats))): + if i < len(self._start_stats): + stats = self._start_stats[i] + rows.append(StartRecord(start=i + 1, objective=stats['objective'], + iterations=None, evaluations=stats['evaluations'], + stop_reason=stats['stop_reason'])) + else: + rows.append(StartRecord(start=i + 1, stop_reason=NOT_STARTED)) + return rows # --- name-space boundary ------------------------------------------------ # def _prefix(self): diff --git a/tests/test_multistart_summary.py b/tests/test_multistart_summary.py new file mode 100644 index 00000000..67353697 --- /dev/null +++ b/tests/test_multistart_summary.py @@ -0,0 +1,314 @@ +"""The per-start summary a multi-start fit writes at the end of a run (#658). + +Several fit types run more than one search from different starting points and report the +best of them. That one number cannot be checked: a run whose twenty starts all reached +the same objective value and a run whose twenty starts all landed somewhere different +used to print exactly the same thing, so nobody could tell whether the answer was worth +believing. The fix is ``Results/multistart_summary.txt``, one row per start sorted by +final objective value, plus a short version on the screen. + +Three families of fit type produce those rows, and each keeps the numbers somewhere +different, so each is checked here: + +* the concurrent local optimizers (``powell`` / ``sim``, and the gradient methods + ``trf`` / ``lbfgs`` / ``gntr``, which share the same base) read them off their per-start + runner objects; +* the metaheuristics (``de`` / ``ade`` / ``ss`` / ``pso``) tally them in the mixin that + drives their restarts, because those methods keep no per-start final value of their own; +* multiple shooting (``ms``) reads them off its per-start homotopy results. Its rows are + checked white-box here rather than end to end, since that fit type needs a + sensitivity-capable simulation backend (``tests/test_shooting_sbml.py`` owns those). + +The formatting itself is pure and file-free (``pybnf.algorithms.multistart_report``), so +the first group of tests exercises it directly. +""" +import numpy as np +import pytest + +from pybnf.algorithms import multistart_report as R +from pybnf.algorithms.multistart_report import StartRecord + +from . import integration_harness as H +from .context import algorithms + + +@pytest.fixture(autouse=True) +def _fakes(monkeypatch): + H.install(monkeypatch) + + +# --------------------------------------------------------------------------- # +# The table itself +# --------------------------------------------------------------------------- # +def _rows(lines): + """The data rows of a written summary, split into columns.""" + return [line.split('\t') for line in lines if line and not line.startswith('#') + and line.split('\t')[0].isdigit()] + + +def test_rows_are_sorted_best_objective_first(): + """The whole point of the table is the shape of the sorted objective column, so the + rows come out best first however the starts happened to be numbered.""" + records = [StartRecord(1, 9.0), StartRecord(2, 1.0), StartRecord(3, 5.0)] + assert [r.start for r in R.sorted_records(records)] == [2, 3, 1] + rows = _rows(R.summary_lines(records)) + assert [row[0] for row in rows] == ['1', '2', '3'] # rank column + assert [row[1] for row in rows] == ['2', '3', '1'] # start numbers + assert [row[2] for row in rows] == ['1', '5', '9'] # objective values + + +def test_a_start_with_no_objective_sorts_last(): + """A start that never ran has nothing to compare, so it goes to the bottom rather than + sorting as if it had scored zero.""" + records = [StartRecord(1, None, stop_reason=R.NOT_STARTED), StartRecord(2, 4.0)] + rows = _rows(R.summary_lines(records)) + assert [row[1] for row in rows] == ['2', '1'] + assert rows[1][2] == 'none' + assert rows[1][5] == R.NOT_STARTED + + +def test_a_start_that_never_stopped_says_so(): + """A start still running when the fit ended (a wall-time budget, usually) is listed + with what it had reached. Dropping it would make a cut-short fit look like a complete + set of starts that agreed with each other.""" + rows = _rows(R.summary_lines([StartRecord(1, 2.0, stop_reason='converged'), + StartRecord(2, 3.0, stop_reason=None)])) + assert rows[0][5] == 'converged' + assert rows[1][5] == R.UNFINISHED + + +def test_missing_counts_are_reported_as_such(): + """A method that keeps no iteration count says so rather than printing a zero, which + would read as a start that took no steps.""" + rows = _rows(R.summary_lines([StartRecord(1, 2.0, iterations=None, evaluations=40), + StartRecord(2, 3.0, iterations=7, evaluations=12)])) + assert rows[0][3] == 'n/a' and rows[0][4] == '40' + assert rows[1][3] == '7' and rows[1][4] == '12' + + +def test_the_plateau_count_separates_agreement_from_disagreement(): + """The number a reader acts on. Starts bunched at the same low value count together + (run more starts, probably pointless); starts spread out do not (run more starts).""" + agreed = [StartRecord(i + 1, 10.0 + i * 1e-6) for i in range(5)] + spread = [StartRecord(i + 1, 10.0 * (i + 1)) for i in range(5)] + assert R.plateau_count(agreed) == 5 + assert R.plateau_count(spread) == 1 + + +def test_the_plateau_count_works_on_a_negative_objective(): + """Log-likelihood objectives are negative, so "within a fraction of the best" has to be + measured against the size of the best value, not against its sign.""" + records = [StartRecord(1, -200.0), StartRecord(2, -199.99), StartRecord(3, -20.0)] + assert R.plateau_count(records) == 2 + + +def test_starts_that_all_failed_are_reported_rather_than_hidden(): + """Every start failing to simulate is the single most useful thing this table can say, + so it survives to both the file and the screen.""" + records = [StartRecord(i + 1, float('inf'), stop_reason='start point failed to simulate') + for i in range(3)] + rows = _rows(R.summary_lines(records)) + assert [row[2] for row in rows] == ['inf', 'inf', 'inf'] + assert R.plateau_count(records) == 0 + console = ' '.join(R.console_lines(records, '/tmp/x.txt')) + assert 'none of the 3 starts produced a usable fit' in console + + +def test_a_reason_carrying_a_tab_cannot_break_its_row(): + """The file is tab separated and the reasons come from the methods themselves, so the + text is flattened before it goes in a cell.""" + rows = _rows(R.summary_lines([StartRecord(1, 1.0, stop_reason='a\treason\nwith breaks'), + StartRecord(2, 2.0)])) + assert len(rows[0]) == 6 + assert rows[0][5] == 'a reason with breaks' + + +# --------------------------------------------------------------------------- # +# The concurrent local optimizers (powell / sim, and the gradient methods) +# --------------------------------------------------------------------------- # +_LOCAL = {'powell': algorithms.PowellAlgorithm, 'sim': algorithms.SimplexAlgorithm} + +# A shallow local mode at the box center and a deeper one off center, so the starts +# genuinely disagree and the table has something to report. +_MODES = [(0.5, [0.0, 0.0], [1.0, 1.0]), (0.5, [6.0, 6.0], [4.0, 4.0])] + + +def _local_config(tmp_path, fit_type, n_starts, **overrides): + tgt, exp = H.write_target(tmp_path, H.multimodal_spec(_MODES)) + base = dict(n_params=2, var_type='uniform_var', bounds=(-10.0, 10.0), + population_size=8, max_iterations=40, random_seed=1234, n_starts=n_starts) + base.update(overrides) + return H.make_config(tmp_path, fit_type, tgt, exp, **base) + + +def _read_summary(alg, name='multistart_summary.txt'): + with open('%s/%s' % (alg.res_dir, name)) as f: + return f.read() + + +@pytest.mark.parametrize('fit_type', list(_LOCAL)) +def test_a_local_multistart_fit_writes_one_row_per_start(tmp_path, fit_type): + """The deliverable, end to end through the real run loop: every start of a concurrent + multi-start fit gets a row, with the objective it reached, the steps it took, the + simulations it cost, and why it stopped.""" + alg = _LOCAL[fit_type](_local_config(tmp_path, fit_type, n_starts=4)) + H.drive(alg) + + text = _read_summary(alg) + rows = _rows(text.splitlines()) + assert len(rows) == 4 + assert sorted(int(row[1]) for row in rows) == [1, 2, 3, 4] + assert 'starts\t4' in text + objectives = [float(row[2]) for row in rows] + assert objectives == sorted(objectives) # best first + # Every start reported a real objective, ran at least one step, cost at least one + # simulation, and said why it stopped. + assert all(np.isfinite(v) for v in objectives) + assert all(int(row[3]) >= 1 and int(row[4]) >= 1 for row in rows) + assert all(row[5] and row[5] != R.UNFINISHED for row in rows) + # The run's reported best fit is the best over all the starts, so no start can beat it. + assert alg.trajectory.best_score() <= min(objectives) + 1e-9 + + +@pytest.mark.parametrize('fit_type', list(_LOCAL)) +def test_a_single_start_fit_writes_no_summary(tmp_path, fit_type): + """One start has nothing to be compared against, so the table would only restate the + reported best fit. Not written.""" + alg = _LOCAL[fit_type](_local_config(tmp_path, fit_type, n_starts=1)) + H.drive(alg) + with pytest.raises(FileNotFoundError): + _read_summary(alg) + + +def test_the_evaluation_count_adds_up_to_the_simulations_the_fit_ran(tmp_path): + """The per-start evaluation counts are a breakdown of the run's own work, so they + account for every completed simulation and none twice.""" + alg = _LOCAL['powell'](_local_config(tmp_path, 'powell', n_starts=3)) + H.drive(alg) + records = alg.multistart_records() + assert sum(r.evaluations for r in records) == alg.completed_simulations + + +def test_an_unfinished_start_still_appears(tmp_path): + """White-box, because a wall-time budget expiring mid-start is awkward to arrange: + a runner that never set a stop reason is still a row, so a fit cut short does not + silently report a smaller, better-agreeing set of starts than it ran.""" + alg = _LOCAL['powell'](_local_config(tmp_path, 'powell', n_starts=3)) + alg.start_run() + alg.runners[0].fval, alg.runners[0].iteration = 1.0, 5 + alg.runners[0].stop_reason = 'converged' + alg.runners[1].fval, alg.runners[1].iteration = 2.0, 3 # still running: no stop reason + alg.runners[2].fval, alg.runners[2].iteration = 3.0, 4 + alg.runners[2].stop_reason = 'converged' + + records = alg.multistart_records() + assert [r.start for r in records] == [1, 2, 3] + assert records[1].stop_reason is None + assert R.UNFINISHED in '\n'.join(R.summary_lines(records)) + + +# --------------------------------------------------------------------------- # +# The metaheuristics, whose starts run one after another +# --------------------------------------------------------------------------- # +_META = {'de': algorithms.DifferentialEvolution, + 'ss': algorithms.ScatterSearch, + 'pso': algorithms.ParticleSwarm, + 'ade': algorithms.AsynchronousDifferentialEvolution} + +_META_BUDGET = {'de': dict(population_size=8, max_iterations=15), + 'ss': dict(population_size=5, max_iterations=6), + 'pso': dict(population_size=8, max_iterations=15), + 'ade': dict(population_size=8, max_iterations=15)} + + +def _meta_config(tmp_path, fit_type, n_starts): + tgt, exp = H.write_target(tmp_path, H.multimodal_spec(_MODES)) + base = dict(n_params=2, var_type='uniform_var', bounds=(-10.0, 10.0), + random_seed=1234, n_starts=n_starts) + base.update(_META_BUDGET[fit_type]) + return H.make_config(tmp_path, fit_type, tgt, exp, **base) + + +@pytest.mark.parametrize('fit_type', list(_META)) +def test_a_metaheuristic_multistart_fit_writes_one_row_per_start(tmp_path, fit_type): + """These methods run their starts one after another and keep no final value of their + own, so the mixin takes the best objective each start reached. Every start still gets + a row, and the iteration column says n/a rather than inventing a shared unit.""" + alg = _META[fit_type](_meta_config(tmp_path, fit_type, n_starts=3)) + H.drive(alg) + + rows = _rows(_read_summary(alg).splitlines()) + assert len(rows) == 3 + assert sorted(int(row[1]) for row in rows) == [1, 2, 3] + objectives = [float(row[2]) for row in rows] + assert objectives == sorted(objectives) and all(np.isfinite(v) for v in objectives) + assert all(row[3] == 'n/a' for row in rows) # no shared iteration count + assert all(int(row[4]) >= 1 for row in rows) # but every start cost work + assert alg.trajectory.best_score() <= min(objectives) + 1e-9 + + +def test_a_single_start_metaheuristic_writes_no_summary(tmp_path): + alg = _META['de'](_meta_config(tmp_path, 'de', n_starts=1)) + H.drive(alg) + with pytest.raises(FileNotFoundError): + _read_summary(alg) + + +def test_a_start_the_fit_never_reached_is_listed_as_such(tmp_path): + """These starts run in sequence, so a fit that stops early never reaches the later + ones. They are listed as starts that never ran, so a six-row table for a twenty-start + fit cannot be misread as twenty starts that agreed.""" + alg = _META['de'](_meta_config(tmp_path, 'de', n_starts=5)) + alg._start_stats = [{'objective': 3.0, 'evaluations': 20, 'stop_reason': 'ended'}, + {'objective': 2.0, 'evaluations': 20, 'stop_reason': 'ended'}] + records = alg.multistart_records() + assert len(records) == 5 + assert [r.stop_reason for r in records[2:]] == [R.NOT_STARTED] * 3 + assert all(r.objective is None for r in records[2:]) + + +# --------------------------------------------------------------------------- # +# Multiple shooting +# --------------------------------------------------------------------------- # +class _FakeOuter: + def __init__(self, n_iterates, n_evaluations): + self.iterates = [None] * n_iterates + self.n_evaluations = n_evaluations + + +class _FakeStage: + def __init__(self, n_iterates, n_evaluations): + self.outer = _FakeOuter(n_iterates, n_evaluations) + + +class _FakeHomotopy: + def __init__(self, score, stages, stop_reason): + self.best_score = score + self.stages = stages + self.stop_reason = stop_reason + + @property + def n_evaluations(self): + return sum(s.outer.n_evaluations for s in self.stages) + + +def test_multiple_shooting_reads_its_rows_off_its_ladder_results(): + """``ms`` drives its own search rather than the per-start step machines the shared base + reads, so its rows come from the homotopy result each start produced. Checked against + the shape of those results, so it needs no simulation backend.""" + from pybnf.algorithms.optimizers.multiple_shooting import MultipleShootingAlgorithm + + alg = MultipleShootingAlgorithm.__new__(MultipleShootingAlgorithm) + alg.start_psets = [object(), object(), object()] + alg.homotopies = [ + _FakeHomotopy(5.0, [_FakeStage(3, 30), _FakeStage(2, 20)], 'converged'), + _FakeHomotopy(1.0, [_FakeStage(4, 40)], 'max_outer'), + ] + + records = alg.multistart_records() + assert [r.start for r in records] == [1, 2, 3] + assert records[0].objective == 5.0 + assert records[0].iterations == 5 and records[0].evaluations == 50 + assert records[1].stop_reason == 'max_outer' + # The third start was never reached: the run stopped before the loop got to it. + assert records[2].objective is None and records[2].stop_reason == R.NOT_STARTED diff --git a/tests/test_shooting_sbml.py b/tests/test_shooting_sbml.py index 2c9c0e0f..d19208fd 100644 --- a/tests/test_shooting_sbml.py +++ b/tests/test_shooting_sbml.py @@ -469,6 +469,18 @@ def test_ms_is_a_registered_refiner_that_starts_from_the_injected_point(tmp_path assert alg.start_psets[0]['S'] == pytest.approx(222.0) +def test_a_reset_clears_the_previous_runs_ladder_results(tmp_path): + """A bootstrap replicate reuses the algorithm object, and everything that reports on + the ladder reads the list of per-start results: the best start behind + ``continuity_defects.txt`` and the best stage trace, and now the per-start summary + (#658). Left uncleared, a replicate could report a start belonging to the replicate + before it, fitted to different resampled data.""" + alg = H.build(_ms_config(tmp_path), 'ms') + alg.homotopies = ['a result from the previous replicate'] + alg.reset(None) + assert alg.homotopies == [] + + def test_a_cmaes_fit_may_name_ms_as_its_refiner(tmp_path): """The config half of the same arm: ``refine_method = ms`` passes validation on a fit that is not itself ``ms``, and ``MSConfig``'s keys come along as a coherent group rather