diff --git a/CHANGELOG.md b/CHANGELOG.md index 35e8a62..f5eb291 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index 1c1b91f..1c200df 100644 --- a/README.md +++ b/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 | diff --git a/app_ods.py b/app_ods.py index 43ed5f7..469f83a 100644 --- a/app_ods.py +++ b/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 diff --git a/catlab/kinetics_engine.py b/catlab/kinetics_engine.py index d07d9ab..8ccd0cf 100644 --- a/catlab/kinetics_engine.py +++ b/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", diff --git a/catlab/metrics.py b/catlab/metrics.py index 35e3bcb..51aff86 100644 --- a/catlab/metrics.py +++ b/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": diff --git a/tests/data/golden_fit.json b/tests/data/golden_fit.json index d9ca223..258ad5b 100644 --- a/tests/data/golden_fit.json +++ b/tests/data/golden_fit.json @@ -67,6 +67,39 @@ "n_fit": 6, "n_anchor": 1 }, + "Pseudo-first (initial drop)": { + "params": [ + 0.019858028474994956, + 0.9983140427907604, + 0.007798 + ], + "R2": 0.9995, + "pred": [ + 0.007798, + 0.005779455553785918, + 0.004290653516883525, + 0.0031853705648597647, + 0.0023648112334283505, + 0.0013033753827344358, + 0.0007183606726425031 + ], + "adj_r2": 0.9992, + "aic": -115.9595, + "aicc": -103.9595, + "label": "k = 1.99 \u00d7 10\u207b\u00b2 min\u207b\u00b9, initial drop = 0.2%", + "t_half": 34.8202, + "k": 0.019858028474994956, + "k_se": 0.0002980204435581425, + "A": 0.9983140427907604, + "A_se": 0.010443044632608181, + "initial_drop_pct": 0.16859572092395858, + "col_k": "k_drop (1/min)", + "r0": 0.00015459183067468733, + "r0_se": null, + "at_bound": [], + "n_fit": 6, + "n_anchor": 1 + }, "Pseudo-second-order": { "params": [ 4.390689650513702, @@ -309,6 +342,39 @@ "n_fit": 6, "n_anchor": 1 }, + "Pseudo-first (initial drop)": { + "params": [ + 0.012923766044946832, + 0.7996660408517524, + 0.007798 + ], + "R2": 0.9623, + "pred": [ + 0.007798, + 0.0051368996291144085, + 0.0042316552213689995, + 0.003485936499722233, + 0.002871631228066929, + 0.001948709306081273, + 0.0013224079479606664 + ], + "adj_r2": 0.9372, + "aic": -93.6675, + "aicc": -81.6675, + "label": "k = 1.29 \u00d7 10\u207b\u00b2 min\u207b\u00b9, initial drop = 20.0%", + "t_half": 36.3351, + "k": 0.012923766044946832, + "k_se": 0.001486996468127161, + "A": 0.7996660408517524, + "A_se": 0.05160495214651284, + "initial_drop_pct": 20.03339591482476, + "col_k": "k_drop (1/min)", + "r0": 8.058996584959205e-05, + "r0_se": null, + "at_bound": [], + "n_fit": 6, + "n_anchor": 1 + }, "Pseudo-second-order": { "params": [ 3.9551552169669417, @@ -551,6 +617,41 @@ "n_fit": 6, "n_anchor": 1 }, + "Pseudo-first (initial drop)": { + "params": [ + 0.011786219109468513, + 0.9999995133128385, + 0.007798 + ], + "R2": 0.9754, + "pred": [ + 0.007798, + 0.0065343541788760815, + 0.00547548157418171, + 0.004588196116782911, + 0.003844692621252743, + 0.0026996093679895823, + 0.0018955717550607329 + ], + "adj_r2": 0.959, + "aic": -92.3204, + "aicc": -80.3204, + "label": "k = 1.18 \u00d7 10\u207b\u00b2 min\u207b\u00b9, initial drop = 0.0%", + "t_half": 58.8099, + "k": 0.011786219109468513, + "k_se": NaN, + "A": 0.9999995133128385, + "A_se": NaN, + "initial_drop_pct": 4.866871614606083e-05, + "col_k": "k_drop (1/min)", + "r0": 9.190889188473599e-05, + "r0_se": null, + "at_bound": [ + "A at upper bound 1" + ], + "n_fit": 6, + "n_anchor": 1 + }, "Pseudo-second-order": { "params": [ 2.2221593816454344, @@ -793,6 +894,39 @@ "n_fit": 6, "n_anchor": 1 }, + "Pseudo-first (initial drop)": { + "params": [ + 0.00503490110441002, + 0.8693007994474234, + 0.007798 + ], + "R2": 0.9853, + "pred": [ + 0.007798, + 0.006285703096390955, + 0.005828468006273035, + 0.005404493145667892, + 0.005011359096444348, + 0.004308802924968661, + 0.0037047400293846945 + ], + "adj_r2": 0.9755, + "aic": -103.6845, + "aicc": -91.6845, + "label": "k = 5.03 \u00d7 10\u207b\u00b3 min\u207b\u00b9, initial drop = 13.1%", + "t_half": 109.8494, + "k": 0.00503490110441002, + "k_se": 0.00032655964621225507, + "A": 0.8693007994474234, + "A_se": 0.016441810860043617, + "initial_drop_pct": 13.069920055257656, + "col_k": "k_drop (1/min)", + "r0": 3.4130626043467894e-05, + "r0_se": null, + "at_bound": [], + "n_fit": 6, + "n_anchor": 1 + }, "Pseudo-second-order": { "params": [ 1.2539472062569104, @@ -1035,6 +1169,41 @@ "n_fit": 6, "n_anchor": 1 }, + "Pseudo-first (initial drop)": { + "params": [ + 0.06119166186304211, + 0.9999999999999999, + 0.007798 + ], + "R2": 0.8842, + "pred": [ + 0.007798, + 0.0031142624876403214, + 0.0012437331164303267, + 0.0004967057436695253, + 0.0001983677949353029, + 3.16384452159908e-05, + 5.046137736277814e-06 + ], + "adj_r2": 0.807, + "aic": -85.4228, + "aicc": -73.4228, + "label": "k = 6.12 \u00d7 10\u207b\u00b2 min\u207b\u00b9, initial drop = 0.0%", + "t_half": 11.3275, + "k": 0.06119166186304211, + "k_se": NaN, + "A": 0.9999999999999999, + "A_se": NaN, + "initial_drop_pct": 1.1102230246251565e-14, + "col_k": "k_drop (1/min)", + "r0": 0.00047717257920800235, + "r0_se": null, + "at_bound": [ + "A at upper bound 1" + ], + "n_fit": 6, + "n_anchor": 1 + }, "Pseudo-second-order": { "params": [ 20.04313929577164, @@ -1277,6 +1446,39 @@ "n_fit": 6, "n_anchor": 1 }, + "Pseudo-first (initial drop)": { + "params": [ + 0.011580987580891917, + 0.5024772109957858, + 0.007798 + ], + "R2": 0.8912, + "pred": [ + 0.007798, + 0.0032934889879341995, + 0.002768297947081283, + 0.0023268556694404697, + 0.0019558072902216696, + 0.0013817812351149773, + 0.0009762308338156723 + ], + "adj_r2": 0.8187, + "aic": -93.0703, + "aicc": -81.0703, + "label": "k = 1.16 \u00d7 10\u207b\u00b2 min\u207b\u00b9, initial drop = 49.8%", + "t_half": 0.4268, + "k": 0.011580987580891917, + "k_se": 0.0022543204501218016, + "A": 0.5024772109957858, + "A_se": 0.05149353492571532, + "initial_drop_pct": 49.75227890042142, + "col_k": "k_drop (1/min)", + "r0": 4.537798388906209e-05, + "r0_se": null, + "at_bound": [], + "n_fit": 6, + "n_anchor": 1 + }, "Pseudo-second-order": { "params": [ 8.180347247509392, @@ -1521,6 +1723,42 @@ "n_fit": 7, "n_anchor": 1 }, + "Pseudo-first (initial drop)": { + "params": [ + 0.019995399507290342, + 0.9999999999999999, + 0.007798 + ], + "R2": 0.9916, + "pred": [ + 0.007798, + 0.005777299147459686, + 0.004280223831654067, + 0.003171086623948663, + 0.0023493608680507086, + 0.0012895345443172776, + 0.0007078092444637061, + 0.00021324690189632837 + ], + "adj_r2": 0.9874, + "aic": -115.2364, + "aicc": -107.2364, + "label": "k = 2.00 \u00d7 10\u207b\u00b2 min\u207b\u00b9, initial drop = 0.0%", + "t_half": 34.6653, + "k": 0.019995399507290342, + "k_se": NaN, + "A": 0.9999999999999999, + "A_se": NaN, + "initial_drop_pct": 1.1102230246251565e-14, + "col_k": "k_drop (1/min)", + "r0": 0.00015592412535785007, + "r0_se": null, + "at_bound": [ + "A at upper bound 1" + ], + "n_fit": 7, + "n_anchor": 1 + }, "Pseudo-second-order": { "params": [ 4.584556326059062, diff --git a/tests/synthetic_data.py b/tests/synthetic_data.py new file mode 100644 index 0000000..38f4d5e --- /dev/null +++ b/tests/synthetic_data.py @@ -0,0 +1,15 @@ +# tests/synthetic_data.py +# Synthetic kinetic datasets for tests. These values are computer-generated +# archetypes, NOT real experimental measurements, so they are safe to commit to +# a public repository. + +# Sparse 5-point time grid in minutes (no t=0 measurement). +T_SPARSE = [20, 40, 80, 180, 240] + +# Removal (%) per archetype, keyed by a short descriptive name. +SYNTHETIC_REMOVAL = { + "A_first_order": [21.9, 39.5, 61.4, 91.3, 95.9], # clean PFO + "B_initial_drop": [35.7, 45.7, 61.8, 83.2, 89.8], # ~24 % fast drop, then PFO + "C_plateau": [21.8, 26.0, 34.3, 47.1, 53.1], # slows early, levels off + "D_dark_like": [7.9, 13.2, 19.1, 25.9, 28.1], # slow adsorption control +} diff --git a/tests/test_app_ui.py b/tests/test_app_ui.py index 2694248..29e0625 100644 --- a/tests/test_app_ui.py +++ b/tests/test_app_ui.py @@ -14,6 +14,7 @@ from streamlit.testing.v1 import AppTest from catlab.kinetics_engine import MODEL_NAMES +from tests.synthetic_data import SYNTHETIC_REMOVAL, T_SPARSE pytestmark = pytest.mark.ui @@ -149,6 +150,27 @@ def test_t0_row_does_not_change_fit(self, app, csv_bytes): for cat in r2_with.index: assert np.isclose(r2_with[cat], r2_without[cat], rtol=1e-9) + def test_five_point_dataset_warns_too_few_points(self, app): + rem = SYNTHETIC_REMOVAL["B_initial_drop"] + rows = ["Time (min),CatA Removal (%)"] + rows += [f"{t},{r}" for t, r in zip(T_SPARSE, rem)] + data = ("\n".join(rows) + "\n").encode("utf-8") + _upload_csv(app, data) + _run_analysis(app) + assert any("informative data points" in w.value for w in app.tabs[0].warning) + + def test_seven_point_dataset_has_no_too_few_points_warning(self, app, csv_bytes): + _upload_csv(app, csv_bytes) + _run_analysis(app) + assert not any("informative data points" in w.value for w in app.tabs[0].warning) + + def test_summary_has_n_fitted_column(self, app, csv_bytes): + _upload_csv(app, csv_bytes) + _run_analysis(app) + summary = _summary(app) + assert summary is not None + assert "n (fitted)" in summary.columns + # ================================================================ # STEP 3 — Tabs 2-9 @@ -169,6 +191,15 @@ def test_linearization_reports_r2_for_each_catalyst(self, app, csv_bytes): return pytest.fail("linearization R² pivot table not found") + def test_initial_drop_dataset_warns_pfo_intercept(self, app): + rem = SYNTHETIC_REMOVAL["B_initial_drop"] + rows = ["Time (min),CatA Removal (%)"] + rows += [f"{t},{r}" for t, r in zip(T_SPARSE, rem)] + data = ("\n".join(rows) + "\n").encode("utf-8") + _upload_csv(app, data) + warnings = [w.value for w in app.tabs[1].warning] + assert any("does not pass through the origin" in w for w in warnings) + class TestTabRemoval: def test_removal_renders_efficiency_plots(self, app, csv_bytes): diff --git a/tests/test_golden_master.py b/tests/test_golden_master.py index 8d0b71f..2f85d8b 100644 --- a/tests/test_golden_master.py +++ b/tests/test_golden_master.py @@ -66,6 +66,9 @@ def _compute_results(): STAT_KEYS = {"R2", "adj_r2", "aic", "aicc", "n_fit", "n_anchor"} FRAGILE_KEYS = { "params", + "A", + "A_se", + "initial_drop_pct", "k", "k_se", "k2", diff --git a/tests/test_initial_drop.py b/tests/test_initial_drop.py new file mode 100644 index 0000000..6f060d5 --- /dev/null +++ b/tests/test_initial_drop.py @@ -0,0 +1,93 @@ +# tests/test_initial_drop.py +# Tests for the "Pseudo-first (initial drop)" model: C(t) = A·C0·exp(-k·t). +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +import numpy as np + +from catlab.kinetics_engine import ( + MODEL_NAMES, + _best_model, + _fit_nonlinear, + _first_order_drop, + _post_first_drop, +) + +C0 = 7.798e-3 # mol/L (~250 ppmS) +T = np.array([10, 20, 40, 60, 90, 120, 180, 240.0]) +DROP_MODEL = "Pseudo-first (initial drop)" + + +def _two_stage_data(seed): + rng = np.random.default_rng(seed) + clean = 0.75 * C0 * np.exp(-0.009 * T) + return clean * (1.0 + 0.02 * rng.standard_normal(len(T))) + + +def _pure_pfo_data(seed): + rng = np.random.default_rng(seed) + clean = C0 * np.exp(-0.009 * T) + return clean * (1.0 + 0.02 * rng.standard_normal(len(T))) + + +def test_first_order_drop_shape(): + # At t=0 the model returns A·C0 (value just after the drop), not C0. + assert np.isclose(_first_order_drop(np.array([0.0]), 0.01, 0.75, C0)[0], 0.75 * C0) + # For t > 0 it decays as A·C0·exp(-k t). + t = np.array([60.0]) + assert np.isclose(_first_order_drop(t, 0.01, 0.75, C0)[0], 0.75 * C0 * np.exp(-0.01 * 60.0)) + + +def test_two_stage_data_selects_initial_drop(): + wins = 0 + for seed in range(50): + res = _fit_nonlinear(T, _two_stage_data(seed), C0) + if _best_model(res, MODEL_NAMES) == DROP_MODEL: + wins += 1 + assert wins >= 45, f"initial-drop model selected only {wins}/50 times" + + +def test_pure_pfo_data_rarely_selects_initial_drop(): + wins = 0 + for seed in range(50): + res = _fit_nonlinear(T, _pure_pfo_data(seed), C0) + if _best_model(res, MODEL_NAMES) == DROP_MODEL: + wins += 1 + assert wins <= 5, f"initial-drop model selected {wins}/50 times on pure PFO data" + + +def test_noise_free_pfo_flags_A_at_upper_bound(): + res = _fit_nonlinear(T, C0 * np.exp(-0.009 * T), C0) + mr = res[DROP_MODEL] + assert mr.get("at_bound") + assert any("A at upper bound" in b for b in mr["at_bound"]) + + +def test_fitted_A_recovers_drop_fraction(): + A_values = [] + drop_pcts = [] + for seed in range(50): + res = _fit_nonlinear(T, _two_stage_data(seed), C0) + mr = res[DROP_MODEL] + A_values.append(mr["A"]) + drop_pcts.append(mr["initial_drop_pct"]) + A_values = np.array(A_values) + drop_pcts = np.array(drop_pcts) + # Fitted A is within 0.05 of 0.75 for most seeds. + assert np.sum(np.abs(A_values - 0.75) < 0.05) >= 45 + # initial_drop_pct ≈ 100 * (1 - 0.75) = 25. + assert abs(float(np.mean(drop_pcts)) - 25.0) < 1.0 + + +def test_t_half_low_A_is_zero(): + # A <= 0.5: half already removed by the drop. + assert _post_first_drop([0.01, 0.4], [None, None], [], C0)["t_half"] == 0.0 + + +def test_t_half_high_A_is_ln2A_over_k(): + k, A = 0.01, 0.9 + expected = np.log(2.0 * A) / k + r = _post_first_drop([k, A], [None, None], [], C0) + assert np.isclose(r["t_half"], expected) diff --git a/tests/test_metrics.py b/tests/test_metrics.py index f4bc673..fc81295 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -14,6 +14,7 @@ _initial_tof_mass, _arrhenius_ci, c_to_user, + linear_intercept_check, ) MW_S = 32.06 @@ -110,3 +111,66 @@ def test_c_to_user_mmol_and_gL(): def test_c_to_user_molL_passthrough(): assert c_to_user(0.001, "mol/L", None) == 0.001 + + +# ---- linear_intercept_check --------------------------------------------- +def test_intercept_check_exact_line_no_deviation(): + x = np.array([10.0, 20.0, 40.0, 80.0]) + y = 2.0 * x + 5.0 # exact line, intercept 5 == expected + chk = linear_intercept_check(x, y, expected=5.0) + assert chk["deviates"] is False + assert np.isclose(chk["intercept"], 5.0) + + +def test_intercept_check_offset_intercept_deviates(): + rng = np.random.default_rng(0) + x = np.linspace(1.0, 10.0, 20) + y = 2.0 * x + 0.30 + rng.normal(0, 1e-3, x.size) + chk = linear_intercept_check(x, y, expected=0.0) + assert chk["deviates"] is True + assert chk["intercept"] > 0.0 + + +def test_intercept_check_two_points_no_deviation(): + x = np.array([1.0, 2.0]) + y = np.array([0.3, 0.5]) # any values; n < 3 must give deviates False + chk = linear_intercept_check(x, y, expected=0.0) + assert chk["deviates"] is False + assert np.isnan(chk["intercept"]) + + +def test_intercept_check_real_first_order_not_deviating(): + # Synthetic first-order removal (no fast initial step): the ln(C0/C) line + # passes through the origin, so the PFO intercept must not deviate. + t = np.array([20.0, 40.0, 80.0, 180.0, 240.0]) + C0 = 7.798e-3 + k = 0.008 + Ct = C0 * np.exp(-k * t) + y = np.log(C0 / Ct) + chk = linear_intercept_check(t, y, expected=0.0) + assert chk["deviates"] is False + + +def test_intercept_check_initial_drop_deviates(): + # Synthetic data with a fast initial step (~24 % removed instantly) followed + # by first-order decay: the ln(C0/C) line does not pass through the origin. + t = np.array([20.0, 40.0, 80.0, 180.0, 240.0]) + C0 = 7.798e-3 + A = 0.76 + k = 0.006 + Ct = A * C0 * np.exp(-k * t) + y = np.log(C0 / Ct) + chk = linear_intercept_check(t, y, expected=0.0) + assert chk["deviates"] is True + assert chk["intercept"] > 0.0 + + +def test_intercept_check_reports_r2(): + from catlab.metrics import linear_intercept_check + + x = np.array([10.0, 20.0, 40.0, 80.0, 160.0]) + straight = linear_intercept_check(x, 0.01 * x, 0.0) + curved = linear_intercept_check(x, np.sqrt(x), 0.0) + assert straight["r2"] > 0.999 + assert curved["r2"] < straight["r2"] + assert np.isnan(linear_intercept_check(x[:2], x[:2], 0.0)["r2"])