{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "daeed3d3",
   "metadata": {},
   "source": [
    "# The Linear-Softplus and Double-Softplus Growth Models\n",
    "\n",
    "The phenotypic library ships two complementary kinetic fits for\n",
    "time-course colony sizing on solid media:\n",
    "\n",
    "- :class:`~phenotypic.analysis.LinearSoftplus` — 4 parameters, no\n",
    "  saturation ceiling. Use when colonies are still in the linear-\n",
    "  growth regime, or when you want the saturation tail discarded as\n",
    "  observation noise. Pruning is ON by default.\n",
    "- :class:`~phenotypic.analysis.DoubleSoftplus` — 4 or 5 parameters,\n",
    "  with a softplus saturation ceiling. Use when the data shows a\n",
    "  clear carrying-capacity plateau and you want to recover ``smax``\n",
    "  (and optionally fit β).\n",
    "\n",
    "Both classes share the same lag-phase math, the same σ-resolution\n",
    "machinery, the same Gaussian-prior-on-``s0`` infrastructure, and the\n",
    "same robust-loss handling. ``DoubleSoftplus`` adds a softplus ceiling\n",
    "on top of the shared core and a per-group binary mode dispatch\n",
    "(``fitted_beta`` vs ``fixed_beta``). This page is a conceptual\n",
    "walkthrough covering:\n",
    "\n",
    "1. The model equations and what each parameter means biologically.\n",
    "2. The two-way mode dispatch on ``DoubleSoftplus`` — how the fitter\n",
    "   picks between a 4-parameter clamped fit (fixed β) and a\n",
    "   5-parameter fit with saturation sharpness as a free parameter.\n",
    "3. Every outlier-robustness mitigation wired into the fits:\n",
    "   saturation pruning (``LinearSoftplus`` only), robust loss,\n",
    "   inverse-variance weighting, the replicate-SEM quantile floor,\n",
    "   the $n=1$ timepoint pooled-std fallback, and the optional\n",
    "   Gaussian prior on the inoculum size $s_0$.\n",
    "\n",
    "The growth-curve literature has pulled in this direction for years —\n",
    "robust M-estimators for outliers\n",
    "([Huber, 1964](https://doi.org/10.1214/aoms/1177703732)),\n",
    "inverse-variance weighting for heteroskedastic timepoints, and\n",
    "weakly-informative priors for ill-posed initial conditions\n",
    "([Gelman et al., 2008](https://doi.org/10.1214/08-AOAS191)).\n",
    "The softplus-lag representation itself and the saturation-clamp form\n",
    "extend the classical 3-phase bacterial-growth framing of\n",
    "[Zwietering et al., 1990](https://doi.org/10.1128/aem.56.6.1875-1881.1990)\n",
    "with smoothly-differentiable transitions that play nicely with\n",
    "gradient-based optimizers.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cf98d7f",
   "metadata": {},
   "source": [
    "## 1. The Model Equation\n",
    "\n",
    "The fit combines a **softplus lag** with an optional **softplus\n",
    "saturation ceiling**. The unclamped (growth-only) form is\n",
    "\n",
    "$$\n",
    "s_{\\mathrm{lin}}(t) \\;=\\; \\frac{v}{\\alpha}\\,\\ln\\!\\bigl(1 + e^{\\alpha(t - \\lambda)}\\bigr) + s_0 ,\n",
    "$$\n",
    "\n",
    "and the clamped form folds this through a second softplus to impose\n",
    "the carrying-capacity ceiling:\n",
    "\n",
    "$$\n",
    "s(t) \\;=\\; s_{\\max} \\;-\\; \\frac{1}{\\beta}\\,\\ln\\!\\bigl(1 + e^{\\beta(s_{\\max} - s_{\\mathrm{lin}}(t))}\\bigr) .\n",
    "$$\n",
    "\n",
    "### Why softplus rather than a piecewise-linear lag?\n",
    "\n",
    "The softplus function\n",
    "\n",
    "$$\n",
    "\\mathrm{sp}_\\alpha(x) \\;=\\; \\frac{1}{\\alpha}\\ln\\!\\bigl(1 + e^{\\alpha x}\\bigr)\n",
    "$$\n",
    "\n",
    "is the canonical $C^\\infty$ approximation of $\\max(x, 0)$: far below\n",
    "$x=0$ it decays exponentially to $0$, far above it grows linearly\n",
    "with slope $1$, and the transition sharpness is set by $\\alpha$. The\n",
    "lag term in the model is a $v$-scaled softplus of $(t - \\lambda)$, so\n",
    "the curve transitions smoothly from \"no growth\" to \"linear growth at\n",
    "rate $v$\" without any discontinuous derivatives. That smoothness\n",
    "matters for gradient-based least-squares solvers: a kink in the\n",
    "model creates a kink in the residual surface, and\n",
    "Levenberg–Marquardt / trust-region solvers do not recover gracefully\n",
    "from kinks.\n",
    "\n",
    "The saturation clamp plays the same role at the ceiling. Subtracting\n",
    "a softplus of the excess $(s_{\\max} - s_{\\mathrm{lin}})$ from\n",
    "$s_{\\max}$ flattens the curve smoothly once it approaches the\n",
    "carrying capacity. When $s_{\\mathrm{lin}}(t)$ is well below\n",
    "$s_{\\max}$ the clamp is inactive (``softplus(large positive)`` is\n",
    "essentially the identity); when it is well above, the clamp pins the\n",
    "output to $s_{\\max}$ (``softplus(large negative)`` is nearly zero).\n",
    "The transition sharpness is set by $\\beta$.\n",
    "\n",
    "### Numerical detail\n",
    "\n",
    "The implementation uses :func:`numpy.logaddexp`, the numerically\n",
    "stable form of $\\ln(1 + e^x)$. Naïve evaluation of\n",
    "$\\ln(1 + e^{\\alpha(t-\\lambda)})$ overflows when\n",
    "$\\alpha(t-\\lambda) \\gtrsim 700$; ``logaddexp(0, x)`` computes the same\n",
    "quantity as $\\max(0, x) + \\ln(1 + e^{-|x|})$ and stays well-behaved\n",
    "across the full range the optimizer visits.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a96fd04c",
   "metadata": {},
   "source": [
    "## 2. Parameter Cheat-Sheet\n",
    "\n",
    "| Symbol | Role | Biological interpretation | Typical scale |\n",
    "|--------|------|---------------------------|---------------|\n",
    "| $v$    | Post-lag linear growth rate | Maximum specific growth rate reached during the exponential/linear phase | Units of $s$ per unit $t$ (e.g. pixels²/hour) |\n",
    "| $s_0$  | Intercept of the post-lag line | Effective inoculum size, back-extrapolated from the linear phase | Same units as $s$ |\n",
    "| $\\lambda$ | Lag duration | Time at which the curve enters the linear phase | Same units as $t$ |\n",
    "| $\\alpha$ | Lag transition sharpness | $\\alpha \\to \\infty$ gives a hard switch-on; $\\alpha \\approx 2$ is very gentle | Dimensionally $1/t$; fit bounds are $[2, 50]$ |\n",
    "| $s_{\\max}$ | Carrying capacity | Plateau imposed by nutrient limitation, space, or toxic by-product accumulation | Same units as $s$ |\n",
    "| $\\beta$ | Saturation transition sharpness | Same role as $\\alpha$, but for the ceiling | Dimensionally $1/s$; fit bounds are $[2, 50]$ in ``fitted_beta`` mode; default scalar is ``10.0`` |\n",
    "\n",
    "Two things worth highlighting:\n",
    "\n",
    "- **$s_0$ is an *intercept*, not the $t=0$ observation.** The\n",
    "  softplus lag is not identically zero for $t < \\lambda$ — it is\n",
    "  *exponentially small but non-zero*. So $s(0) = s_0 + (v/\\alpha)\\,\n",
    "  \\ln(1 + e^{-\\alpha\\lambda})$, which is numerically very close to\n",
    "  $s_0$ for typical $\\alpha\\lambda \\gg 1$ but not exactly equal. In\n",
    "  the limit of a perfectly sharp switch-on $(\\alpha \\to \\infty)$, the\n",
    "  softplus collapses to $\\max(t - \\lambda, 0)$ and $s_0$ is recovered\n",
    "  exactly as the pre-lag baseline.\n",
    "- **$v$ and $\\lambda$ are only separately identifiable with enough\n",
    "  sub-saturation points.** On a fully-clamped curve, many $(v,\n",
    "  \\lambda)$ pairs produce the same plateau, and the fit can drift\n",
    "  into implausible regions. Section 3 discusses how pruning and\n",
    "  priors address this.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb308327",
   "metadata": {},
   "source": [
    "## 3. DoubleSoftplus: Per-Group Mode Dispatch\n",
    "\n",
    "``LinearSoftplus`` is single-mode by construction — it always fits\n",
    "the 4-parameter unclamped form. ``DoubleSoftplus`` is the saturating\n",
    "variant, and within it the fitter picks one of two variants per\n",
    "group, recorded in the ``DoubleSoftplus_mode`` column of the results.\n",
    "\n",
    "| Mode | Trigger | Free parameters |\n",
    "|------|---------|-----------------|\n",
    "| ``\"fixed_beta\"`` | ``beta`` user-supplied, **or** no shoulder detected | $v, s_0, \\lambda, \\alpha$ (clamped, $\\beta$ held constant) |\n",
    "| ``\"fitted_beta\"`` | ``beta is None`` and a shoulder is detected | $v, s_0, \\lambda, \\alpha, \\beta$ |\n",
    "\n",
    "When a curve has no clear plateau, it likely belongs to\n",
    "``LinearSoftplus`` instead — using ``DoubleSoftplus`` on\n",
    "non-saturating data still produces a fit (with $\\beta$ pinned at the\n",
    "module default and ``smax`` set to the per-group observed max), but\n",
    "the inferred carrying capacity will be sampling-window-dependent\n",
    "rather than biologically meaningful.\n",
    "\n",
    "### Shoulder detection\n",
    "\n",
    "Shoulder detection lives in ``_has_saturation_shoulder`` and is\n",
    "deliberately conservative:\n",
    "\n",
    "1. **Dynamic-range gate.** If the finite-$y$ range is $\\leq 5\\%$ of\n",
    "   $|\\mathrm{median}(y)|$, the group is treated as flat noise and no\n",
    "   shoulder is declared. This guards against coincidental ratio tests\n",
    "   firing on curves that never really grew.\n",
    "2. **Peak-slope estimation.** The gradient $dy/dt$ is smoothed with a\n",
    "   small box filter whose width scales with the series length\n",
    "   (``min(5, max(3, len(y)//4))``); the maximum of the smoothed\n",
    "   gradient becomes the reference \"peak slope\" — robust to\n",
    "   single-point noise spikes that would otherwise inflate a raw max.\n",
    "3. **Terminal-slope collapse.** The minimum of the last few *raw*\n",
    "   gradients (``tail_k = max(2, min(3, len(dy_dt) // 5))``) is\n",
    "   compared against ``shoulder_slope_ratio × peak_slope`` (default\n",
    "   $0.05 \\times$ peak). Using the min (not the mean) of the tail\n",
    "   keeps the criterion sensitive to short plateaus where just the\n",
    "   final $2$–$3$ points have collapsed to near-zero slope.\n",
    "\n",
    "If the terminal-slope collapse condition fires, the group has a\n",
    "shoulder; otherwise it does not.\n",
    "\n",
    "### When to pick ``LinearSoftplus`` vs ``DoubleSoftplus``\n",
    "\n",
    "Picking the right class up front is more honest than relying on\n",
    "runtime dispatch to silently switch models. If a fit group has not\n",
    "plateaued within the observed window and you fit a saturated model,\n",
    "pinning $s_{\\max}$ to the observed maximum silently encodes *\"this\n",
    "is the carrying capacity\"* — even though the data only supports\n",
    "*\"growth has not plateaued yet\"*. Those are very different claims.\n",
    "The split into two classes makes that choice explicit:\n",
    "\n",
    "- Use ``LinearSoftplus`` when growth is still in the linear regime.\n",
    "  $s_{\\max}$ and $\\beta$ never appear in the output.\n",
    "- Use ``DoubleSoftplus`` when there is a real plateau worth\n",
    "  recovering. The mode column records whether $\\beta$ was fit or\n",
    "  held at the default.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c25c6e6",
   "metadata": {},
   "source": [
    "## 4. Outlier Mitigations\n",
    "\n",
    "Field colony data is rarely well-behaved. A typical run has:\n",
    "\n",
    "- A handful of catastrophically wrong measurements — bubble\n",
    "  artefacts, segmentation errors, contamination spikes,\n",
    "  mis-registered timepoints.\n",
    "- Heteroskedastic noise — early timepoints often have a factor-of-10\n",
    "  smaller absolute $\\sigma$ than late ones, yet contribute the same\n",
    "  number of residuals in an unweighted fit.\n",
    "- Replicate coverage that varies timepoint-to-timepoint — some\n",
    "  timepoints have many plates, others have only one.\n",
    "- An inoculum whose magnitude is known *a priori* to a few percent\n",
    "  but that the optimizer has no access to — so $s_0$ routinely drifts\n",
    "  into biologically implausible regions when the early-lag fit is\n",
    "  under-determined.\n",
    "\n",
    "Both ``LinearSoftplus`` and ``DoubleSoftplus`` ship the same set of\n",
    "mitigations, with one exception: saturation pruning lives only on\n",
    "``LinearSoftplus`` (in ``DoubleSoftplus`` the saturation IS the\n",
    "model — pruning it would defeat the fit). The mitigations are\n",
    "orthogonal axes of the fit, composable, and the library defaults\n",
    "already sit at the robust end of each dial. Several power-user knobs\n",
    "that used to be public ``__init__`` kwargs (``stderr_floor_quantile``,\n",
    "``v_upper``, the pruning slope-ratio and buffer) are now class\n",
    "constants on the shared base — set them by subclass override if you\n",
    "need to deviate.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14104dd9",
   "metadata": {},
   "source": [
    "### 4.1 Robust Loss (``loss``, ``f_scale``)\n",
    "\n",
    ":func:`scipy.optimize.least_squares` exposes a family of robust\n",
    "$\\rho$-functions that down-weight large residuals\n",
    "([Huber, 1964](https://doi.org/10.1214/aoms/1177703732); the\n",
    "[scipy reference](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html)\n",
    "lists the exact forms). Both ``LinearSoftplus`` and ``DoubleSoftplus``\n",
    "default to ``loss=\"huber\"`` with ``f_scale=1.0``.\n",
    "\n",
    "The optimizer minimises $\\tfrac{1}{2}\\,\\mathrm{cost}$ where\n",
    "\n",
    "$$\n",
    "\\mathrm{cost}(x) \\;=\\; \\sum_i \\rho\\!\\Bigl(r_i(x)^2 / f_{\\mathrm{scale}}^2\\Bigr) ,\n",
    "$$\n",
    "\n",
    "and $\\rho$ is the chosen loss function applied to the squared,\n",
    "rescaled residuals. With $\\rho(z) = z$ (``loss=\"linear\"``) this is\n",
    "the classical least-squares objective; with the Huber $\\rho$,\n",
    "residuals well below the $f_\\mathrm{scale}$ knee contribute like\n",
    "$z$ (i.e. squared error), while residuals well above it contribute\n",
    "proportionally to $\\sqrt{z}$ (i.e. absolute error, reducing their\n",
    "influence).\n",
    "\n",
    "``f_scale`` is the knee of the loss. Its units depend on whether the\n",
    "fit is weighted:\n",
    "\n",
    "- **Unweighted fit.** $f_\\mathrm{scale}$ is measured in the same units\n",
    "  as the residual — so $f_\\mathrm{scale} = 1$ means \"any residual\n",
    "  larger than $1$ unit of $s$ starts getting downweighted.\"\n",
    "- **Weighted fit.** Residuals are pre-divided by $\\sigma$ before\n",
    "  entering the loss, so $f_\\mathrm{scale}$ is measured in multiples\n",
    "  of $\\sigma$ — $f_\\mathrm{scale} = 3$ treats points more than\n",
    "  $\\sim\\!3\\sigma$ off the curve as outliers.\n",
    "\n",
    "Choice of ``loss``:\n",
    "\n",
    "- ``\"linear\"`` — classical least-squares, no outlier protection.\n",
    "- ``\"huber\"`` (default) — quadratic inside the knee, linear outside.\n",
    "  The least aggressive robust option; recovers least-squares\n",
    "  efficiency on clean data while limiting the influence of true\n",
    "  outliers.\n",
    "- ``\"soft_l1\"`` — smooth approximation to $\\ell_1$ loss; similar to\n",
    "  Huber but with a smoother transition at the knee.\n",
    "- ``\"cauchy\"`` — heavy-tailed; strongly suppresses outliers but can\n",
    "  introduce multiple local minima.\n",
    "- ``\"arctan\"`` — bounded influence function (asymptotes to a\n",
    "  constant); the most aggressive outlier suppression available in\n",
    "  scipy.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8622d29e",
   "metadata": {},
   "source": [
    "### 4.2 Inverse-Variance Weighting (``stderr_label``)\n",
    "\n",
    "When replicates are available, the fitter divides each residual by a\n",
    "per-timepoint standard error,\n",
    "\n",
    "$$\n",
    "r_i \\;=\\; \\frac{y_i - s(t_i)}{\\sigma_i} ,\n",
    "$$\n",
    "\n",
    "so tight-replicate timepoints pull the fit harder than noisy ones.\n",
    "This addresses heteroskedasticity directly: at late timepoints\n",
    "colony variance is often 5–10× the early-timepoint variance, and\n",
    "without weighting those residuals dominate the loss purely because\n",
    "they are larger in absolute terms, not because they are more\n",
    "informative.\n",
    "\n",
    "Priority order inside ``_resolve_y_stderr``:\n",
    "\n",
    "1. **User-supplied ``stderr_label``.** If the user passes a column\n",
    "   name and that column exists in the aggregated group, those\n",
    "   $\\sigma$ values are used verbatim. No pool fallback — the user\n",
    "   has opted into their own $\\sigma$ semantics.\n",
    "2. **Auto-derived replicate SEM** (when ``stderr_label is None``).\n",
    "   During ``analyze`` the fitter computes\n",
    "   ``data.groupby(groupby + [time_label])[on].transform('sem')`` and\n",
    "   stores it in a column named ``{on}_stderr``. That column is then\n",
    "   read as $\\sigma$ for each aggregated row.\n",
    "3. **No weights** (return ``None``) if neither path is available.\n",
    "\n",
    "The auto-SEM path is the common case: the user does not need to\n",
    "pre-compute $\\sigma$, and the per-timepoint standard error naturally\n",
    "captures both the measurement noise and any per-plate variability\n",
    "contributing to that timepoint.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14422d2a",
   "metadata": {},
   "source": [
    "### 4.3 Replicate-SEM Quantile Floor (class constant ``_STDERR_FLOOR_QUANTILE``)\n",
    "\n",
    "Pure inverse-variance weighting has a well-known failure mode. When\n",
    "replicates happen to agree by chance — or when the measurement\n",
    "process is discretised so a handful of reps land on identical\n",
    "values — the SEM collapses toward zero and the corresponding\n",
    "$1/\\sigma^2$ weight blows up. A few \"lucky\" timepoints can then pin\n",
    "the entire fit, because the optimizer will sacrifice a huge amount\n",
    "of fit quality elsewhere to gain a tiny improvement on the\n",
    "$1/\\sigma^2 \\to \\infty$ residuals.\n",
    "\n",
    "The shared base class sets ``_STDERR_FLOOR_QUANTILE = 0.25`` (a\n",
    "class constant, not an ``__init__`` kwarg). The procedure is:\n",
    "\n",
    "1. Collect all finite, positive $\\sigma$ values in the group.\n",
    "2. Compute the $q$-th quantile (default $q = 0.25$) of that set.\n",
    "3. Replace any $\\sigma_i$ below the quantile with the quantile value:\n",
    "   $\\sigma_i \\leftarrow \\max(\\sigma_i, \\sigma_q)$.\n",
    "\n",
    "The ratio between the stiffest and softest point in a group is\n",
    "therefore capped at $\\sigma_{\\max} / \\sigma_q$, which at $q = 0.25$\n",
    "is at most a factor of about $16$ under typical empirical $\\sigma$\n",
    "distributions. That is large enough to preserve legitimate\n",
    "heteroskedasticity (orders-of-magnitude differences across time) but\n",
    "small enough that a coincidentally-agreeing replicate cannot dominate\n",
    "the fit.\n",
    "\n",
    "Two notes:\n",
    "\n",
    "- **Disabling.** Subclass ``LinearSoftplus`` (or ``DoubleSoftplus``)\n",
    "  and set ``_STDERR_FLOOR_QUANTILE = None`` to recover raw inverse-\n",
    "  variance weighting. This is appropriate when the ``stderr_label``\n",
    "  column is trusted across its full range, or when you are\n",
    "  deliberately probing the effect of the floor.\n",
    "- **All-zero / all-NaN fallback.** When a fit group has *no* finite\n",
    "  positive $\\sigma$ at all (every timepoint has $n=1$, or noise-free\n",
    "  synthetic replicates that all coincide), the model bypasses the\n",
    "  weighted pathway entirely and falls back to unweighted residuals.\n",
    "  The floor is inapplicable in that case by construction.\n",
    "\n",
    "The floor applies uniformly to both user-supplied ``stderr_label``\n",
    "columns and the auto-derived replicate-SEM column, so the same\n",
    "calibration works across both paths.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6c0d3f9",
   "metadata": {},
   "source": [
    "### 4.4 Pooled Point-Level Std for $n=1$ Timepoints\n",
    "\n",
    "A related pathology: if a fit group has a mix of multi-replicate and\n",
    "singleton ($n=1$) timepoints, the singletons get $\\sigma =\n",
    "\\mathrm{NaN}$ from the SEM calculation. A naïve fix — filling the\n",
    "NaN with a tiny $\\varepsilon$ — would let the singleton point\n",
    "dominate the fit with $1/\\varepsilon^2$ weight, the same pathology\n",
    "the quantile floor addresses but arriving via a different route.\n",
    "\n",
    "The fitter's auto-$\\sigma$ path pre-computes a *pooled point-level\n",
    "std* alongside the SEM column during ``analyze``. Concretely:\n",
    "\n",
    "1. ``{on}_stderr`` — per-timepoint standard error, via\n",
    "   ``groupby(groupby + [time_label])[on].transform('sem')``.\n",
    "2. ``{on}_std_pool`` — per-timepoint *std* broadcast and then\n",
    "   median-reduced across the fit group:\n",
    "\n",
    "   ```python\n",
    "   std_pt  = groupby(groupby + [time_label])[on].transform('std')\n",
    "   pool    = std_pt.groupby(groupby).transform('median')\n",
    "   ```\n",
    "\n",
    "   So each fit group gets one scalar: the median of its\n",
    "   multi-replicate timepoints' per-timepoint stds.\n",
    "\n",
    "Inside ``_resolve_y_stderr``, NaN / zero entries in the SEM vector\n",
    "are filled with the pool value rather than $\\varepsilon$. The $n=1$\n",
    "rows end up with $\\sigma \\approx$ the typical per-timepoint noise\n",
    "level, so their weight is commensurate with — rather than dominant\n",
    "over — the multi-replicate rows.\n",
    "\n",
    "Accepted trade-off: median-of-transform ensures that timepoints with\n",
    "*more* replicates contribute *proportionally more rows* (and hence\n",
    "more weight) to the pool. This tends to enlarge the pool-$\\sigma$\n",
    "relative to a simple unweighted median of per-timepoint stds, which\n",
    "is the safer direction — pool-filled $n=1$ rows will never dominate\n",
    "over their better-supported $n \\geq 2$ neighbours.\n",
    "\n",
    "Three edge cases:\n",
    "\n",
    "- **User-supplied ``stderr_label`` path.** The pool is *not* used\n",
    "  here; NaN / zero entries fall back to the $\\varepsilon$-scaled\n",
    "  fill. The user has committed to their own $\\sigma$ semantics, so\n",
    "  the library doesn't silently cross-contaminate them with\n",
    "  replicate-derived pool values.\n",
    "- **Fully-singleton fit group.** If every timepoint has $n = 1$, the\n",
    "  pool is itself NaN (std needs $n \\geq 2$). ``_resolve_y_stderr``\n",
    "  returns ``None`` and the fit runs unweighted — the same\n",
    "  fallback as the \"no finite positive $\\sigma$\" case.\n",
    "- **Quantile floor interaction.** Pool-filled entries are *not*\n",
    "  included in the set used to compute the floor quantile — the\n",
    "  floor is derived from finite, positive, non-filled $\\sigma$ values\n",
    "  only — but they *are* subsequently lifted to the floor via the\n",
    "  final ``np.maximum``. So pool-filled rows inherit both the pool\n",
    "  baseline and any floor that is larger than the pool.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e07912b2",
   "metadata": {},
   "source": [
    "### 4.5 Saturation Pruning (``LinearSoftplus.prune_saturated``)\n",
    "\n",
    "Pruning lives only on ``LinearSoftplus`` — ``DoubleSoftplus`` models\n",
    "the saturation explicitly and trimming the plateau would defeat the\n",
    "fit.\n",
    "\n",
    "Post-saturation points cause two problems for a non-saturating model.\n",
    "First, they are often where contamination overruns a colony or where\n",
    "segmentation confuses neighbour-colony merges, pulling the late-time\n",
    "values upward or downward in biologically meaningless ways. Second,\n",
    "even perfectly clean plateau points are **uninformative** about the\n",
    "lag-phase parameters $(\\lambda, \\alpha, s_0)$ — yet they contribute\n",
    "an equal share of residuals unless something trims them off.\n",
    "\n",
    "``_prepare_group`` implements a hybrid heuristic. A point is flagged\n",
    "as saturated when **both** of the following are true:\n",
    "\n",
    "1. **Amplitude gate.**\n",
    "   $y \\;\\geq\\; \\min(y) + 0.90\\,\\bigl(\\max(y) - \\min(y)\\bigr)$.\n",
    "   The point has climbed into the upper $10\\%$ of the observed\n",
    "   dynamic range.\n",
    "2. **Sustained derivative run.** Three consecutive smoothed-gradient\n",
    "   values below ``_PRUNE_SLOPE_RATIO × peak_slope`` (default\n",
    "   $0.05 \\times$ peak). A point is part of a *sustained* low-slope\n",
    "   run, not just an isolated dip.\n",
    "\n",
    "Both gates must agree before any trimming happens. Only rows past\n",
    "the latest index satisfying both — plus a small ``_PRUNE_BUFFER``\n",
    "class constant (default $2$) of rows kept after the cut so the fit\n",
    "still sees some plateau evidence — are removed.\n",
    "\n",
    "Why both gates?\n",
    "\n",
    "- The **amplitude gate** is structurally immune to lag-phase noise.\n",
    "  During the lag phase $y \\approx s_0$, far below the\n",
    "  $0.9 \\times$ amplitude target, so pruning can never fire on\n",
    "  lag-phase points no matter how flat the derivative is.\n",
    "- The **sustained-run gate** rejects transient mid-growth dips. A\n",
    "  single low-slope point in the middle of the linear phase (e.g. a\n",
    "  measurement where two consecutive timepoints happened to land at\n",
    "  similar values) does not satisfy the \"three in a row\" requirement.\n",
    "\n",
    "Two additional guards:\n",
    "\n",
    "- **Tail-growth short-circuit.** If the mean of the tail-window\n",
    "  smoothed slope is *still* above the threshold, pruning returns\n",
    "  the group unchanged. Curves that never plateau within the\n",
    "  observation window are not trimmed.\n",
    "- **Degenerate-input guard.** If the smoothed peak slope is\n",
    "  non-finite or non-positive (all-NaN $y$, or a constant series),\n",
    "  pruning returns the group unchanged.\n",
    "\n",
    "Set ``prune_saturated=False`` to disable pruning entirely. Subclass\n",
    "``LinearSoftplus`` and override ``_PRUNE_SLOPE_RATIO`` /\n",
    "``_PRUNE_BUFFER`` for tuning; these are class constants, not\n",
    "``__init__`` kwargs.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30b4358d",
   "metadata": {},
   "source": [
    "### 4.6 Gaussian Prior on $s_0$ (``s0_prior``, ``s0_prior_cv``, ``s0_prior_sigma``)\n",
    "\n",
    "When the earliest observed timepoint is already well past the lag\n",
    "(e.g. the first capture is at $t = 12$ hr and the lag was\n",
    "$\\lambda = 4$ hr), the fit is under-determined: $s_0$ and $\\lambda$\n",
    "trade off freely along an approximately linear manifold and the\n",
    "optimizer has no principled basis for choosing one point on the\n",
    "manifold over another. A weakly-informative Gaussian prior on $s_0$\n",
    "breaks the symmetry without driving the fit toward implausible\n",
    "regions — the standard weakly-informative-prior argument from\n",
    "[Gelman et al., 2008](https://doi.org/10.1214/08-AOAS191).\n",
    "\n",
    "### Implementation as a virtual residual\n",
    "\n",
    "The prior is added to the least-squares residual vector as a single\n",
    "extra row:\n",
    "\n",
    "$$\n",
    "r_{\\mathrm{prior}} \\;=\\; \\frac{s_0 - \\mu_{\\mathrm{prior}}}{\\sigma_{\\mathrm{prior}}} .\n",
    "$$\n",
    "\n",
    "Because least-squares with Gaussian noise is equivalent to maximum-\n",
    "likelihood under a Gaussian likelihood, appending a Gaussian-shaped\n",
    "virtual residual is equivalent to maximum *a posteriori* (MAP)\n",
    "estimation with a Gaussian prior on $s_0$. The optimizer sees no\n",
    "distinction between the prior residual and the data residuals; the\n",
    "prior simply adds $(s_0 - \\mu)^2 / \\sigma^2$ to the objective,\n",
    "which is exactly the negative log-density of a Gaussian prior up to\n",
    "additive constants.\n",
    "\n",
    "### ``s0_prior`` dispatch (by type)\n",
    "\n",
    "The ``s0_prior`` argument is polymorphic — the type dispatches which\n",
    "$\\mu$-estimation strategy is used.\n",
    "\n",
    "| Value | Behaviour |\n",
    "|-------|-----------|\n",
    "| ``None`` or ``False`` | No prior (default). Residual vector is unchanged. |\n",
    "| ``True`` | Data-ground on ``self.on``. $\\mu$ is the median of the target column at the earliest observed timepoint within each effective group. |\n",
    "| ``str`` (column name) | Data-ground on the named column. $\\mu$ is the median of ``data[s0_prior]`` at the earliest observed timepoint within each effective group. Useful when the measurement column (e.g. ``Shape_Area``) is not a faithful inoculum estimate, but another column (e.g. a calibrated ``Inoculum_Size`` observation from a separate imaging step) is. |\n",
    "| positive ``int`` / ``float`` | Scalar $\\mu$, applied uniformly to every fit group. |\n",
    "| anything else | Raises :class:`TypeError`. |\n",
    "| non-positive or non-finite scalar | Raises :class:`ValueError`. |\n",
    "\n",
    "The ``bool`` branch comes *before* the numeric branch in the\n",
    "underlying helper because ``bool`` is a subclass of ``int`` in\n",
    "Python — so ``isinstance(True, int)`` is ``True`` and an unguarded\n",
    "numeric branch would consume the boolean inputs.\n",
    "\n",
    "### Choosing $\\sigma$: ``s0_prior_cv`` vs ``s0_prior_sigma``\n",
    "\n",
    "The prior $\\sigma$ is set by exactly one of two mutually exclusive\n",
    "knobs:\n",
    "\n",
    "| Argument | Behaviour |\n",
    "|----------|-----------|\n",
    "| ``s0_prior_cv`` | Coefficient of variation: $\\sigma = \\mathrm{cv} \\times \\mu$ |\n",
    "| ``s0_prior_sigma`` | Absolute σ: $\\sigma = \\mathrm{sigma}$, independent of $\\mu$ |\n",
    "\n",
    "If both are set, construction raises :class:`ValueError`. If neither\n",
    "is set and the prior is engaged, the helper applies CV $= 0.05$ as a\n",
    "moderately informative default — the same effective default as the\n",
    "legacy ``s0_prior_factor=0.05``. Pick CV when the prior strength\n",
    "should scale with magnitude (typical for biologically positive\n",
    "quantities); pick absolute σ when you have a hard tolerance budget.\n",
    "\n",
    "### Empirical-Bayes pooling (``s0_prior_groupby``)\n",
    "\n",
    "For column-backed priors (``s0_prior = True`` or a string),\n",
    "``s0_prior_groupby`` lets you pool $\\mu$ estimation over a *coarser*\n",
    "grouping than ``groupby``. For example, if the fit grouping is\n",
    "``[\"Media\", \"Strain\", \"Plate\"]`` but inoculation spread varies\n",
    "per-media and replicates within a media/strain combination are\n",
    "exchangeable, setting ``s0_prior_groupby=[\"Media\"]`` computes the\n",
    "prior $\\mu$ as the median of the earliest-timepoint observations\n",
    "*within each media*, then broadcasts that one $\\mu$ to every\n",
    "strain/plate within the media.\n",
    "\n",
    "This is the classical empirical-Bayes move —\n",
    "[Efron & Morris, 1975](https://doi.org/10.1080/01621459.1975.10479864)\n",
    "and the shrinkage literature that followed — for situations where\n",
    "per-fit-group $\\mu$ estimates would be noisy but the group-level\n",
    "average is stable.\n",
    "\n",
    "Constraints, verified at ``analyze`` time:\n",
    "\n",
    "- ``s0_prior_groupby`` must be a strict subset of ``groupby``\n",
    "  (coarser, not different). Violating this raises\n",
    "  :class:`ValueError`.\n",
    "- Columns referenced in ``s0_prior_groupby`` must exist in the data.\n",
    "- ``s0_prior_groupby`` must be non-empty if supplied. Pass ``None``\n",
    "  (the default) to fall back to using the full fit ``groupby`` for\n",
    "  prior-$\\mu$ estimation.\n",
    "\n",
    "### When the prior does (and does not) engage\n",
    "\n",
    "The helper's ``is_configured`` property summarises the engagement\n",
    "logic: the prior attempts to apply when either a column label is set\n",
    "(``s0_prior=True`` or ``s0_prior=<str>``) or a scalar ``direct_mean``\n",
    "is set (``s0_prior`` was a positive number). If neither applies, the\n",
    "prior is disabled.\n",
    "\n",
    "Even when configured, the prior can still skip a group when:\n",
    "\n",
    "- The broadcast $\\mu$ column is missing for the group (the earliest\n",
    "  timepoint had no usable column values — e.g. all NaN).\n",
    "- The resulting $\\mu$ or $\\sigma$ is non-finite or non-positive.\n",
    "\n",
    "In either case, ``_inoc_stats(group)`` returns ``None`` and the fit\n",
    "runs without the prior residual for that group. Other groups in the\n",
    "same ``analyze`` call are unaffected.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0679f6ce",
   "metadata": {},
   "source": [
    "## 5. Putting It Together\n",
    "\n",
    "The mitigations are orthogonal axes of the fit. The library defaults\n",
    "already sit at the robust end of every dial:\n",
    "\n",
    "| Argument | Lives on | Default | Addresses |\n",
    "|----------|----------|---------|-----------|\n",
    "| ``loss`` | both | ``\"huber\"`` | Catastrophic outliers (spikes, mis-segmentation) |\n",
    "| ``f_scale`` | both | ``1.0`` | Knee of the robust loss |\n",
    "| ``prune_saturated`` | ``LinearSoftplus`` | ``True`` | Uninformative plateau points, tail contamination |\n",
    "| ``stderr_label`` | both | ``None`` | Auto-derives replicate SEM during aggregation |\n",
    "| ``s0_prior`` | both | ``None`` | Opt-in; off by default because it changes the objective |\n",
    "| ``s0_prior_cv`` | both | ``None`` (defaults to ``0.05`` if prior engaged) | Relative prior strength |\n",
    "| ``s0_prior_sigma`` | both | ``None`` | Absolute prior σ (mutually exclusive with cv) |\n",
    "| ``s0_prior_groupby`` | both | ``None`` | Empirical-Bayes pooling (opt-in) |\n",
    "| ``smax`` | ``DoubleSoftplus`` | ``None`` (per-group max) | Carrying capacity |\n",
    "| ``beta`` | ``DoubleSoftplus`` | ``None`` (mode dispatch) | Saturation transition sharpness |\n",
    "| ``shoulder_slope_ratio`` | ``DoubleSoftplus`` | ``0.05`` | Shoulder-detection threshold for fitted_beta vs fixed_beta |\n",
    "\n",
    "Class constants (subclass override; not ``__init__`` kwargs):\n",
    "\n",
    "| Constant | Default | Lives on |\n",
    "|----------|---------|----------|\n",
    "| ``_V_UPPER`` | ``50.0`` | shared base |\n",
    "| ``_STDERR_FLOOR_QUANTILE`` | ``0.25`` | shared base |\n",
    "| ``_PRUNE_SLOPE_RATIO`` | ``0.05`` | ``LinearSoftplus`` |\n",
    "| ``_PRUNE_BUFFER`` | ``2`` | ``LinearSoftplus`` |\n",
    "\n",
    "The per-group mode dispatch on ``DoubleSoftplus``\n",
    "(``fitted_beta`` / ``fixed_beta``) runs on top of these and needs no\n",
    "configuration in the common case — $\\beta$ is fit when a shoulder is\n",
    "present, held at the default when no shoulder is detected.\n",
    "\n",
    "## References\n",
    "\n",
    "- Efron, B., & Morris, C. (1975). Data analysis using Stein's estimator and its generalizations. *Journal of the American Statistical Association*, 70(350), 311–319. https://doi.org/10.1080/01621459.1975.10479864\n",
    "- Gelman, A., Jakulin, A., Pittau, M. G., & Su, Y.-S. (2008). A weakly informative default prior distribution for logistic and other regression models. *The Annals of Applied Statistics*, 2(4), 1360–1383. https://doi.org/10.1214/08-AOAS191\n",
    "- Huber, P. J. (1964). Robust estimation of a location parameter. *The Annals of Mathematical Statistics*, 35(1), 73–101. https://doi.org/10.1214/aoms/1177703732\n",
    "- Zwietering, M. H., Jongenburger, I., Rombouts, F. M., & van 't Riet, K. (1990). Modeling of the bacterial growth curve. *Applied and Environmental Microbiology*, 56(6), 1875–1881. https://doi.org/10.1128/aem.56.6.1875-1881.1990\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}