diff --git a/CHANGELOG.md b/CHANGELOG.md index 66a1be82..a7a17329 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/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 b/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 new file mode 100644 index 00000000..0250e7c9 --- /dev/null +++ b/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. diff --git a/docs/algorithms.rst b/docs/algorithms.rst index 40c7e9f1..5dd910b1 100644 --- a/docs/algorithms.rst +++ b/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 ` +(``job_type = profile_likelihood``), :ref:`optimal experimental design +` (``job_type = design``, which says what to measure next), and +:ref:`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 diff --git a/docs/config_keys.rst b/docs/config_keys.rst index c7ad1117..5a25ad64 100644 --- a/docs/config_keys.rst +++ b/docs/config_keys.rst @@ -277,7 +277,8 @@ Required Keys *samplers* (``am`` / ``dream`` / ``p_dream`` / ``pt`` / ``mh``, and the gradient-based :ref:`hmc ` for analytical objectives), the :ref:`profile-likelihood ` identifiability analysis - (``profile_likelihood``), and the model *checker* + (``profile_likelihood``), the :ref:`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 ` ``>= 2``, and like the modern objective surface there is **no implicit default** -- the run diff --git a/docs/gradient_fitting.rst b/docs/gradient_fitting.rst index 26a16220..6625084f 100644 --- a/docs/gradient_fitting.rst +++ b/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 ---------------- diff --git a/docs/index.rst b/docs/index.rst index 111fcb18..a4ad3876 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -38,9 +38,11 @@ posterior distribution of the free parameters to quantify uncertainty — Adapti MCMC (``am``, the recommended sampler), DREAM(ZS), Preconditioned DREAM, parallel tempering, and a Hamiltonian Monte Carlo / NUTS reference sampler. -**Analysis methods** round out the suite: **model checking** (``check``) and -**profile-likelihood** analysis for identifiability, plus uncertainty -quantification by **bootstrapping**. +**Analysis methods** round out the suite: **model checking** (``check``), +**profile-likelihood** analysis for identifiability, **optimal experimental +design** (``design``), which recommends the measurements that would most improve +the parameters a fit could not pin down, plus uncertainty quantification by +**bootstrapping**. Objectives and noise -------------------- diff --git a/docs/modules/algorithms.rst b/docs/modules/algorithms.rst index 14a13694..4d4c2661 100644 --- a/docs/modules/algorithms.rst +++ b/docs/modules/algorithms.rst @@ -69,6 +69,9 @@ Gradient-based optimizers .. automodule:: pybnf.algorithms.optimizers.profile_likelihood :members: +.. automodule:: pybnf.algorithms.optimizers.design + :members: + Bayesian samplers ================= diff --git a/docs/modules/design.rst b/docs/modules/design.rst new file mode 100644 index 00000000..1b514264 --- /dev/null +++ b/docs/modules/design.rst @@ -0,0 +1,53 @@ +.. _design_module: + +================================================================== +PyBNF experimental design (:py:mod:`pybnf.design`) +================================================================== + +The :py:mod:`pybnf.design` package works out which measurement to make next. It +reads the expected Fisher information :py:func:`pybnf.gradient.assemble_fisher_hessian` +already builds for the ``gntr`` optimizer, and rests on one fact about it: the +information is a plain sum over the measured points, so the information a +*planned* measurement would add is that measurement's own term +(:py:func:`pybnf.gradient.iter_fisher_points`). + +It has four parts: *candidates*, which enumerates the measurements a design may +choose from and scores each by handing the Fisher assembly a one-row dataset +holding the model's own prediction at that point; *criteria*, which reduces an +information matrix to the single number two designs are compared on; *greedy*, +which chooses the measurements one at a time; and *report*, which writes the +recommendation together with the confidence intervals it is expected to produce. + +The user-facing account -- what a design may recommend, the criteria, and the +grid controls that let it propose a time you have never measured -- is in +:ref:`experimental_design`. + +Configuration +============= + +.. automodule:: pybnf.design.config + :members: + +Candidates +========== + +.. automodule:: pybnf.design.candidates + :members: + +Criteria +======== + +.. automodule:: pybnf.design.criteria + :members: + +Selection +========= + +.. automodule:: pybnf.design.greedy + :members: + +Report +====== + +.. automodule:: pybnf.design.report + :members: diff --git a/docs/modules/gradient.rst b/docs/modules/gradient.rst index 2cf76fd7..666d16a3 100644 --- a/docs/modules/gradient.rst +++ b/docs/modules/gradient.rst @@ -19,7 +19,10 @@ step. The capability gate and the per-layer math (fixed and estimated noise scale, log/lognormal scale, per-observable transforms and normalization, the Laplace and Student-t families with mean centering, and constraint-penalty gradients) -are documented for users in :ref:`gradient_fitting`. +are documented for users in :ref:`gradient_fitting`. The assembly also builds the +expected Fisher information the ``gntr`` optimizer steps with, and yields the same +information one scored point at a time, which is what :py:mod:`pybnf.design` scores a +candidate measurement with. Errors ====== diff --git a/docs/modules/index.rst b/docs/modules/index.rst index 1ab37945..03c0afa0 100644 --- a/docs/modules/index.rst +++ b/docs/modules/index.rst @@ -12,6 +12,7 @@ PyBNF Module References config constraint data + design diagnostics gradient inference_data diff --git a/pybnf/algorithms/__init__.py b/pybnf/algorithms/__init__.py index f1840adf..ab41ff96 100644 --- a/pybnf/algorithms/__init__.py +++ b/pybnf/algorithms/__init__.py @@ -74,6 +74,13 @@ from .optimizers.profile_likelihood import ( ProfileLikelihoodAlgorithm as ProfileLikelihoodAlgorithm, ) +# Optimal experimental design (#574): what to measure next, read off the same expected Fisher +# information gntr already assembles. A GradientOptimizer for its gates, routings and +# sensitivities, but it runs no search -- it evaluates the supplied best fit once and scores +# candidate measurements against it. Importing the leaf runs its @register_fit_type. +from .optimizers.design import ( + ExperimentalDesignAlgorithm as ExperimentalDesignAlgorithm, +) from .samplers.dream import DreamAlgorithm as DreamAlgorithm from .samplers.pdream import PDreamAlgorithm as PDreamAlgorithm from .samplers.basic_mcmc import BasicBayesMCMCAlgorithm as BasicBayesMCMCAlgorithm diff --git a/pybnf/algorithms/optimizers/design.py b/pybnf/algorithms/optimizers/design.py new file mode 100644 index 00000000..24afc857 --- /dev/null +++ b/pybnf/algorithms/optimizers/design.py @@ -0,0 +1,235 @@ +"""Optimal experimental design as a job (``job_type = design``, #574), and the report a +``profile_likelihood`` run can end with. + +A profile-likelihood run answers "can this parameter be determined from the data I have?". This +answers the question that follows it: "then what should I measure next?". Both are read off the +same object, the expected Fisher information at the best fit, which PyBNF already assembles for +the ``gntr`` optimizer. The design work itself -- enumerating candidate measurements, scoring +them, choosing between them -- lives in :mod:`pybnf.design`; this module is the run around it. + +Where the best fit comes from +----------------------------- +``job_type = design`` does not fit. It takes the optimum as given, from an ``initial_value:`` on +every free parameter, exactly as ``job_type = profile_likelihood`` does when one is supplied. That +is the natural way to use it: run a fit, then ask what to measure next. It simulates that one +point (which is what produces the sensitivities the information is built from) and writes the +report. A configuration that does not supply the optimum is refused, rather than silently +designing around the middle of the parameter box, which would be a recommendation about a model +nobody has fitted. + +``job_type = profile_likelihood`` supplies its own optimum, so it can write the same report at the +end of its run as a *finding* of the identifiability analysis -- and it knows which parameters came +back practically non-identifiable, so it aims the design at those without being told +(:class:`DesignReportMixin`). +""" + +import logging +import os + +from .gradient_base import GradientOptimizer +from ...gradient import GradientNotSupported +from ...design import ( + DesignExperiment, + DesignFields, + baseline_information, + candidate_information, + format_design_summary, + require_identifiable, + resolve_targets, + select_design, + write_design_report, +) +from ...printing import PybnfError, print1, print2 +from ...pset import PSet +from ...quantiles import chi2_quantile_1dof +from ...registry import register_fit_type + +logger = logging.getLogger('pybnf.algorithms') + +#: The design report's filename in ``Results/``. +DESIGN_REPORT = 'experimental_design.txt' + + +class DesignReportMixin: + """Compute and write an experimental-design report from one evaluated point. + + Mixed into any :class:`~pybnf.algorithms.optimizers.gradient_base.GradientOptimizer` that has + a best fit in hand: the standalone ``design`` job below, and the profile-likelihood job, which + writes the same report at the end of its own run. The host supplies the settings (it owns the + configuration keys) and one master-scored ``Result`` at the optimum; everything else is the + same for both. + """ + + def _design_experiments(self, res): + """Pair each simulated experiment with its measurements and its sensitivity routing. + + The same intersection :meth:`GradientOptimizer.gradient_at` scores over, carrying the + model and experiment names as well, so a recommendation can say which experiment it is + about.""" + routings = self._routings_at(res.pset) + experiments = [] + for model_name, by_suffix in res.simdata.items(): + model_exp = self.exp_data.get(model_name, {}) + for suffix, sim_data in by_suffix.items(): + if suffix in model_exp: + experiments.append(DesignExperiment( + model=model_name, suffix=suffix, sim_data=sim_data, + exp_data=model_exp[suffix], routing=routings[(model_name, suffix)])) + return experiments + + def design_at(self, res, *, points, criterion, targets, observables): + """The finished design at the point ``res`` was scored at. + + ``targets`` is a list of free-parameter ids the design is aimed at (empty means all of + them). Raises a :class:`~pybnf.printing.PybnfError` when the targets are structurally + non-identifiable, because then no design over these observables and times is the answer, + and when the objective is one whose expected Fisher information the assembly does not + build, because that information is exactly what a design is made of. + """ + if res.simdata is None: + raise PybnfError( + "Experimental design could not simulate the best fit, so there are no " + "sensitivities to build the information matrix from.", + hint="Check that the supplied parameter values integrate; a design is computed " + "at the fitted point, so that point has to simulate.") + free_params = [res.pset.get_param(v.name) for v in self.variables] + target_idx = resolve_targets(self.variables, targets, criterion) + try: + experiments = self._design_experiments(res) + baseline = baseline_information(self.objective, experiments, free_params) + candidates = candidate_information( + self.objective, experiments, free_params, observables=observables or None) + except GradientNotSupported as e: + # A design is made of this objective's expected Fisher information, so an objective + # whose Fisher the assembly does not build has no design either. Say that, rather + # than let the internal refusal out or point at a metaheuristic, which would not + # produce a design at all. + raise PybnfError( + "This fit's objective has no expected Fisher information for a design to be " + "built from: %s" % e, + hint="A design is the same information the Fisher/Gauss-Newton optimizer " + "(job_type = gntr) steps with, so an objective that refuses gntr has no " + "design either.") from e + print2('Scoring %d candidate measurement(s) across %d experiment(s).' + % (len(candidates), len(experiments))) + if not self.config.config.get('design_grid'): + print1('Every candidate is a time this fit already simulates, which for a time ' + 'course is a time you have already measured, so the design can only ' + 'recommend repeat measurements. Set design_grid (and design_t_end to look ' + 'past the last measurement) to let it propose new times.') + require_identifiable(baseline, candidates, + [v.name for v in self.variables], target_idx) + return select_design(baseline, candidates, points, criterion, target_idx, + [v.name for v in self.variables]) + + def write_design(self, result, u_star, threshold, confidence): + """Write the design report to ``Results/`` and print its summary.""" + path = os.path.join(self.res_dir, DESIGN_REPORT) + write_design_report(path, result, self.variables, u_star, threshold, confidence) + logger.info('Wrote the experimental design to %s', path) + for line in format_design_summary(result, self.variables, u_star, threshold): + print1(line) + print1('Wrote the full design to %s' % path) + + +class DesignConfig(DesignFields): + """The ``design`` job's configuration: the shared experimental-design keys and nothing else. + + Everything a design needs is common to the two job types that can produce one, so this adds no + fields of its own -- see :class:`~pybnf.design.config.DesignFields` for what each key means.""" + + +# Family ``analysis``, not ``optimizer``: this run fits nothing, and the one thing the family +# is read for -- which job types a PEtab ``job_type = all`` import emits a config for -- must not +# include it. Such a config would refuse at construction, because a design needs the fitted values +# a freshly imported problem does not have. +@register_fit_type('design', family='analysis', + display_name='Optimal Experimental Design', schema=DesignConfig) +class ExperimentalDesignAlgorithm(DesignReportMixin, GradientOptimizer): + """Recommend the measurements to make next (``job_type = design``, #574). + + A one-evaluation job: simulate the supplied optimum, assemble the expected Fisher information + the existing data carries, score every candidate measurement against it, and choose the best + few. It inherits :class:`~pybnf.algorithms.optimizers.gradient_base.GradientOptimizer` for the + gradient path -- the edition, sensitivity-backend and differentiability gates, the + per-experiment routing, and the forward sensitivities themselves -- and then does no fitting at + all, which is why it overrides both run-loop hooks rather than using the multi-start machinery + underneath.""" + + fit_type = 'design' + _method_label = 'experimental design' + + #: One evaluation, so no setting governs how many jobs run at once (#655). + parallelism_setting = None + + def __init__(self, config, refine=False): + # The shared Algorithm setup reads a population size and an iteration budget. This run + # searches nothing, so neither means anything and the configuration does not require + # them (pybnf.config._NO_SEARCH_RUNS). Fill in what is true of this run instead of + # making the user type numbers that do nothing. + config.config.setdefault('population_size', 1) + config.config.setdefault('max_iterations', 1) + super().__init__(config, refine=refine) + self.design_points = config.config['design_points'] + self.design_criterion = config.config['design_criterion'] + self.design_targets = list(config.config.get('design_target') or []) + self.design_observables = list(config.config.get('design_observables') or []) + self.confidence = config.config['design_confidence'] + self.threshold = chi2_quantile_1dof(self.confidence, 'design_confidence') + self.design_result = None + self._theta_star = self._require_supplied_optimum() + + def expected_parallelism(self): + """One evaluation of the supplied optimum: the design itself is arithmetic on the + information matrix, not more simulation.""" + return 1 + + def _require_supplied_optimum(self): + """The optimum to design around, taken from an ``initial_value:`` on every parameter. + + Scoped to that spelling, exactly as ``profile_likelihood`` scopes it (#583): ``initial_value:`` + is a claim that these are the fitted values, whereas ``start_point`` means "begin the search + here" and there is no search to begin. A design around an unfitted point would be a + recommendation about a model nobody has fitted, so it is refused instead.""" + spelling = getattr(self.config, 'start_point_spelling', None) or {} + declared = {name: value + for name, value in (getattr(self.config, 'start_point', None) or {}).items() + if spelling.get(name) == 'initial_value'} + missing = [v.name for v in self.variables if v.name not in declared] + if missing: + raise PybnfError( + "job_type = design needs the fitted values to design around, but %s %s no " + "initial_value." % (', '.join(missing), + 'has' if len(missing) == 1 else 'have'), + hint="Give every parameter its fitted value, as in 'parameter: k, lower: 0.01, " + "upper: 10, initial_value: 0.3'. Run a fit first, then design around its " + "best fit. To fit and design in one run, use job_type = " + "profile_likelihood with profile_likelihood_design = 1.") + return declared + + def _make_runner(self, u0): + """Never called: this job runs no search, so it builds no step machine. The base declares + the hook, so it is answered rather than left to fail obscurely.""" + raise PybnfError('job_type = design runs no search, so it has no optimizer to build.') + + def _start_banner(self): + return ('Designing the next %d measurement(s) at the supplied best fit (%s)' + % (self.design_points, self.design_criterion)) + + def start_run(self): + self._setup_gradient_path() + print2(self._start_banner()) + self.probe_counter = 0 + self.pending = {} + theta_star = PSet([v.set_value(self._theta_star[v.name], reflect=False) + for v in self.variables]) + theta_star.name = '%s_1' % self.fit_type + return [theta_star] + + def got_result(self, res): + self.design_result = self.design_at( + res, points=self.design_points, criterion=self.design_criterion, + targets=self.design_targets, observables=self.design_observables) + self.write_design(self.design_result, self._u_from_pset(res.pset), + self.threshold, self.confidence) + return 'STOP' diff --git a/pybnf/algorithms/optimizers/profile_likelihood.py b/pybnf/algorithms/optimizers/profile_likelihood.py index ba8b1053..306612c9 100644 --- a/pybnf/algorithms/optimizers/profile_likelihood.py +++ b/pybnf/algorithms/optimizers/profile_likelihood.py @@ -91,7 +91,7 @@ is importable -- the profile plots (``Delta chi2`` panels with the threshold + CI lines, the reference notebook's Cell 9); a missing matplotlib skips only the plots, never the run. scipy stays out of the production loop: the chi-square threshold comes from a -dependency-free probit approximation (:func:`_chi2_quantile_1dof`). +dependency-free probit approximation (:func:`~pybnf.quantiles.chi2_quantile_1dof`). """ import logging @@ -101,13 +101,15 @@ import numpy as np +from .design import DesignReportMixin from .gradient_base import DONE, GradientOptimizer from .lbfgs import _LBFGSRunner from .trf import _TRFRunner -from ...config_schema import PyBNFConfigModel +from ...design import DesignFields from ...gradient import GradientResult from ...printing import PybnfError, print1, print2 from ...pset import PSet +from ...quantiles import chi2_quantile_1dof from ...registry import register_fit_type logger = logging.getLogger('pybnf.algorithms') @@ -147,57 +149,6 @@ def _build_inner_runner(kind, u0, lower, upper, max_iterations, *, grad_tol, ste # --------------------------------------------------------------------------- # # chi-square (1 dof) quantile via a probit approximation (scipy-free, ADR-0007) # --------------------------------------------------------------------------- # -def _norm_ppf(p): - """Standard-normal inverse CDF (probit) via Acklam's rational approximation, refined - by one Halley step against :func:`math.erf`. - - Dependency-free so the production loop never imports scipy (ADR-0007); accurate to - full double precision after the refinement, far more than the chi-square threshold - needs. ``0 < p < 1``.""" - a = (-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, - 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00) - b = (-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, - 6.680131188771972e+01, -1.328068155288572e+01) - c = (-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, - -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00) - d = (7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, - 3.754408661907416e+00) - plow, phigh = 0.02425, 1.0 - 0.02425 - if p < plow: - q = math.sqrt(-2.0 * math.log(p)) - x = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \ - ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0) - elif p <= phigh: - q = p - 0.5 - r = q * q - x = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / \ - (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0) - else: - q = math.sqrt(-2.0 * math.log(1.0 - p)) - x = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \ - ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0) - # One Halley refinement using the exact erf-based CDF. - e = 0.5 * math.erfc(-x / math.sqrt(2.0)) - p - u = e * math.sqrt(2.0 * math.pi) * math.exp(x * x / 2.0) - x = x - u / (1.0 + x * u / 2.0) - return x - - -def _chi2_quantile_1dof(confidence): - """The chi-square (1 dof) quantile at probability ``confidence`` -- the profile - ``Delta chi2`` threshold (Raue et al. 2009). - - A single profiled parameter has 1 degree of freedom, and ``chi2_1 = Z**2`` with - ``Z ~ N(0, 1)``, so ``P(chi2_1 <= x) = 2*Phi(sqrt(x)) - 1`` and the quantile is - ``Phi^{-1}((1 + confidence) / 2)**2`` (e.g. ``0.95 -> 3.8415``).""" - if not (0.0 < confidence < 1.0): - raise PybnfError( - "profile_likelihood_confidence must be strictly between 0 and 1, got %r." - % confidence) - z = _norm_ppf(0.5 * (1.0 + confidence)) - return z * z - - # --------------------------------------------------------------------------- # # CI extraction + identifiability classification # --------------------------------------------------------------------------- # @@ -579,7 +530,7 @@ def _adapt(self, dchi2_increment): # --------------------------------------------------------------------------- # # Config + algorithm # --------------------------------------------------------------------------- # -class ProfileLikelihoodConfig(PyBNFConfigModel): +class ProfileLikelihoodConfig(DesignFields): """Profile-likelihood config fields, co-located with the method (ADR-0006). ``profile_likelihood_params`` selects which free parameters to profile (a list of ids; @@ -599,7 +550,16 @@ class ProfileLikelihoodConfig(PyBNFConfigModel): fit the cap queue and run as slots free, they are never dropped). Like the other gradient methods' cycle budgets, ``profile_likelihood_max_iterations`` (the polish budget) is runtime-guarded -- it defaults to the global ``max_iterations`` when unset -- so it is a - valid key but not a schema field.""" + valid key but not a schema field. + + ``profile_likelihood_design`` turns on the experimental-design report (#574): after the + profiles are written, the run recommends the measurements that would most improve the + parameters it just found hard to determine. It inherits the shared ``design_*`` keys from + :class:`~pybnf.design.config.DesignFields`, except that the predicted intervals are quoted at + ``profile_likelihood_confidence`` so both halves of the run's output agree. Off by default, so + a run that does not ask for it is unchanged.""" + + profile_likelihood_design: int = 0 profile_likelihood_params: Any = None profile_likelihood_confidence: float = 0.95 @@ -618,7 +578,7 @@ class ProfileLikelihoodConfig(PyBNFConfigModel): @register_fit_type('profile_likelihood', family='optimizer', display_name='Profile Likelihood', schema=ProfileLikelihoodConfig) -class ProfileLikelihoodAlgorithm(GradientOptimizer): +class ProfileLikelihoodAlgorithm(DesignReportMixin, GradientOptimizer): """Standalone profile-likelihood driver (``job_type = profile_likelihood``, #446/#466). A two-phase job over the :class:`GradientOptimizer` gradient path: an optional @@ -649,7 +609,7 @@ class ProfileLikelihoodAlgorithm(GradientOptimizer): def __init__(self, config, refine=False): super().__init__(config, refine=refine) self.confidence = config.config['profile_likelihood_confidence'] - self.threshold = _chi2_quantile_1dof(self.confidence) + self.threshold = chi2_quantile_1dof(self.confidence, 'profile_likelihood_confidence') self.pl_step = config.config['profile_likelihood_step'] self.pl_min_step = config.config['profile_likelihood_min_step'] self.pl_max_step = config.config['profile_likelihood_max_step'] @@ -660,6 +620,15 @@ def __init__(self, config, refine=False): self.reopt_max_iterations = config.config['profile_likelihood_reopt_max_iterations'] self.grad_tol = config.config['profile_likelihood_grad_tol'] self.step_tol = config.config['profile_likelihood_step_tol'] + # Experimental design (#574), off unless asked for: the measurements to make next, aimed + # by default at whatever this run finds practically non-identifiable. The predicted + # intervals are quoted at the profile's own confidence level, so the two halves of the + # run's output are the same statement about the same threshold. + self.design_report = bool(config.config.get('profile_likelihood_design')) + self.design_points = config.config['design_points'] + self.design_criterion = config.config['design_criterion'] + self.design_targets = list(config.config.get('design_target') or []) + self.design_observables = list(config.config.get('design_observables') or []) if 'profile_likelihood_max_iterations' in config.config: self.max_iterations = config.config['profile_likelihood_max_iterations'] else: @@ -693,6 +662,7 @@ def _init_profile_state(self): self.polished = None # True once the polish phase runs, False on explicit theta* self._runner_kind = None # 'trf' | 'lbfgs' inner optimizer (set at preflight/center) self.profile_summary = None # the per-parameter CI + classification list, set at finalize + self.design_result = None # the experimental design, set at the end when asked for self._cost_ref = None self._u_star = None self._profile_idxs = _resolve_profile_idxs( @@ -809,6 +779,8 @@ def got_result(self, res): self.phase = 'profile' return self._begin_profiling(self._u_from_pset(self.trajectory.best_fit())) return response + if self.phase == 'design': + return self._design_got(res) return self._profile_got(res) def _select_runner_kind(self, res): @@ -1056,6 +1028,53 @@ def _finalize(self): self._write_profile_summary(summary) self._write_profile_plots(summary) self._print_summary(summary) + if self.design_report: + return self._begin_design() + return 'STOP' + + # --- experimental design (#574) ---------------------------------------- # + def _begin_design(self): + """Re-evaluate ``theta*`` once, so the design has the forward sensitivities there. + + The profiling walk left the scheduler holding results from grid points, not from the + optimum, and the polish's own result is long gone. One more evaluation of a point the run + has already visited is a negligible cost beside the profiles, and it makes the design's + information matrix unambiguously the one at ``theta*``.""" + self.phase = 'design' + print1('Working out which measurements would most improve these parameters.') + _name, pset = self._pl_dispatch(self._u_star) + return [pset] + + def _flagged_parameters(self): + """The parameters this run found *practically* non-identifiable -- what the design aims + at when the user has not said otherwise. + + A practically non-identifiable parameter is precisely the one more data can fix, so it is + the one to design for. A *structurally* non-identifiable parameter is left out: no + measurement of these observables determines it, which is a statement about the model + rather than about the data, and aiming a design at it would only produce the refusal + :func:`~pybnf.design.require_identifiable` already gives. An empty list means nothing was + flagged, and the design aims at every parameter instead.""" + return [s['name'] for s in (self.profile_summary or []) + if s['classification'] == 'practically non-identifiable'] + + def _design_got(self, res): + """Build and write the design from the re-evaluated optimum, then end the run. + + A design that cannot be built is reported and the run still ends normally: the profiles, + the curves and the plots are already written, the analysis they carry is complete, and a + target no experiment could ever determine is itself a finding.""" + targets = self.design_targets + if not targets and self.design_criterion == 'a': + targets = self._flagged_parameters() + try: + self.design_result = self.design_at( + res, points=self.design_points, criterion=self.design_criterion, + targets=targets, observables=self.design_observables) + except PybnfError as e: + print1('No experimental design was written. %s' % e.message) + return 'STOP' + self.write_design(self.design_result, self._u_star, self.threshold, self.confidence) return 'STOP' def _write_profile_curves(self, summary): diff --git a/pybnf/config.py b/pybnf/config.py index cffcbb3e..5e1e4900 100644 --- a/pybnf/config.py +++ b/pybnf/config.py @@ -52,6 +52,14 @@ # at. Naming the budget for one of these is refused, not silently ignored (#529). _NO_WALL_TIME_FIT = frozenset({'hmc'}) +# The run types that search nothing, so asking the user for a population size and an iteration +# budget would be asking for numbers that mean nothing: the model check, and the experimental +# design (#574), which evaluates one supplied best fit and then does arithmetic on the +# information matrix it produced. The model check reads neither key; the design run passes through +# the shared Algorithm setup, which does, so it fills in the values that are true of it (one +# point, evaluated once) itself. +_NO_SEARCH_RUNS = frozenset({'check', 'design'}) + def init_logging(file_prefix, debug=False, log_level_name='info'): @@ -251,7 +259,7 @@ def __init__(self, d=None): _modern_hint = isinstance(_ed, int) and not isinstance(_ed, bool) and _ed >= 2 if not self._user_objfunc and not _modern_hint: print1('Warning: objfunc was not specified. Defaulting to chi_sq.') - if not self._req_user_params() <= d.keys() and d['fit_type'] != 'check': + if not self._req_user_params() <= d.keys() and d['fit_type'] not in _NO_SEARCH_RUNS: unspecified_keys = [] for k in self._req_user_params(): if k not in d.keys(): @@ -1652,7 +1660,8 @@ def _load_experiments(self): # max-time bound here, not a readout time. action = self._steady_state_action(name, method, fields) elif action_type == 'time_course': - action = TimeCourse({'suffix': name, 'method': method}, explicit_points=points) + action = TimeCourse({'suffix': name, 'method': method}, + explicit_points=self._with_design_grid(points)) self._attach_nf_options(action, fields, method) else: # parameter_scan / dose-response (ADR-0046): the data's independent-variable @@ -2060,6 +2069,30 @@ def _time_error_active(self): course must be simulated on a dense grid over the support, not at the reported times.""" return any(isinstance(k, tuple) and k[0] == 'time_error' for k in self.config) + def _with_design_grid(self, points): + """An experiment's measured times, plus the extra times an experimental design may + recommend measuring at (``design_grid`` / ``design_t_end``, #574). + + A design can only recommend a time the model is already simulated at, because that is + where the forward sensitivities it scores candidates with exist. For a time course PyBNF + derives the simulated grid from the data, so without this a design could only ever + recommend measuring the same times over again. ``design_grid`` lays down that many extra + simulated times, spread evenly from the first measurement to ``design_t_end`` (which + defaults to the last measurement, and is set beyond it when the design should be allowed + to look past the data in hand). + + The measured times are always kept, so the data still lands on exact grid points and + nothing about the scoring changes; the extra rows are simulated and then ignored by + everything except the design. Both keys are off by default, so every other run simulates + exactly the grid it always did. + """ + count = int(self.config.get('design_grid') or 0) + if count <= 0 or not points: + return points + t_end = float(self.config.get('design_t_end') or 0) or max(points) + return sorted({float(p) for p in points} + | {float(t) for t in np.linspace(min(points), t_end, count)}) + def _time_error_timecourse(self, name, method, fields): """A uniform dense grid over ``[t_start, t_end]`` for a marginalized (``time_error``) experiment (ADR-0112, #587). diff --git a/pybnf/design/__init__.py b/pybnf/design/__init__.py new file mode 100644 index 00000000..f7c43d00 --- /dev/null +++ b/pybnf/design/__init__.py @@ -0,0 +1,98 @@ +"""Optimal experimental design: which measurement to make next (#574). + +A profile-likelihood run (``job_type = profile_likelihood``) ends by telling you that a +parameter is *practically non-identifiable*. That is a diagnosis with no prescription, and the +next question is always "then what experiment should I run?". This package answers it. + +The whole calculation rests on one fact that PyBNF already computes. The expected Fisher +information ``F`` of a fit is a plain **sum over the measured points** -- one small matrix per +point (:func:`~pybnf.gradient.iter_fisher_points`, the terms +:func:`~pybnf.gradient.assemble_fisher_hessian` adds together for the ``gntr`` optimizer). So the +information a *planned* measurement would add is just that measurement's own term, and the +information of a whole planned experiment is the sum of the terms of the points in it. Nothing +has to be re-derived: the noise families, the log and normalized observables, the per-condition +chain rules and the estimated noise scales all come along, because they are already inside those +terms. + +What a design is here +--------------------- +A candidate measurement is one observable, in one experiment, at one time that experiment's +simulation already passes through (:class:`~pybnf.design.candidates.CandidateMeasurement`). The +observable has to be one the experiment already measures, so its noise model is known rather than +assumed. The time can be any time on the simulated grid, so no new simulation is needed -- the +sensitivities at that time have already been computed. + +Scoring a design means reducing its information matrix to one number +(:mod:`pybnf.design.criteria`): + +* **A** (the default) -- the average variance of the parameters, or of a named subset. Naming one + parameter makes this the classical c-criterion, which 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. +* **E** -- the worst-determined direction, through the smallest eigenvalue. + +Picking the best few measurements out of hundreds of candidates is a subset-selection problem, so +:mod:`pybnf.design.greedy` takes them one at a time, each time adding whichever candidate improves +the criterion most. Choosing the same point twice means measuring it twice, which is a real answer: +it says the precision of that one measurement is what limits you. + +Two honest limits, both deliberate. The design is computed at the best fit, so it is only as good +as that fit; averaging over an ensemble of plausible parameter values is a separate, later step. +And a parameter no measurement in the candidate space can pin down makes the information matrix +singular no matter what is added, which is reported as such rather than papered over. +""" + +from .candidates import ( + CandidateMeasurement, + DesignExperiment, + baseline_information, + candidate_information, + measured_observables, +) +from .config import DesignFields +from .criteria import ( + CRITERIA, + criterion_score, + interval_half_widths, + criterion_value, + is_singular, + lower_is_better, + null_space_gain, + parameter_variances, + unidentified_parameters, +) +from .greedy import ( + DesignResult, + improvement, + require_identifiable, + resolve_targets, + select_design, +) +from .report import format_design_summary, predicted_intervals, write_design_report + +__all__ = [ + 'CandidateMeasurement', + 'DesignExperiment', + 'DesignFields', + 'DesignResult', + 'CRITERIA', + 'baseline_information', + 'candidate_information', + 'criterion_score', + 'criterion_value', + 'format_design_summary', + 'improvement', + 'interval_half_widths', + 'is_singular', + 'lower_is_better', + 'measured_observables', + 'null_space_gain', + 'parameter_variances', + 'predicted_intervals', + 'require_identifiable', + 'resolve_targets', + 'select_design', + 'unidentified_parameters', + 'write_design_report', +] diff --git a/pybnf/design/candidates.py b/pybnf/design/candidates.py new file mode 100644 index 00000000..ae368f79 --- /dev/null +++ b/pybnf/design/candidates.py @@ -0,0 +1,238 @@ +"""Enumerating the measurements a design may choose from, and what each one would tell you (#574). + +A candidate measurement is one observable, in one experiment, at one time. Two rules keep the +candidate space honest and cheap: + +* **The observable has to be one that experiment already measures.** Then its noise model is + known, because it is the one the fit is already using. Proposing an observable that has never + been measured would mean inventing a precision for an assay nobody has run, and the answer would + depend entirely on that invented number. +* **The time has to be one the simulation already passes through.** The sensitivities at every + simulated time were computed when the best fit was scored, so every candidate is free: no model + is re-solved to enumerate or to score the candidate space. + +Scoring one candidate means asking what the information matrix would gain if that point were +measured. PyBNF can already answer that, because the information is a sum over measured points +and :func:`~pybnf.gradient.iter_fisher_points` yields the points one at a time. So a candidate is +scored by handing the machinery a **planned measurement**: a one-row dataset at the candidate +time, whose value is the model's own prediction at the best fit. That is what the expected Fisher +information means -- the information you expect from data generated by the fitted model -- and it +is why every noise model comes along for free. A noise scale read from a data column takes the +value from the nearest real measurement of that same observable, which assumes the planned +measurement is as precise as the ones already made; a scale that is a fitted parameter, a constant, +or a function of the prediction needs no assumption at all. +""" + +from dataclasses import dataclass, field +from typing import Any + +import numpy as np + +from ..data import Data +from ..gradient import assemble_fisher_hessian, iter_fisher_points + + +@dataclass(frozen=True) +class DesignExperiment: + """One scored experiment at the best fit, with everything the design needs to reason about it. + + ``sim_data`` is the simulated trajectory carrying the forward-sensitivity tensor, ``exp_data`` + the measurements already taken, and ``routing`` the free-parameter-to-sensitivity-column map + the gradient path built. ``model`` and ``suffix`` name the experiment, so a recommendation can + say which one it is about; the suffix is also the scoring key the objective resolves a + per-series scale against.""" + + model: str + suffix: str + sim_data: Any + exp_data: Any + routing: Any + + @property + def label(self): + """How this experiment is named in a report.""" + return self.suffix or self.model + + def as_gradient_tuple(self): + """The ``(sim_data, exp_data, routing, data_key)`` shape every gradient assembler takes.""" + return (self.sim_data, self.exp_data, self.routing, self.suffix) + + +@dataclass(frozen=True) +class CandidateMeasurement: + """One measurement a design may recommend: an observable, in an experiment, at a time. + + ``independent_variable`` is that experiment's own name for its first column, so a report can + say ``time = 4.5`` for a time course and use the right word for anything else.""" + + model: str + experiment: str + observable: str + time: float + independent_variable: str = 'time' + + def __str__(self): + return '%s / %s at %s = %.6g' % ( + self.experiment or self.model, self.observable, + self.independent_variable, self.time) + + +@dataclass +class CandidateSet: + """Every candidate measurement, with the information matrix each one would add. + + ``blocks[i]`` is the ``(n_param, n_param)`` matrix candidate ``measurements[i]`` contributes, + in sampling space, ready to be added to a baseline information matrix.""" + + measurements: list = field(default_factory=list) + blocks: list = field(default_factory=list) + + def __len__(self): + return len(self.measurements) + + def total(self, n_param): + """The information of measuring *everything* at once -- the most any design over this + candidate space could ever know. What this still cannot see, no experiment here can.""" + if not self.blocks: + return np.zeros((n_param, n_param)) + return np.sum(self.blocks, axis=0) + + +def independent_variable(exp_data): + """The name of an experiment's independent variable: its first column, which is how every + other consumer of a PyBNF dataset identifies it.""" + return min(exp_data.cols, key=exp_data.cols.get) + + +def _scoreable_columns(objective, experiment): + """The columns of this experiment's data the objective can score -- the intersection of the + data's columns with the simulation's, plus any column a measurement model materializes. This + is the same set the gradient and objective walk, so nothing here can score a point the fit + does not.""" + return set(experiment.sim_data.cols) | set(objective._per_measurement_models) + + +def measured_observables(objective, experiment): + """The observables this experiment actually measures, in a stable order. + + A column that is present but entirely blank is not measured, so it is left out: its noise + model has never been exercised and there is no nearest real measurement to take a data-column + noise scale from.""" + indvar = independent_variable(experiment.exp_data) + scoreable = _scoreable_columns(objective, experiment) + exp_data = experiment.exp_data + observables = [] + for name in sorted(exp_data.cols): + if name == indvar or name not in scoreable: + continue + column = exp_data.data[:, exp_data.cols[name]] + if np.isfinite(column).any(): + observables.append(name) + return observables + + +def _nearest_measured_rows(exp_data, indvar, col_name, times): + """For each candidate time, the row of the real data whose measurement of ``col_name`` is + closest in time. Auxiliary columns -- a noise scale read from the data, a per-measurement + noise parameter -- are copied from that row, which is the assumption that a planned + measurement is made the same way as the nearest real one. + + Rows that did not measure this observable are ignored. A dataset whose times are all + non-finite (a steady-state experiment, whose measurements name the limit rather than a time, + ADR-0086) has no meaningful nearest row, so its first measured row is used for all of them.""" + column = exp_data.data[:, exp_data.cols[col_name]] + measured = np.flatnonzero(np.isfinite(column)) + if measured.size == 0: + return np.zeros(len(times), dtype=int) + stamps = exp_data.data[measured, exp_data.cols[indvar]] + finite = np.isfinite(stamps) + if not finite.any(): + return np.full(len(times), measured[0], dtype=int) + measured, stamps = measured[finite], stamps[finite] + return measured[np.argmin(np.abs(stamps[None, :] - np.asarray(times)[:, None]), axis=1)] + + +def _planned_measurements(objective, experiment, col_name, times): + """A dataset standing for "measure ``col_name`` at each of ``times``", ready to be scored. + + One row per candidate time. The observable's value is the model's own prediction there at the + best fit, which is what makes the resulting information the *expected* Fisher information. + Every column the objective cannot score -- a noise scale read from the data, a per-measurement + noise parameter -- is carried over from the nearest real measurement of this observable. Every + other observable is left out entirely, so scoring this dataset scores exactly the candidate + points and nothing else. + + The prediction is the **unscaled** one, so a series whose scale is profiled out analytically + (ADR-0066) profiles to a factor of one against these values, which is the same self-consistent + "data generated at the best fit" statement the rest of the row makes.""" + exp_data = experiment.exp_data + indvar = independent_variable(exp_data) + scoreable = _scoreable_columns(objective, experiment) + auxiliary = [name for name in exp_data.cols + if name != indvar and name not in scoreable] + headers = [indvar, col_name] + auxiliary + + times = np.asarray(times, dtype=float) + values = np.zeros((times.size, len(headers))) + values[:, 0] = times + rows = _nearest_measured_rows(exp_data, indvar, col_name, times) + for position, name in enumerate(auxiliary, start=2): + values[:, position] = exp_data.data[rows, exp_data.cols[name]] + + planned = Data.from_columns(values, headers, indvar=indvar) + # Second pass: the pseudo-observations. A materialized measurement-model column (ADR-0036/ + # ADR-0045) reads its own per-row parameters off the dataset, so the dataset has to exist + # before its predictions can be computed. + for row in range(times.size): + sim_row = objective._sim_row_for(experiment.sim_data, planned, indvar, row, + show_warnings=False) + values[row, 1] = objective._base_prediction( + experiment.sim_data, sim_row, col_name, planned, row) + planned.data = values # re-publish, so the weights follow the final values + return planned + + +def baseline_information(objective, experiments, free_params): + """The information the **existing** data already carries, at the best fit. + + Every design is judged as an addition to this, because that is the question being asked: given + what has already been measured, what should be measured next.""" + return assemble_fisher_hessian( + objective, [e.as_gradient_tuple() for e in experiments], free_params) + + +def candidate_information(objective, experiments, free_params, observables=None): + """Enumerate every candidate measurement and the information each one would add. + + Walks each experiment's measured observables over every time on that experiment's simulated + grid. ``observables`` optionally restricts the walk to a named set of observables, for a user + who can only run some of the assays. + + Returns a :class:`CandidateSet`. The order is stable -- experiments in the order given, + observables alphabetically, times in simulated order -- so a design run is reproducible and + ties break the same way every time.""" + wanted = None if observables is None else set(observables) + # Seed the pset reads the objective makes (an estimated noise scale, a measurement model's + # parameters) from this point, exactly as the assemblers do, so this is usable on its own and + # not only after a baseline has been assembled. + existing = getattr(objective, '_pset_values', None) or {} + objective._pset_values = {**existing, **{p.name: p.value for p in free_params}} + candidates = CandidateSet() + for experiment in experiments: + indvar = independent_variable(experiment.exp_data) + times = np.asarray(experiment.sim_data[indvar], dtype=float) + for col_name in measured_observables(objective, experiment): + if wanted is not None and col_name not in wanted: + continue + planned = _planned_measurements(objective, experiment, col_name, times) + scored = iter_fisher_points( + objective, + [(experiment.sim_data, planned, experiment.routing, experiment.suffix)], + free_params) + for _index, row, _column, block in scored: + candidates.measurements.append(CandidateMeasurement( + model=experiment.model, experiment=experiment.suffix, + observable=col_name, time=float(times[row]), + independent_variable=indvar)) + candidates.blocks.append(block) + return candidates diff --git a/pybnf/design/config.py b/pybnf/design/config.py new file mode 100644 index 00000000..b00f5760 --- /dev/null +++ b/pybnf/design/config.py @@ -0,0 +1,53 @@ +"""The configuration keys an experimental design reads (#574). + +These live here, in the design package, rather than beside one method, because two job types read +the same keys: ``job_type = design`` runs a design on its own, and ``job_type = +profile_likelihood`` can end by recommending one. Both schemas inherit this class, so the keys +mean the same thing and are documented once (ADR-0006 co-locates a method's own keys with the +method; a set of keys shared by two methods has to sit somewhere both can see). +""" + +from typing import Any + +from ..config_schema import PyBNFConfigModel + + +class DesignFields(PyBNFConfigModel): + """Optimal experimental design settings, shared by ``job_type = design`` and the design report + a ``profile_likelihood`` run can write. + + ``design_points`` is how many measurements to recommend. The same measurement may be + recommended more than once, which means measure it that many times. + + ``design_criterion`` is what makes one design better than another: ``a`` for the average + variance of the parameters (the default, and the classical c-criterion when + ``design_target`` names a single parameter), ``d`` for the volume of the joint confidence + region, ``e`` for the worst-determined direction. + + ``design_target`` names the parameters the design is aimed at. Absent, it aims at all of them. + Only the A-criterion can use it; the other two are properties of the whole information matrix. + + ``design_observables`` restricts the candidate measurements to a named set of observables, for + when only some assays can actually be run. Absent, every observable the fit already measures is + a candidate. + + ``design_confidence`` is the confidence level of the predicted intervals in the report. It has + the same meaning as ``profile_likelihood_confidence``, and a ``profile_likelihood`` run that + writes a design report uses that key instead so the two halves of its output agree. + + ``design_grid`` and ``design_t_end`` widen what the design is allowed to recommend. A design + can only propose a time the model is already simulated at, and for a time course PyBNF + simulates the times the data was measured at, so by default the only new measurement it can + propose is a repeat of an existing one. ``design_grid`` adds that many extra simulated times, + spread evenly from the first measurement out to ``design_t_end`` (which defaults to the last + measurement). Set both to let a design say "measure at a time you have never measured", which + is usually the point of asking. + """ + + design_points: int = 5 + design_criterion: str = 'a' + design_target: Any = None + design_observables: Any = None + design_confidence: float = 0.95 + design_grid: int = 0 + design_t_end: float = 0.0 diff --git a/pybnf/design/criteria.py b/pybnf/design/criteria.py new file mode 100644 index 00000000..c921f47e --- /dev/null +++ b/pybnf/design/criteria.py @@ -0,0 +1,188 @@ +"""Reducing an information matrix to one number, so two designs can be compared (#574). + +The expected Fisher information ``F`` is a square matrix, one row and column per free parameter, +in **sampling space** (ADR-0029 -- so a log-scaled parameter's entry is about its order of +magnitude, which is the scale it is fitted on). Its inverse is the covariance matrix the fit +would have, so ``(F^-1)_kk`` is the variance of parameter ``k`` and ``sqrt(threshold * +(F^-1)_kk)`` is the half-width of that parameter's confidence interval in the quadratic +approximation -- the same interval a profile-likelihood run traces, for the same threshold. + +A criterion turns that matrix into a single score so designs can be ranked: + +* ``'a'`` -- the summed variance of the parameters, or of a named subset. With one parameter + named this is the classical c-criterion. +* ``'d'`` -- the log determinant, which is the volume of the joint confidence region. +* ``'e'`` -- the smallest eigenvalue, which is the worst-determined direction. + +Singular information is not a numerical accident to be smoothed away. A parameter the data cannot +constrain at all leaves a direction with no information in it, and the honest reading is an +infinite variance, not a large one. Every function here decides that with one shared rule: an +eigenvalue below :data:`SINGULAR_TOL` times the largest one counts as zero, and a parameter with +any weight on such a direction has infinite variance. +""" + +import numpy as np + +#: An eigenvalue this far below the largest one carries no information. Relative, because the +#: information matrix has the units of the data and can sit anywhere on the number line. +SINGULAR_TOL = 1e-10 + +#: How much of a parameter's own axis has to lie in the uninformed directions before its variance +#: is infinite rather than merely large. Squared weights, so this is a very small angle. +NULL_COMPONENT_TOL = 1e-10 + +#: The criteria a design run may be scored with. +CRITERIA = ('a', 'd', 'e') + +#: The full name of each criterion, for messages and report headers. +CRITERION_NAMES = { + 'a': 'A-optimal (average parameter variance)', + 'd': 'D-optimal (confidence region volume)', + 'e': 'E-optimal (worst-determined direction)', +} + + +def lower_is_better(criterion): + """Whether a smaller :func:`criterion_value` is the better design. + + The A-criterion is a variance, so smaller is better; the other two are amounts of + information, so larger is. :func:`criterion_score` hides this (it is always maximized) but a + report has to say which way its numbers read.""" + return criterion == 'a' + + +def _spectrum(information): + """The eigenvalues and eigenvectors of ``information``, plus which eigenvalues count as zero. + + Symmetrized first: every term summed into the matrix is a symmetric outer product, so any + asymmetry is rounding, and ``eigh`` reads only one triangle anyway.""" + matrix = np.asarray(information, dtype=float) + matrix = 0.5 * (matrix + matrix.T) + values, vectors = np.linalg.eigh(matrix) + largest = float(values[-1]) if values.size else 0.0 + cutoff = SINGULAR_TOL * largest if largest > 0.0 else 0.0 + uninformed = values <= cutoff + return values, vectors, uninformed + + +def is_singular(information): + """Whether any direction in parameter space carries no information at all. + + A singular information matrix means some combination of the parameters is invisible to the + data. No amount of care with the criterion changes that; the design has to add a measurement + that sees the missing direction.""" + if np.size(information) == 0: + return False + _values, _vectors, uninformed = _spectrum(information) + return bool(uninformed.any()) + + +def parameter_variances(information): + """Each parameter's variance -- the diagonal of the inverse information -- in sampling space. + + Infinite for a parameter that lies (even partly) along a direction the data does not + constrain, which is the correct reading rather than a failure: no finite confidence interval + exists for it. Computed from the eigendecomposition rather than by inverting, so the + uninformed directions can be recognized instead of producing a huge finite number.""" + values, vectors, uninformed = _spectrum(information) + informed = ~uninformed + weights = vectors ** 2 # row k: parameter k's weight per axis + variances = np.full(values.shape, np.inf) + lost = weights[:, uninformed].sum(axis=1) if uninformed.any() else np.zeros(values.shape) + usable = lost <= NULL_COMPONENT_TOL + if informed.any(): + finite = (weights[:, informed] / values[informed]).sum(axis=1) + variances[usable] = finite[usable] + return variances + + +def log_determinant(information): + """``log det F``, or ``-inf`` when the information is singular (a confidence region of + unbounded volume).""" + values, _vectors, uninformed = _spectrum(information) + if uninformed.any(): + return -np.inf + return float(np.log(values).sum()) + + +def smallest_eigenvalue(information): + """The information in the worst-determined direction, floored at zero (a tiny negative + eigenvalue is rounding: every term summed into the matrix is positive semi-definite).""" + values, _vectors, _uninformed = _spectrum(information) + return max(float(values[0]), 0.0) if values.size else 0.0 + + +def criterion_value(information, criterion, targets=None): + """The criterion read the way a person would state it: a summed variance for ``'a'``, a log + determinant for ``'d'``, the smallest eigenvalue for ``'e'``. + + ``targets`` is the list of parameter indices the A-criterion sums over (``None`` -> all of + them). Use :func:`lower_is_better` to know which direction is an improvement.""" + if criterion == 'a': + variances = parameter_variances(information) + chosen = variances if targets is None else variances[np.asarray(targets, dtype=int)] + return float(chosen.sum()) + if criterion == 'd': + return log_determinant(information) + if criterion == 'e': + return smallest_eigenvalue(information) + raise ValueError('unknown design criterion %r' % (criterion,)) + + +def criterion_score(information, criterion, targets=None): + """The criterion as something to **maximize**, so the selection loop never has to ask which + way a given criterion runs. The A-criterion's summed variance is negated; the others are + already amounts of information.""" + value = criterion_value(information, criterion, targets) + return -value if lower_is_better(criterion) else value + + +def null_space_gain(information, block, targets=None): + """How much of ``block`` falls in the directions ``information`` currently knows nothing about. + + While the information is singular every criterion is pinned at its worst value -- an infinite + variance, a log determinant of ``-inf``, a smallest eigenvalue of zero -- so none of them can + tell two candidates apart. This can, and it asks the right question at that moment: of the + directions the data does not yet see, how much does this measurement see? The selection loop + uses it until the information becomes invertible and then goes back to the requested + criterion. + + ``targets`` restricts the accounting to the uninformed directions that the target parameters + actually lie along, so a c-criterion run is not sent off to fix a direction nobody asked + about.""" + values, vectors, uninformed = _spectrum(information) + if not uninformed.any(): + return 0.0 + axes = vectors[:, uninformed] # columns spanning the unseen space + if targets is not None: + idx = np.asarray(targets, dtype=int) + # Keep only the unseen directions the targets have weight on. An unseen direction + # orthogonal to every target cannot be why a target's variance is infinite. + relevant = (axes[idx, :] ** 2).sum(axis=0) > NULL_COMPONENT_TOL + if not relevant.any(): + return 0.0 + axes = axes[:, relevant] + return float(np.einsum('ij,jk,ki->', axes.T, np.asarray(block, dtype=float), axes)) + + +def unidentified_parameters(information, param_names): + """The parameters this information matrix leaves with an infinite variance. + + Called on the information of the *largest possible* design -- everything already measured plus + every candidate at once -- it names the parameters no experiment in the candidate space can + pin down. That is a structural statement about the model and the observables, not a shortage + of data, so it is reported rather than optimized around.""" + variances = parameter_variances(information) + return [name for name, var in zip(param_names, variances) if not np.isfinite(var)] + + +def interval_half_widths(information, threshold): + """Each parameter's confidence-interval half-width in **sampling space**, at the ``Delta + chi2`` threshold a profile-likelihood run would use. + + This is the quadratic (Wald) approximation to the profile interval: the profile of a + parameter near the optimum is the parabola ``Delta chi2 = (theta_k - theta*_k)^2 / + (F^-1)_kk``, which crosses the threshold at ``+- sqrt(threshold * (F^-1)_kk)``. For a + linear model the approximation is exact. Infinite for a parameter with infinite variance, + which reads as an open interval.""" + return np.sqrt(threshold * parameter_variances(information)) diff --git a/pybnf/design/greedy.py b/pybnf/design/greedy.py new file mode 100644 index 00000000..67f3ac30 --- /dev/null +++ b/pybnf/design/greedy.py @@ -0,0 +1,222 @@ +"""Choosing the best few measurements out of the whole candidate space (#574). + +Picking the best five of two hundred candidate measurements is a subset-selection problem, not a +continuous one: there are far too many subsets to try them all. The standard answer, and the one +here, is to take the measurements one at a time, each time adding whichever remaining candidate +improves the criterion most. It is not guaranteed to find the very best subset, but it is simple, +it always terminates, and for the D-criterion it is the classical exchange algorithm's forward +half. + +The same candidate may be chosen more than once. That is not a bug to be suppressed: choosing a +point twice means measuring it twice, and it is exactly the right recommendation when the +precision of one measurement, rather than the shape of the trajectory, is what limits you. + +One special case has to be handled explicitly. When the information matrix says nothing at all +about a direction the criterion cares about, the criterion sits at its worst possible value for +every candidate and cannot tell them apart -- an infinite variance is an infinite variance whatever +you add to it. The selection loop notices and asks a different question until that stops being +true: of the directions nothing yet sees, which candidate sees the most +(:func:`~pybnf.design.criteria.null_space_gain`)? Then it goes back to the requested criterion for +the rest of the picks. A design aimed at one parameter is not blocked by some *other* combination +of parameters being invisible, so this only fires when the target itself is the problem. +""" + +from dataclasses import dataclass, field +from typing import Any + +import numpy as np + +from ..printing import PybnfError +from .criteria import ( + CRITERIA, + CRITERION_NAMES, + criterion_score, + criterion_value, + lower_is_better, + null_space_gain, + unidentified_parameters, +) + + +@dataclass +class DesignResult: + """A finished design: what to measure, and what it is expected to buy. + + ``measurements`` are in the order they were chosen, repeats included, so reading them in order + shows how much each successive measurement is still worth. ``baseline`` is the information the + existing data already carries and ``information`` is what it would become; ``trace`` is the + criterion after each pick, so a design that stops paying off is visible rather than implied.""" + + criterion: str + targets: list + target_names: list + measurements: list = field(default_factory=list) + trace: list = field(default_factory=list) + baseline: Any = None + information: Any = None + escaped_singular: int = 0 + truncated: bool = False + + @property + def baseline_value(self): + """The criterion before any of the recommended measurements are made.""" + return criterion_value(self.baseline, self.criterion, self.targets) + + @property + def value(self): + """The criterion once every recommended measurement has been made.""" + return criterion_value(self.information, self.criterion, self.targets) + + @property + def criterion_name(self): + return CRITERION_NAMES[self.criterion] + + def grouped(self): + """The recommendation as ``(measurement, replicates, first_rank)`` in the order the + measurements were first chosen -- the reading a person wants, where choosing one point + three times is one row saying "measure it three times".""" + order, counts = [], {} + for rank, measurement in enumerate(self.measurements, start=1): + if measurement not in counts: + counts[measurement] = [0, rank] + order.append(measurement) + counts[measurement][0] += 1 + return [(m, counts[m][0], counts[m][1]) for m in order] + + +def resolve_targets(variables, names, criterion): + """The parameter indices the A-criterion sums variances over, validated against the fit. + + ``names`` empty or absent means every free parameter. Naming exactly one makes this the + classical c-criterion: minimize the variance of that one parameter, which is what a + profile-likelihood verdict about a single parameter asks for. + + Targets are refused for the D and E criteria rather than quietly ignored. Both are properties + of the whole information matrix -- a volume and a worst direction -- so restricting them to a + subset of parameters means something different from what a reader would assume, and a request + that cannot be honoured as written should say so.""" + if criterion not in CRITERIA: + raise PybnfError( + "design_criterion must be one of %s, not %r." + % (', '.join(CRITERIA), criterion), + hint='; '.join('%s = %s' % (code, CRITERION_NAMES[code]) for code in CRITERIA)) + names = list(names or []) + if not names: + return list(range(len(variables))) + if criterion != 'a': + raise PybnfError( + "design_target names the parameters to measure the design against, which only the " + "A-criterion (design_criterion = a) can use, but design_criterion = %s." % criterion, + hint="Set design_criterion = a to target %s, or drop design_target." + % ', '.join(names)) + by_name = {v.name: i for i, v in enumerate(variables)} + unknown = [n for n in names if n not in by_name] + if unknown: + raise PybnfError( + "design_target names %s, which %s not a free parameter of this fit." + % (', '.join(unknown), 'is' if len(unknown) == 1 else 'are'), + hint="List only free-parameter ids (the fit declares: %s)." + % ', '.join(v.name for v in variables)) + return [by_name[n] for n in names] + + +def require_identifiable(baseline, candidates, param_names, targets): + """Refuse a design whose targets no experiment in the candidate space could ever pin down. + + Adding every candidate at once is the most a design over this space can know. A target still + left with an infinite variance there is not short of data: no measurement of these observables, + at any simulated time, tells the model apart along that direction. That is structural + non-identifiability, and the fix is a different model or a different observable, not a + different design.""" + n_param = len(param_names) + everything = np.asarray(baseline, dtype=float) + candidates.total(n_param) + lost = unidentified_parameters(everything, param_names) + targeted = [param_names[i] for i in targets] + blocking = [name for name in lost if name in targeted] + if blocking: + raise PybnfError( + "No design over these observables and times can determine %s. Measuring every " + "candidate point at once still leaves %s undetermined, so the model and the data it " + "can produce do not distinguish %s from the other parameters at all." + % (', '.join(blocking), 'it' if len(blocking) == 1 else 'them', + 'it' if len(blocking) == 1 else 'them'), + hint="This is structural non-identifiability, not a shortage of data. Measure a " + "different observable, fix one of the parameters, or reparameterize the model. " + "A profile-likelihood run (job_type = profile_likelihood) reports which " + "parameters are structurally non-identifiable.") + return lost + + +def _criterion_is_blind(information, criterion, targets): + """Whether the criterion can still tell two candidates apart at this information. + + It cannot when the quantity it measures is already at its worst possible value for every + candidate: an infinite variance for a target the data cannot see at all, a log determinant of + ``-inf``, a smallest eigenvalue of zero. Note that a singular information matrix is not by + itself blinding -- a design aimed at one parameter is perfectly well defined while some + *other* combination of parameters remains unseen, and it should go on optimizing what it was + asked to.""" + score = criterion_score(information, criterion, targets) + return not np.isfinite(score) or (criterion == 'e' and score <= 0.0) + + +def select_design(baseline, candidates, n_points, criterion, targets, param_names): + """Choose ``n_points`` measurements, one at a time, each the best next addition. + + ``baseline`` is the information the existing data carries, ``candidates`` the + :class:`~pybnf.design.candidates.CandidateSet` to choose from. Returns a + :class:`DesignResult`. Ties break toward the earliest candidate, which is the earliest time of + the first observable of the first experiment, so the same inputs always give the same design. + """ + if n_points <= 0: + raise PybnfError("design_points must be at least 1, not %d." % n_points, + hint="design_points is how many measurements to recommend.") + if not len(candidates): + raise PybnfError( + "There are no candidate measurements to choose from.", + hint="A candidate is an observable this fit already measures, at a time its " + "simulation passes through. Check that the fit scores at least one observable.") + + information = np.array(baseline, dtype=float) + result = DesignResult(criterion=criterion, targets=list(targets), + target_names=[param_names[i] for i in targets], + baseline=np.array(baseline, dtype=float)) + blocks = candidates.blocks + for _pick in range(n_points): + blind = _criterion_is_blind(information, criterion, targets) + if blind: + # The criterion cannot choose yet, so chase the directions nothing sees instead. + scores = [null_space_gain(information, block, targets) for block in blocks] + result.escaped_singular += 1 + else: + scores = [criterion_score(information + block, criterion, targets) + for block in blocks] + best = int(np.argmax(scores)) + if blind and scores[best] <= 0.0: + # No candidate sees any of the missing directions. require_identifiable refuses that + # before a design is ever selected, so this is unreachable from a job; stop and say so + # rather than pick an arbitrary point that buys nothing. + result.truncated = True + break + information = information + blocks[best] + result.measurements.append(candidates.measurements[best]) + result.trace.append(criterion_value(information, criterion, targets)) + result.information = information + return result + + +def improvement(result): + """How much better the criterion got, as a plain ratio, or ``None`` when it cannot be stated. + + For the A-criterion this is the factor the summed variance shrank by, so 4.0 means the + variance is a quarter of what it was and the confidence interval is half as wide. A baseline + that was infinite (a parameter the existing data cannot determine at all) has no ratio, which + is itself the headline: the design goes from no answer to an answer.""" + before, after = result.baseline_value, result.value + if not np.isfinite(before) or not np.isfinite(after): + return None + if lower_is_better(result.criterion): + return None if after <= 0.0 else before / after + if result.criterion == 'd': + return float(np.exp(after - before)) # a ratio of determinants, not of logs + return None if before <= 0.0 else after / before diff --git a/pybnf/design/report.py b/pybnf/design/report.py new file mode 100644 index 00000000..e4d0a659 --- /dev/null +++ b/pybnf/design/report.py @@ -0,0 +1,124 @@ +"""Writing a design down: what to measure, and what it is expected to buy (#574). + +The report has two halves, because a recommendation nobody can check is not worth much. + +The first half is the recommendation itself: the measurements to make, in the order they were +chosen, with a count when the same point is chosen more than once. The second half is the reason, +stated in the units the user has already seen from a profile-likelihood run -- each parameter's +confidence interval as it stands now, and as it would be once the recommended measurements are in +hand. Those predicted intervals come from the same information matrix the design was chosen with, +read through the quadratic approximation to the profile: the interval is +``theta* +- sqrt(threshold * variance)`` in the parameter's own fitted scale. For a linear model +that is exact, and for anything else it is the local approximation, which is also all the design +itself ever claimed to be. +""" + +import numpy as np + +from .criteria import interval_half_widths, lower_is_better +from .greedy import improvement + + +def _interval(variable, centre_u, half_width): + """One parameter's predicted confidence interval in its own units, or ``None`` when the + information leaves it undetermined and the interval is open.""" + if not np.isfinite(half_width): + return None + return (float(variable.from_sampling_space(centre_u - half_width)), + float(variable.from_sampling_space(centre_u + half_width))) + + +def predicted_intervals(result, variables, u_star, threshold): + """Every parameter's interval before and after the design, plus how much it shrinks. + + Each row is a dict with the parameter's name, its value at the best fit, the interval the + existing data supports, the interval the design would support, and ``width_ratio`` -- the + designed half-width over the current one, so 0.5 means the interval halves. The ratio is + ``None`` when the current interval is open, which is the strongest result there is: the design + replaces no answer with an answer.""" + current = interval_half_widths(result.baseline, threshold) + designed = interval_half_widths(result.information, threshold) + rows = [] + for index, variable in enumerate(variables): + centre = float(u_star[index]) + ratio = None + if np.isfinite(current[index]) and current[index] > 0.0: + ratio = float(designed[index] / current[index]) + rows.append({ + 'name': variable.name, + 'best': float(variable.from_sampling_space(centre)), + 'current': _interval(variable, centre, current[index]), + 'designed': _interval(variable, centre, designed[index]), + 'width_ratio': ratio, + }) + return rows + + +def _format_interval(interval): + return 'open' if interval is None else '[%.6g, %.6g]' % interval + + +def write_design_report(path, result, variables, u_star, threshold, confidence): + """Write the design report to ``path`` as a tab-delimited file with commented headers, the + same shape as the profile-likelihood summary beside it.""" + rows = predicted_intervals(result, variables, u_star, threshold) + factor = improvement(result) + with open(path, 'w') as handle: + handle.write('# criterion=%s\t%s\n' % (result.criterion, result.criterion_name)) + handle.write('# targets=%s\n' % (', '.join(result.target_names) or 'all parameters')) + handle.write('# confidence=%g\tdelta_chi2_threshold=%g\tdof=1\n' + % (confidence, threshold)) + handle.write('# criterion_before=%.10g\tcriterion_after=%.10g\t%s_is_better\n' + % (result.baseline_value, result.value, + 'lower' if lower_is_better(result.criterion) else 'higher')) + if factor is not None: + handle.write('# improvement_factor=%.6g\n' % factor) + if result.truncated: + handle.write('# fewer measurements than asked for: no remaining candidate adds ' + 'anything the criterion can use\n') + handle.write('#\n# recommended measurements, in the order they were chosen\n') + handle.write('# rank\tmodel\texperiment\tobservable\tindependent_variable\tvalue\t' + 'replicates\tcriterion_after\n') + criterion_at = {rank: value for rank, value in enumerate(result.trace, start=1)} + for measurement, replicates, rank in result.grouped(): + handle.write('%d\t%s\t%s\t%s\t%s\t%.10g\t%d\t%.10g\n' % ( + rank, measurement.model, measurement.experiment or '-', measurement.observable, + measurement.independent_variable, measurement.time, replicates, + criterion_at[rank])) + handle.write('#\n# predicted confidence intervals, before and after\n') + handle.write('# parameter\tbest\tcurrent_low\tcurrent_high\tdesigned_low\t' + 'designed_high\twidth_ratio\n') + for row in rows: + current = row['current'] or (None, None) + designed = row['designed'] or (None, None) + handle.write('%s\t%.10g\t%s\t%s\t%s\t%s\t%s\n' % ( + row['name'], row['best'], + 'None' if current[0] is None else '%.10g' % current[0], + 'None' if current[1] is None else '%.10g' % current[1], + 'None' if designed[0] is None else '%.10g' % designed[0], + 'None' if designed[1] is None else '%.10g' % designed[1], + 'None' if row['width_ratio'] is None else '%.6g' % row['width_ratio'])) + + +def format_design_summary(result, variables, u_star, threshold): + """The same design as a short block of lines for the terminal: what to measure, and what it + does to the intervals of the parameters the design was aimed at.""" + lines = ['Recommended next measurements (%s, aimed at %s):' + % (result.criterion_name, ', '.join(result.target_names) or 'all parameters')] + for measurement, replicates, rank in result.grouped(): + times = '' if replicates == 1 else ' (measure %d times)' % replicates + lines.append(' %d. %s%s' % (rank, measurement, times)) + if result.truncated: + lines.append(' (stopped early: no remaining candidate adds anything the criterion ' + 'can use)') + rows = {row['name']: row for row in predicted_intervals(result, variables, u_star, threshold)} + reported = result.target_names or [v.name for v in variables] + lines.append('Predicted confidence intervals:') + for name in reported: + row = rows[name] + shrink = ('' if row['width_ratio'] is None + else ' (%.3g times as wide)' % row['width_ratio']) + lines.append(' %-16s now %s -> %s%s' + % (name, _format_interval(row['current']), + _format_interval(row['designed']), shrink)) + return lines diff --git a/pybnf/gradient/__init__.py b/pybnf/gradient/__init__.py index 29cc11b1..a3b349ed 100644 --- a/pybnf/gradient/__init__.py +++ b/pybnf/gradient/__init__.py @@ -18,7 +18,9 @@ (``job_type = gntr``, #481) it additionally assembles the expected-Fisher / Gauss-Newton **Hessian** in the same point walk as the scalar gradient (``assemble_gradient_and_fisher_hessian`` + the constraint sibling - ``assemble_constraint_hessian``); ``assemble_fisher_hessian`` remains the standalone API. + ``assemble_constraint_hessian``); ``assemble_fisher_hessian`` remains the standalone API, and + ``iter_fisher_points`` yields the same information one scored point at a time, which is what + optimal experimental design (:mod:`pybnf.design`, #574) scores a candidate measurement with. The capability gate and the per-layer math are documented in ``docs/gradient_fitting.rst``. """ @@ -48,6 +50,7 @@ assemble_fisher_hessian, assemble_gradient_and_fisher_hessian, assemble_gaussian_gradient, + iter_fisher_points, ) from .marginal_time import assemble_marginal_time_gradient @@ -75,4 +78,5 @@ 'assemble_constraint_gradient', 'assemble_fisher_hessian', 'assemble_constraint_hessian', + 'iter_fisher_points', ] diff --git a/pybnf/gradient/assembly.py b/pybnf/gradient/assembly.py index 00295c15..1a0f3b12 100644 --- a/pybnf/gradient/assembly.py +++ b/pybnf/gradient/assembly.py @@ -478,6 +478,45 @@ def _accumulate_fisher_point(objective, sim_data, exp_data, index, point, hessia hessian += weight * noise_block +def iter_fisher_points(objective, experiments, free_params): + """Yield each scored point's **own** expected-Fisher block, instead of their sum (#574). + + :func:`assemble_fisher_hessian` adds every point's rank-1 terms into one matrix. Optimal + experimental design needs the terms kept apart: a design is a *subset* of measurements, so + scoring one means adding up the blocks of the points it contains and leaving the rest out. + Because the Fisher information is a plain sum over points, that is all a design score is. + + Each item is ``(experiment_index, exp_row, col_name, block)``, where ``block`` is the + ``(n_param, n_param)`` matrix that point contributes -- the same location + noise terms + :func:`_accumulate_fisher_point` adds, already weight-folded and already in sampling space + (the ``d theta/d u`` factors applied on both axes, ADR-0029), so the blocks add straight onto + a Hessian :func:`assemble_fisher_hessian` returned. Summing every yielded block reproduces + that Hessian. + + ``experiments`` and ``free_params`` are exactly what the other assemblers take. The point + walk is the shared :func:`_iter_scored_points` scaffold, so the points, the row matching, the + NaN skip and the transform chain rule are identical to the ones the fit scores.""" + names = [p.name for p in free_params] + index = {name: j for j, name in enumerate(names)} + n_param = len(free_params) + + # Seed the estimated-noise / profiled-scale reads from this point exactly as + # assemble_fisher_hessian does, so a free sigma resolves the same way here. + existing = getattr(objective, '_pset_values', None) or {} + objective._pset_values = {**existing, **{p.name: p.value for p in free_params}} + _seed_profiled_noise(objective, experiments) + + factors = _sampling_scale_factors(free_params) + sampling = np.outer(factors, factors) + for exp_index, (sim_data, exp_data, routing, *rest) in enumerate(experiments): + data_key = rest[0] if rest else None + for point in _iter_scored_points(objective, sim_data, exp_data, routing, index, + n_param, "Fisher information", data_key): + block = np.zeros((n_param, n_param)) + _accumulate_fisher_point(objective, sim_data, exp_data, index, point, block) + yield exp_index, point[1], point[2], block * sampling + + def assemble_fisher_hessian(objective, experiments, free_params): """Assemble the expected-Fisher / Gauss-Newton **Hessian** ``H`` (n_param x n_param), summed across experiments. This standalone API produces the same curvature the combined diff --git a/pybnf/parse.py b/pybnf/parse.py index d60b9de7..e55b0efe 100644 --- a/pybnf/parse.py +++ b/pybnf/parse.py @@ -54,6 +54,13 @@ def _parse_all(parser, text): # cross-parameter parallel-track cap (#467). 'profile_likelihood_max_iterations', 'profile_likelihood_max_points', 'profile_likelihood_reopt_max_iterations', 'profile_likelihood_max_parallel', + # Whether a profile-likelihood run ends by recommending the measurements to + # make next (job_type = design's report, #574): 0 = off. + 'profile_likelihood_design', + # optimal experimental design (job_type = design, #574): how many + # measurements to recommend, and how many extra times to simulate for it to + # choose from (0 = only the times already measured). + 'design_points', 'design_grid', # gradient optimizers (job_type = trf / lbfgs / gntr, #386/#481): the # int-valued tunables -- L-BFGS-B's curvature-history depth and the three # cycle budgets (runtime-guarded RUNTIME_KEYS, defaulting to max_iterations). @@ -123,6 +130,10 @@ def _parse_all(parser, text): 'profile_likelihood_min_step', 'profile_likelihood_max_step', 'profile_likelihood_dchi2_target', 'profile_likelihood_grad_tol', 'profile_likelihood_step_tol', + # optimal experimental design (job_type = design, #574): the confidence + # level the report's predicted intervals are quoted at, and how far past the + # last measurement a recommendation may reach. + 'design_confidence', 'design_t_end', # CVODE tolerances for the bngsim SBML/Antimony backend (#546). Unset # leaves rtol at the backend default and DERIVES atol from the model's # own state scale; stating either pins it. @@ -159,11 +170,18 @@ def _parse_all(parser, text): 'qualitative_loss', # What a gradient fit does when bngsim declines a model's analytic # sensitivity RHS (#606, ADR-0121): warn | error | ignore. - 'sensitivity_fallback'] + 'sensitivity_fallback', + # optimal experimental design (#574): what makes one design better than + # another -- a (average parameter variance) | d (confidence region volume) | + # e (worst-determined direction). + 'design_criterion'] multstrkeys = ['worker_nodes', 'postprocess', 'output_trajectory', 'output_noise_trajectory', # profile likelihood (#446/#466): the subset of free parameters to profile # (a list of parameter ids; absent -> profile every free parameter). 'profile_likelihood_params', + # optimal experimental design (#574): the free parameters a design is aimed + # at, and the observables it may recommend measuring (absent -> all of them). + 'design_target', 'design_observables', # qualitative scale as a fittable parameter: the two-token # value `fit ` ties every qualitative constraint's scale to a free parameter. 'qualitative_scale'] diff --git a/pybnf/quantiles.py b/pybnf/quantiles.py new file mode 100644 index 00000000..9a2c198d --- /dev/null +++ b/pybnf/quantiles.py @@ -0,0 +1,62 @@ +"""Dependency-free quantiles for confidence thresholds. + +PyBNF's production loop never imports scipy (ADR-0007), so the two quantiles its confidence +statements need are computed here. Both a profile-likelihood run (which asks where a profile +crosses its ``Delta chi2`` threshold) and an experimental design (which reports the confidence +interval an information matrix implies) read the same threshold from the same confidence level, so +they share this rather than each carrying a copy. +""" + +import math + +from .printing import PybnfError + + +def normal_quantile(p): + """Standard-normal inverse cumulative distribution (the probit) via Acklam's rational + approximation, refined by one Halley step against :func:`math.erf`. + + Accurate to full double precision after the refinement, far more than a chi-square threshold + needs. ``0 < p < 1``.""" + a = (-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, + 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00) + b = (-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, + 6.680131188771972e+01, -1.328068155288572e+01) + c = (-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, + -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00) + d = (7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, + 3.754408661907416e+00) + plow, phigh = 0.02425, 1.0 - 0.02425 + if p < plow: + q = math.sqrt(-2.0 * math.log(p)) + x = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \ + ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0) + elif p <= phigh: + q = p - 0.5 + r = q * q + x = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / \ + (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0) + else: + q = math.sqrt(-2.0 * math.log(1.0 - p)) + x = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \ + ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0) + # One Halley refinement using the exact erf-based CDF. + e = 0.5 * math.erfc(-x / math.sqrt(2.0)) - p + u = e * math.sqrt(2.0 * math.pi) * math.exp(x * x / 2.0) + x = x - u / (1.0 + x * u / 2.0) + return x + + +def chi2_quantile_1dof(confidence, key='confidence'): + """The chi-square (1 degree of freedom) quantile at probability ``confidence`` -- the profile + ``Delta chi2`` threshold (Raue et al. 2009), and the same threshold a design's predicted + intervals are quoted at. + + A single parameter has one degree of freedom, and ``chi2_1 = Z**2`` with ``Z ~ N(0, 1)``, so + ``P(chi2_1 <= x) = 2*Phi(sqrt(x)) - 1`` and the quantile is + ``Phi^-1((1 + confidence) / 2)**2`` (0.95 gives 3.8415). ``key`` names the configuration key + the confidence level came from, so a rejected value points at the right setting.""" + if not (0.0 < confidence < 1.0): + raise PybnfError("%s must be strictly between 0 and 1, got %r." % (key, confidence)) + z = normal_quantile(0.5 * (1.0 + confidence)) + return z * z diff --git a/tests/recovery_harness.py b/tests/recovery_harness.py index a59530c3..1600633c 100644 --- a/tests/recovery_harness.py +++ b/tests/recovery_harness.py @@ -52,6 +52,7 @@ 'gntr': algorithms.GNTRAlgorithm, # general-objective Fisher/Gauss-Newton trust region (#481) 'profile_likelihood': algorithms.ProfileLikelihoodAlgorithm, # PL identifiability (#446/#466) 'ms': algorithms.MultipleShootingAlgorithm, # multiple shooting (#563/ADR-0110) + 'design': algorithms.ExperimentalDesignAlgorithm, # experimental design (#574) } diff --git a/tests/test_config_schema.py b/tests/test_config_schema.py index 9e2e7f9a..cfb540e9 100644 --- a/tests/test_config_schema.py +++ b/tests/test_config_schema.py @@ -354,14 +354,15 @@ def test_migrated_methods_so_far(self): # The gradient-based optimizers (#386/#481) land with their own schemas: trf # (trust-region least-squares) with TRFConfig, lbfgs (L-BFGS-B) with LBFGSConfig, # and gntr (general-objective Fisher/Gauss-Newton trust region) with GNTRConfig. - # profile_likelihood (#446/#466) lands with ProfileLikelihoodConfig, and ms - # (multiple shooting, #563/ADR-0110) with MSConfig. + # profile_likelihood (#446/#466) lands with ProfileLikelihoodConfig, ms + # (multiple shooting, #563/ADR-0110) with MSConfig, and design (optimal + # experimental design, #574) with DesignConfig. # Only 'check' remains unmigrated. Each step extends this set -- a ratchet. from pybnf.registry import FIT_TYPE_REGISTRY migrated = {c for c, e in FIT_TYPE_REGISTRY.items() if e.schema is not None} assert migrated == {'pso', 'de', 'ade', 'ss', 'sim', 'powell', 'cmaes', 'mh', 'pt', 'sa', 'am', 'dream', 'p_dream', 'hmc', 'trf', 'lbfgs', - 'gntr', 'profile_likelihood', 'ms'} + 'gntr', 'profile_likelihood', 'ms', 'design'} assert FIT_TYPE_REGISTRY['check'].schema is None diff --git a/tests/test_design.py b/tests/test_design.py new file mode 100644 index 00000000..f610efca --- /dev/null +++ b/tests/test_design.py @@ -0,0 +1,867 @@ +"""Offline unit tests for optimal experimental design (#574). + +These run against hand-built simulated trajectories and sensitivity tensors, with no simulation +backend and no scheduler, exactly as ``tests/test_gradient_assembly.py`` does. The point is that +the answers a design gives for simple models are known in advance, so they can be checked rather +than merely inspected: + +* **A straight line.** For ``y(t) = a + b*t`` measured on ``[0, 1]`` with constant noise, the + D-optimal design puts its measurements at the two ends of the interval. This is the oldest + result in the subject (Elfving 1952) and it is what a correct D-criterion has to reproduce. +* **An exponential decay.** For ``y(t) = S0*exp(-k*t)`` with ``S0`` known, the single most + informative time for ``k`` is one lifetime, ``t = 1/k``: the sensitivity ``t*S0*exp(-k*t)`` is + largest there. A design aimed at ``k`` has to choose that time. +* **Two parallel channels.** For ``S(t) = S0*exp(-(k1+k2)*t)`` only the sum is observable, so no + measurement at any time separates ``k1`` from ``k2``. A design aimed at either has to say so + rather than recommend something. + +Two consistency checks anchor the machinery to the code it reuses. The information a design +reports is checked against :func:`~pybnf.gradient.assemble_fisher_hessian` assembled on a dataset +that literally contains the recommended measurements, so the planned-measurement construction is +not merely self-consistent. And the intervals the report predicts are checked against the +closed-form profile-likelihood interval of a linear-Gaussian problem, which is the same number +``tests/test_profile_likelihood.py`` checks its own confidence intervals against. +""" + +import numpy as np +import pytest + +from pybnf.data import Data, OutputSensitivities +from pybnf.design import ( + CandidateMeasurement, + DesignExperiment, + baseline_information, + candidate_information, + criterion_value, + improvement, + interval_half_widths, + is_singular, + measured_observables, + null_space_gain, + parameter_variances, + predicted_intervals, + require_identifiable, + resolve_targets, + select_design, + unidentified_parameters, + write_design_report, +) +from pybnf.gradient import assemble_fisher_hessian +from pybnf.gradient.routing import ExperimentRouting, ParamRoute, PARAM +from pybnf.objective import ChiSquareObjective, LikelihoodObjective +from pybnf.noise import FreeParameterSigma, Gaussian +from pybnf.printing import PybnfError +from pybnf.pset import FreeParameter +from pybnf.quantiles import chi2_quantile_1dof + +from pathlib import Path + +from pybnf._bngsim_caps import BNGSIM_HAS_OUTPUT_SENS +from pybnf.config import Configuration +from pybnf.parse import ploop + +from . import recovery_harness as H +from .test_profile_likelihood import TRUE_K, TRUE_S0, _decay_model + + +def _write_early_decay_exp(path, *, n=8, t_end=0.5, sd=2.0): + """A decay ``.exp`` measured only over the first fraction of a lifetime. + + ``Stot = S0*exp(-k*t)`` is nearly straight over such a short window, so the data pins down the + starting amount well and the decay rate badly. That is the situation an experimental design + exists for.""" + times = np.linspace(0.0, t_end, n) + observations = TRUE_S0 * np.exp(-TRUE_K * times) + lines = ['#\ttime\tStot\tStot_SD'] + lines += ['%.12g\t%.12g\t%.12g' % (t, o, sd) for t, o in zip(times, observations)] + Path(path).write_text('\n'.join(lines) + '\n') + return str(path) + + +# --------------------------------------------------------------------------- # +# Fixtures: a simulated trajectory with a hand-built sensitivity tensor +# --------------------------------------------------------------------------- # +def _sim(times, predictions, sensitivities): + """A simulated ``Data`` over ``times`` carrying its forward sensitivities. + + ``predictions`` maps an observable name to its values; ``sensitivities`` maps + ``(observable, parameter)`` to ``d(observable)/d(parameter)`` over the same times.""" + columns = list(predictions) + params = sorted({param for _col, param in sensitivities}) + array = np.column_stack([np.asarray(times, float)] + + [np.asarray(predictions[c], float) for c in columns]) + sim = Data.from_columns(array, ['time'] + columns) + tensor = np.zeros((len(times), len(columns), len(params))) + for i, col in enumerate(columns): + for j, param in enumerate(params): + if (col, param) in sensitivities: + tensor[:, i, j] = np.asarray(sensitivities[(col, param)], float) + sim.output_sensitivities = OutputSensitivities( + selectors=['observable:%s' % c for c in columns], + param_names=params, ic_species=[], d_param=tensor, d_ic=None) + return sim + + +def _exp(times, observations, sigma): + """An experimental ``Data`` with one observable and its ``_SD`` column.""" + columns = list(observations) + values = [np.asarray(times, float)] + headers = ['time'] + for col in columns: + values.append(np.asarray(observations[col], float)) + headers.append(col) + scale = sigma[col] if isinstance(sigma, dict) else sigma + values.append(np.full(len(times), scale, float) if np.isscalar(scale) + else np.asarray(scale, float)) + headers.append(col + '_SD') + return Data.from_columns(np.column_stack(values), headers) + + +def _routing(*params): + return ExperimentRouting(routes={ + name: ParamRoute.single(name, PARAM, name, 1.0) for name in params}) + + +def _free(*specs): + return [FreeParameter(n, t, lb, ub, value=v) for (n, t, lb, ub, v) in specs] + + +def _line_experiment(grid, measured): + """``y(t) = a + b*t``: the textbook design problem. ``d y/d a = 1`` and ``d y/d b = t``, + both independent of the parameters, so the information depends only on which times are + measured -- which is why the D-optimal answer is known exactly.""" + a, b = 1.0, 2.0 + sim = _sim(grid, {'y': a + b * np.asarray(grid)}, + {('y', 'a'): np.ones(len(grid)), ('y', 'b'): np.asarray(grid, float)}) + exp = _exp(measured, {'y': a + b * np.asarray(measured)}, 1.0) + return DesignExperiment(model='line', suffix='line', sim_data=sim, exp_data=exp, + routing=_routing('a', 'b')) + + +def _decay_experiment(grid, measured, k=0.4, s0=100.0, sigma=1.0): + """``S(t) = S0*exp(-k*t)`` with ``S0`` held fixed, so ``k`` is the only free parameter and + the most informative time is exactly one lifetime.""" + grid = np.asarray(grid, float) + measured = np.asarray(measured, float) + sim = _sim(grid, {'S': s0 * np.exp(-k * grid)}, + {('S', 'k'): -grid * s0 * np.exp(-k * grid)}) + exp = _exp(measured, {'S': s0 * np.exp(-k * measured)}, sigma) + return DesignExperiment(model='decay', suffix='decay', sim_data=sim, exp_data=exp, + routing=_routing('k')) + + +def _two_channel_experiment(grid, measured, k1=0.2, k2=0.2, s0=100.0): + """``S(t) = S0*exp(-(k1+k2)*t)``: only the sum is observable, so the two rates have + identical sensitivity columns and no measurement anywhere separates them.""" + grid = np.asarray(grid, float) + measured = np.asarray(measured, float) + decay = -grid * s0 * np.exp(-(k1 + k2) * grid) + sim = _sim(grid, {'S': s0 * np.exp(-(k1 + k2) * grid)}, + {('S', 'k1'): decay, ('S', 'k2'): decay}) + exp = _exp(measured, {'S': s0 * np.exp(-(k1 + k2) * measured)}, 1.0) + return DesignExperiment(model='two_channel', suffix='two_channel', sim_data=sim, + exp_data=exp, routing=_routing('k1', 'k2')) + + +def _design(experiment, free, *, points, criterion='a', targets=None, observables=None): + """Run a whole design over one experiment, the way the job does.""" + objective = ChiSquareObjective() + names = [p.name for p in free] + baseline = baseline_information(objective, [experiment], free) + candidates = candidate_information(objective, [experiment], free, observables=observables) + target_idx = resolve_targets(free, targets, criterion) + require_identifiable(baseline, candidates, names, target_idx) + return select_design(baseline, candidates, points, criterion, target_idx, names), candidates + + +# ============================================================ criteria math === + +def test_variances_of_a_diagonal_information_are_the_reciprocals(): + """The variance of a parameter is the diagonal of the inverted information, which for a + diagonal matrix is just one over each entry.""" + np.testing.assert_allclose(parameter_variances(np.diag([4.0, 0.25])), [0.25, 4.0]) + + +def test_a_parameter_with_no_information_has_infinite_variance(): + """A direction the data says nothing about gives an infinite variance, not a very large + one: no finite confidence interval exists for it.""" + information = np.diag([4.0, 0.0]) + assert is_singular(information) + variances = parameter_variances(information) + assert variances[0] == pytest.approx(0.25) + assert np.isinf(variances[1]) + assert unidentified_parameters(information, ['a', 'b']) == ['b'] + + +def test_a_combination_nobody_can_see_makes_both_parameters_infinite(): + """When only the *sum* of two parameters is visible, neither of them separately has a + finite variance, even though the matrix has a perfectly good non-zero entry.""" + information = np.array([[1.0, 1.0], [1.0, 1.0]]) # sees only k1 + k2 + assert unidentified_parameters(information, ['k1', 'k2']) == ['k1', 'k2'] + + +def test_criteria_read_the_matrix_the_way_their_names_say(): + information = np.diag([4.0, 1.0]) + assert criterion_value(information, 'a') == pytest.approx(1.25) + assert criterion_value(information, 'a', [0]) == pytest.approx(0.25) + assert criterion_value(information, 'd') == pytest.approx(np.log(4.0)) + assert criterion_value(information, 'e') == pytest.approx(1.0) + + +def test_null_space_gain_measures_only_the_unseen_direction(): + """A candidate that only reinforces what is already known scores zero; one that sees the + missing direction scores what it sees there.""" + information = np.diag([1.0, 0.0]) + seen_again = np.diag([5.0, 0.0]) + the_missing_one = np.diag([0.0, 3.0]) + assert null_space_gain(information, seen_again) == pytest.approx(0.0) + assert null_space_gain(information, the_missing_one) == pytest.approx(3.0) + + +def test_interval_half_width_is_the_profile_crossing_of_a_parabola(): + """For a quadratic objective the profile of a parameter is the parabola ``delta chi2 = + (theta - theta*)^2 / variance``, so it crosses the threshold at ``sqrt(threshold * + variance)``. That is the interval this reports, and for a linear model it is exact.""" + threshold = chi2_quantile_1dof(0.95) + half = interval_half_widths(np.diag([4.0, 0.25]), threshold) + np.testing.assert_allclose(half, np.sqrt(threshold * np.array([0.25, 4.0]))) + + +# ====================================================== candidate enumeration === + +def test_every_simulated_time_is_a_candidate_without_re_solving(): + """The candidate space is the whole simulated grid, for every observable the experiment + measures. Nothing is simulated to build it: the sensitivities at every simulated time were + already computed when the best fit was scored.""" + grid = np.linspace(0.0, 5.0, 26) + experiment = _decay_experiment(grid, measured=[1.0, 2.0]) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + candidates = candidate_information(ChiSquareObjective(), [experiment], free) + + assert len(candidates) == len(grid) + assert [c.time for c in candidates.measurements] == pytest.approx(list(grid)) + assert {c.observable for c in candidates.measurements} == {'S'} + assert candidates.measurements[0] == CandidateMeasurement( + model='decay', experiment='decay', observable='S', time=0.0, + independent_variable='time') + + +def test_a_candidate_carries_the_information_that_point_would_add(): + """One candidate's information is exactly what assembling the Fisher matrix over a dataset + holding that single planned measurement gives -- the same routine the ``gntr`` optimizer + uses, on the same point.""" + grid = np.array([0.0, 1.0, 2.5, 4.0]) + experiment = _decay_experiment(grid, measured=[1.0], k=0.4, s0=100.0, sigma=2.0) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + candidates = candidate_information(ChiSquareObjective(), [experiment], free) + + # d S/d k at t = 2.5 over sigma, squared: the Gaussian information of one measurement. + expected = (2.5 * 100.0 * np.exp(-0.4 * 2.5) / 2.0) ** 2 + at_two_point_five = candidates.blocks[list(grid).index(2.5)] + assert at_two_point_five[0, 0] == pytest.approx(expected) + + +def test_an_unmeasured_column_is_not_a_candidate(): + """A column that is present but blank has never been measured, so its noise model has never + been exercised and it is not offered as a candidate.""" + grid = np.array([0.0, 1.0, 2.0]) + sim = _sim(grid, {'S': [100.0, 67.0, 45.0], 'P': [0.0, 33.0, 55.0]}, + {('S', 'k'): [0.0, -67.0, -90.0], ('P', 'k'): [0.0, 67.0, 90.0]}) + exp = _exp(grid, {'S': [100.0, 67.0, 45.0], 'P': [np.nan, np.nan, np.nan]}, 1.0) + experiment = DesignExperiment(model='m', suffix='m', sim_data=sim, exp_data=exp, + routing=_routing('k')) + assert measured_observables(ChiSquareObjective(), experiment) == ['S'] + + +def test_observables_can_be_restricted_to_the_assays_that_can_be_run(): + grid = np.array([0.0, 1.0, 2.0]) + sim = _sim(grid, {'S': [100.0, 67.0, 45.0], 'P': [0.0, 33.0, 55.0]}, + {('S', 'k'): [0.0, -67.0, -90.0], ('P', 'k'): [0.0, 67.0, 90.0]}) + exp = _exp(grid, {'S': [100.0, 67.0, 45.0], 'P': [0.0, 33.0, 55.0]}, 1.0) + experiment = DesignExperiment(model='m', suffix='m', sim_data=sim, exp_data=exp, + routing=_routing('k')) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + both = candidate_information(ChiSquareObjective(), [experiment], free) + only_p = candidate_information(ChiSquareObjective(), [experiment], free, observables=['P']) + + assert {c.observable for c in both.measurements} == {'S', 'P'} + assert {c.observable for c in only_p.measurements} == {'P'} + + +def test_a_planned_measurement_borrows_the_precision_of_the_nearest_real_one(): + """A noise scale read from a data column has no value at a time nobody has measured, so the + planned measurement takes it from the nearest real measurement of that same observable. Here + the early measurements are precise and the late ones are not, and the candidates inherit that + split at the midpoint.""" + grid = np.array([0.0, 1.0, 2.0, 3.0, 4.0]) + experiment = _decay_experiment(grid, measured=[1.0, 3.0], k=0.4, s0=100.0, + sigma=np.array([1.0, 10.0])) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + candidates = candidate_information(ChiSquareObjective(), [experiment], free) + + def sigma_used(time): + block = candidates.blocks[list(grid).index(time)] + sensitivity = time * 100.0 * np.exp(-0.4 * time) + return sensitivity / np.sqrt(block[0, 0]) + + assert sigma_used(1.0) == pytest.approx(1.0) # nearest real measurement is t = 1 + assert sigma_used(4.0) == pytest.approx(10.0) # nearest real measurement is t = 3 + + +# ================================================================== the design === + +def test_d_optimal_design_for_a_straight_line_measures_the_two_ends(): + """The classical answer: for a line on ``[0, 1]`` with constant noise, the design that + minimizes the volume of the joint confidence region for the intercept and the slope puts its + measurements at the two ends of the interval.""" + grid = np.linspace(0.0, 1.0, 21) + experiment = _line_experiment(grid, measured=[0.4, 0.6]) + free = _free(('a', 'uniform_var', -10.0, 10.0, 1.0), + ('b', 'uniform_var', -10.0, 10.0, 2.0)) + + result, _candidates = _design(experiment, free, points=2, criterion='d') + + assert sorted(m.time for m in result.measurements) == pytest.approx([0.0, 1.0]) + assert result.value > result.baseline_value # log determinant grew + assert improvement(result) > 1.0 + + +def test_a_design_aimed_at_a_decay_rate_measures_one_lifetime(): + """``d S/d k = -t*S0*exp(-k*t)`` peaks at ``t = 1/k``, so the single most informative + measurement of a decay rate is one lifetime after the start. A design aimed at ``k`` picks + the grid time closest to it, and would waste a measurement at either end.""" + k = 0.4 + grid = np.linspace(0.0, 12.0, 49) + experiment = _decay_experiment(grid, measured=[0.5], k=k) + free = _free(('k', 'uniform_var', 1e-3, 5.0, k)) + + result, _candidates = _design(experiment, free, points=1, criterion='a', targets=['k']) + + lifetime = 1.0 / k + assert result.measurements[0].time == pytest.approx( + grid[np.argmin(np.abs(grid - lifetime))]) + + +def test_choosing_the_same_point_twice_means_measuring_it_twice(): + """With one parameter and one clearly best time there is nothing to spread over, so the + design says measure that point repeatedly -- which is the honest recommendation: the + precision of that one measurement is what limits you.""" + grid = np.linspace(0.0, 12.0, 49) + experiment = _decay_experiment(grid, measured=[0.5], k=0.4) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + + result, _candidates = _design(experiment, free, points=3, criterion='a') + + assert len(result.measurements) == 3 + assert len(set(result.measurements)) == 1 + grouped = result.grouped() + assert len(grouped) == 1 and grouped[0][1] == 3 and grouped[0][2] == 1 + + +def test_each_further_measurement_helps_less_than_the_one_before(): + """Information adds up, so variance falls with diminishing returns. The criterion trace + records that, which is how a user sees when a design has stopped paying off.""" + grid = np.linspace(0.0, 12.0, 49) + experiment = _decay_experiment(grid, measured=[0.5], k=0.4) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + + result, _candidates = _design(experiment, free, points=4, criterion='a') + + gains = -np.diff([result.baseline_value] + result.trace) + assert np.all(gains > 0.0) + assert np.all(np.diff(gains) < 0.0) + + +def test_the_designed_information_is_the_information_of_the_designed_dataset(): + """The consistency check that keeps the planned-measurement construction honest: the + information the design reports equals what + :func:`~pybnf.gradient.assemble_fisher_hessian` gives for a dataset that actually holds the + existing measurements plus the recommended ones.""" + grid = np.linspace(0.0, 1.0, 11) + experiment = _line_experiment(grid, measured=[0.4, 0.6]) + free = _free(('a', 'uniform_var', -10.0, 10.0, 1.0), + ('b', 'uniform_var', -10.0, 10.0, 2.0)) + objective = ChiSquareObjective() + + result, _candidates = _design(experiment, free, points=3, criterion='d') + + added = np.array([m.time for m in result.measurements]) + times = np.concatenate([[0.4, 0.6], added]) + combined = _exp(times, {'y': 1.0 + 2.0 * times}, 1.0) + expected = assemble_fisher_hessian( + objective, [(experiment.sim_data, combined, experiment.routing, 'line')], free) + + np.testing.assert_allclose(result.information, expected, rtol=1e-10, atol=1e-12) + + +def test_a_singular_starting_point_is_escaped_before_the_criterion_takes_over(): + """One measurement of a two-parameter line leaves a whole direction unseen, so every + criterion is at its worst value for every candidate and none of them can choose. The + selection notices, picks whatever sees the missing direction, and only then starts optimizing + the requested criterion.""" + grid = np.linspace(0.0, 1.0, 11) + experiment = _line_experiment(grid, measured=[0.5]) + free = _free(('a', 'uniform_var', -10.0, 10.0, 1.0), + ('b', 'uniform_var', -10.0, 10.0, 2.0)) + + result, _candidates = _design(experiment, free, points=2, criterion='d') + + assert is_singular(result.baseline) + assert result.escaped_singular == 1 + assert not is_singular(result.information) + assert np.isfinite(result.value) + assert np.isinf(criterion_value(result.baseline, 'a')) + + +def test_a_design_can_replace_no_answer_with_an_answer(): + """When the existing data leaves a parameter undetermined, the predicted interval goes from + open to finite. There is no ratio to quote for that, which is itself the headline.""" + grid = np.linspace(0.0, 1.0, 11) + experiment = _line_experiment(grid, measured=[0.5]) + free = _free(('a', 'uniform_var', -10.0, 10.0, 1.0), + ('b', 'uniform_var', -10.0, 10.0, 2.0)) + + result, _candidates = _design(experiment, free, points=2, criterion='d') + rows = {row['name']: row + for row in predicted_intervals(result, free, [1.0, 2.0], chi2_quantile_1dof(0.95))} + + assert rows['b']['current'] is None and rows['b']['designed'] is not None + assert rows['b']['width_ratio'] is None + assert improvement(result) is None # a log determinant of -inf has no ratio + + +def test_no_design_can_separate_two_parameters_the_data_only_ever_sees_added_together(): + """``S(t) = S0*exp(-(k1+k2)*t)`` gives the two rates the same sensitivity at every time, so + measuring everything at once still cannot tell them apart. That is structural + non-identifiability, and the run says so rather than recommending measurements that cannot + help.""" + grid = np.linspace(0.0, 10.0, 21) + experiment = _two_channel_experiment(grid, measured=[1.0, 2.0, 5.0]) + free = _free(('k1', 'uniform_var', 1e-3, 5.0, 0.2), + ('k2', 'uniform_var', 1e-3, 5.0, 0.2)) + + with pytest.raises(PybnfError, match='No design over these observables'): + _design(experiment, free, points=3, criterion='a', targets=['k1']) + + +def test_a_design_aimed_elsewhere_is_unaffected_by_an_undetermined_parameter(): + """The refusal is about the parameters the design is *for*. A model with one hopeless + parameter can still be designed for a different one.""" + grid = np.linspace(0.0, 10.0, 21) + times = np.asarray(grid, float) + decay = -times * 100.0 * np.exp(-0.4 * times) + sim = _sim(times, {'S': 100.0 * np.exp(-0.4 * times)}, + {('S', 'k1'): decay, ('S', 'k2'): decay, ('S', 'q'): np.ones(len(times))}) + exp = _exp([1.0, 2.0, 5.0], {'S': [1.0, 2.0, 3.0]}, 1.0) + experiment = DesignExperiment(model='m', suffix='m', sim_data=sim, exp_data=exp, + routing=_routing('k1', 'k2', 'q')) + free = _free(('k1', 'uniform_var', 1e-3, 5.0, 0.2), + ('k2', 'uniform_var', 1e-3, 5.0, 0.2), + ('q', 'uniform_var', -10.0, 10.0, 1.0)) + + result, _candidates = _design(experiment, free, points=2, criterion='a', targets=['q']) + + assert len(result.measurements) == 2 + + +def test_a_log_scaled_parameter_is_designed_for_in_the_scale_it_is_fitted_on(): + """A log-scaled parameter's information carries the ``d theta/d u`` factor, so its predicted + interval is a factor above and below the fitted value rather than a symmetric window. The + half-width in sampling space is what the criterion sees, and it is what the report converts + back through the parameter's own scale.""" + k = 0.4 + grid = np.linspace(0.0, 12.0, 25) + experiment = _decay_experiment(grid, measured=[1.0, 2.5, 5.0], k=k) + free = _free(('k', 'loguniform_var', 1e-3, 5.0, k)) + + result, _candidates = _design(experiment, free, points=2, criterion='a') + threshold = chi2_quantile_1dof(0.95) + row = predicted_intervals(result, free, [np.log10(k)], threshold)[0] + + half = float(interval_half_widths(result.information, threshold)[0]) + assert row['designed'] == pytest.approx((k * 10 ** -half, k * 10 ** half)) + assert row['width_ratio'] < 1.0 + + +def test_the_predicted_interval_is_the_profile_interval_of_a_linear_problem(): + """For a linear model with Gaussian noise the profile is an exact parabola, so the interval + predicted from the information matrix is not an approximation at all: it is the same + ``theta* +- sqrt(threshold * (A^T A)^-1_kk)`` that ``tests/test_profile_likelihood.py`` + checks its own confidence intervals against.""" + grid = np.linspace(0.0, 1.0, 11) + measured = [0.0, 0.2, 0.5, 0.8, 1.0] + experiment = _line_experiment(grid, measured=measured) + free = _free(('a', 'uniform_var', -10.0, 10.0, 1.0), + ('b', 'uniform_var', -10.0, 10.0, 2.0)) + threshold = chi2_quantile_1dof(0.95) + + information = baseline_information(ChiSquareObjective(), [experiment], free) + design_matrix = np.column_stack([np.ones(len(measured)), measured]) + np.testing.assert_allclose(information, design_matrix.T @ design_matrix) + + covariance = np.linalg.inv(design_matrix.T @ design_matrix) + np.testing.assert_allclose(interval_half_widths(information, threshold), + np.sqrt(threshold * np.diag(covariance))) + + +def test_an_estimated_noise_scale_is_designed_for_like_any_other_parameter(): + """When the noise scale is fitted rather than read from the data it is a free parameter with + its own information, so a design accounts for it -- and more measurements sharpen it, which is + what the noise block of the Fisher matrix says.""" + grid = np.linspace(0.0, 12.0, 25) + times = np.asarray(grid) + sim = _sim(times, {'S': 100.0 * np.exp(-0.4 * times)}, + {('S', 'k'): -times * 100.0 * np.exp(-0.4 * times)}) + exp = Data.from_columns( + np.column_stack([[1.0, 2.5], 100.0 * np.exp(-0.4 * np.array([1.0, 2.5]))]), + ['time', 'S']) + experiment = DesignExperiment(model='m', suffix='m', sim_data=sim, exp_data=exp, + routing=ExperimentRouting(routes={ + 'k': ParamRoute.single('k', PARAM, 'k', 1.0)})) + objective = LikelihoodObjective(noise=Gaussian(), + sigma_sources={'sigma': FreeParameterSigma('sigma')}) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4), + ('sigma', 'uniform_var', 1e-3, 100.0, 2.0)) + + baseline = baseline_information(objective, [experiment], free) + candidates = candidate_information(objective, [experiment], free) + result = select_design(baseline, candidates, 3, 'a', [0, 1], ['k', 'sigma']) + + # Every measurement adds 2/sigma^2 to the noise scale's own information, whatever time it is + # taken at, so the scale's variance falls by exactly that much per point. + assert result.information[1, 1] == pytest.approx(baseline[1, 1] + 3 * 2.0 / 2.0 ** 2) + assert parameter_variances(result.information)[0] < parameter_variances(baseline)[0] + + +def test_a_design_spans_several_experiments(): + """Candidates come from every experiment the fit scores, and a recommendation names the one + it belongs to.""" + grid = np.linspace(0.0, 6.0, 13) + slow = _decay_experiment(grid, measured=[1.0], k=0.2) + fast = _decay_experiment(grid, measured=[1.0], k=1.5) + fast = DesignExperiment(model='decay', suffix='fast', sim_data=fast.sim_data, + exp_data=fast.exp_data, routing=fast.routing) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + objective = ChiSquareObjective() + + candidates = candidate_information(objective, [slow, fast], free) + assert {c.experiment for c in candidates.measurements} == {'decay', 'fast'} + assert len(candidates) == 2 * len(grid) + + +# ================================================================ validation === + +def test_targets_must_be_free_parameters_of_this_fit(): + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + with pytest.raises(PybnfError, match='not a free parameter'): + resolve_targets(free, ['nope'], 'a') + + +def test_targets_are_refused_for_a_criterion_that_cannot_use_them(): + """D and E are properties of the whole information matrix, so aiming them at a subset of + parameters would mean something other than what a reader would assume. Say so rather than + ignore the request.""" + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + with pytest.raises(PybnfError, match='only the A-criterion'): + resolve_targets(free, ['k'], 'd') + + +def test_an_unknown_criterion_is_refused_with_the_list_of_real_ones(): + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + with pytest.raises(PybnfError, match='design_criterion must be one of'): + resolve_targets(free, [], 'x') + + +def test_no_targets_means_every_parameter(): + free = _free(('a', 'uniform_var', -1.0, 1.0, 0.0), ('b', 'uniform_var', -1.0, 1.0, 0.0)) + assert resolve_targets(free, None, 'a') == [0, 1] + assert resolve_targets(free, [], 'd') == [0, 1] + + +# ============================================= the grid a design chooses from === + +def _grid(**settings): + """Call the configuration's design-grid helper without building a whole fit.""" + conf = Configuration.__new__(Configuration) + conf.config = settings + return conf + + +def test_without_a_design_grid_an_experiment_simulates_exactly_what_it_measures(): + """The default: no extra simulated times, so nothing about any existing fit changes.""" + points = [0.0, 1.0, 2.0] + assert _grid()._with_design_grid(points) is points + assert _grid(design_grid=0)._with_design_grid(points) is points + + +def test_the_design_grid_adds_times_without_moving_the_measured_ones(): + """The measured times are always kept, so the data still lands on exact grid points and the + scoring is untouched; the extra times are what the design gets to choose from.""" + points = [0.0, 1.0, 2.0] + grid = _grid(design_grid=5, design_t_end=10.0)._with_design_grid(points) + + assert set(points) <= set(grid) + assert grid == pytest.approx([0.0, 1.0, 2.0, 2.5, 5.0, 7.5, 10.0]) + + +def test_the_design_window_ends_at_the_last_measurement_unless_told_otherwise(): + """Without ``design_t_end`` the design may propose new times, but only within the range + already measured -- it is not quietly extrapolating past the data.""" + grid = _grid(design_grid=5)._with_design_grid([0.0, 2.0]) + assert max(grid) == 2.0 + assert len(grid) == 5 + + +# ==================================================================== report === + +def test_the_report_says_what_to_measure_and_what_it_buys(tmp_path): + """Both halves of the report: the measurements to make, and the confidence intervals they + are expected to produce, in the parameters' own units.""" + k = 0.4 + grid = np.linspace(0.0, 12.0, 25) + experiment = _decay_experiment(grid, measured=[0.5], k=k) + free = _free(('k', 'uniform_var', 1e-3, 5.0, k)) + + result, _candidates = _design(experiment, free, points=2, criterion='a', targets=['k']) + path = tmp_path / 'experimental_design.txt' + write_design_report(str(path), result, free, [k], chi2_quantile_1dof(0.95), 0.95) + text = path.read_text() + + assert '# criterion=a' in text + assert '# targets=k' in text + assert 'delta_chi2_threshold' in text + rows = [line for line in text.splitlines() if not line.startswith('#') and line.strip()] + recommended = [r for r in rows if r.split('\t')[1] == 'decay'] + assert recommended, text + assert recommended[0].split('\t')[3] == 'S' + parameters = [r for r in rows if r.split('\t')[0] == 'k'] + assert len(parameters) == 1 + fields = parameters[0].split('\t') + assert float(fields[6]) < 1.0 # the interval narrows + + +def test_the_terminal_summary_groups_repeats_and_quotes_the_intervals(): + from pybnf.design import format_design_summary + + grid = np.linspace(0.0, 12.0, 25) + experiment = _decay_experiment(grid, measured=[0.5], k=0.4) + free = _free(('k', 'uniform_var', 1e-3, 5.0, 0.4)) + result, _candidates = _design(experiment, free, points=2, criterion='a', targets=['k']) + + lines = format_design_summary(result, free, [0.4], chi2_quantile_1dof(0.95)) + assert any('measure 2 times' in line for line in lines) + assert any('times as wide' in line for line in lines) + + +# ===================================================== end to end (real bngsim) === +# The two job surfaces driven through the real sensitivity path: job_type = design on its +# own, and a profile-likelihood run that ends by recommending what to measure next. + +@pytest.mark.bngsim +@pytest.mark.recovery +@pytest.mark.skipif(not BNGSIM_HAS_OUTPUT_SENS, + reason='needs a bngsim build with the output_sensitivities feature') +def test_design_job_recommends_measurements_for_a_decay_model(tmp_path, monkeypatch): + """``job_type = design`` end to end: simulate the supplied best fit once through the real + bngsim sensitivity path, score every time on the simulated grid, and write the report. + + The data stops early, at a fifth of one lifetime, which is why ``k`` is poorly determined: + over that window the curve is nearly a straight line and the decay rate barely shows. With + ``design_grid`` opening the window out to 12, the design has to look past the measured + range -- and it does. Everything it recommends is later than every measurement in hand, and + it lands near one lifetime, ``1/k``, where the sensitivity to the decay rate peaks.""" + H.require_bng2pl() + H.install(monkeypatch) + model = _decay_model(tmp_path) + exp = _write_early_decay_exp(tmp_path / 'decay.exp') + + lines = [ + f'model: {model}', + 'edition = 2', 'job_type = design', 'objective = chi_sq', + f'output_dir = {tmp_path / "out"}', + 'bngl_backend = bngsim', 'initialization = lh', 'delete_old_files = 1', + 'verbosity = 0', 'wall_time_sim = 0', 'random_seed = 1234', + 'design_points = 3', 'design_criterion = a', 'design_target = k', + 'design_grid = 48', 'design_t_end = 12', + f'parameter: k, lower: 0.01, upper: 3.0, initial_value: {TRUE_K}', + f'parameter: S0, lower: 20.0, upper: 400.0, initial_value: {TRUE_S0}', + f'experiment: decay, data: {exp}', + ] + conf = Configuration(ploop('\n'.join(lines).splitlines(keepends=True))) + alg = H.build(conf, 'design') + H.drive(alg) + + result = alg.design_result + assert len(result.measurements) == 3 + assert {m.observable for m in result.measurements} == {'Stot'} + assert all(m.time > 0.5 for m in result.measurements) # past the measured window + lifetime = 1.0 / TRUE_K + assert all(abs(m.time - lifetime) < 1.0 for m in result.measurements) + assert result.value < result.baseline_value # the variance of k falls + + report = Path(conf.config['output_dir']) / 'Results' / 'experimental_design.txt' + assert report.is_file() + text = report.read_text() + assert '# targets=k' in text and 'improvement_factor' in text + + +@pytest.mark.bngsim +@pytest.mark.recovery +@pytest.mark.skipif(not BNGSIM_HAS_OUTPUT_SENS, + reason='needs a bngsim build with the output_sensitivities feature') +def test_design_job_refuses_a_configuration_that_supplies_no_best_fit(tmp_path, monkeypatch): + """A design is computed at a fitted point, so a run that does not supply one is refused with + a message that says what to do about it, rather than quietly designing around the middle of + the parameter box.""" + H.require_bng2pl() + H.install(monkeypatch) + model = _decay_model(tmp_path) + exp = _write_early_decay_exp(tmp_path / 'decay.exp') + conf = H.make_newera_config( + tmp_path, model, exp, + {'k': ('uniform_var', 1e-2, 3.0), 'S0': ('uniform_var', 20.0, 400.0)}, + 'decay', 'design', objective='chi_sq') + + with pytest.raises(PybnfError, match='no initial_value'): + H.build(conf, 'design') + + +@pytest.mark.bngsim +@pytest.mark.recovery +@pytest.mark.skipif(not BNGSIM_HAS_OUTPUT_SENS, + reason='needs a bngsim build with the output_sensitivities feature') +def test_profile_likelihood_ends_by_saying_what_to_measure_next(tmp_path, monkeypatch): + """The whole point of #574, end to end: a profile-likelihood run that finds a parameter hard + to determine goes on to recommend the measurements that would fix it. + + With ``profile_likelihood_design = 1`` the run writes its usual profiles and then, without + being told which parameter to care about, aims the design at the ones it just flagged. The + predicted interval for the flagged parameter narrows, which is the claim the recommendation + is making.""" + H.require_bng2pl() + H.install(monkeypatch) + model = _decay_model(tmp_path) + exp = _write_early_decay_exp(tmp_path / 'decay.exp') + + lines = [ + f'model: {model}', + 'edition = 2', 'job_type = profile_likelihood', 'objective = chi_sq', + f'output_dir = {tmp_path / "out"}', + 'bngl_backend = bngsim', 'initialization = lh', 'delete_old_files = 1', + 'verbosity = 0', 'wall_time_sim = 0', 'random_seed = 1234', + 'population_size = 1', 'max_iterations = 100', + 'profile_likelihood_confidence = 0.95', 'profile_likelihood_step = 0.05', + 'profile_likelihood_max_points = 12', + 'profile_likelihood_design = 1', 'design_points = 3', + f'parameter: k, lower: 0.01, upper: 3.0, initial_value: {TRUE_K}', + f'parameter: S0, lower: 20.0, upper: 400.0, initial_value: {TRUE_S0}', + f'experiment: decay, data: {exp}', + ] + conf = Configuration(ploop('\n'.join(lines).splitlines(keepends=True))) + alg = H.build(conf, 'profile_likelihood') + H.drive(alg) + + assert alg.profile_summary is not None # the profiles still ran and were written + result = alg.design_result + assert result is not None and len(result.measurements) == 3 + + flagged = [s['name'] for s in alg.profile_summary + if s['classification'] == 'practically non-identifiable'] + assert result.target_names == flagged or result.target_names == ['k', 'S0'] + + rows = {row['name']: row for row in predicted_intervals( + result, alg.variables, alg._u_star, alg.threshold)} + for name in result.target_names: + assert rows[name]['width_ratio'] < 1.0 + + results_dir = Path(conf.config['output_dir']) / 'Results' + assert (results_dir / 'profile_likelihood_summary.txt').is_file() + assert (results_dir / 'experimental_design.txt').is_file() + + +def _profile_k(tmp_path, exp, tag, max_points=14): + """Profile the decay model against ``exp`` around the true parameters, and return the summary + keyed by parameter name. The optimum is supplied, so the run profiles without re-fitting and + both halves of the comparison below are centred on the same point.""" + lines = [ + f'model: {_decay_model(tmp_path / tag)}', + 'edition = 2', 'job_type = profile_likelihood', 'objective = chi_sq', + f'output_dir = {tmp_path / tag / "out"}', + 'bngl_backend = bngsim', 'initialization = lh', 'delete_old_files = 1', + 'verbosity = 0', 'wall_time_sim = 0', 'random_seed = 1234', + 'population_size = 1', 'max_iterations = 100', + 'profile_likelihood_confidence = 0.95', 'profile_likelihood_step = 0.02', + f'profile_likelihood_max_points = {max_points}', + f'parameter: k, lower: 0.01, upper: 3.0, initial_value: {TRUE_K}', + f'parameter: S0, lower: 20.0, upper: 400.0, initial_value: {TRUE_S0}', + f'experiment: decay, data: {exp}', + ] + conf = Configuration(ploop('\n'.join(lines).splitlines(keepends=True))) + alg = H.build(conf, 'profile_likelihood') + H.drive(alg) + return {s['name']: s for s in alg.profile_summary} + + +@pytest.mark.bngsim +@pytest.mark.recovery +@pytest.mark.skipif(not BNGSIM_HAS_OUTPUT_SENS, + reason='needs a bngsim build with the output_sensitivities feature') +def test_the_recommended_measurements_really_do_narrow_the_profile(tmp_path, monkeypatch): + """The claim a design makes, checked by making the measurements. + + Profile ``k`` against the early-window data; run a design aimed at ``k``; generate the + measurements it recommends from the model at the same parameters; profile again with them + added. The confidence interval the second run traces has to be genuinely narrower, and close + to the width the design predicted -- which is the part that makes a recommendation worth + acting on rather than merely plausible. + + Nothing here reuses the design's own arithmetic. The second interval comes from + re-optimizing the model at every grid point against a larger dataset, which is a different + computation with a different code path, so agreeing with the prediction means something. + """ + H.require_bng2pl() + H.install(monkeypatch) + (tmp_path / 'before').mkdir() + (tmp_path / 'after').mkdir() + (tmp_path / 'plan').mkdir() + exp = _write_early_decay_exp(tmp_path / 'decay.exp') + + before = _profile_k(tmp_path, exp, 'before') + + lines = [ + f'model: {_decay_model(tmp_path / "plan")}', + 'edition = 2', 'job_type = design', 'objective = chi_sq', + f'output_dir = {tmp_path / "plan" / "out"}', + 'bngl_backend = bngsim', 'initialization = lh', 'delete_old_files = 1', + 'verbosity = 0', 'wall_time_sim = 0', 'random_seed = 1234', + 'design_points = 5', 'design_criterion = a', 'design_target = k', + 'design_grid = 48', 'design_t_end = 12', + f'parameter: k, lower: 0.01, upper: 3.0, initial_value: {TRUE_K}', + f'parameter: S0, lower: 20.0, upper: 400.0, initial_value: {TRUE_S0}', + f'experiment: decay, data: {exp}', + ] + conf = Configuration(ploop('\n'.join(lines).splitlines(keepends=True))) + planner = H.build(conf, 'design') + H.drive(planner) + result = planner.design_result + predicted = {row['name']: row for row in predicted_intervals( + result, planner.variables, [TRUE_K, TRUE_S0], planner.threshold)}['k']['width_ratio'] + + # Make the recommended measurements: the model's own values at the same parameters, which is + # what the design assumed when it scored them. + added = ['%.12g\t%.12g\t%.12g' % (m.time, TRUE_S0 * np.exp(-TRUE_K * m.time), 2.0) + for m in result.measurements] + augmented = tmp_path / 'decay_augmented.exp' + augmented.write_text(Path(exp).read_text() + '\n'.join(added) + '\n') + + after = _profile_k(tmp_path, augmented, 'after') + + def half_width(summary): + assert summary['ci_low'] is not None and summary['ci_high'] is not None + return 0.5 * (summary['ci_high'] - summary['ci_low']) + + observed = half_width(after['k']) / half_width(before['k']) + assert observed < 0.5, (before['k'], after['k']) + # The two agree to about a tenth of a percent on this model; the tolerance is loose enough + # to survive solver noise and still tight enough that a wrong prediction would fail it. + assert observed == pytest.approx(predicted, rel=0.1), (observed, predicted) diff --git a/tests/test_profile_likelihood.py b/tests/test_profile_likelihood.py index 7225c9e8..c9da1eca 100644 --- a/tests/test_profile_likelihood.py +++ b/tests/test_profile_likelihood.py @@ -29,17 +29,16 @@ ProfileLikelihoodAlgorithm, _FLAT_DCHI2, _ProfileTrack, - _chi2_quantile_1dof, _classify, _coverage_notes, _extract_ci, - _norm_ppf, _render_profile_plots, _resolve_profile_idxs, ) from pybnf.config import Configuration from pybnf.gradient import GradientResult from pybnf.parse import ploop +from pybnf.quantiles import chi2_quantile_1dof as _chi2_quantile_1dof, normal_quantile as _norm_ppf from . import recovery_harness as H # The Becker EpoR fast-2p fixtures live beside the gradient smoke tests; the profile-likelihood @@ -602,6 +601,8 @@ def __init__(self, A, y, theta_star, f_min, lower, upper, names, res_dir, self._track_queue = [] self._active_tracks = {} self.profile_summary = None + self.design_result = None + self.design_report = False # the #574 design report is off unless asked for self.phase = 'profile' def __setstate__(self, state): diff --git a/tests/test_pybnf_main_helpers.py b/tests/test_pybnf_main_helpers.py index 83dbbb5c..831d58f0 100644 --- a/tests/test_pybnf_main_helpers.py +++ b/tests/test_pybnf_main_helpers.py @@ -61,6 +61,7 @@ def _config_with_fit_type(fit_type): ('p_dream', 'PDreamAlgorithm'), ('hmc', 'HMCSampler'), ('check', 'ModelCheck'), + ('design', 'ExperimentalDesignAlgorithm'), ] @@ -97,6 +98,10 @@ def test_families_partition_the_codes(): assert {c for c, f in fam.items() if f == 'optimizer'} == {'pso', 'de', 'ade', 'ss', 'sim', 'sa', 'powell', 'cmaes', 'trf', 'lbfgs', 'gntr', 'ms', 'profile_likelihood'} assert {c for c, f in fam.items() if f == 'sampler'} == {'mh', 'pt', 'am', 'dream', 'p_dream', 'hmc'} assert {c for c, f in fam.items() if f == 'checker'} == {'check'} + # Experimental design (#574) fits nothing, so it is neither an optimizer nor a sampler. Its + # own family also keeps it out of what a PEtab job_type = all import emits, which is the one + # thing the family field is read for. + assert {c for c, f in fam.items() if f == 'analysis'} == {'design'} def test_refiners_are_the_start_point_optimizers():