feat(kinetics): too-few-points warning, intercept checks, 'Pseudo-first (initial drop)' model

CHANGELOG.md

@@ -2,6 +2,22 @@
 
 ## Unreleased
 
+### Added
+- **Tab 1 warns when there are too few points for model selection.** When a
+  catalyst has fewer than `MIN_FIT_POINTS` informative data points (the t=0
+  anchor excluded), the summary shows a warning that AICc over-penalises
+  two-parameter models and cannot evaluate three-parameter ones. The summary
+  table now also reports `n (fitted)` per catalyst.
+- **Tab 2 flags linearisation intercepts that contradict the model.** The
+  Pseudo-first and Pseudo-second-order linearisations are checked against their
+  theoretical intercepts (0 and 1/C₀, 95 % CI, t = 0 anchor excluded). A PFO
+  line that does not pass through the origin warns about a fast initial removal
+  step; a PSO mismatch is noted only when the straight line fits (R² ≥ 0.95).
+- **New kinetic model: "Pseudo-first (initial drop)".** Fits
+  C(t) = A·C₀·e^(−kt) to model a fast initial removal (rapid adsorption or a
+  fast initial reaction) followed by first-order decay. A = 1 recovers plain
+  PFO; the reported "initial drop" is 100·(1 − A).
+
 ### Fixed
 - **CSV uploads work in every tab.** All tabs read the same uploaded file and
   `read_csv` left the cursor at its end, so on a CSV upload every tab after the

README.md

@@ -48,7 +48,7 @@ pip install -r requirements-lock.txt
 
 | Tab | Module | Description |
 |-----|--------|-------------|
-| 1 | **Kinetic Fitting** | Fit 8 kinetic models, AICc model selection, k±SE, r₀, t½ |
+| 1 | **Kinetic Fitting** | Fit 9 kinetic models, AICc model selection, k±SE, r₀, t½ |
 | 2 | **Linearization** | Linear transforms (1/C vs t, ln(C₀/C) vs t) with best-model summary |
 | 3 | **Removal Efficiency** | Desulfurization efficiency (%) vs time + bar chart |
 | 4 | **TON / TOF** | Option A: site-based (metal catalysts) · Option B: mass-normalized (carbon-based) |
@@ -69,6 +69,7 @@ All models fitted by nonlinear least squares with C₀ locked.
 |-------|--------------------|----|-------|
 | Zero-order | $C_t = C_0 - k_0 t$ | $C_0 / (2k_0)$ | Simplified |
 | Pseudo-first-order | $C_t = C_0\, e^{-k_{app}t}$ | $\ln 2 / k_{app}$ | Simplified |
+| Pseudo-first (initial drop) | $C_t = A\,C_0\, e^{-k t}$ | $\ln(2A)/k$ | Semi-empirical |
 | Pseudo-second-order | $C_t = C_0 / (1 + k_2 C_0 t)$ | $1 / (k_2 C_0)$ | Simplified |
 | Elovich | $C_t = C_0 - \frac{1}{\beta}\ln(1+\alpha\beta t)$ | $(e^{C_0\beta/2}-1)/(\alpha\beta)$ | Phenomenological |
 | Langmuir-Hinshelwood | $dC/dt = -k_{LH} K C / (1+KC)$ (ODE) | $\ln2/(k_{LH}K) + C_0/(2k_{LH})$ | Mechanistic |

app_ods.py

@@ -45,6 +45,7 @@
     _power_law,
     _avrami,
     _double_exponential,
+    _first_order_drop,
     _r2,
     _fmt_sci,
     _fmt_thalf,
@@ -55,7 +56,14 @@
     _anchor_mask,
 )
 from catlab import __version__
-from catlab.metrics import _C0_both, _initial_tof_site, _initial_tof_mass, _arrhenius_ci, c_to_user
+from catlab.metrics import (
+    _C0_both,
+    _initial_tof_site,
+    _initial_tof_mass,
+    _arrhenius_ci,
+    c_to_user,
+    linear_intercept_check,
+)
 
 # Re-exported for backward compatibility (tests/test_app_ods.py).
 from catlab.kinetics_engine import _get_valid_models  # noqa: F401
@@ -283,6 +291,7 @@ def _sidebar_settings():
 **Model classes:**
 - *Mechanistic*: L-H (surface-reaction based)
 - *Simplified mechanistic*: Zero-, Pseudo-first-, Pseudo-second-order, Power-Law
