Report each start's final objective value at the end of a multi-start fit (#658)

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

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 <multistart_summary>`.
+
 .. _alg-ss:
 
 Scatter Search

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 <multistart_summary>`.
+
   Example:
 
     * ``n_starts = 10``

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? <multistart_summary>` 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.

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. 

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).