Source code for pybnf.noise.source
"""The ``SigmaSource`` abstraction (ADR-0021): where a noise model's noise
parameter comes from.
ADR-0011 made the ``NoiseModel`` a pure per-point kernel and named the *noise
parameter source* as the objective wrapper's job, but left it hard-coded -- the
``_SD`` data column, the magic free parameters ``sigma__FREE`` / ``r__FREE``, and
the ``neg_bin_r`` constant each lived in a different objfunc subclass. This module
lifts those three into one first-class abstraction so they can be selected **per
observable** (the #410 engine) and named freely (dissolving the magic strings, as
M2.3 did for priors).
A ``SigmaSource`` answers two questions the per-point engine needs: its
``value(...)`` for one observation, and whether it is ``estimated``. The
**estimated** flag is load-bearing -- it is what decides whether the family's
likelihood normalizer is summed: a fixed source (a data column, a constant)
contributes only the family's ``data_fit``, while an estimated source (a free
parameter) contributes the full ``nll`` including the normalizer. This is exactly
ADR-0011's "normalizer retained iff the noise parameter is estimated", now keyed
off the source rather than hard-coded per objfunc.
"""
import re
from abc import ABC, abstractmethod
import numpy as np
from ..printing import PybnfError
# A PEtab per-measurement placeholder symbol (``observableParameter1_*`` /
# ``noiseParameter1_*``): in a row-varying noiseFormula it stands for the *row's* token
# (a number or an estimated id), read from the binding table rather than substituted once
# (ADR-0045). Mirrors ``petab.formula._PLACEHOLDER`` (kept local so the noise engine carries
# no petab dependency on this path).
_PLACEHOLDER = re.compile(r'(?:observable|noise)Parameter\d')
[docs]
class SigmaSource(ABC):
"""Where a noise model reads its noise parameter for one observation.
``estimated`` decides normalizer inclusion (see the module docstring); it is
``False`` by default and overridden to ``True`` by the free-parameter source.
"""
estimated = False
[docs]
@abstractmethod
def value(self, owner, exp_data, exp_row, col_name):
"""The noise-parameter value for one observation. ``owner`` is the objective
(used by the free-parameter source to read the resolved pset values);
``exp_data``/``exp_row``/``col_name`` locate the point (used by the data-
column source)."""
[docs]
def required_free_param(self):
"""The ``__FREE`` parameter name this source requires the fit to declare,
or ``None`` if it sources no free parameter (used by ``_load_variables``
validation, ADR-0021)."""
return None
[docs]
def required_free_params(self):
"""The **set** of free-parameter names this source requires the fit to
declare. Defaults to the single :meth:`required_free_param` (or empty);
:class:`FormulaSigma` overrides it because an expression sigma references
several (ADR-0044). The objective unions these across its sources
(``required_free_noise_params``)."""
name = self.required_free_param()
return {name} if name is not None else set()
[docs]
def exp_column(self, col_name):
"""The experimental-data column this source consumes, or ``None`` if it
reads no data column -- so ``_check_columns`` can exempt it from the
unused-column error."""
return None
[docs]
class DataColumnSigma(SigmaSource):
"""The noise parameter read per point from an experimental-data column named
``<observable><suffix>`` (``chi_sq`` / ``lognormal``: the ``_SD`` column). It is
a *fixed* source -- the value is data, not estimated -- so the caller drops the
likelihood normalizer. The suffix is explicit (default ``_SD``) so a non-Gaussian
family can read a differently-named column without the "standard deviation"
misnomer (ADR-0021)."""
estimated = False
def __init__(self, suffix='_SD'):
self.suffix = suffix
[docs]
def value(self, owner, exp_data, exp_row, col_name):
column = col_name + self.suffix
try:
idx = exp_data.cols[column]
except KeyError:
# Todo: Check for this and throw the error before all the workers get created.
raise PybnfError(f'Column {column} not found',
f"Column {column} was not found in the experimental data. A noise model that reads its "
f"scale from the data requires a {column} column corresponding to {col_name}, giving the "
"per-point noise scale (e.g. the standard deviation) of that variable. ")
return exp_data.data[exp_row, idx]
[docs]
class FreeParameterSigma(SigmaSource):
"""The noise parameter estimated as a free parameter, resolved by name from the
pset (``chi_sq_dynamic``'s ``sigma__FREE``, ``neg_bin_dynamic``'s ``r__FREE``,
or a per-observable ``__FREE`` parameter). It is *estimated*, so the caller
keeps the likelihood normalizer."""
estimated = True
def __init__(self, name):
self.name = name
[docs]
def value(self, owner, exp_data, exp_row, col_name):
# ``owner._pset_values`` is the {name: value} map the objective's
# evaluate_multiple builds once per evaluation from the pset (ADR-0021).
return owner._pset_values[self.name]
[docs]
class ConstantSigma(SigmaSource):
"""The noise parameter held at a fixed configuration constant (``neg_bin``'s
``neg_bin_r``; the native ``fix_at`` source). Fixed, so the caller drops the
likelihood normalizer."""
estimated = False
def __init__(self, value):
self.const = value
[docs]
class RelativeSigma(SigmaSource):
"""The noise parameter proportional to the measurement: constant-coefficient-of-
variation noise, ``sigma = cv * |observation|`` (the native ``relative`` source,
ADR-0031). This is the honest heteroscedastic noise model the legacy ``norm_sos``
objfunc fits -- a Gaussian whose standard deviation scales with the data, so the
squared residual is normalized per point by the measurement (``cv = 1`` reproduces
``((sim - exp) / exp)**2`` up to the proper ``1/2``). The ``cv`` is a fixed
constant, so this source is *fixed* and the caller drops the likelihood
normalizer."""
estimated = False
def __init__(self, cv=1.0):
self.cv = cv
[docs]
def value(self, owner, exp_data, exp_row, col_name):
observation = exp_data.data[exp_row, exp_data.cols[col_name]]
# abs() keeps sigma positive for a negative measurement; the Gaussian uses
# sigma**2, so the sign never mattered (norm_sos squared the raw ratio).
return self.cv * abs(observation)
[docs]
class FormulaSigma(SigmaSource):
"""The noise parameter given by an **expression** over free parameters (+ constants),
evaluated against the current PSet at each point (the native ``formula`` source,
ADR-0044).
PEtab lets the ``noiseFormula`` be an arithmetic expression -- e.g. ``0.1 +
0.05*scaling`` after a per-observable ``observableParameter`` placeholder is substituted
in -- a per-observable sigma that is neither a data column (``DataColumnSigma``) nor a
single free parameter (``FreeParameterSigma``). This source closes that gap: it holds a
PEtab-math expression string and, lazily on first use, compiles it to a vectorized
``numpy`` callable over its free symbols (the same compile-once-per-worker pattern as
:class:`~pybnf.measurement.MeasurementModel`, ADR-0036 -- a ``lambdify``\\ d callable is
not picklable, so it is excluded from ``__getstate__`` and rebuilt worker-side).
It is *estimated* (it reads estimated parameters), so the caller keeps the family's
likelihood normalizer (ADR-0011). Its free symbols are the nuisance free parameters the
fit must declare (:meth:`required_free_params`); they resolve from ``owner._pset_values``
exactly as :class:`FreeParameterSigma`'s single name does."""
estimated = True
def __init__(self, formula, names=None):
#: The PEtab-math expression (over free-parameter ids + numeric constants).
self.formula = formula
#: The ordered free-symbol names (PSet keys at eval time); ``None`` until first
#: resolved (a parse), then the compiler's canonical sorted order. Picklable.
self.names = list(names) if names is not None else None
self._func = None # lambdify callable; not pickled (rebuilt lazily worker-side)
def __getstate__(self):
state = self.__dict__.copy()
state['_func'] = None
return state
def _ensure_names(self):
if self.names is None:
from ..petab.formula import formula_free_symbols
self.names = formula_free_symbols(self.formula)
return self.names
def _callable(self):
if self._func is None:
from ..petab.formula import compile_petab_formula
# allowed_symbols = the formula's own free symbols, so validation is a no-op
# here (the symbols were validated as declared free parameters at config load);
# adopt the compiler's canonical name ordering for positional calling.
func, names = compile_petab_formula(self.formula, self._ensure_names())
self._func, self.names = func, names
return self._func, self.names
[docs]
def value(self, owner, exp_data, exp_row, col_name):
func, names = self._callable()
return float(func(*[owner._pset_values[n] for n in names]))
[docs]
class PerMeasurementFormulaSigma(SigmaSource):
"""The noise parameter given by an expression that references a **row-varying** PEtab
per-measurement placeholder, bound per data point from the experimental data's
binding table (the native ``formula`` source over a placeholder, ADR-0045).
ADR-0044's :class:`FormulaSigma` covers a noiseFormula whose placeholder is *constant
across an observable's rows* (substituted away to a per-observable σ). When the
placeholder's ``noiseParameters`` token instead **differs row to row** -- a per-condition
or per-timepoint estimated σ -- it cannot be substituted to one symbol; it must be bound
at scoring time from the row's token. This source keeps the placeholder as a free symbol
of the (cached, lambdified) expression and, at ``value(owner, exp_data, exp_row,
col_name)``, reads ``exp_data.measurement_params[col_name][placeholder][exp_row]`` for the
row's token -- a **number** binds as itself, a **parameter id** resolves from
``owner._pset_values`` (a per-row estimated nuisance, ADR-0034). Any non-placeholder free
symbol resolves from the PSet exactly as :class:`FormulaSigma`'s do.
It is *estimated* (a per-row token id is an estimated parameter), lazy-compiled, and not
pickled (the same compile-once-per-worker pattern as :class:`FormulaSigma`). The binding
table lives on the experimental :class:`~pybnf.data.Data` (carried there by ``config.py``
from the importer's sidecar), not on the source, so the source survives dask scatter
independently of the data."""
estimated = True
def __init__(self, formula, names=None):
#: The PEtab-math expression, with its per-measurement placeholder(s) kept as symbols.
self.formula = formula
#: The ordered free-symbol names (placeholders + fixed PSet names); ``None`` until a
#: parse resolves them, then the compiler's canonical sorted order. Picklable.
self.names = list(names) if names is not None else None
self._func = None # lambdify callable; not pickled (rebuilt lazily worker-side)
def __getstate__(self):
state = self.__dict__.copy()
state['_func'] = None
return state
def _ensure_names(self):
if self.names is None:
from ..petab.formula import formula_free_symbols
self.names = formula_free_symbols(self.formula)
return self.names
def _callable(self):
if self._func is None:
from ..petab.formula import compile_petab_formula
func, names = compile_petab_formula(self.formula, self._ensure_names())
self._func, self.names = func, names
return self._func, self.names
[docs]
def required_free_params(self):
# Only the FIXED (non-placeholder) symbols are always-required PSet names; the
# per-row placeholder tokens are validated against the binding table by config.py
# (they vary by row, so they are not a fixed required-name set here).
return {n for n in self._ensure_names() if not _PLACEHOLDER.match(n)}
[docs]
def value(self, owner, exp_data, exp_row, col_name):
func, names = self._callable()
bindings = self._row_bindings(exp_data, col_name)
args = []
for name in names:
if _PLACEHOLDER.match(name):
args.append(self._resolve_token(owner, bindings, name, exp_row, col_name))
else:
args.append(owner._pset_values[name])
return float(func(*args))
@staticmethod
def _row_bindings(exp_data, col_name):
table = getattr(exp_data, 'measurement_params', None)
if not table or col_name not in table:
raise PybnfError(
f"A row-varying noise formula for '{col_name}' needs a per-measurement "
f"binding table, but none is attached to the experimental data. The "
f"experiment must declare 'measurement_params: <file>.tsv' (ADR-0045).")
return table[col_name]
@staticmethod
def _resolve_token(owner, bindings, placeholder, exp_row, col_name):
try:
token = bindings[placeholder][exp_row]
except (KeyError, IndexError):
raise PybnfError(
f"The per-measurement binding table for '{col_name}' has no value for "
f"placeholder '{placeholder}' at data row {exp_row} (ADR-0045).")
try:
return float(token) # a per-row numeric token, inlined
except (TypeError, ValueError):
return owner._pset_values[token] # a per-row parameter id, from the PSet
[docs]
class ColumnMeanSigma(SigmaSource):
"""The noise parameter set to the observable's experimental column mean -- one
scalar shared by every point of the column (the native ``column_mean`` source,
ADR-0031). This is the heteroscedastic-across-observables model the legacy
``ave_norm_sos`` objfunc fits -- each variable's residuals are normalized by that
variable's own scale (its mean), so a large-magnitude observable does not dominate
a small one (``sigma = ybar`` reproduces ``((sim - exp) / ybar)**2`` up to the
proper ``1/2``). Fixed (it is data, not estimated), so the caller drops the
likelihood normalizer."""
estimated = False
[docs]
def value(self, owner, exp_data, exp_row, col_name):
# The mean is a per-column constant; recomputing it per point is O(1) amortized
# over the small data PyBNF fits and keeps the source stateless (no per-exp_data
# cache to invalidate across models/suffixes). The legacy ave_norm_sos
# precomputes it once in evaluate(); the result is identical.
return np.average(exp_data[col_name])