+- *Semi-empirical*: Pseudo-first (initial drop) — two-stage (fast initial step + first-order)
 - *Phenomenological / empirical*: Elovich (chemisorption heterogeneity),
   Avrami (nucleation/growth — uncommon in ODS; use with caution),
   Double-Exponential (two-site parallel decay — high overfitting risk with < 10 points)
@@ -443,6 +452,8 @@ def _fit_curve(model, params, t_fine, C0):
         return _zero_order(t_fine, params[0], C0)
     elif model == "Pseudo-first":
         return _first_order(t_fine, params[0], C0)
+    elif model == "Pseudo-first (initial drop)":
+        return _first_order_drop(t_fine, params[0], params[1], C0)
     elif model == "Pseudo-second-order":
         return _second_order(t_fine, params[0], C0)
     elif model == "Elovich":
@@ -692,6 +703,21 @@ def _tab_kinetics(cfg, uploaded):
             )
         all_results[col] = _fit_nonlinear(t_fit, Ct_fit, C0)
 
+        # Warn once per catalyst when there are too few informative points for
+        # AICc model selection to be trustworthy (the t=0 anchor is not counted).
+        n_fit = next(
+            (all_results[col][m]["n_fit"] for m in model_names if "n_fit" in all_results[col][m]),
+            None,
+        )
+        if n_fit is not None and n_fit < MIN_FIT_POINTS:
+            cat_label = col.replace(" Removal (%)", "").strip()
+            st.warning(
+                f"⚠️ {cat_label}: only {n_fit} informative data points (t = 0 not counted). "
+                f"With fewer than {MIN_FIT_POINTS} points, AICc penalises every 2-parameter "
+                f"model heavily and cannot evaluate 3-parameter models, so the selected model "
+                f"is not reliable. Add measurements (especially early times, e.g. 5-10 min)."
+            )
+
     if not all_results:
         return
 
@@ -839,6 +865,7 @@ def _tab_kinetics(cfg, uploaded):
             title_lin = f"{cat_label} — Pseudo-second-order  |  1/C vs t"
         elif best in (
             "Pseudo-first",
+            "Pseudo-first (initial drop)",
             "Zero-order",
             "Power-Law",
             "Avrami",
@@ -919,6 +946,7 @@ def _tab_kinetics(cfg, uploaded):
                 "Adj-R²": br.get("adj_r2", "N/A"),
                 "AICc": br.get("aicc", "N/A"),
                 "AIC": br.get("aic", "N/A"),
+                "n (fitted)": br.get("n_fit", ""),
                 "t½ (min)": _fmt_thalf(t_half),
                 "r₀ (mol/L/min)": _fmt_sci(r0),
                 "r₀/m (mol/g/min)": _fmt_sci(r0_m),
@@ -1005,6 +1033,11 @@ def _tab_kinetics(cfg, uploaded):
                             "AICc": mr.get("aicc", "N/A"),
                             "AIC": mr.get("aic", "N/A"),
                             "t½ (min)": _fmt_thalf(mr.get("t_half", float("nan"))),
+                            "Initial drop (%)": (
+                                f"{mr['initial_drop_pct']:.1f}"
+                                if mr.get("initial_drop_pct") is not None
+                                else "–"
+                            ),
                             "Note": note,
                         }
                     )
@@ -1018,6 +1051,7 @@ def _tab_kinetics(cfg, uploaded):
                             "AICc": "–",
                             "AIC": "–",
                             "t½ (min)": "–",
+                            "Initial drop (%)": "–",
                             "Note": note,
                         }
                     )
@@ -1137,6 +1171,58 @@ def _tab_linearization(cfg, uploaded):
         st.pyplot(fig)
         plt.close(fig)
 
