from __future__ import annotations
from typing import Any, Dict, List, Tuple
import numpy as np
import pandas as pd
from pydantic import field_validator
from phenotypic.analysis.abc_._linear_softplus_base import (
_DEFAULT_BETA,
_LinearSoftplusBase,
)
from phenotypic.schema import (
LINEAR_CAP_AND_LAG_MODEL,
MODEL_METRICS,
)
_MODE_FIXED_BETA = "fixed_beta"
_MODE_FITTED_BETA = "fitted_beta"
[docs]
class LinearCapAndLagModel(_LinearSoftplusBase):
r"""Linear-softplus growth fitter with a softplus saturation ceiling.
Fits a linear post-lag growth phase with a softplus lag transition
*and* a softplus saturation ceiling:
.. math::
s_{\text{unclamped}}(t) =
\frac{v}{\alpha}\, \ln\!\bigl(1 + e^{\alpha(t-\lambda)}\bigr) + s_0
.. math::
s(t) = s_{\max}
- \frac{1}{\beta}\,\ln\!\bigl(1 + e^{\beta(s_{\max} - s_{\text{unclamped}}(t))}\bigr)
Use this class when colonies show a clear carrying-capacity plateau
in the observation window. For pre-saturation linear growth, use
:class:`LinearLagModel` instead.
Per-group mode dispatch:
The fit picks one of two variants per fit group, recorded in
``LINEAR_CAP_AND_LAG_MODEL.mode``:
- ``"fitted_beta"`` — 5-parameter fit. Triggered when ``beta``
is ``None`` *and* a saturation shoulder is detected in the
group (smoothed tail slope flattens below
``shoulder_slope_ratio`` times the peak slope).
- ``"fixed_beta"`` — 4-parameter fit with ``beta`` held
constant. Triggered when the user supplied an explicit
scalar ``beta``, *or* when no shoulder is detected. The
effective ``beta`` is ``self.beta`` when set, else the
module default (``10.0``).
Pruning is intentionally not exposed on this class — the saturation
plateau IS the model, so dropping the tail would defeat the fit.
Attributes:
smax (float | None): Fixed carrying capacity. ``None`` falls back
to the per-group observed max.
beta (float | None): Saturation transition sharpness. ``None``
(default) opts into per-group mode dispatch — fit when a
shoulder is present, otherwise held at the module default.
Set a positive scalar to force ``"fixed_beta"`` mode
unconditionally.
stderr_label (str | None): Column providing per-timepoint standard
errors used as weights. Same semantics as :class:`LinearLagModel`.
s0_prior (bool | float | str | None): Unified Gaussian-prior
source for ``s0``. Same dispatch as :class:`LinearLagModel`.
s0_prior_cv (float | None): CV coefficient for the prior σ
(``σ = cv × µ``). Mutually exclusive with
``s0_prior_sigma``. Defaults to ``None``; if neither knob
is set and the prior is engaged, CV=0.05 is applied.
s0_prior_sigma (float | None): Absolute σ for the prior.
Mutually exclusive with ``s0_prior_cv``.
s0_prior_groupby (List[str] | None): Optional coarser grouping
for empirical-Bayes pooling of the per-group prior ``µ``.
shoulder_slope_ratio (float): Fraction of peak ``ds/dt`` below
which the tail slope counts as a saturation shoulder for
mode dispatch. Defaults to ``0.05``.
Output column naming:
``analyze`` emits metric-qualified columns
``LinearCapAndLagModel_<metric>_<parameter>`` (e.g. fitting
``on="Shape_Area"`` yields ``LinearCapAndLagModel_Area_v``), plus
qualified fit-quality columns ``ModelMetrics_<metric>_<label>``. The
``<metric>`` segment is ``self.on`` with a recognized
measurement-category prefix stripped.
.. note::
**``f_scale`` is unit-sensitive only on the unweighted fit
path.** The inherited ``f_scale`` (see :class:`ModelFitter`) is
the Huber/robust inlier–outlier threshold expressed in *residual
units*, and those units depend on whether the fit is weighted:
- **Weighted** (``stderr_label`` set, or the default
auto-derived replicate SEM when timepoints carry ≥2
replicates): residuals are divided by σ and are therefore
dimensionless, so ``f_scale=1.0`` means "one standard error"
and is invariant to the units of ``on``. No retuning is
needed when the measurement scale changes.
- **Unweighted** (no σ — e.g. single-replicate timepoints):
residuals are in the native units of ``on``, so ``f_scale``
is an absolute size threshold. If those units change (e.g.
radius in px → mm, which shrinks residuals ~50×) ``f_scale``
must be rescaled to match, or the default robust
``loss="huber"`` never reaches its linear arm and silently
degrades to ordinary least squares — losing all outlier
suppression. ``loss="linear"`` ignores ``f_scale`` and is
unaffected.
"""
_measurement_infoclass = LINEAR_CAP_AND_LAG_MODEL
smax: float | None = None
beta: float | None = None
shoulder_slope_ratio: float = 0.05
@field_validator("beta")
@classmethod
def _validate_beta(cls, value: float | None) -> float | None:
"""Reject a non-positive or non-finite scalar ``beta``.
Args:
value: The candidate ``beta``. ``None`` opts into per-group
mode dispatch and is left untouched.
Returns:
The validated ``beta`` value.
Raises:
ValueError: If ``beta`` is not ``None`` and is not a
positive finite number.
"""
if value is not None and (not np.isfinite(value) or value <= 0):
raise ValueError(
f"beta must be None or a positive finite number, "
f"got {value!r}."
)
return value
@field_validator("shoulder_slope_ratio")
@classmethod
def _validate_shoulder_slope_ratio(cls, value: float) -> float:
"""Reject a ``shoulder_slope_ratio`` outside the open ``(0, 1)``.
Args:
value: The candidate ``shoulder_slope_ratio``.
Returns:
The validated ratio.
Raises:
ValueError: If the ratio is not finite or not in ``(0, 1)``.
"""
if not np.isfinite(value) or value <= 0 or value >= 1:
raise ValueError(
f"shoulder_slope_ratio must be in (0, 1), "
f"got {value!r}."
)
return value
# ------------------------------------------------------------------ #
# Model math
# ------------------------------------------------------------------ #
[docs]
@staticmethod
def model_func(
t: np.ndarray | float,
v: float,
s0: float,
lam: float,
alpha: float,
smax: float,
beta: float = _DEFAULT_BETA,
) -> float | np.ndarray:
r"""Linear-softplus growth curve with softplus saturation ceiling.
Args:
t: Time (scalar or array).
v: Post-lag growth rate.
s0: Initial size.
lam: Lag duration.
alpha: Lag transition sharpness.
smax: Carrying capacity (saturation ceiling).
beta: Saturation transition sharpness.
Returns:
Predicted size at ``t``; scalar when ``t`` is scalar,
otherwise an array.
"""
unclamped = _LinearSoftplusBase._lag_softplus(
t=t, v=v, s0=s0, lam=lam, alpha=alpha
)
softplus_sat = np.logaddexp(0.0, beta * (smax - unclamped)) / beta
return smax - softplus_sat
# ------------------------------------------------------------------ #
# Saturation-shoulder detection and per-group mode dispatch
# ------------------------------------------------------------------ #
def _has_saturation_shoulder(self, group: pd.DataFrame) -> bool:
"""Return ``True`` iff ``group`` shows a saturation shoulder.
Compares a robust peak slope (smoothed ``dy/dt`` maximum) to
the minimum of the last few raw gradients. A shoulder is
declared when the curve's terminal slope fell to
``self.shoulder_slope_ratio`` times the peak or below — a
criterion that still fires on curves with a short
post-saturation tail (as few as 2 plateau points), where
tail-mean statistics would be diluted by transitional points.
A dynamic-range guard rejects flat-noise groups that would
otherwise trivially satisfy the ratio test.
"""
if len(group) < 6:
return False
g = group.sort_values(self.time_label).reset_index(drop=True)
y = g[self.on].to_numpy(dtype=float)
t = g[self.time_label].to_numpy(dtype=float)
y_finite = y[np.isfinite(y)]
if y_finite.size == 0:
return False
med_abs = abs(float(np.median(y_finite))) + 1e-9
if float(y_finite.max() - y_finite.min()) <= 0.05 * med_abs:
return False
dy_dt = np.gradient(y, t)
# Peak slope: smoothed max — robust to single-point noise spikes.
window = min(5, max(3, len(y) // 4))
smoothed = np.convolve(dy_dt, np.ones(window) / window, mode="same")
peak_slope = float(np.nanmax(smoothed))
if not np.isfinite(peak_slope) or peak_slope <= 0:
return False
# Terminal slope: min of the last few raw gradients. Using the
# min (not mean) keeps the criterion sensitive to short plateaus
# where the final 2-3 points have collapsed to near-zero slope
# but earlier tail points are still in the transition.
tail_k = max(2, min(3, len(dy_dt) // 5))
tail_min = float(np.nanmin(dy_dt[-tail_k:]))
if not np.isfinite(tail_min):
return False
return tail_min <= self.shoulder_slope_ratio * peak_slope
def _mode_for(self, group: pd.DataFrame) -> str:
"""Pick the fit variant for ``group``.
Two-way dispatch:
1. User-explicit scalar ``beta`` → ``"fixed_beta"``.
2. Shoulder detected (``beta is None``) → ``"fitted_beta"``.
3. Otherwise (``beta is None``, no shoulder) → ``"fixed_beta"``
with the module-default ``β``.
"""
if self.beta is not None:
return _MODE_FIXED_BETA
if self._has_saturation_shoulder(group):
return _MODE_FITTED_BETA
return _MODE_FIXED_BETA
# ------------------------------------------------------------------ #
# Per-group loss kwargs (smax, beta on top of the base σ + prior).
# ------------------------------------------------------------------ #
def _extra_loss_kwargs(self, group: pd.DataFrame) -> Dict[str, Any]:
kw = super()._extra_loss_kwargs(group)
kw["smax"] = self._smax_for(group)
if self._mode_for(group) == _MODE_FIXED_BETA:
kw["beta"] = (
float(self.beta) if self.beta is not None else _DEFAULT_BETA
)
# In fitted-beta mode the 5th optimizer entry supplies beta
# directly; no kwarg needed.
return kw
def _smax_for(self, group: pd.DataFrame) -> float:
if self.smax is not None:
return float(self.smax)
return float(group[self.on].max())
# ------------------------------------------------------------------ #
# Required fit hooks — 4-or-5 vector depending on mode.
# ------------------------------------------------------------------ #
def _initial_guess(self, group: pd.DataFrame) -> list[float]:
"""Heuristic initial guess for the per-group optimizer vector.
Returns ``[v, s0, lam, alpha]`` in fixed-beta mode, or
``[v, s0, lam, alpha, beta]`` in fitted-beta mode. The fifth
entry seeds at the module-default ``beta``.
"""
g = group.sort_values(self.time_label)
y = g[self.on].to_numpy(dtype=float)
t = g[self.time_label].to_numpy(dtype=float)
stats = self._inoc_stats(group)
if stats is not None:
s0_init = stats[0]
else:
s0_init = float(np.median(y[: max(2, len(y) // 4)]))
cut = max(2, int(0.4 * len(y)))
if len(y) - cut >= 2 and np.ptp(t[cut:]) > 0:
slope = float(np.polyfit(t[cut:], y[cut:], 1)[0])
else:
slope = 1.0
v_init = float(np.clip(slope, 1e-4, self._V_UPPER))
y_range = max(y.max() - s0_init, 1e-6)
crossing_mask = y > s0_init + 0.1 * y_range
if crossing_mask.any():
lam_init = float(t[np.argmax(crossing_mask)])
else:
lam_init = float(t[0])
alpha_init = 10.0
guess = [v_init, s0_init, lam_init, alpha_init]
if self._mode_for(group) == _MODE_FITTED_BETA:
guess.append(_DEFAULT_BETA)
return guess
def _bounds(self, group: pd.DataFrame) -> Tuple[List[float], List[float]]:
"""Parameter bounds ``(lower, upper)``.
Four entries in fixed-beta mode; a fifth ``(2.0, 50.0)`` pair
appended in fitted-beta mode, matching the ``alpha`` range
convention (softplus is ~linear below 2 and effectively a
hard step above 50).
"""
t_max = float(group[self.time_label].max())
if t_max <= 0:
t_max = 1.0
stats = self._inoc_stats(group)
if stats is not None:
s0_upper = float(max(3.0 * stats[0], 1e-6))
else:
s0_upper = float(group[self.on].max()) or 1.0
if s0_upper <= 0:
s0_upper = 1.0
lower = [0.0, 0.0, 0.0, 2.0]
upper = [self._V_UPPER, s0_upper, t_max, 50.0]
if self._mode_for(group) == _MODE_FITTED_BETA:
lower.append(2.0)
upper.append(50.0)
return lower, upper
def _unpack_params(
self, x: np.ndarray, group: pd.DataFrame
) -> Dict[Any, Any]:
v, s0, lam, alpha = (float(x[i]) for i in range(4))
smax_val = self._smax_for(group)
# Derive mode from the optimizer-vector length. The vector
# length was set by ``_initial_guess`` (which called
# ``_mode_for``); reading it back here is both faster than
# re-running shoulder detection and structurally consistent
# with whatever mode set the bounds and x0 — no risk of the
# detector flipping between calls if a subclass mutates the
# group in-place.
if len(x) >= 5:
mode = _MODE_FITTED_BETA
beta_val = float(x[4])
else:
mode = _MODE_FIXED_BETA
beta_val = (
float(self.beta) if self.beta is not None else _DEFAULT_BETA
)
return {
LINEAR_CAP_AND_LAG_MODEL.v : v,
LINEAR_CAP_AND_LAG_MODEL.s0 : s0,
LINEAR_CAP_AND_LAG_MODEL.lam : lam,
LINEAR_CAP_AND_LAG_MODEL.alpha: alpha,
LINEAR_CAP_AND_LAG_MODEL.smax : smax_val,
LINEAR_CAP_AND_LAG_MODEL.beta : beta_val,
LINEAR_CAP_AND_LAG_MODEL.mode : mode,
}
def _predict_kwargs(self, row) -> Dict[str, Any]:
smax_val = row[LINEAR_CAP_AND_LAG_MODEL.smax]
beta_val = row[LINEAR_CAP_AND_LAG_MODEL.beta]
return {
"v" : float(row[LINEAR_CAP_AND_LAG_MODEL.v]),
"s0" : float(row[LINEAR_CAP_AND_LAG_MODEL.s0]),
"lam" : float(row[LINEAR_CAP_AND_LAG_MODEL.lam]),
"alpha": float(row[LINEAR_CAP_AND_LAG_MODEL.alpha]),
"smax" : float(smax_val),
"beta" : _DEFAULT_BETA if pd.isna(beta_val) else float(beta_val),
}
def _hover_fields(self) -> List[Tuple[str, Any, str]]:
return [
("v", LINEAR_CAP_AND_LAG_MODEL.v, ".4f"),
("s0", LINEAR_CAP_AND_LAG_MODEL.s0, ".3f"),
("lambda", LINEAR_CAP_AND_LAG_MODEL.lam, ".3f"),
("alpha", LINEAR_CAP_AND_LAG_MODEL.alpha, ".2f"),
("smax", LINEAR_CAP_AND_LAG_MODEL.smax, ".3f"),
("beta", LINEAR_CAP_AND_LAG_MODEL.beta, ".2f"),
("mode", LINEAR_CAP_AND_LAG_MODEL.mode, ""),
("RMSE", MODEL_METRICS.RMSE, ".4f"),
]
LinearCapAndLagModel.__doc__ = LINEAR_CAP_AND_LAG_MODEL.append_rst_to_doc(LinearCapAndLagModel)
