"""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
[docs]
@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
[docs]
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)
[docs]
@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)
[docs]
@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)
[docs]
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)
[docs]
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)
[docs]
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