+    # ── Intercept consistency checks ─────────────────────────────
+    # Each linearised form has a theoretical intercept (Pseudo-first -> 0,
+    # Pseudo-second-order -> 1/C0).  A measured intercept far
+    # from that value signals a model/mechanism mismatch.  The check runs on the
+    # non-anchor points only (t = 0 with C = C0 is fitted exactly by construction).
+    # Zero-order (C vs t -> C0) is not checked: on any decelerating run the
+    # straight line misses C0, which repeats the Pseudo-first warning below.
+    intercept_checks = [
+        (
+            "Pseudo-first",
+            lambda t_, Ct: (t_[Ct > 0], np.log(C0 / np.maximum(Ct[Ct > 0], 1e-15))),
+            0.0,
+        ),
+        ("Pseudo-second-order", lambda t_, Ct: (t_[Ct > 0], 1.0 / Ct[Ct > 0]), 1.0 / C0),
+    ]
+    for col in removal_cols:
+        removal = df[col].dropna().values[: len(t)].astype(float)
+        Ct = C0 * (1 - removal / 100.0)
+        cat_label = col.replace(" Removal (%)", "").strip()
+        keep = ~_anchor_mask(t, Ct, C0)
+        t_nl = t[keep]
+        Ct_nl = Ct[keep]
+        if len(t_nl) < 3:
+            continue
+        for mname, transform, expected in intercept_checks:
+            x_vals, y_vals = transform(t_nl, Ct_nl)
+            if len(x_vals) < 3:
+                continue
+            chk = linear_intercept_check(x_vals, y_vals, expected)
+            if not chk["deviates"]:
+                continue
+            a = chk["intercept"]
+            ci = chk["ci_half"]
+            if mname == "Pseudo-first" and a > 0:
+                x_pct = 100.0 * (1.0 - np.exp(-a))
+                st.warning(
+                    f"⚠️ {cat_label}: the ln(C₀/C) line does not pass through the origin "
+                    f"(intercept = {a:.4g} ± {ci:.4g}). This corresponds to ≈ {x_pct:.0f} % "
+                    f"removed almost instantly at the start (X = 100·(1 − e^(−a))) — a fast "
+                    f"initial step such as rapid adsorption or a fast initial reaction. A "
+                    f"pseudo-first-order curve forced through C₀ will underfit; see the "
+                    f"'Pseudo-first (initial drop)' model in Tab 1."
+                )
+            elif chk["r2"] >= 0.95:
+                # A missed intercept only means something when the straight line
+                # otherwise fits; on curved data (e.g. Zero-order on any
+                # decelerating run) it fires every time and is noise.
+                st.info(
+                    f"ℹ️ {cat_label}: the {mname} linearisation intercept "
+                    f"({a:.4g} ± {ci:.4g}) deviates from its theoretical value {expected:.4g}."
+                )
+
     if not summary_rows:
         return
 

catlab/kinetics_engine.py

@@ -31,6 +31,7 @@
 N_PARAMS = {
     "Zero-order": 1,
     "Pseudo-first": 1,
+    "Pseudo-first (initial drop)": 2,
     "Pseudo-second-order": 1,
     "Elovich": 2,
     "L-H": 2,
@@ -83,6 +84,7 @@
     "Power-Law": [(1, "n", 0.1, 5.0)],
     "Avrami": [(1, "n", 0.1, 4.0)],
     "Double-Exponential": [(2, "A", 0.0, 1.0)],
+    "Pseudo-first (initial drop)": [(1, "A", 0.0, 1.0)],
 }
 BOUND_NEAR = 1e-3  # window: within 1e-3 * max(|bound|, 1) of the bound
 BOUND_STEP = 1e-6  # inward probe step, as a fraction of max(|bound|, 1)
