python/dpt/validation/identifiability.py

1"""Bounded offline sensitivity diagnostics with explicit parameter and noise units.23The caller supplies a small Jacobian from a declared diagnostic experiment. This4module does not download images or retain Jacobians inside a production solve.5NumPy, supplied by the optional scientific environment, is loaded only on use.6"""78from __future__ import annotations910import importlib11import math12from collections.abc import Sequence13from dataclasses import dataclass14from typing import Any1516from dpt.contracts import ContractError, NumericalError, finite_scalar171819@dataclass(frozen=True, slots=True)20class SensitivitySpectrum:21    singular_values: tuple[float, ...]22    right_directions: tuple[tuple[float, ...], ...]23    rank: int24    threshold: float25    condition_number: float | None26    observations: int27    parameters: int282930def _scaled_entry(value: float, unit: float, weight: float) -> float:31    value = finite_scalar(value, "Jacobian entry")32    root_weight = math.sqrt(finite_scalar(weight, "precision weight", minimum=0))33    if value == 0 or root_weight == 0:34        return 0.035    mantissa, exponent = 1.0, 036    for factor in (value, unit, root_weight):37        coefficient, power = math.frexp(factor)38        mantissa *= coefficient39        exponent += power40    try:41        return math.ldexp(mantissa, exponent)42    except OverflowError as error:43        raise NumericalError("scaled sensitivity exceeds binary64 range") from error444546# region book:scaled-sensitivity-diagnostics47def sensitivity_spectrum(48    jacobian_rows: Sequence[Sequence[float]],49    parameter_scales: Sequence[float],50    *,51    precision_weights: Sequence[float] | None = None,52    relative_threshold: float = 1e-8,53    absolute_threshold: float = 0.0,54) -> SensitivitySpectrum:55    """SVD of sqrt(W) J S; right directions use dimensionless chart coordinates.5657    Singular values and an explicitly chosen threshold describe local sensitivity,58    not global identifiability or clinical recovery. Duplicate views add no new59    independent directions, although weighting can change threshold-defined rank.60    Noise correlations require an explicitly prewhitened Jacobian; diagonal61    precision weights must not stand in for an unknown covariance model.62    """63    rows, columns = len(jacobian_rows), len(parameter_scales)64    if not 1 <= rows <= 4096 or not 1 <= columns <= 32:65        raise ContractError("offline diagnostic is bounded to 4096 observations and 32 parameters")66    if any(len(row) != columns for row in jacobian_rows):67        raise ContractError("every Jacobian row must contain each active parameter")68    units = tuple(finite_scalar(value, "parameter scale", minimum=0) for value in parameter_scales)69    if min(units) == 0:70        raise ContractError("parameter scales must be positive")71    relative = finite_scalar(relative_threshold, "relative threshold", minimum=0)72    absolute = finite_scalar(absolute_threshold, "absolute threshold", minimum=0)73    if relative >= 1:74        raise ContractError("relative rank threshold must be less than one")75    weights = tuple(precision_weights) if precision_weights is not None else (1.0,) * rows76    if len(weights) != rows:77        raise ContractError("one fixed precision weight is required per observation")78    scaled = [79        [_scaled_entry(value, units[column], weights[index]) for column, value in enumerate(row)]80        for index, row in enumerate(jacobian_rows)81    ]82    if not all(math.isfinite(value) for row in scaled for value in row):83        raise NumericalError("scaled sensitivity exceeds binary64 range")84    # Zero rows preserve the missing directions when observations < parameters,85    # without requesting the large square left-singular-vector matrix.86    scaled.extend([[0.0] * columns for _ in range(max(0, columns - rows))])87    np: Any = importlib.import_module("numpy")88    _, values, right = np.linalg.svd(np.asarray(scaled, dtype=np.float64), full_matrices=False)89    singular = tuple(float(value) for value in values)90    if not all(map(math.isfinite, singular)):91        raise NumericalError("sensitivity decomposition produced nonfinite singular values")92    threshold = max(absolute, relative * singular[0])93    rank = sum(value > threshold for value in singular)94    condition = singular[0] / singular[-1] if rank == columns else None95    if condition is not None and not math.isfinite(condition):96        condition = None97    return SensitivitySpectrum(98        singular,99        tuple(tuple(float(value) for value in row) for row in right),100        rank,101        threshold,102        condition,103        rows,104        columns,105    )106107108# endregion book:scaled-sensitivity-diagnostics109