Recommend the measurements to make next (#574)

CHANGELOG.md

@@ -5,6 +5,37 @@ All notable changes to PyBNF are documented below. This project adheres to
 
 ## [Unreleased]
 
+### Added
+- **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
 - **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

docs/adr/0129-a-planned-measurement-is-a-one-row-dataset-at-the-fitted-point-so-experimental-design-is-a-sum-over-fisher-terms-pybnf-already-assembles-rather-than-new-sensitivity-math.md

@@ -0,0 +1,157 @@
+# A planned measurement is a one-row dataset at the fitted point, so experimental design is a sum over Fisher terms PyBNF already assembles rather than new sensitivity math (issue #574)
+
+## Status
+
+Accepted. Builds directly on ADR-0080 (the expected-Fisher noise block) and the `gntr`
+work behind it, and answers the question a `profile_likelihood` run (#446/#466) leaves
+open.
+
+## The gap
+
+`job_type = profile_likelihood` ends by classifying each parameter. A **practically
+non-identifiable** verdict says the data does not pin the parameter down. The user's next
+question is "then what should I measure?", and PyBNF had no answer.
+
+## Why this was cheap
+
+`assemble_fisher_hessian` already builds the expected Fisher information at a point, for
+the whole objective surface PyBNF supports: estimated noise scales, log scales,
+normalization and trajectory transforms, the Laplace and Student-t families, count
+families, constraint penalties. It was built for the EFIM trust-region optimizer.
+
+That matrix is what optimal experimental design optimizes. And it is a **plain sum over
+scored points** — one small positive semi-definite matrix per measurement. So:
+
+* the information a planned measurement would add is that measurement's own term;
+* the information of a planned experiment is the sum of the terms in it;
+* choosing a design is choosing which terms to add.
+
+None of that needs new sensitivity math. What it needs is the terms kept apart instead of
+summed, which is one new generator, `iter_fisher_points`, walking the same
+`_iter_scored_points` scaffold the existing assemblers walk. Summing everything it yields
+reproduces `assemble_fisher_hessian` exactly, and that shipped path is untouched.
+
+## The decision
+
+**A planned measurement is a one-row dataset at the fitted point, and its information is
+whatever the existing Fisher assembly returns for it.**
+
+The dataset holds the candidate time, the model's own prediction there as the
+pseudo-observation, and every auxiliary column carried over from the nearest real
+measurement of that observable.
+
+The pseudo-observation is the point of the construction. The *expected* Fisher information
+is by definition the information expected from data generated by the fitted model, so
+setting the observation to the prediction is not a convenience, it is the definition. It
+also makes every noise model work without a special case:
+
+| noise scale | what the planned row gives it |
+|---|---|
+| a data column (`chi_sq`'s `_SD`) | the nearest real measurement's value |
+| relative to the observation | `cv` times the prediction, which is right |
+| a fitted parameter | read from the parameter set, unchanged |
+| a function of the prediction | the prediction at the candidate time |
+| the column mean | the mean of the planned column |
+
+An analytically profiled per-series scale (ADR-0066) profiles to a factor of one against
+these values, which is the same self-consistent statement, so the unscaled prediction is
+what is written.
+
+## The candidate space, and its one real limit
+
+A candidate is an observable the experiment **already measures**, at a time the model is
+**already simulated at**.
+
+The first rule keeps the noise model honest: proposing an observable nobody has measured
+means inventing a precision for an assay nobody has run, and the answer would be a
+function of that invented number.
+
+The second rule is what makes the whole thing free — the sensitivities at every simulated
+time were computed when the best fit was scored, so enumerating and scoring hundreds of
+candidates costs no simulation. But it has a consequence that the first implementation
+made obvious and that is worth stating plainly: **for a time course, PyBNF derives the
+simulated grid from the data**, so with nothing else in place a design could only ever
+recommend measuring the same times over again.
+
+Hence `design_grid` and `design_t_end`, which add extra simulated times, out past the last
+measurement when asked. The measured times are always kept in the grid, so the data still
+lands on exact grid points and the scoring is bit-for-bit what it was; the extra rows are
+simulated and ignored by everything except the design. Both default to off. This follows
+the precedent ADR-0112 set for `time_error`, which also needed a data-bearing time course
+simulated on a grid the data did not dictate.
+
+## Criteria, and one non-obvious case
+
+Three, on the information matrix in **sampling space** (so a log-scaled parameter is judged
+on its order of magnitude, which is the scale it is fitted on): `a` (summed parameter
+variance, restricted to `design_target` — the classical c-criterion when one parameter is
+named), `d` (log determinant), `e` (smallest eigenvalue).
+
+`design_target` is refused for `d` and `e` rather than ignored. Both are properties of the
+whole matrix, so a restricted version would mean something other than what the name says.
+
+The non-obvious case is a criterion that **cannot tell candidates apart**. When a target
+parameter has infinite variance, every candidate leaves it infinite; when the log
+determinant is `-inf`, adding one term rarely rescues it; the smallest eigenvalue is zero
+for every candidate until the matrix has full rank. In each case the ranking is flat and
+the selection would pick arbitrarily. So while that holds, the selection maximizes a
+different quantity — how much of a candidate's information falls in the directions nothing
+currently sees — and switches back to the requested criterion the moment it can
+discriminate.
+
+The refinement that matters: this is decided from **the criterion**, not from whether the
+matrix is singular. A design aimed at one parameter is perfectly well posed while some
+*other* combination of parameters stays invisible, and it should go on optimizing what it
+was asked to. Getting this wrong is not a rounding error; it returns an empty design.
+
+## Structural non-identifiability is reported, not optimized around
+
+Adding every candidate at once is the most any design over the space could know. A target
+still left with an infinite variance there cannot be fixed by any experiment in the space,
+and the run says so and stops. This is the design half of the same fact a profile
+likelihood reports as a flat profile. A `profile_likelihood` run that hits it still ends
+normally: the profiles are already written, and an impossible design is a finding.
+
+## Two surfaces, one calculation
+
+* `job_type = design` takes the optimum as given (`initial_value:` on every parameter,
+  scoped to that spelling exactly as profile likelihood scopes it, #583), evaluates it
+  once, and writes the report. It runs no search, which is why it is exempt from having to
+  declare `population_size` and `max_iterations`.
+* `profile_likelihood_design = 1` makes a profile-likelihood run end with the same report,
+  around the optimum it just found and aimed at the parameters it just flagged as
+  practically non-identifiable. Structurally non-identifiable parameters are deliberately
+  left out of that aim: no design fixes them, so aiming at them would only produce the
+  refusal above.
+
+Both are off or absent by default, so no existing run changes.
+
+## The report is in the units of the question
+
+The design is quoted as each parameter's confidence interval now and after — the same
+quantity `profile_likelihood_summary.txt` reports, at the same threshold. That is
+`theta* ± sqrt(threshold · (F^-1)_kk)` in the parameter's own scale: exact for a linear
+model, and local for anything else, which is also all a design computed at one point ever
+claimed to be. It makes the recommendation checkable rather than merely plausible, which
+is the property a prescription needs and a diagnosis does not.
+
+## Consequences
+
+* `assemble_fisher_hessian` and every `gntr` fit are untouched; the new generator is a
+  second walk over the same scaffold.
+* `iter_fisher_points` is public, so the per-point information is available to anything
+  else that wants it.
+* The quantile helpers moved out of `profile_likelihood.py` into `pybnf/quantiles.py`,
+  because a confidence threshold is now read by two features rather than one.
+* `design_grid` changes an experiment's simulated grid when set. A normalization computed
+  over the trajectory (ADR-0053) would see the denser grid, which is a better estimate of
+  the same thing but not the identical number. Off by default.
+
+## Deliberately not in this cut
+
+* **Robust design over an ensemble.** Averaging the criterion over posterior draws, or over
+  the points a profile-likelihood run already computed, needs a simulation per ensemble
+  member and needs its own evidence that it beats the local design. Filed separately.
+* **Proposing a new condition.** A perturbation is a model modification, and PyBNF has no
+  vocabulary for proposing one.
+* **Closed-loop design.** This emits a ranked recommendation. A person runs the assay.

docs/algorithms.rst

@@ -830,7 +830,9 @@ optimizers driven by exact forward parameter sensitivities:
 
 All three converge far faster than the metaheuristics near a good fit, and the same
 sensitivity machinery drives profile-likelihood identifiability analysis
-(``job_type = profile_likelihood``) and :ref:`multiple shooting <multiple_shooting>`
+(``job_type = profile_likelihood``), :ref:`optimal experimental design
+<experimental_design>` (``job_type = design``, which says what to measure next), and
+:ref:`multiple shooting <multiple_shooting>`
 (``job_type = ms``) — which changes the fit's *transcription* rather than its search,
 cutting each time course into segments joined by continuity constraints so that a long
 horizon stops hiding the answer from every optimizer. These methods, the noise families and

docs/config_keys.rst

@@ -277,7 +277,8 @@ Required Keys
   *samplers* (``am`` / ``dream`` / ``p_dream`` / ``pt`` / ``mh``, and the
   gradient-based :ref:`hmc <alg-hmc>` for analytical objectives), the
   :ref:`profile-likelihood <gradient_fitting>` identifiability analysis
-  (``profile_likelihood``), and the model *checker*
+  (``profile_likelihood``), the :ref:`experimental design <experimental_design>` that
+  says what to measure next (``design``), and the model *checker*
   (``check``), not just fitting. The value names the specific
   procedure; the key names the kind of job. Requires :ref:`edition <edition>` ``>= 2``,
   and like the modern objective surface there is **no implicit default** -- the run

docs/gradient_fitting.rst

@@ -424,6 +424,114 @@ subset to profile; default all), ``profile_likelihood_step`` / ``profile_likelih
 concurrent directional walks; ``0`` = all of them).
 
 
+.. _experimental_design:
+
+
+Experimental design (what to measure next)
+-------------------------------------------
+
+A profile-likelihood run tells you that a parameter is *practically non-identifiable*. The next
+question is always "then what experiment should I run?", and ``job_type = design`` answers it.
+
+The answer comes from the same object the ``gntr`` optimizer already builds: the expected Fisher
+information at the best fit. That matrix is a **sum over the measured points** — one small
+contribution per measurement — so the information a *planned* measurement would add is just that
+measurement's own contribution, and the information of a whole planned experiment is the sum of the
+ones in it. Its inverse is the covariance the fit would have, which is why a design can be quoted
+in the units you already read off a profile-likelihood run: the confidence interval each parameter
+would end up with.
+
+**What a design may recommend.** One observable, in one experiment, at one time. The observable has
+to be one that experiment already measures, so its noise model is the one your fit is already
+using rather than an invented precision for an assay nobody has run. The time has to be one the
+model is already simulated at, because that is where the forward sensitivities exist — which also
+means enumerating and scoring the whole candidate space costs no simulation at all. The same
+measurement may be recommended more than once; that means measure it that many times, and it is the
+right answer when the precision of one measurement, rather than the shape of the trajectory, is
+what limits you.
+
+.. important::
+
+   For a time course PyBNF simulates the times your data was measured at, so **by default a design
+   can only recommend repeating a measurement you have already made**. ``design_grid = N`` adds
+   ``N`` extra simulated times for it to choose from, and ``design_t_end`` moves the far end of
+   that window past your last measurement. Set both when you want to be told to measure at a time
+   you never have. Your measured times are always kept in the grid, so the scoring of your data is
+   unchanged; the extra rows are simulated and then ignored by everything except the design.
+
+**Choosing between designs.** ``design_criterion`` says what makes one design better than another:
+
+* ``a`` (the default) — the summed variance of the parameters, or of the ones named by
+  ``design_target``. Naming a single parameter is the classical **c-criterion**, and it is the one
+  that answers a profile-likelihood verdict: "``k_deg`` came back practically non-identifiable"
+  becomes "measure whatever pins down ``k_deg``".
+* ``d`` — the volume of the joint confidence region, through the log determinant of the
+  information. The all-round choice when no one parameter is the problem.
+* ``e`` — the worst-determined direction, through the smallest eigenvalue.
+
+``design_target`` applies only to ``a``. The other two are properties of the whole information
+matrix, so naming targets for them is refused rather than quietly ignored.
+
+The measurements are chosen one at a time, each time adding whichever candidate improves the
+criterion most. That is the standard treatment of a subset-selection problem and it is not
+guaranteed to find the single best set of ``design_points`` measurements, but it always terminates
+and the report shows the criterion after each pick, so a design that has stopped paying off is
+visible rather than implied.
+
+**Running it.** ``job_type = design`` does not fit. Give it the fitted values as an
+``initial_value:`` on every parameter, exactly as a profile-likelihood run takes its optimum; it
+simulates that one point and writes the report::
+
+    edition = 2
+    model: model.bngl
+    experiment: myexp, data: mydata.exp
+    output_dir = output/design
+
+    bngl_backend = bngsim
+    job_type = design
+    objective = chi_sq
+
+    design_points = 5                  # how many measurements to recommend
+    design_criterion = a
+    design_target = k_deg              # the parameter the design is aimed at
+    design_grid = 50                   # 50 extra candidate times to choose from
+    design_t_end = 120                 # ... out to t = 120, past the last measurement
+
+    parameter: k_deg, lower: 1e-4, upper: 1e2, initial_value: 0.017
+    parameter: k_syn, lower: 1e-4, upper: 1e2, initial_value: 3.1
+
+A configuration that supplies no fitted values is refused: a design computed at an unfitted point
+is a recommendation about a model nobody has fitted.
+
+**Getting it from a profile-likelihood run instead.** Set ``profile_likelihood_design = 1`` and the
+identifiability run ends by writing the same report, around the optimum it just found and aimed at
+the parameters it just flagged as practically non-identifiable. It reads the same ``design_*`` keys,
+except that the predicted intervals are quoted at ``profile_likelihood_confidence`` so both halves
+of the output are the same statement. Off by default, so a run that does not ask for it is
+unchanged.
+
+**Reading the report.** ``Results/experimental_design.txt`` has two halves. The first is the
+recommendation: the measurements to make, in the order they were chosen, with a replicate count
+and the criterion after each pick. The second is what they are expected to buy: every parameter's
+confidence interval as it stands now and as it would be once the recommended measurements are in
+hand, with the ratio of the two widths. Those predicted intervals are
+:math:`\theta^\* \pm \sqrt{\Delta\chi^2\,(F^{-1})_{kk}}` in the parameter's own fitted scale
+— the quadratic approximation to the profile, which is exact for a linear model and local for
+anything else, as is the design itself.
+
+**When no design can help.** If measuring every candidate at once still leaves a target parameter
+undetermined, the run says so instead of recommending measurements. That is structural
+non-identifiability: the model and the data it can produce do not distinguish that parameter at
+all, and the fix is a different observable, a fixed parameter, or a reparameterization — not more
+data. A profile-likelihood run reports which parameters are in that position.
+
+**What this cut does not do.** The design is computed at the best fit, so it is only as good as
+that fit; averaging the criterion over an ensemble of plausible parameter values (posterior draws,
+or the points a profile-likelihood run already computed) is a separate step and is not implemented.
+Neither is proposing a new experimental *condition*, which would mean proposing a model
+modification. The design is recorded in ADR-0129.
+
+
 What it computes
 ----------------