diff --git a/CHANGELOG.md b/CHANGELOG.md index 9b042ddbb..56d5dec3b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -119,6 +119,23 @@ ### Added +- **Wald SPRT classification for CAT** (`fast_mlsirm.sprt_classify`; in + `mlsirm_core::exposure`). Single-cut binary-response sequential probability + ratio test: point hypotheses at `theta_cut -/+ delta`, cumulative binary + log-likelihood ratio under the D=1 logistic 3PL, and inclusive + first-crossing decisions against the log Wald boundaries + `A = ln((1-beta)/alpha)`, `B = ln(beta/(1-alpha))` -> `"above"`/`"below"`/ + `"continue"` with 1-based `n_used`; the full `llr_trace` is returned as an + offline diagnostic (entries past `n_used` are counterfactual replay + values). Verified against R catIrt `termSPRT.R`/`logLik.brm.R`/`p.brm.R` + and Thompson (2007, doi:10.7275/fq3r-zz60); Reckase (1983), Eggen (1999), + and Wald (1947) are cited as historical origins via Thompson (not directly + read). Log-likelihood ratios are computed in stable log space (softplus / + log-sigmoid), so extreme-but-valid parameters that saturate the response + probability to numerical 0/1 yield finite LLRs instead of errors. Pinned + 17-digit interior-crossing oracle, error-path and 500-rep Monte-Carlo + structural-invariant tests; 4 executed mutation kills (swapped boundaries, + dropped guessing floor, collapsed null hypothesis, off-by-one `n_used`). - **Owen-approximate posterior-predictive EPV item selection** (`fast_mlsirm.epv_select`; in `mlsirm_core::exposure`). Deliberately reduced scope of van der Linden's (1998, doi:10.1007/BF02294775) minimum diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 1120fee07..b3dd20d34 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -42,7 +42,8 @@ use mlsirm_core::exposure::{ a_stratified as core_a_stratified, ccat_select as core_ccat_select, epv_select as core_epv_select, kl_information as core_kl_information, kl_select as core_kl_select, owen_cat as core_owen_cat, owen_update as core_owen_update, - sympson_hetter as core_sympson_hetter, AStratifiedConfig, SympsonHetterConfig, + sprt_classify as core_sprt_classify, sympson_hetter as core_sympson_hetter, AStratifiedConfig, + SympsonHetterConfig, }; use mlsirm_core::facets::fit_facets as core_fit_facets; use mlsirm_core::factor::{ @@ -2430,6 +2431,57 @@ fn py_epv_select( Ok(out.into()) } +/// Single-cut binary-response Wald SPRT classification for CAT +/// (`mlsirm_core::exposure::sprt_classify`). D = 1 logistic 3PL; point +/// hypotheses at `theta_cut -/+ delta`; log Wald boundaries +/// A = ln((1-beta)/alpha), B = ln(beta/(1-alpha)) with inclusive +/// first-crossing decisions ("above"/"below"/"continue"). `llr_trace` +/// entries past `n_used` are offline counterfactual replay values. +/// +/// References (APA 7th; see the core module comment for read/not-read +/// source status): +/// Thompson, N. A. (2007). A practitioner's guide for variable-length +/// computerized classification testing. Practical Assessment, Research & +/// Evaluation, 12(1). https://doi.org/10.7275/fq3r-zz60 (READ) +/// Nydick, S. W. (2014). catIrt (R package). (READ: termSPRT.R, +/// logLik.brm.R, p.brm.R) +/// Eggen, T. J. H. M. (1999). Applied Psychological Measurement, 23(3), +/// 249-261. (NOT read; historical citation via Thompson) +/// Reckase, M. D. (1983). A procedure for decision making using tailored +/// testing. (NOT read; historical citation via Thompson) +/// Wald, A. (1947). Sequential analysis. Wiley. (NOT read; boundary forms +/// verified through the READ sources above) +#[pyfunction] +fn py_sprt_classify( + py: Python<'_>, + a: PyReadonlyArray1<'_, f64>, + b: PyReadonlyArray1<'_, f64>, + c: PyReadonlyArray1<'_, f64>, + responses: PyReadonlyArray1<'_, u8>, + theta_cut: f64, + delta: f64, + alpha: f64, + beta: f64, +) -> PyResult> { + let res = core_sprt_classify( + a.as_slice()?, + b.as_slice()?, + c.as_slice()?, + responses.as_slice()?, + theta_cut, + delta, + alpha, + beta, + ) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("decision", res.decision)?; + out.set_item("n_used", res.n_used)?; + out.set_item("llr", res.llr)?; + out.set_item("llr_trace", numpy::PyArray1::from_slice(py, &res.llr_trace))?; + Ok(out.into()) +} + /// Horn's parallel analysis for principal-component retention /// (`mlsirm_core::parallel`; oracle: CRAN paran 1.5.6, PCA path). `data` is /// a flattened row-major `n_persons * n_items` matrix; `centile` is 0 for @@ -6250,6 +6302,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(py_owen_cat, m)?)?; m.add_function(wrap_pyfunction!(py_ccat_select, m)?)?; m.add_function(wrap_pyfunction!(py_epv_select, m)?)?; + m.add_function(wrap_pyfunction!(py_sprt_classify, m)?)?; m.add_function(wrap_pyfunction!(guttman_lambdas, m)?)?; m.add_function(wrap_pyfunction!(tenberge_mu, m)?)?; m.add_function(wrap_pyfunction!(cronbach_alpha, m)?)?; diff --git a/crates/mlsirm-core/src/exposure.rs b/crates/mlsirm-core/src/exposure.rs index c4e0de42f..e0a0b4cdc 100644 --- a/crates/mlsirm-core/src/exposure.rs +++ b/crates/mlsirm-core/src/exposure.rs @@ -1403,3 +1403,192 @@ pub fn epv_select( predictive, }) } + +// ===================== Wald SPRT classification for CAT ===================== +// +// `sprt_classify` implements single-cut, binary-response SPRT classification +// (Wald's sequential probability ratio test applied to IRT classification +// testing). Two point hypotheses around the cut score, +// theta0 = theta_cut - delta, theta1 = theta_cut + delta, +// are compared through the cumulative binary log-likelihood ratio under the +// D = 1 logistic 3PL +// P_i(theta) = c_i + (1 - c_i) / (1 + exp(-a_i (theta - b_i))), +// LLR_k = sum_{i<=k} [ u_i ln(P_i(theta1)/P_i(theta0)) +// + (1 - u_i) ln((1 - P_i(theta1))/(1 - P_i(theta0))) ], +// against the log Wald boundaries +// A = ln((1 - beta) / alpha), B = ln(beta / (1 - alpha)). +// Responses are walked in order and the FIRST crossing decides (inclusive +// comparisons, matching catIrt): LLR_k >= A -> "above" with n_used = k; +// LLR_k <= B -> "below" with n_used = k; no crossing -> "continue" with +// n_used = len(responses). +// +// CITATION GOVERNANCE / SCOPE (adversarial spec review, sprt_spec_review.md): +// boundaries and the binary log-likelihood-ratio form were verified against +// READ sources: catIrt R/termSPRT.R + R/logLik.brm.R + R/p.brm.R (GitHub +// swnydick/catIrt) and Thompson (2007), p. 7. Reckase (1983) and Eggen +// (1999) are historical citations via Thompson and were NOT directly read. +// This function implements only a single-cut binary 3PL SPRT with D = 1 +// logistic-scale item parameters; it is not a multi-cut, polytomous, or +// D = 1.7 compatibility layer (parameters calibrated on the D = 1.7 metric +// must be rescaled by the caller, a_D1 = 1.7 * a_D17, before use). +// +// The returned decision/n_used are first-crossing SPRT results. llr_trace is +// computed for ALL supplied responses as an offline diagnostic; entries after +// n_used are counterfactual replay values - live CAT would terminate at +// n_used and would not administer later items. +// +// References (APA 7th): +// Wald, A. (1947). Sequential analysis. Wiley. (NOT read; boundary forms +// verified through the sources below) +// Thompson, N. A. (2007). A practitioner's guide for variable-length +// computerized classification testing. Practical Assessment, Research & +// Evaluation, 12(1). https://doi.org/10.7275/fq3r-zz60 (READ: p. 7 +// likelihood-ratio form and Wald decision points) +// Nydick, S. W. (2014). catIrt: An R package for simulating IRT-based +// computerized adaptive tests. (READ: R/termSPRT.R boundary and +// inclusive-comparison conventions; R/logLik.brm.R binary log +// likelihood; R/p.brm.R D = 1 3PL) +// Eggen, T. J. H. M. (1999). Item selection in adaptive testing with the +// sequential probability ratio test. Applied Psychological Measurement, +// 23(3), 249-261. (NOT read; historical citation via Thompson) +// Reckase, M. D. (1983). A procedure for decision making using tailored +// testing. (NOT read; historical citation via Thompson) + +/// Result of [`sprt_classify`]. `decision` is `"above"`, `"below"`, or +/// `"continue"`; `n_used` is the 1-based count of responses consumed by the +/// first boundary crossing (or all responses when no crossing occurs); +/// `llr_trace` holds the cumulative log-likelihood ratio after every supplied +/// response (entries past `n_used` are offline counterfactuals); `llr` is the +/// final trace entry. +#[derive(Debug, Clone)] +pub struct SprtResult { + pub decision: &'static str, + pub n_used: usize, + pub llr_trace: Vec, + pub llr: f64, +} + +/// Single-cut binary-response Wald SPRT classification (see module comment +/// above for the exact verified contract and source status). +pub fn sprt_classify( + a: &[f64], + b: &[f64], + c: &[f64], + responses: &[u8], + theta_cut: f64, + delta: f64, + alpha: f64, + beta: f64, +) -> Result { + let n = a.len(); + if n == 0 { + return Err("sprt_classify: item pool is empty".into()); + } + if b.len() != n || c.len() != n || responses.len() != n { + return Err(format!( + "sprt_classify: length mismatch (a: {}, b: {}, c: {}, responses: {})", + n, + b.len(), + c.len(), + responses.len() + )); + } + for i in 0..n { + if !a[i].is_finite() || a[i] <= 0.0 { + return Err(format!("sprt_classify: a[{i}] must be finite and > 0")); + } + if !b[i].is_finite() { + return Err(format!("sprt_classify: b[{i}] must be finite")); + } + if !c[i].is_finite() || !(0.0..1.0).contains(&c[i]) { + return Err(format!( + "sprt_classify: c[{i}] must be finite and in [0, 1)" + )); + } + if responses[i] > 1 { + return Err(format!("sprt_classify: responses[{i}] must be 0 or 1")); + } + } + if !theta_cut.is_finite() { + return Err("sprt_classify: theta_cut must be finite".into()); + } + if !delta.is_finite() || delta <= 0.0 { + return Err("sprt_classify: delta must be finite and > 0".into()); + } + for (name, v) in [("alpha", alpha), ("beta", beta)] { + if !v.is_finite() || v <= 0.0 || v >= 1.0 { + return Err(format!( + "sprt_classify: {name} must be finite and in (0, 1)" + )); + } + } + if alpha + beta >= 1.0 { + return Err("sprt_classify: alpha + beta must be < 1".into()); + } + + let upper = ((1.0 - beta) / alpha).ln(); + let lower = (beta / (1.0 - alpha)).ln(); + let theta0 = theta_cut - delta; + let theta1 = theta_cut + delta; + // Stable softplus ln(1 + e^z): shift by max(z, 0) so exp never overflows. + let softplus = |z: f64| -> f64 { + if z > 0.0 { + z + (-z).exp().ln_1p() + } else { + z.exp().ln_1p() + } + }; + // Stable log-probabilities under the D = 1 logistic 3PL + // P = c + (1 - c) sigmoid(z), z = a (theta - b) (crate CAT convention; + // catIrt p.brm.R). ln(1 - P) = ln(1 - c) - softplus(z) always; ln(P) + // needs the log-sigmoid branch -softplus(-z) only when c = 0 (for c > 0 + // the direct form is bounded below by c and stays finite). + let ln_p = |z: f64, ci: f64| -> f64 { + if ci > 0.0 { + (ci + (1.0 - ci) / (1.0 + (-z).exp())).ln() + } else { + -softplus(-z) + } + }; + + let mut llr_trace = Vec::with_capacity(n); + let mut cum = 0.0_f64; + let mut decision = "continue"; + let mut n_used = n; + for i in 0..n { + let z0 = a[i] * (theta0 - b[i]); + let z1 = a[i] * (theta1 - b[i]); + let inc = if responses[i] == 1 { + // ln(P(theta1)) - ln(P(theta0)), each log computed stably. + ln_p(z1, c[i]) - ln_p(z0, c[i]) + } else { + // ln(1-P(theta1)) - ln(1-P(theta0)); the ln(1-c) terms cancel. + softplus(z0) - softplus(z1) + }; + // Defensive: unreachable for validated inputs with the stable forms + // above (kept as a hard failure rather than silently propagating). + if !inc.is_finite() { + return Err(format!( + "sprt_classify: non-finite log-likelihood-ratio increment at item {i}" + )); + } + cum += inc; + llr_trace.push(cum); + // First crossing decides; inclusive comparisons (catIrt termSPRT.R). + if decision == "continue" { + if cum >= upper { + decision = "above"; + n_used = i + 1; + } else if cum <= lower { + decision = "below"; + n_used = i + 1; + } + } + } + Ok(SprtResult { + decision, + n_used, + llr: *llr_trace.last().unwrap(), + llr_trace, + }) +} diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index a0d606f57..e42a2e96a 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -97,6 +97,7 @@ owen_cat as owen_cat, ccat_select as ccat_select, epv_select as epv_select, + sprt_classify as sprt_classify, ) from .rsm import fit_rsm as fit_rsm, RsmFit as RsmFit from .mixed import fit_mixed_items as fit_mixed_items, MixedFormatFit as MixedFormatFit, MixedItemParameters as MixedItemParameters @@ -254,6 +255,7 @@ "owen_cat", "ccat_select", "epv_select", + "sprt_classify", "AStratifiedResult", "omega_total_1f", "omega_total_1f_from_data", diff --git a/python/fast_mlsirm/exposure.py b/python/fast_mlsirm/exposure.py index c685aee00..1065d1d90 100644 --- a/python/fast_mlsirm/exposure.py +++ b/python/fast_mlsirm/exposure.py @@ -619,4 +619,98 @@ def epv_select( "selected": int(r["selected"]), "epv": np.asarray(r["epv"]), "predictive": np.asarray(r["predictive"]), + } + +def sprt_classify( + a: np.ndarray, + b: np.ndarray, + c: np.ndarray | None = None, + *, + responses: np.ndarray, + theta_cut: float, + delta: float, + alpha: float = 0.05, + beta: float = 0.05, +) -> dict: + """Single-cut binary-response Wald SPRT classification for CAT. + + Compares the point hypotheses ``theta0 = theta_cut - delta`` and + ``theta1 = theta_cut + delta`` through the cumulative binary + log-likelihood ratio under the D = 1 logistic 3PL + ``P_i(theta) = c_i + (1 - c_i) / (1 + exp(-a_i (theta - b_i)))`` + against the log Wald boundaries ``A = ln((1 - beta) / alpha)`` and + ``B = ln(beta / (1 - alpha))``. Responses are walked in order and the + FIRST inclusive crossing decides: ``LLR_k >= A`` -> ``"above"``, + ``LLR_k <= B`` -> ``"below"`` (``n_used = k``, 1-based); no crossing -> + ``"continue"`` with ``n_used = len(responses)``. All numerics run in the + Rust core (``mlsirm_core::exposure::sprt_classify``). + + ``llr_trace`` is returned for ALL supplied responses as an offline + diagnostic; entries past ``n_used`` are counterfactual replay values (a + live CAT would stop at ``n_used`` and never administer later items). + Parameters calibrated on the D = 1.7 normal-ogive metric must be + rescaled by the caller (``a_D1 = 1.7 * a_D17``) before use. + + Source status: the boundary and likelihood-ratio forms were verified + against R catIrt ``termSPRT.R``/``logLik.brm.R``/``p.brm.R`` (READ) and + Thompson (2007, p. 7; READ). Reckase (1983), Eggen (1999), and Wald + (1947) were NOT read and are cited as historical origins via Thompson. + + References (APA 7th ed.): + Thompson, N. A. (2007). A practitioner's guide for variable-length + computerized classification testing. *Practical Assessment, + Research & Evaluation, 12*(1). + https://doi.org/10.7275/fq3r-zz60 + Eggen, T. J. H. M. (1999). Item selection in adaptive testing with + the sequential probability ratio test. *Applied Psychological + Measurement, 23*(3), 249-261. + https://doi.org/10.1177/01466219922031365 + Reckase, M. D. (1983). A procedure for decision making using + tailored testing. In D. J. Weiss (Ed.), *New horizons in + testing* (pp. 237-255). Academic Press. + Wald, A. (1947). *Sequential analysis*. Wiley. (NOT read; boundary + forms verified through the READ sources above.) + """ + from . import _core + + # Reject complex input BEFORE the dtype casts: the casts would silently + # discard imaginary parts (complex laundering). + for name, arr in (("a", a), ("b", b), ("c", c), ("responses", responses)): + if arr is not None and np.iscomplexobj(np.asarray(arr)): + raise ValueError(f"{name} must be real-valued") + a = np.asarray(a, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + if a.ndim != 1 or b.ndim != 1: + raise ValueError("a and b must be 1-D arrays") + if c is None: + c = np.zeros_like(a) + c = np.asarray(c, dtype=np.float64) + if c.ndim != 1: + raise ValueError("c must be a 1-D array") + # Validate responses BEFORE the uint8 cast (casts truncate/wrap). + resp = np.asarray(responses) + if resp.ndim != 1: + raise ValueError("responses must be a 1-D array") + if resp.dtype == np.bool_: + resp = resp.astype(np.uint8) + else: + resp_f = np.asarray(resp, dtype=np.float64) + if not np.all(np.isin(resp_f, (0.0, 1.0))): + raise ValueError("responses must contain only 0 and 1") + resp = resp_f.astype(np.uint8) + r = _core.py_sprt_classify( + np.ascontiguousarray(a), + np.ascontiguousarray(b), + np.ascontiguousarray(c), + np.ascontiguousarray(resp), + float(theta_cut), + float(delta), + float(alpha), + float(beta), + ) + return { + "decision": str(r["decision"]), + "n_used": int(r["n_used"]), + "llr": float(r["llr"]), + "llr_trace": np.asarray(r["llr_trace"]), } \ No newline at end of file diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 0f0b57a78..76b61c776 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -6609,4 +6609,95 @@ def test_error_paths(self): with pytest.raises(ValueError): epv_select(a, np.array([0.0 + 1j, 0.5]), c, administered=adm, mu=0.0, sig2=1.0) with pytest.raises(ValueError): - epv_select(a, b, np.array([0.0, 0.1 + 1j]), administered=adm, mu=0.0, sig2=1.0) \ No newline at end of file + epv_select(a, b, np.array([0.0, 0.1 + 1j]), administered=adm, mu=0.0, sig2=1.0) + +class TestSprtClassify: + """sprt_classify wrapper: every assert reads crate-returned dict values.""" + + def test_pinned_oracle(self): + from fast_mlsirm import sprt_classify + + r = sprt_classify( + np.full(5, 2.0), + np.zeros(5), + np.full(5, 0.1), + responses=np.array([1, 1, 1, 1, 0]), + theta_cut=0.0, + delta=0.5, + alpha=0.05, + beta=0.05, + ) + assert r["decision"] == "above" + assert r["n_used"] == 4 + trace = np.array( + [ + 0.79567203915954553, + 1.5913440783190911, + 2.3870161174786366, + 3.1826881566381821, + 2.1826881566381821, + ] + ) + np.testing.assert_allclose(r["llr_trace"], trace, rtol=0, atol=5e-15) + assert abs(r["llr"] - 2.1826881566381821) < 5e-15 + + def test_default_c_and_bool_responses(self): + from fast_mlsirm import sprt_classify + + r = sprt_classify( + np.array([1.0, 1.1]), + np.array([0.0, 0.2]), + responses=np.array([True, False]), + theta_cut=0.0, + delta=0.3, + ) + assert r["decision"] == "continue" + assert r["n_used"] == 2 + assert r["llr_trace"].shape == (2,) + assert r["llr"] == r["llr_trace"][-1] + + def test_rejects_bad_responses(self): + from fast_mlsirm import sprt_classify + + with pytest.raises(ValueError, match="0 and 1"): + sprt_classify( + np.array([1.0, 1.0]), + np.zeros(2), + responses=np.array([1, 2]), + theta_cut=0.0, + delta=0.5, + ) + with pytest.raises(ValueError, match="0 and 1"): + sprt_classify( + np.array([1.0, 1.0]), + np.zeros(2), + responses=np.array([1.0, 0.5]), + theta_cut=0.0, + delta=0.5, + ) + + def test_rejects_complex_input(self): + from fast_mlsirm import sprt_classify + + with pytest.raises(ValueError, match="real-valued"): + sprt_classify( + np.array([1.0 + 1j, 1.0]), + np.zeros(2), + responses=np.array([1, 0]), + theta_cut=0.0, + delta=0.5, + ) + + def test_core_validation_propagates(self): + from fast_mlsirm import sprt_classify + + with pytest.raises(ValueError, match="alpha \\+ beta"): + sprt_classify( + np.array([1.0, 1.0]), + np.zeros(2), + responses=np.array([1, 0]), + theta_cut=0.0, + delta=0.5, + alpha=0.6, + beta=0.5, + ) \ No newline at end of file diff --git a/tests/unit/exposure_tests.rs b/tests/unit/exposure_tests.rs index 408f80eaa..7fe620f2d 100644 --- a/tests/unit/exposure_tests.rs +++ b/tests/unit/exposure_tests.rs @@ -29,7 +29,7 @@ use crate::exposure::{ a_stratified, ccat_select, eap_interim, epv_select, kl_information, kl_select, owen_cat, - owen_update, p3pl, sympson_hetter, AStratifiedConfig, Lcg, SympsonHetterConfig, + owen_update, p3pl, sprt_classify, sympson_hetter, AStratifiedConfig, Lcg, SympsonHetterConfig, }; fn pool30() -> (Vec, Vec, Vec) { @@ -1633,3 +1633,209 @@ fn epv_mc500_invariants() { } } } + +/// Spec-review pinned oracle (sprt_spec_review.md): homogeneous 5-item pool, +/// responses [1,1,1,1,0], theta_cut = 0, delta = 0.5, alpha = beta = 0.05. +/// Every assert reads crate outputs (r.decision, r.n_used, r.llr_trace, +/// r.llr). Killing mutants: SPRT-M1 (swapped boundaries classify at k = 1), +/// SPRT-M2 (dropping c gives trace [1,2,3,4,3] and crossing at k = 3), +/// SPRT-M3 (theta_cut in one likelihood: no crossing), SPRT-M4 (off-by-one +/// n_used), SPRT-M5 (final-LLR-only decision returns "continue" because the +/// counterfactual final LLR 2.1826... is inside the band). +#[test] +fn sprt_pinned_oracle_interior_crossing() { + let a = [2.0; 5]; + let b = [0.0; 5]; + let c = [0.1; 5]; + let r = sprt_classify(&a, &b, &c, &[1, 1, 1, 1, 0], 0.0, 0.5, 0.05, 0.05).unwrap(); + assert_eq!(r.decision, "above"); + assert_eq!(r.n_used, 4); + let trace_oracle = [ + 0.79567203915954553, + 1.5913440783190911, + 2.3870161174786366, + 3.1826881566381821, + 2.1826881566381821, + ]; + assert_eq!(r.llr_trace.len(), 5); + for k in 0..5 { + assert!( + (r.llr_trace[k] - trace_oracle[k]).abs() < 5e-15, + "llr_trace[{k}] = {}, oracle {}", + r.llr_trace[k], + trace_oracle[k] + ); + } + assert!( + (r.llr - 2.1826881566381821).abs() < 5e-15, + "final llr = {}", + r.llr + ); + // The final LLR sits strictly inside (B, A): a final-LLR-only mutant + // (SPRT-M5) would say "continue"; first-crossing SPRT says "above". + let upper = (0.95_f64 / 0.05).ln(); + assert!(r.llr < upper && r.llr > -upper); + assert!(r.llr_trace[3] >= upper); +} + +/// Below-decision symmetry and heterogeneous-pool behavior, reading only +/// crate outputs. All-wrong responses on an informative pool must cross the +/// lower boundary before the pool is exhausted, and n_used must mark the +/// first crossing (every earlier trace entry strictly inside the band). +#[test] +fn sprt_below_decision_first_crossing() { + let a = [1.8, 1.2, 2.0, 0.9, 1.5, 1.6, 1.1, 1.9]; + let b = [-0.3, 0.4, 0.1, -0.8, 0.6, -0.1, 0.9, 0.2]; + let c = [0.0, 0.05, 0.2, 0.1, 0.0, 0.15, 0.0, 0.25]; + let r = sprt_classify(&a, &b, &c, &[0; 8], 0.0, 0.5, 0.05, 0.05).unwrap(); + assert_eq!(r.decision, "below"); + let lower = (0.05_f64 / 0.95).ln(); + let upper = (0.95_f64 / 0.05).ln(); + assert!(r.n_used >= 1 && r.n_used < 8, "n_used = {}", r.n_used); + assert!(r.llr_trace[r.n_used - 1] <= lower); + for k in 0..r.n_used - 1 { + assert!( + r.llr_trace[k] > lower && r.llr_trace[k] < upper, + "premature crossing at {k}" + ); + } +} + +/// No-crossing path: wide error rates push the boundaries out so a short, +/// mixed response set stays inside the band -> "continue" with n_used = n. +#[test] +fn sprt_continue_when_no_crossing() { + let a = [1.0, 1.1]; + let b = [0.0, 0.2]; + let c = [0.0, 0.0]; + let r = sprt_classify(&a, &b, &c, &[1, 0], 0.0, 0.3, 0.05, 0.05).unwrap(); + assert_eq!(r.decision, "continue"); + assert_eq!(r.n_used, 2); + assert_eq!(r.llr_trace.len(), 2); + assert_eq!(r.llr, r.llr_trace[1]); +} + +#[test] +fn sprt_error_paths() { + let a = [1.0, 1.2]; + let b = [0.0, 0.5]; + let c = [0.0, 0.1]; + let u = [1u8, 0u8]; + assert!(sprt_classify(&[], &[], &[], &[], 0.0, 0.5, 0.05, 0.05) + .unwrap_err() + .contains("empty")); + assert!(sprt_classify(&a, &b[..1], &c, &u, 0.0, 0.5, 0.05, 0.05) + .unwrap_err() + .contains("length mismatch")); + assert!( + sprt_classify(&[1.0, -0.5], &b, &c, &u, 0.0, 0.5, 0.05, 0.05) + .unwrap_err() + .contains("a[1]") + ); + assert!( + sprt_classify(&a, &[0.0, f64::NAN], &c, &u, 0.0, 0.5, 0.05, 0.05) + .unwrap_err() + .contains("b[1]") + ); + assert!(sprt_classify(&a, &b, &[0.0, 1.0], &u, 0.0, 0.5, 0.05, 0.05) + .unwrap_err() + .contains("c[1]")); + assert!(sprt_classify(&a, &b, &c, &[1, 2], 0.0, 0.5, 0.05, 0.05) + .unwrap_err() + .contains("responses[1]")); + assert!( + sprt_classify(&a, &b, &c, &u, f64::INFINITY, 0.5, 0.05, 0.05) + .unwrap_err() + .contains("theta_cut") + ); + assert!(sprt_classify(&a, &b, &c, &u, 0.0, 0.0, 0.05, 0.05) + .unwrap_err() + .contains("delta")); + assert!(sprt_classify(&a, &b, &c, &u, 0.0, 0.5, 0.0, 0.05) + .unwrap_err() + .contains("alpha")); + assert!(sprt_classify(&a, &b, &c, &u, 0.0, 0.5, 0.05, 1.0) + .unwrap_err() + .contains("beta")); + assert!(sprt_classify(&a, &b, &c, &u, 0.0, 0.5, 0.6, 0.5) + .unwrap_err() + .contains("alpha + beta")); +} + +/// 500-rep Monte-Carlo structural invariants, reading only crate outputs: +/// trace finite and same length as the pool; decision consistent with the +/// first crossing of the inclusive Wald boundaries (no earlier crossing +/// before n_used; crossing entry beyond the matching boundary; "continue" +/// iff no entry ever leaves the open band, with n_used = n). +#[test] +#[ignore = "slow Monte-Carlo suite; run explicitly with --ignored"] +fn sprt_mc500_invariants() { + let mut lcg = Lcg(20260725); + for rep in 0..500 { + let n = 5 + (lcg.next_f64() * 26.0) as usize; + let mut a = Vec::with_capacity(n); + let mut b = Vec::with_capacity(n); + let mut c = Vec::with_capacity(n); + let mut u = Vec::with_capacity(n); + for _ in 0..n { + a.push(0.5 + 2.0 * lcg.next_f64()); + b.push(-2.0 + 4.0 * lcg.next_f64()); + c.push(0.3 * lcg.next_f64()); + u.push((lcg.next_f64() < 0.5) as u8); + } + let delta = 0.2 + 0.6 * lcg.next_f64(); + let alpha = 0.01 + 0.3 * lcg.next_f64(); + let beta = 0.01 + 0.3 * lcg.next_f64(); + let cut = -1.0 + 2.0 * lcg.next_f64(); + let r = sprt_classify(&a, &b, &c, &u, cut, delta, alpha, beta) + .unwrap_or_else(|e| panic!("rep {rep}: {e}")); + let upper = ((1.0 - beta) / alpha).ln(); + let lower = (beta / (1.0 - alpha)).ln(); + assert_eq!(r.llr_trace.len(), n, "rep {rep}"); + assert!(r.llr_trace.iter().all(|v| v.is_finite()), "rep {rep}"); + assert!((r.llr - r.llr_trace[n - 1]).abs() == 0.0, "rep {rep}"); + for k in 0..r.n_used.saturating_sub(1) { + assert!( + r.llr_trace[k] > lower && r.llr_trace[k] < upper, + "rep {rep}: crossing before n_used at {k}" + ); + } + match r.decision { + "above" => assert!(r.llr_trace[r.n_used - 1] >= upper, "rep {rep}"), + "below" => assert!(r.llr_trace[r.n_used - 1] <= lower, "rep {rep}"), + "continue" => { + assert_eq!(r.n_used, n, "rep {rep}"); + assert!( + r.llr_trace.iter().all(|v| *v > lower && *v < upper), + "rep {rep}" + ); + } + other => panic!("rep {rep}: unexpected decision {other}"), + } + } +} + +/// Regression for the saturation defect (impl review): extreme but valid +/// parameters (a = 50, delta = 20 -> z = -/+1000) drive the 2PL response +/// probability to numerical 0/1. The stable log-space forms must return a +/// finite LLR (+1000 for a correct response, -1000 for a wrong one, since +/// ln P -> -softplus(-z) and ln(1-P) -> -softplus(z)), not Err. Asserts +/// read crate outputs (decision, n_used, llr, llr_trace). +#[test] +fn sprt_extreme_parameters_stay_finite() { + let r = sprt_classify(&[50.0], &[0.0], &[0.0], &[1], 0.0, 20.0, 0.05, 0.05).unwrap(); + assert_eq!(r.decision, "above"); + assert_eq!(r.n_used, 1); + assert!((r.llr - 1000.0).abs() < 1e-9, "llr = {}", r.llr); + assert_eq!(r.llr_trace.len(), 1); + let r0 = sprt_classify(&[50.0], &[0.0], &[0.0], &[0], 0.0, 20.0, 0.05, 0.05).unwrap(); + assert_eq!(r0.decision, "below"); + assert_eq!(r0.n_used, 1); + assert!((r0.llr + 1000.0).abs() < 1e-9, "llr = {}", r0.llr); + // Nonzero guessing floor at the same extremity: ln P is bounded below + // by ln(c), so the correct-response increment is ln(~1) - ln(c) = + // -ln(0.2) exactly at saturation. + let rc = sprt_classify(&[50.0], &[0.0], &[0.2], &[1], 0.0, 20.0, 0.05, 0.05).unwrap(); + assert!((rc.llr - (-0.2_f64.ln())).abs() < 1e-12, "llr = {}", rc.llr); + assert!(rc.llr.is_finite()); +}