@@ -129,6 +131,7 @@ def _params_at_bound(f, t, Ct, popt, specs):
 MODEL_NAMES = [
     "Zero-order",
     "Pseudo-first",
+    "Pseudo-first (initial drop)",
     "Pseudo-second-order",
     "Elovich",
     "L-H",
@@ -199,6 +202,19 @@ def _first_order(t, k, C0):
     return C0 * np.exp(-k * t)
 
 
+def _first_order_drop(t, k, A, C0):
+    """Pseudo-first-order with a fast initial step: C(t) = A·C0·exp(−k·t).
+
+    (1 − A) is the fraction removed almost instantly at the start (rapid
+    adsorption or a fast initial reaction).  A = 1 is plain PFO.  At t = 0 the
+    function returns A·C0 (the value just AFTER the drop); the (t=0, C0) anchor
+    is excluded from fitting and its re-expanded "pred" stays C0, the value just
+    BEFORE the drop (see _anchor_mask).
+    """
+    t = np.asarray(t, dtype=float)
+    return A * C0 * np.exp(-k * t)
+
+
 def _second_order(t, k, C0):
     return C0 / (1 + k * C0 * t)
 
@@ -422,6 +438,12 @@ def _anchor_mask(time, Ct, C0):
     points are therefore left out of fitting and of every statistic.
     A t = 0 point whose concentration differs from C0 (e.g. after a dark
     adsorption step) is informative and is kept.
+
+    The exception is "Pseudo-first (initial drop)", which returns A·C0 (not C0)
+    at t = 0 because of its fast initial step.  Its anchor is still excluded
+    from fitting — otherwise it would force A = 1 and negate the model's
+    purpose — and its re-expanded "pred" at the anchor stays C0, the value just
+    before the drop.
     """
     t = np.asarray(time, dtype=float)
     Ct = np.asarray(Ct, dtype=float)
@@ -458,6 +480,31 @@ def _post_first(p, se, at_bound, C0):
     }
 
 
+def _post_first_drop(p, se, at_bound, C0):
+    k, A = p
+    initial_drop_pct = 100.0 * (1.0 - A)
+    r0 = k * A * C0
+    if A > 0.5 and k > 0:
+        t_half = round(np.log(2.0 * A) / k, 4)
+    elif A <= 0.5:
+        t_half = 0.0  # half already removed by the drop
+    else:
+        t_half = float("nan")  # k == 0
+    return {
+        "label": f"k = {_fmt_sci(k)} min⁻¹, initial drop = {initial_drop_pct:.1f}%",
+        "t_half": t_half,
+        "k": k,
+        "k_se": se[0],
+        "A": A,
+        "A_se": se[1],
+        "initial_drop_pct": initial_drop_pct,
+        "col_k": "k_drop (1/min)",
+        "r0": r0,
+        "r0_se": None,
+        "at_bound": at_bound,
+    }
+
+
 def _post_second(p, se, at_bound, C0):
     k2 = p[0]
     r0 = k2 * C0**2
@@ -565,6 +612,16 @@ def _de_prepare(p, se):
 _MODEL_SPECS = [
     ("Zero-order", _zero_order, [1e-6], ([0], [np.inf]), 5000, False, None, _post_zero),
     ("Pseudo-first", _first_order, [0.01], ([0], [np.inf]), 5000, False, None, _post_first),
+    (
+        "Pseudo-first (initial drop)",
+        _first_order_drop,
+        [0.01, 0.9],
+        ([0, 0], [np.inf, 1.0]),
+        10000,
+        True,
+        None,
+        _post_first_drop,
+    ),
     ("Pseudo-second-order", _second_order, [1.0], ([0], [np.inf]), 5000, False, None, _post_second),
     (
         "Elovich",

catlab/metrics.py

@@ -59,6 +59,48 @@ def _arrhenius_ci(cov, n_T):
     return Ea_ci, lnA_ci, df
 
 
+def linear_intercept_check(x, y, expected, conf=0.95):
+    """Test whether a linearisation passes through its theoretical intercept.
+
+    Fits y = slope * x + intercept with scipy.stats.linregress and flags the
+    intercept as "deviating" when |intercept - expected| exceeds the confidence
+    half-width t(0.975, n-2) * intercept_stderr.  With fewer than three points
+    there is no valid intercept uncertainty, so deviates is returned False.
+
+    Returns {"intercept", "ci_half", "expected", "deviates", "r2"}; r2 is the
+    R² of the linear fit, so callers can ignore an intercept mismatch on a line
+    that does not fit the data anyway.
+    """
+    x = np.asarray(x, dtype=float)
+    y = np.asarray(y, dtype=float)
+    n = len(x)
+    if n < 3:
+        return {
+            "intercept": float("nan"),
+            "ci_half": float("nan"),
+            "expected": float(expected),
+            "deviates": False,
+            "r2": float("nan"),
+        }
+    reg = scipy_stats.linregress(x, y)
+    df = n - 2
+    t_crit = scipy_stats.t.ppf(1.0 - (1.0 - conf) / 2.0, df)
+    ci_half = t_crit * reg.intercept_stderr
+    if ci_half == 0 or not np.isfinite(ci_half):
+        # Zero residual: the intercept is known exactly, so a machine-epsilon
+        # offset must not count as a real deviation.
+        deviates = not np.isclose(reg.intercept, expected, rtol=1e-12, atol=1e-15)
+    else:
+        deviates = abs(reg.intercept - expected) > ci_half
+    return {
+        "intercept": float(reg.intercept),
+        "ci_half": float(ci_half),
+        "expected": float(expected),
+        "deviates": bool(deviates),
+        "r2": float(reg.rvalue**2),
+    }
+
+
 def c_to_user(Ct_mol, unit, mw):
     """Convert a concentration from mol/L to the user's display unit."""
     if unit == "ppmS":