From ddf3cd33df8f2606882619c14cc1c5b57b332fdc Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 19:14:33 +0900 Subject: [PATCH 01/59] feat(irt): expose Rust polytomous predictions --- crates/fast-mlsirm-py/src/lib.rs | 28 ++++++++ crates/mlsirm-core/src/poly.rs | 64 +++++++++++++++++++ ...6-public-polytomous-prediction-boundary.md | 46 +++++++++++++ python/fast_mlsirm/_legacy_init.py | 4 +- python/fast_mlsirm/polytomous.py | 60 +++++++++++++++++ tests/test_polytomous_predictions.py | 36 +++++++++++ 6 files changed, 237 insertions(+), 1 deletion(-) create mode 100644 docs/adr/0016-public-polytomous-prediction-boundary.md create mode 100644 tests/test_polytomous_predictions.py diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index d324e8688..a77f50663 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -6299,6 +6299,33 @@ fn grm_cell_logprobs(base: f64, thresholds: PyReadonlyArray1<'_, f64>) -> PyResu Ok(core_grm_logprobs(base, thresholds.as_slice()?)) } +/// Batched public GRM/GPCM category probabilities and expected category scores. +#[pyfunction] +#[pyo3(signature = (theta, slope, cat_params, n_items, n_cat, model))] +fn polytomous_predictions( + py: Python<'_>, + theta: PyReadonlyArray1<'_, f64>, + slope: PyReadonlyArray1<'_, f64>, + cat_params: PyReadonlyArray1<'_, f64>, + n_items: usize, + n_cat: usize, + model: &str, +) -> PyResult> { + let result = mlsirm_core::poly::polytomous_predictions( + theta.as_slice()?, + slope.as_slice()?, + cat_params.as_slice()?, + n_items, + n_cat, + parse_poly_model(model)?, + ) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("probabilities", PyArray1::from_vec(py, result.probabilities))?; + out.set_item("expected", PyArray1::from_vec(py, result.expected))?; + Ok(out.into()) +} + /// Unidimensional polytomous marginal-EM fit (Rust compute path). `model` is /// "grm" (default) or "gpcm"; `y` holds integer categories `0..n_cat-1`. #[pyfunction] @@ -9321,6 +9348,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(empirical_reliability, m)?)?; m.add_function(wrap_pyfunction!(gpcm_cell_logprobs, m)?)?; m.add_function(wrap_pyfunction!(grm_cell_logprobs, m)?)?; + m.add_function(wrap_pyfunction!(polytomous_predictions, m)?)?; m.add_function(wrap_pyfunction!(fit_poly_unidim, m)?)?; m.add_function(wrap_pyfunction!(fit_nominal, m)?)?; m.add_function(wrap_pyfunction!(poly_person_fit, m)?)?; diff --git a/crates/mlsirm-core/src/poly.rs b/crates/mlsirm-core/src/poly.rs index e0c48d9ed..b85520dc5 100644 --- a/crates/mlsirm-core/src/poly.rs +++ b/crates/mlsirm-core/src/poly.rs @@ -196,6 +196,70 @@ pub enum PolyModel { Gpcm, } +/// Batched category probabilities and expected integer category scores for a +/// fitted unidimensional GRM or GPCM item bank. +#[derive(Debug, Clone, PartialEq)] +pub struct PolytomousPredictions { + pub probabilities: Vec, + pub expected: Vec, +} + +pub fn polytomous_predictions( + theta: &[f64], + slope: &[f64], + cat_params: &[f64], + n_items: usize, + n_cat: usize, + model: PolyModel, +) -> Result { + if theta.is_empty() || theta.iter().any(|value| !value.is_finite()) { + return Err("theta must be a non-empty finite vector".into()); + } + if n_items == 0 || n_cat < 2 { + return Err("n_items must be positive and n_cat must be at least 2".into()); + } + validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?; + let cells = crate::checked_mul_usize( + theta.len(), + n_items, + "theta count * n_items exceeds the prediction buffer size", + )?; + let probability_cells = crate::checked_mul_usize( + cells, + n_cat, + "prediction cells * n_cat exceeds the probability buffer size", + )?; + let mut probabilities = Vec::with_capacity(probability_cells); + let mut expected = Vec::with_capacity(cells); + let scores: Vec = (0..n_cat).map(|category| category as f64).collect(); + for &person_theta in theta { + for item in 0..n_items { + let base = slope[item] * person_theta; + let params = &cat_params[item * (n_cat - 1)..(item + 1) * (n_cat - 1)]; + let logprobs = match model { + PolyModel::Grm => grm_logprobs(base, params), + PolyModel::Gpcm => { + let mut intercepts = Vec::with_capacity(n_cat); + intercepts.push(0.0); + intercepts.extend_from_slice(params); + gpcm_logprobs(base, &scores, &intercepts) + } + }; + let mut mean = 0.0; + for (category, logprob) in logprobs.into_iter().enumerate() { + let probability = logprob.exp(); + probabilities.push(probability); + mean += category as f64 * probability; + } + expected.push(mean); + } + } + Ok(PolytomousPredictions { + probabilities, + expected, + }) +} + /// Result of [`fit_poly_unidim`]. `slope[i]` is item `i`'s discrimination `a_i`; /// `cat_params[i]` holds the `K-1` free category parameters (GPCM additive /// intercepts, or GRM cumulative thresholds). `n_iter` counts completed M-steps; diff --git a/docs/adr/0016-public-polytomous-prediction-boundary.md b/docs/adr/0016-public-polytomous-prediction-boundary.md new file mode 100644 index 000000000..d196ee5ed --- /dev/null +++ b/docs/adr/0016-public-polytomous-prediction-boundary.md @@ -0,0 +1,46 @@ +# ADR-0016: Public polytomous prediction boundary + +Status: **Proposed** +Date: 2026-08-23 + +## Context + +Downstream consumers need GRM/GPCM category probabilities and expected category +scores from a fitted item bank. The package already owns the stable Samejima +(1969) and Muraki (1992) category kernels in Rust, but only private `_core` +cell helpers expose them. Consumers therefore risk reimplementing parameter +conversion and probability arithmetic. + +## Decision + +Expose `polytomous_category_probabilities(fit, theta)` and +`polytomous_expected_response(fit, theta)` from the Python package. Both accept +the existing `PolytomousFit` and a finite one-dimensional trait vector. A +single Rust batch kernel validates item parameters, computes probabilities, +and derives expected integer category scores. PyO3 remains an implementation +boundary; consumers never import `_core`. + +The output shapes are persons x items x categories and persons x items, +respectively. GRM uses the fitted decreasing cumulative-boundary intercepts; +GPCM uses the fitted additive category intercepts. No alternate threshold or +step convention is inferred. + +## Evidence and consequences + +Tests cover GRM/GPCM normalization, expected-score identity, extreme finite +traits, invalid parameters, and delegation through the compiled Rust module. +This adds no dependency and no second numerical implementation. Samejima +(1969) grounds GRM cumulative probabilities; Muraki (1992, 1993) grounds GPCM +category probabilities and expected integer scores. + +## References + +- Samejima, F. (1969). *Estimation of latent ability using a response pattern + of graded scores*. Psychometrika, 34(S1), 1-97. + https://doi.org/10.1007/BF03372160 +- Muraki, E. (1992). A generalized partial credit model: Application of an EM + algorithm. *Applied Psychological Measurement, 16*(2), 159-176. + https://doi.org/10.1177/014662169201600206 +- Muraki, E. (1993). Information functions of the generalized partial credit + model. *Applied Psychological Measurement, 17*(4), 351-363. + https://doi.org/10.1177/014662169301700403 diff --git a/python/fast_mlsirm/_legacy_init.py b/python/fast_mlsirm/_legacy_init.py index b4f72ca79..0031c6d7d 100644 --- a/python/fast_mlsirm/_legacy_init.py +++ b/python/fast_mlsirm/_legacy_init.py @@ -218,7 +218,7 @@ load_serving_bundle as load_serving_bundle, score_respondents as score_respondents) from .preprocessing import irtree_expand as irtree_expand -from .polytomous import fit_polytomous as fit_polytomous, PolytomousFit as PolytomousFit, score_polytomous as score_polytomous, information_polytomous as information_polytomous, fit_lsirm_polytomous as fit_lsirm_polytomous, PolyLsirmFit as PolyLsirmFit, polytomous_information_criteria as polytomous_information_criteria, item_fit_polytomous as item_fit_polytomous, m2_polytomous as m2_polytomous, local_dependence_polytomous as local_dependence_polytomous, fit_nominal_polytomous as fit_nominal_polytomous, NominalFit as NominalFit, person_fit_polytomous as person_fit_polytomous, cat_simulate_polytomous as cat_simulate_polytomous, dif_polytomous as dif_polytomous, u3_person_fit_polytomous as u3_person_fit_polytomous, u3_cutoff_polytomous as u3_cutoff_polytomous +from .polytomous import fit_polytomous as fit_polytomous, PolytomousFit as PolytomousFit, score_polytomous as score_polytomous, information_polytomous as information_polytomous, polytomous_category_probabilities as polytomous_category_probabilities, polytomous_expected_response as polytomous_expected_response, fit_lsirm_polytomous as fit_lsirm_polytomous, PolyLsirmFit as PolyLsirmFit, polytomous_information_criteria as polytomous_information_criteria, item_fit_polytomous as item_fit_polytomous, m2_polytomous as m2_polytomous, local_dependence_polytomous as local_dependence_polytomous, fit_nominal_polytomous as fit_nominal_polytomous, NominalFit as NominalFit, person_fit_polytomous as person_fit_polytomous, cat_simulate_polytomous as cat_simulate_polytomous, dif_polytomous as dif_polytomous, u3_person_fit_polytomous as u3_person_fit_polytomous, u3_cutoff_polytomous as u3_cutoff_polytomous from .dif import (mantel_haenszel_dif as mantel_haenszel_dif, logistic_dif as logistic_dif, mantel_haenszel_dif_purified as mantel_haenszel_dif_purified, logistic_dif_purified as logistic_dif_purified, @@ -515,6 +515,8 @@ "fit_polytomous", "score_polytomous", "information_polytomous", + "polytomous_category_probabilities", + "polytomous_expected_response", "fit_lsirm_polytomous", "PolyLsirmFit", "polytomous_information_criteria", diff --git a/python/fast_mlsirm/polytomous.py b/python/fast_mlsirm/polytomous.py index 0e8fed84a..9db232033 100644 --- a/python/fast_mlsirm/polytomous.py +++ b/python/fast_mlsirm/polytomous.py @@ -30,6 +30,8 @@ "fit_polytomous", "score_polytomous", "information_polytomous", + "polytomous_category_probabilities", + "polytomous_expected_response", "PolyLsirmFit", "fit_lsirm_polytomous", "polytomous_information_criteria", @@ -83,6 +85,64 @@ class PolytomousFit: thresholds: np.ndarray | None = None +def _polytomous_predictions( + fit: PolytomousFit, + theta: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """Validate a fitted bank and delegate its prediction grid to Rust.""" + if not isinstance(fit, PolytomousFit): + raise TypeError("fit must be a PolytomousFit") + th = np.asarray(theta, dtype=np.float64) + if th.ndim != 1 or th.size == 0 or not np.all(np.isfinite(th)): + raise ValueError("theta must be a non-empty finite 1-D array") + slope = np.asarray(fit.slope, dtype=np.float64) + cat_params = np.asarray(fit.cat_params, dtype=np.float64) + if slope.ndim != 1 or slope.size == 0: + raise ValueError("fit.slope must be a non-empty 1-D array") + if ( + cat_params.ndim != 2 + or cat_params.shape[0] != slope.size + or cat_params.shape[1] < 1 + ): + raise ValueError("fit.cat_params must be n_items x (n_cat - 1)") + if not np.all(np.isfinite(slope)) or not np.all(np.isfinite(cat_params)): + raise ValueError("fit item parameters must be finite") + model = fit.model.lower() if type(fit.model) is str else "" + if model not in VALID_POLY_MODELS: + raise ValueError(f"fit.model must be one of {sorted(VALID_POLY_MODELS)}") + core = _core_module() + if core is None or not hasattr(core, "polytomous_predictions"): + raise RuntimeError("polytomous predictions require the compiled Rust core") + result = core.polytomous_predictions( + th, + slope, + cat_params.reshape(-1), + int(slope.size), + int(cat_params.shape[1] + 1), + model, + ) + probabilities = np.asarray(result["probabilities"], dtype=np.float64).reshape( + th.size, slope.size, cat_params.shape[1] + 1 + ) + expected = np.asarray(result["expected"], dtype=np.float64).reshape( + th.size, slope.size + ) + return probabilities, expected + + +def polytomous_category_probabilities( + fit: PolytomousFit, + theta: np.ndarray, +) -> np.ndarray: + """Return ``P(Y=k | theta, item)`` as persons x items x categories.""" + return _polytomous_predictions(fit, theta)[0] + + +def polytomous_expected_response(fit: PolytomousFit, theta: np.ndarray) -> np.ndarray: + """Return ``E[Y | theta, item]`` as a persons x items matrix.""" + return _polytomous_predictions(fit, theta)[1] + + def _core_module(): """Return the compiled Rust core module, or ``None`` if it is unavailable.""" try: diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py new file mode 100644 index 000000000..abdda9831 --- /dev/null +++ b/tests/test_polytomous_predictions.py @@ -0,0 +1,36 @@ +"""Public Rust-owned GRM/GPCM prediction contract.""" + +import numpy as np +import pytest + +from fast_mlsirm import ( + PolytomousFit, + polytomous_category_probabilities, + polytomous_expected_response, +) + + +@pytest.mark.parametrize( + ("model", "cat_params"), + [ + ("grm", np.array([[1.0, -1.0], [0.5, -0.5]])), + ("gpcm", np.array([[0.2, -0.4], [-0.1, 0.3]])), + ], +) +def test_public_polytomous_predictions_are_normalized_and_consistent(model, cat_params): + fit = PolytomousFit(model, np.array([1.0, 0.8]), cat_params, 0.0, 0) + theta = np.array([-1000.0, 0.0, 1000.0]) + + probabilities = polytomous_category_probabilities(fit, theta) + expected = polytomous_expected_response(fit, theta) + + assert probabilities.shape == (3, 2, 3) + assert np.all(np.isfinite(probabilities)) + assert np.allclose(probabilities.sum(axis=2), 1.0) + assert np.allclose(expected, probabilities @ np.arange(3, dtype=np.float64)) + + +def test_public_polytomous_predictions_reject_invalid_grm_thresholds(): + fit = PolytomousFit("grm", np.array([1.0]), np.array([[-1.0, 1.0]]), 0.0, 0) + with pytest.raises(ValueError, match="strictly decreasing"): + polytomous_category_probabilities(fit, np.array([0.0])) From a6cfbe67be331bcb0d1b3c29e51ba65ce203ad0d Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 03:32:32 -0700 Subject: [PATCH 02/59] test(irt): bound polytomous prediction output --- tests/test_polytomous_predictions.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index abdda9831..0565f2379 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -3,6 +3,7 @@ import numpy as np import pytest +import fast_mlsirm.polytomous as polytomous_module from fast_mlsirm import ( PolytomousFit, polytomous_category_probabilities, @@ -34,3 +35,17 @@ def test_public_polytomous_predictions_reject_invalid_grm_thresholds(): fit = PolytomousFit("grm", np.array([1.0]), np.array([[-1.0, 1.0]]), 0.0, 0) with pytest.raises(ValueError, match="strictly decreasing"): polytomous_category_probabilities(fit, np.array([0.0])) + + +def test_public_polytomous_predictions_bound_output_before_native_dispatch(monkeypatch): + """Reject oversized logical prediction grids before Rust/output allocation.""" + + fit = PolytomousFit("grm", np.array([1.0]), np.array([[0.0]]), 0.0, 0) + theta = np.broadcast_to(np.array([0.0]), (10_000_001,)) + + def unexpected_core_discovery(): + raise AssertionError("oversized prediction grid reached compiled-core discovery") + + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match="20,000,000.*prediction"): + polytomous_category_probabilities(fit, theta) From 45c9049afb7dfbd90d7c8fd0a533d30aaff882ec Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 03:36:05 -0700 Subject: [PATCH 03/59] fix(irt): guard polytomous prediction resources --- .../_polytomous_prediction_safety.py | 74 +++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 python/fast_mlsirm/_polytomous_prediction_safety.py diff --git a/python/fast_mlsirm/_polytomous_prediction_safety.py b/python/fast_mlsirm/_polytomous_prediction_safety.py new file mode 100644 index 000000000..418c4e761 --- /dev/null +++ b/python/fast_mlsirm/_polytomous_prediction_safety.py @@ -0,0 +1,74 @@ +"""Resource guard for public polytomous prediction grids. + +The public GRM/GPCM prediction surface returns dense persons x items x categories +probabilities plus persons x items expected scores. This adapter rejects an +oversized logical output grid before compiled-core discovery/allocation while +leaving all probability arithmetic in the Rust implementation. +""" + +from __future__ import annotations + +from functools import wraps +from types import ModuleType +from typing import Any + +import numpy as np + +MAX_POLYTOMOUS_PREDICTION_CELLS = 20_000_000 +_GUARD_MARKER = "_fast_mlsirm_polytomous_prediction_resource_guard" + + +def _trusted_length(value: object) -> int | None: + """Return a callback-free logical vector length when the container is trusted.""" + + value_type = type(value) + if value_type is np.ndarray: + return int(value.size) if value.ndim == 1 else None + if value_type is list or value_type is tuple: + return len(value) + return None + + +def _prediction_cells(module: ModuleType, fit: object, theta: object) -> int | None: + """Return the dense probability-cell count when shape evidence is inert.""" + + if type(fit) is not module.PolytomousFit: + return None + theta_count = _trusted_length(theta) + if theta_count is None or theta_count == 0: + return None + + slope = fit.slope + cat_params = fit.cat_params + if type(slope) is not np.ndarray or type(cat_params) is not np.ndarray: + return None + if slope.ndim != 1 or slope.size == 0: + return None + if ( + cat_params.ndim != 2 + or cat_params.shape[0] != slope.size + or cat_params.shape[1] < 1 + ): + return None + + return theta_count * int(slope.size) * int(cat_params.shape[1] + 1) + + +def install(module: ModuleType) -> None: + """Install an idempotent pre-dispatch resource guard on prediction calls.""" + + original = module._polytomous_predictions + if getattr(original, _GUARD_MARKER, False): + return + + @wraps(original) + def guarded_predictions(fit: Any, theta: Any): + cells = _prediction_cells(module, fit, theta) + if cells is not None and cells > MAX_POLYTOMOUS_PREDICTION_CELLS: + raise ValueError( + f"{MAX_POLYTOMOUS_PREDICTION_CELLS:,}-cell polytomous prediction limit exceeded" + ) + return original(fit, theta) + + setattr(guarded_predictions, _GUARD_MARKER, True) + module._polytomous_predictions = guarded_predictions From b1e906c17cde20a95e1cc86d48601f412b54790b Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 03:36:30 -0700 Subject: [PATCH 04/59] fix(irt): install polytomous prediction resource guard --- python/fast_mlsirm/__init__.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index b918201c2..13e4c4dc4 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -6,11 +6,13 @@ from importlib.metadata import version as _distribution_version from . import _legacy_init as _legacy_init +from . import polytomous as _polytomous from . import reliability as _reliability from . import scaling as _scaling from . import validation as _validation from ._fleiss_control_safety import install as _install_fleiss_control_safety from ._icc_control_safety import install as _install_icc_control_safety +from ._polytomous_prediction_safety import install as _install_polytomous_prediction_safety from ._scaling_control_safety import install as _install_scaling_control_safety # Harden historical public adapters before copying legacy exports. These @@ -19,6 +21,7 @@ _install_icc_control_safety(_reliability) _install_scaling_control_safety(_scaling) _install_fleiss_control_safety(_validation) +_install_polytomous_prediction_safety(_polytomous) _legacy_init.icc = _reliability.icc _legacy_init.bradley_terry_mm = _scaling.bradley_terry_mm _legacy_init.fleiss_kappa = _validation.fleiss_kappa @@ -26,7 +29,9 @@ del ( _install_fleiss_control_safety, _install_icc_control_safety, + _install_polytomous_prediction_safety, _install_scaling_control_safety, + _polytomous, _reliability, _scaling, ) From d841bbb2a1ee1bad8bef5246fea7b634238df863 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 03:37:19 -0700 Subject: [PATCH 05/59] docs: record polytomous prediction resource guard --- docs/changelog.d/1280-polytomous-prediction-resource.md | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 docs/changelog.d/1280-polytomous-prediction-resource.md diff --git a/docs/changelog.d/1280-polytomous-prediction-resource.md b/docs/changelog.d/1280-polytomous-prediction-resource.md new file mode 100644 index 000000000..50c8fc9a5 --- /dev/null +++ b/docs/changelog.d/1280-polytomous-prediction-resource.md @@ -0,0 +1,7 @@ +# Bound public polytomous prediction resources + +## Fixed + +- Reject public GRM/GPCM prediction grids above 20,000,000 dense probability cells before compiled-core discovery or output allocation. +- Keep GRM/GPCM category-probability and expected-score arithmetic in the existing Rust implementation; the Python guard is resource admission only. +- Track the remaining direct private-PyO3/core allocation bound in issue #1280 rather than treating the public guard as full native completion. From d025b7d237d8db7ca97a5611606c6285d5870895 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 19:46:01 +0900 Subject: [PATCH 06/59] fix(irt): enforce prediction limit at owner --- python/fast_mlsirm/__init__.py | 5 -- .../_polytomous_prediction_safety.py | 74 ------------------- python/fast_mlsirm/polytomous.py | 6 ++ 3 files changed, 6 insertions(+), 79 deletions(-) delete mode 100644 python/fast_mlsirm/_polytomous_prediction_safety.py diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index 13e4c4dc4..b918201c2 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -6,13 +6,11 @@ from importlib.metadata import version as _distribution_version from . import _legacy_init as _legacy_init -from . import polytomous as _polytomous from . import reliability as _reliability from . import scaling as _scaling from . import validation as _validation from ._fleiss_control_safety import install as _install_fleiss_control_safety from ._icc_control_safety import install as _install_icc_control_safety -from ._polytomous_prediction_safety import install as _install_polytomous_prediction_safety from ._scaling_control_safety import install as _install_scaling_control_safety # Harden historical public adapters before copying legacy exports. These @@ -21,7 +19,6 @@ _install_icc_control_safety(_reliability) _install_scaling_control_safety(_scaling) _install_fleiss_control_safety(_validation) -_install_polytomous_prediction_safety(_polytomous) _legacy_init.icc = _reliability.icc _legacy_init.bradley_terry_mm = _scaling.bradley_terry_mm _legacy_init.fleiss_kappa = _validation.fleiss_kappa @@ -29,9 +26,7 @@ del ( _install_fleiss_control_safety, _install_icc_control_safety, - _install_polytomous_prediction_safety, _install_scaling_control_safety, - _polytomous, _reliability, _scaling, ) diff --git a/python/fast_mlsirm/_polytomous_prediction_safety.py b/python/fast_mlsirm/_polytomous_prediction_safety.py deleted file mode 100644 index 418c4e761..000000000 --- a/python/fast_mlsirm/_polytomous_prediction_safety.py +++ /dev/null @@ -1,74 +0,0 @@ -"""Resource guard for public polytomous prediction grids. - -The public GRM/GPCM prediction surface returns dense persons x items x categories -probabilities plus persons x items expected scores. This adapter rejects an -oversized logical output grid before compiled-core discovery/allocation while -leaving all probability arithmetic in the Rust implementation. -""" - -from __future__ import annotations - -from functools import wraps -from types import ModuleType -from typing import Any - -import numpy as np - -MAX_POLYTOMOUS_PREDICTION_CELLS = 20_000_000 -_GUARD_MARKER = "_fast_mlsirm_polytomous_prediction_resource_guard" - - -def _trusted_length(value: object) -> int | None: - """Return a callback-free logical vector length when the container is trusted.""" - - value_type = type(value) - if value_type is np.ndarray: - return int(value.size) if value.ndim == 1 else None - if value_type is list or value_type is tuple: - return len(value) - return None - - -def _prediction_cells(module: ModuleType, fit: object, theta: object) -> int | None: - """Return the dense probability-cell count when shape evidence is inert.""" - - if type(fit) is not module.PolytomousFit: - return None - theta_count = _trusted_length(theta) - if theta_count is None or theta_count == 0: - return None - - slope = fit.slope - cat_params = fit.cat_params - if type(slope) is not np.ndarray or type(cat_params) is not np.ndarray: - return None - if slope.ndim != 1 or slope.size == 0: - return None - if ( - cat_params.ndim != 2 - or cat_params.shape[0] != slope.size - or cat_params.shape[1] < 1 - ): - return None - - return theta_count * int(slope.size) * int(cat_params.shape[1] + 1) - - -def install(module: ModuleType) -> None: - """Install an idempotent pre-dispatch resource guard on prediction calls.""" - - original = module._polytomous_predictions - if getattr(original, _GUARD_MARKER, False): - return - - @wraps(original) - def guarded_predictions(fit: Any, theta: Any): - cells = _prediction_cells(module, fit, theta) - if cells is not None and cells > MAX_POLYTOMOUS_PREDICTION_CELLS: - raise ValueError( - f"{MAX_POLYTOMOUS_PREDICTION_CELLS:,}-cell polytomous prediction limit exceeded" - ) - return original(fit, theta) - - setattr(guarded_predictions, _GUARD_MARKER, True) - module._polytomous_predictions = guarded_predictions diff --git a/python/fast_mlsirm/polytomous.py b/python/fast_mlsirm/polytomous.py index 9db232033..05e9e413d 100644 --- a/python/fast_mlsirm/polytomous.py +++ b/python/fast_mlsirm/polytomous.py @@ -110,6 +110,12 @@ def _polytomous_predictions( model = fit.model.lower() if type(fit.model) is str else "" if model not in VALID_POLY_MODELS: raise ValueError(f"fit.model must be one of {sorted(VALID_POLY_MODELS)}") + prediction_cells = int(th.size) * int(slope.size) * int(cat_params.shape[1] + 1) + if prediction_cells > 20_000_000: + raise ValueError( + f"prediction grid of {prediction_cells:,} cells exceeds the " + "20,000,000 prediction-cell limit" + ) core = _core_module() if core is None or not hasattr(core, "polytomous_predictions"): raise RuntimeError("polytomous predictions require the compiled Rust core") From 8d6b62abb10ea5981561bd4be84caa3077ddb48e Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:11:43 -0700 Subject: [PATCH 07/59] test(irt): reject callback-bearing prediction evidence --- tests/test_polytomous_predictions.py | 43 ++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index 0565f2379..d3871fe2b 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -11,6 +11,11 @@ ) +class _ArrayProvider: + def __array__(self, *args, **kwargs): + raise AssertionError("caller __array__ callback executed") + + @pytest.mark.parametrize( ("model", "cat_params"), [ @@ -49,3 +54,41 @@ def unexpected_core_discovery(): monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) with pytest.raises(ValueError, match="20,000,000.*prediction"): polytomous_category_probabilities(fit, theta) + + +@pytest.mark.parametrize("field", ["theta", "slope", "cat_params"]) +def test_public_polytomous_predictions_reject_array_providers_before_callbacks( + monkeypatch, field +): + """Caller array protocols must not run while prediction evidence is admitted.""" + + theta = np.array([0.0]) + slope = np.array([1.0]) + cat_params = np.array([[0.0]]) + if field == "theta": + theta = _ArrayProvider() + elif field == "slope": + slope = _ArrayProvider() + else: + cat_params = _ArrayProvider() + fit = PolytomousFit("grm", slope, cat_params, 0.0, 0) + + def unexpected_core_discovery(): + raise AssertionError("untrusted evidence reached compiled-core discovery") + + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match="trusted NumPy array or built-in sequence"): + polytomous_category_probabilities(fit, theta) + + +def test_public_polytomous_predictions_reject_complex_before_native_dispatch(monkeypatch): + """Do not silently project complex theta evidence onto the real line.""" + + fit = PolytomousFit("gpcm", np.array([1.0]), np.array([[0.0]]), 0.0, 0) + + def unexpected_core_discovery(): + raise AssertionError("complex evidence reached compiled-core discovery") + + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match="real-valued"): + polytomous_expected_response(fit, np.array([0.0 + 1.0j])) From 2a7cf5565a8ac3b468a3a8f575599760454b8260 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:13:03 -0700 Subject: [PATCH 08/59] fix(irt): admit prediction evidence before NumPy coercion --- .../_polytomous_prediction_admission.py | 166 ++++++++++++++++++ 1 file changed, 166 insertions(+) create mode 100644 python/fast_mlsirm/_polytomous_prediction_admission.py diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py new file mode 100644 index 000000000..0c4882df7 --- /dev/null +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -0,0 +1,166 @@ +"""Callback-free public evidence admission for GRM/GPCM predictions.""" + +from __future__ import annotations + +from types import ModuleType + +import numpy as np + + +_NUMPY_INTEGER_SCALAR_TYPES = ( + np.int8, + np.int16, + np.int32, + np.int64, + np.intp, + np.longlong, + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.uintp, + np.ulonglong, +) +_NUMPY_FLOAT_SCALAR_TYPES = (np.float16, np.float32, np.float64, np.longdouble) +_NUMPY_COMPLEX_SCALAR_TYPES = (np.complex64, np.complex128, np.clongdouble) +_TRUSTED_SCALAR_TYPES = ( + bool, + int, + float, + complex, + np.bool_, + *_NUMPY_INTEGER_SCALAR_TYPES, + *_NUMPY_FLOAT_SCALAR_TYPES, + *_NUMPY_COMPLEX_SCALAR_TYPES, +) +_MAX_PREDICTION_CELLS = 20_000_000 +_MARKER = "__fast_mlsirm_polytomous_prediction_admission__" + + +def _preflight(value: object, name: str) -> None: + """Reject callback-bearing and over-budget evidence before NumPy conversion.""" + + error = f"{name} must be a trusted NumPy array or built-in sequence" + resource_error = ( + f"{name} exceeds the {_MAX_PREDICTION_CELLS:,}-cell prediction evidence budget" + ) + frames: list[list[object]] = [[value, 0, 0, False]] + active_container_ids: set[int] = set() + subtree_cells: dict[int, int] = {} + + def add_to_parent(count: int) -> None: + if count > _MAX_PREDICTION_CELLS: + raise ValueError(resource_error) + if not frames: + return + parent = frames[-1] + subtotal = int(parent[2]) + if count > _MAX_PREDICTION_CELLS - subtotal: + raise ValueError(resource_error) + parent[2] = subtotal + count + + while frames: + frame = frames[-1] + item = frame[0] + item_type = type(item) + + if item_type is np.ndarray: + frames.pop() + add_to_parent(int(item.size)) + continue + + if any(item_type is scalar_type for scalar_type in _TRUSTED_SCALAR_TYPES): + frames.pop() + add_to_parent(1) + continue + + if item_type is not list and item_type is not tuple: + raise ValueError(error) + + item_id = id(item) + if not bool(frame[3]): + if item_id in subtree_cells: + frames.pop() + add_to_parent(subtree_cells[item_id]) + continue + if item_id in active_container_ids: + raise ValueError(error) + active_container_ids.add(item_id) + frame[3] = True + + child_index = int(frame[1]) + if child_index < len(item): + frame[1] = child_index + 1 + frames.append([item[child_index], 0, 0, False]) + continue + + count = int(frame[2]) + active_container_ids.remove(item_id) + subtree_cells[item_id] = count + frames.pop() + add_to_parent(count) + + +def _real_numeric_array(value: object, name: str) -> np.ndarray: + """Materialize trusted evidence only when float64 conversion is lossless.""" + + _preflight(value, name) + raw = np.asarray(value) + if np.iscomplexobj(raw): + raise ValueError(f"{name} must be real-valued") + if raw.dtype.kind not in ("b", "i", "u", "f"): + raise ValueError(f"{name} must be a numeric array") + with np.errstate(over="ignore", invalid="ignore"): + converted = raw.astype(np.float64, copy=False) + round_tripped = converted.astype(raw.dtype, copy=False) + if raw.dtype.kind == "f": + lossless = np.array_equal(raw, round_tripped, equal_nan=True) + else: + lossless = np.array_equal(raw, round_tripped) + if not lossless: + raise ValueError(f"{name} must be exactly representable as float64") + return np.ascontiguousarray(converted) + + +def install(module: ModuleType) -> None: + """Install trusted prediction evidence admission on the polytomous module.""" + + current = module._polytomous_predictions + if getattr(current, _MARKER, False): + return + + def guarded_predictions(fit, theta): + if type(fit) is not module.PolytomousFit: + raise TypeError("fit must be a PolytomousFit") + + trusted_theta = _real_numeric_array(theta, "theta") + trusted_slope = _real_numeric_array(fit.slope, "fit.slope") + trusted_cat_params = _real_numeric_array(fit.cat_params, "fit.cat_params") + + if ( + trusted_theta.ndim == 1 + and trusted_slope.ndim == 1 + and trusted_cat_params.ndim == 2 + and trusted_slope.size > 0 + and trusted_cat_params.shape[0] == trusted_slope.size + and trusted_cat_params.shape[1] >= 1 + ): + n_cat = int(trusted_cat_params.shape[1]) + 1 + cells = int(trusted_theta.size) * int(trusted_slope.size) * n_cat + if cells > _MAX_PREDICTION_CELLS: + raise ValueError( + f"prediction grid of {cells:,} cells exceeds the " + f"{_MAX_PREDICTION_CELLS:,} prediction-cell limit" + ) + + trusted_fit = module.PolytomousFit( + fit.model, + trusted_slope, + trusted_cat_params, + fit.loglik, + fit.n_iter, + ) + return current(trusted_fit, trusted_theta) + + setattr(guarded_predictions, _MARKER, True) + module._polytomous_predictions = guarded_predictions From 941843a9db775f5e18030c76d73669ef86686afc Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:13:31 -0700 Subject: [PATCH 09/59] fix(irt): install trusted prediction evidence admission --- python/fast_mlsirm/__init__.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index b918201c2..aa4650bd0 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -6,19 +6,24 @@ from importlib.metadata import version as _distribution_version from . import _legacy_init as _legacy_init +from . import polytomous as _polytomous from . import reliability as _reliability from . import scaling as _scaling from . import validation as _validation from ._fleiss_control_safety import install as _install_fleiss_control_safety from ._icc_control_safety import install as _install_icc_control_safety +from ._polytomous_prediction_admission import ( + install as _install_polytomous_prediction_admission, +) from ._scaling_control_safety import install as _install_scaling_control_safety # Harden historical public adapters before copying legacy exports. These -# wrappers validate and normalize semantic controls only; result arithmetic -# remains in the existing Rust-backed implementations. +# wrappers validate and normalize semantic controls/evidence only; result +# arithmetic remains in the existing Rust-backed implementations. _install_icc_control_safety(_reliability) _install_scaling_control_safety(_scaling) _install_fleiss_control_safety(_validation) +_install_polytomous_prediction_admission(_polytomous) _legacy_init.icc = _reliability.icc _legacy_init.bradley_terry_mm = _scaling.bradley_terry_mm _legacy_init.fleiss_kappa = _validation.fleiss_kappa @@ -26,7 +31,9 @@ del ( _install_fleiss_control_safety, _install_icc_control_safety, + _install_polytomous_prediction_admission, _install_scaling_control_safety, + _polytomous, _reliability, _scaling, ) From 10dc8015a20239ea21db4647be05078ea8a00dab Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:14:09 -0700 Subject: [PATCH 10/59] test(irt): bound prediction grid before dense copy --- tests/test_polytomous_predictions.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index d3871fe2b..89116610b 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -3,6 +3,7 @@ import numpy as np import pytest +import fast_mlsirm._polytomous_prediction_admission as prediction_admission import fast_mlsirm.polytomous as polytomous_module from fast_mlsirm import ( PolytomousFit, @@ -56,6 +57,24 @@ def unexpected_core_discovery(): polytomous_category_probabilities(fit, theta) +def test_public_polytomous_predictions_bound_output_before_dense_copy(monkeypatch): + """Output-cell admission must precede contiguous materialization of broadcast theta.""" + + fit = PolytomousFit("grm", np.array([1.0]), np.array([[0.0]]), 0.0, 0) + theta = np.broadcast_to(np.array([0.0]), (10_000_001,)) + + def unexpected_copy(*args, **kwargs): + raise AssertionError("oversized prediction grid reached contiguous materialization") + + def unexpected_core_discovery(): + raise AssertionError("oversized prediction grid reached compiled-core discovery") + + monkeypatch.setattr(prediction_admission.np, "ascontiguousarray", unexpected_copy) + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match="20,000,000.*prediction"): + polytomous_category_probabilities(fit, theta) + + @pytest.mark.parametrize("field", ["theta", "slope", "cat_params"]) def test_public_polytomous_predictions_reject_array_providers_before_callbacks( monkeypatch, field From 192eb7d0e2b388d8c4e78c8095b1befb5bb195a8 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:14:47 -0700 Subject: [PATCH 11/59] fix(irt): bound prediction grid before dense materialization --- python/fast_mlsirm/_polytomous_prediction_admission.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py index 0c4882df7..0decfff02 100644 --- a/python/fast_mlsirm/_polytomous_prediction_admission.py +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -102,7 +102,7 @@ def add_to_parent(count: int) -> None: def _real_numeric_array(value: object, name: str) -> np.ndarray: - """Materialize trusted evidence only when float64 conversion is lossless.""" + """Admit trusted evidence and verify lossless float64 representation.""" _preflight(value, name) raw = np.asarray(value) @@ -119,7 +119,7 @@ def _real_numeric_array(value: object, name: str) -> np.ndarray: lossless = np.array_equal(raw, round_tripped) if not lossless: raise ValueError(f"{name} must be exactly representable as float64") - return np.ascontiguousarray(converted) + return converted def install(module: ModuleType) -> None: @@ -153,6 +153,9 @@ def guarded_predictions(fit, theta): f"{_MAX_PREDICTION_CELLS:,} prediction-cell limit" ) + trusted_theta = np.ascontiguousarray(trusted_theta) + trusted_slope = np.ascontiguousarray(trusted_slope) + trusted_cat_params = np.ascontiguousarray(trusted_cat_params) trusted_fit = module.PolytomousFit( fit.model, trusted_slope, From 49e201a43a4d04eeeaae8e2100b08416c1cbc861 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:14:59 -0700 Subject: [PATCH 12/59] docs(changelog): record prediction evidence admission --- docs/changelog.d/1281-polytomous-prediction-admission.md | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 docs/changelog.d/1281-polytomous-prediction-admission.md diff --git a/docs/changelog.d/1281-polytomous-prediction-admission.md b/docs/changelog.d/1281-polytomous-prediction-admission.md new file mode 100644 index 000000000..fc395aca0 --- /dev/null +++ b/docs/changelog.d/1281-polytomous-prediction-admission.md @@ -0,0 +1,7 @@ +# Polytomous prediction evidence admission + +## Fixed + +- Reject callback-bearing, complex, non-numeric, and lossy GRM/GPCM prediction evidence before NumPy narrowing or compiled-core discovery. +- Preserve trusted exact NumPy and built-in sequence inputs while enforcing the 20,000,000-cell output ceiling before contiguous dense materialization. +- Keep category-probability and expected-score arithmetic in the Rust prediction kernel; Python performs validation, bounded materialization, and marshalling only. From 425e310b5e185fabdffd9a743d4016ca579c9876 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:33:50 -0700 Subject: [PATCH 13/59] test(poly): require native prediction allocation ceiling --- .../tests/polytomous_predictions_resource.rs | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 crates/mlsirm-core/tests/polytomous_predictions_resource.rs diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs new file mode 100644 index 000000000..065f68bbd --- /dev/null +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -0,0 +1,26 @@ +use mlsirm_core::poly::{polytomous_predictions, PolyModel}; + +#[test] +fn native_predictions_reject_oversized_probability_grid_before_parameter_validation() { + // 312_501 * 64 = 20_000_064 probability cells. Deliberately keep the + // parameter slices empty: the resource envelope must be decided from the + // request shape before any parameter-dependent allocation or validation. + let error = polytomous_predictions(&[0.0], &[], &[], 312_501, 64, PolyModel::Gpcm) + .expect_err("oversized native prediction grid must fail closed"); + + assert!( + error.contains("20,000,000") && error.contains("prediction"), + "unexpected native prediction resource error: {error}" + ); +} + +#[test] +fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { + let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) + .expect("small valid prediction should remain accepted"); + + assert_eq!(prediction.probabilities.len(), 2); + assert_eq!(prediction.expected.len(), 1); + assert!((prediction.probabilities.iter().sum::() - 1.0).abs() < 1e-12); + assert!((prediction.expected[0] - 0.5).abs() < 1e-12); +} From ae48db32c6b6c2ed00d0f14afd2b1fc4fa805fce Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:38:17 -0700 Subject: [PATCH 14/59] fix(poly): bound native prediction output allocation --- crates/mlsirm-core/src/poly.rs | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/crates/mlsirm-core/src/poly.rs b/crates/mlsirm-core/src/poly.rs index b85520dc5..2ab9e5f4b 100644 --- a/crates/mlsirm-core/src/poly.rs +++ b/crates/mlsirm-core/src/poly.rs @@ -17,6 +17,7 @@ pub(crate) const POLY_MAX_CAT: usize = 64; pub(crate) const POLY_MAX_ITER: usize = 100_000; +pub(crate) const POLY_MAX_PREDICTION_CELLS: usize = 20_000_000; fn validate_observed_categories( y: &[usize], @@ -218,7 +219,6 @@ pub fn polytomous_predictions( if n_items == 0 || n_cat < 2 { return Err("n_items must be positive and n_cat must be at least 2".into()); } - validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?; let cells = crate::checked_mul_usize( theta.len(), n_items, @@ -229,6 +229,12 @@ pub fn polytomous_predictions( n_cat, "prediction cells * n_cat exceeds the probability buffer size", )?; + if probability_cells > POLY_MAX_PREDICTION_CELLS { + return Err(format!( + "polytomous prediction grid of {probability_cells} cells exceeds the 20,000,000 prediction-cell limit" + )); + } + validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?; let mut probabilities = Vec::with_capacity(probability_cells); let mut expected = Vec::with_capacity(cells); let scores: Vec = (0..n_cat).map(|category| category as f64).collect(); From 4eda7b25de51ff662ae1f752a3027a4a0e727ca8 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 04:38:55 -0700 Subject: [PATCH 15/59] docs(poly): record native prediction allocation ceiling --- docs/changelog.d/1280-polytomous-prediction-resource.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/changelog.d/1280-polytomous-prediction-resource.md b/docs/changelog.d/1280-polytomous-prediction-resource.md index 50c8fc9a5..af807f8a2 100644 --- a/docs/changelog.d/1280-polytomous-prediction-resource.md +++ b/docs/changelog.d/1280-polytomous-prediction-resource.md @@ -1,7 +1,7 @@ -# Bound public polytomous prediction resources +# Bound polytomous prediction resources ## Fixed - Reject public GRM/GPCM prediction grids above 20,000,000 dense probability cells before compiled-core discovery or output allocation. -- Keep GRM/GPCM category-probability and expected-score arithmetic in the existing Rust implementation; the Python guard is resource admission only. -- Track the remaining direct private-PyO3/core allocation bound in issue #1280 rather than treating the public guard as full native completion. +- Apply the same 20,000,000-cell ceiling inside the Rust `polytomous_predictions` owner before item-parameter validation or `Vec::with_capacity`, so direct core/PyO3 callers cannot bypass the public resource envelope. +- Keep GRM/GPCM category-probability and expected-score arithmetic in the existing Rust implementation; the added checks govern request size and allocation only. From 4cb3117e4dda30d5eec7de326a5aac2f3c486368 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 07:10:34 -0700 Subject: [PATCH 16/59] test(irt): bound prediction evidence rank before NumPy --- tests/test_polytomous_predictions.py | 45 ++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index 89116610b..d7e86b5dd 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -100,6 +100,51 @@ def unexpected_core_discovery(): polytomous_category_probabilities(fit, theta) +@pytest.mark.parametrize( + ("field", "over_rank", "expected_rank"), + [ + ("theta", [[0.0]], 1), + ("slope", [[1.0]], 1), + ("cat_params", [[[0.0]]], 2), + ], +) +def test_public_polytomous_predictions_reject_over_rank_before_numpy_materialization( + monkeypatch, field, over_rank, expected_rank +): + """Known over-rank evidence must fail before NumPy can materialize it.""" + + theta = np.array([0.0]) + slope = np.array([1.0]) + cat_params = np.array([[0.0]]) + if field == "theta": + theta = over_rank + elif field == "slope": + slope = over_rank + else: + cat_params = over_rank + fit = PolytomousFit("grm", slope, cat_params, 0.0, 0) + + def unexpected_numpy_materialization(*args, **kwargs): + raise AssertionError("over-rank evidence reached NumPy materialization") + + def unexpected_core_discovery(): + raise AssertionError("over-rank evidence reached compiled-core discovery") + + monkeypatch.setattr(prediction_admission.np, "asarray", unexpected_numpy_materialization) + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match=rf"{field} must be at most {expected_rank}-D"): + polytomous_category_probabilities(fit, theta) + + +def test_public_polytomous_predictions_preserve_exact_numpy_row_sequences(): + """A trusted NumPy row nested in the 2-D category matrix remains compatible.""" + + fit = PolytomousFit("gpcm", [1.0], [np.array([0.0])], 0.0, 0) + probabilities = polytomous_category_probabilities(fit, [0.0]) + assert probabilities.shape == (1, 1, 2) + assert np.allclose(probabilities.sum(axis=2), 1.0) + + def test_public_polytomous_predictions_reject_complex_before_native_dispatch(monkeypatch): """Do not silently project complex theta evidence onto the real line.""" From 854c3dfd89ac49b714096e7cf393a108b611450f Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 07:11:49 -0700 Subject: [PATCH 17/59] fix(irt): reject over-rank prediction evidence early --- .../_polytomous_prediction_admission.py | 53 ++++++++++++------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py index 0decfff02..1a98b13f7 100644 --- a/python/fast_mlsirm/_polytomous_prediction_admission.py +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -37,18 +37,20 @@ _MARKER = "__fast_mlsirm_polytomous_prediction_admission__" -def _preflight(value: object, name: str) -> None: - """Reject callback-bearing and over-budget evidence before NumPy conversion.""" +def _preflight(value: object, name: str, max_ndim: int) -> None: + """Reject callback-bearing, over-rank, and over-budget evidence before NumPy.""" error = f"{name} must be a trusted NumPy array or built-in sequence" + rank_error = f"{name} must be at most {max_ndim}-D" resource_error = ( f"{name} exceeds the {_MAX_PREDICTION_CELLS:,}-cell prediction evidence budget" ) - frames: list[list[object]] = [[value, 0, 0, False]] + # frame = [item, next_child_index, logical_cells, subtree_rank, entered] + frames: list[list[object]] = [[value, 0, 0, 1, False]] active_container_ids: set[int] = set() - subtree_cells: dict[int, int] = {} + subtree_metrics: dict[int, tuple[int, int]] = {} - def add_to_parent(count: int) -> None: + def add_to_parent(count: int, subtree_rank: int) -> None: if count > _MAX_PREDICTION_CELLS: raise ValueError(resource_error) if not frames: @@ -58,53 +60,66 @@ def add_to_parent(count: int) -> None: if count > _MAX_PREDICTION_CELLS - subtotal: raise ValueError(resource_error) parent[2] = subtotal + count + parent[3] = max(int(parent[3]), 1 + subtree_rank) while frames: frame = frames[-1] item = frame[0] item_type = type(item) + depth = len(frames) - 1 if item_type is np.ndarray: + if depth + int(item.ndim) > max_ndim: + raise ValueError(rank_error) frames.pop() - add_to_parent(int(item.size)) + add_to_parent(int(item.size), int(item.ndim)) continue if any(item_type is scalar_type for scalar_type in _TRUSTED_SCALAR_TYPES): + if depth > max_ndim: + raise ValueError(rank_error) frames.pop() - add_to_parent(1) + add_to_parent(1, 0) continue if item_type is not list and item_type is not tuple: raise ValueError(error) + if depth + 1 > max_ndim: + raise ValueError(rank_error) item_id = id(item) - if not bool(frame[3]): - if item_id in subtree_cells: + if not bool(frame[4]): + cached = subtree_metrics.get(item_id) + if cached is not None: + count, subtree_rank = cached + if depth + subtree_rank > max_ndim: + raise ValueError(rank_error) frames.pop() - add_to_parent(subtree_cells[item_id]) + add_to_parent(count, subtree_rank) continue if item_id in active_container_ids: raise ValueError(error) active_container_ids.add(item_id) - frame[3] = True + frame[4] = True child_index = int(frame[1]) if child_index < len(item): frame[1] = child_index + 1 - frames.append([item[child_index], 0, 0, False]) + frames.append([item[child_index], 0, 0, 1, False]) continue count = int(frame[2]) + subtree_rank = int(frame[3]) active_container_ids.remove(item_id) - subtree_cells[item_id] = count + subtree_metrics[item_id] = (count, subtree_rank) frames.pop() - add_to_parent(count) + add_to_parent(count, subtree_rank) -def _real_numeric_array(value: object, name: str) -> np.ndarray: +def _real_numeric_array(value: object, name: str, max_ndim: int) -> np.ndarray: """Admit trusted evidence and verify lossless float64 representation.""" - _preflight(value, name) + _preflight(value, name, max_ndim) raw = np.asarray(value) if np.iscomplexobj(raw): raise ValueError(f"{name} must be real-valued") @@ -133,9 +148,9 @@ def guarded_predictions(fit, theta): if type(fit) is not module.PolytomousFit: raise TypeError("fit must be a PolytomousFit") - trusted_theta = _real_numeric_array(theta, "theta") - trusted_slope = _real_numeric_array(fit.slope, "fit.slope") - trusted_cat_params = _real_numeric_array(fit.cat_params, "fit.cat_params") + trusted_theta = _real_numeric_array(theta, "theta", 1) + trusted_slope = _real_numeric_array(fit.slope, "fit.slope", 1) + trusted_cat_params = _real_numeric_array(fit.cat_params, "fit.cat_params", 2) if ( trusted_theta.ndim == 1 From 3df104a1a3357b469a87282c85410b57756d310b Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 07:12:27 -0700 Subject: [PATCH 18/59] docs(changelog): record prediction rank admission --- docs/changelog.d/1291-polytomous-prediction-rank.md | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 docs/changelog.d/1291-polytomous-prediction-rank.md diff --git a/docs/changelog.d/1291-polytomous-prediction-rank.md b/docs/changelog.d/1291-polytomous-prediction-rank.md new file mode 100644 index 000000000..2e44773c6 --- /dev/null +++ b/docs/changelog.d/1291-polytomous-prediction-rank.md @@ -0,0 +1,7 @@ +# Polytomous prediction evidence rank admission + +## Fixed + +- Reject trusted-but-over-rank `theta`, slope, and category-parameter evidence before NumPy materialization or compiled-core discovery. +- Preserve exact 1-D theta/slope and 2-D category-parameter inputs, including exact NumPy row arrays nested in trusted built-in category matrices. +- Keep the existing prediction/evidence cell ceilings and Rust-owned GRM/GPCM probability and expected-score arithmetic unchanged. From df2922df968200e6a78eedda211cfb6f91ed444f Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:01:55 -0700 Subject: [PATCH 19/59] fix(predictions): preflight all evidence before NumPy --- .../_polytomous_prediction_admission.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py index 1a98b13f7..7fe5a324a 100644 --- a/python/fast_mlsirm/_polytomous_prediction_admission.py +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -116,10 +116,9 @@ def add_to_parent(count: int, subtree_rank: int) -> None: add_to_parent(count, subtree_rank) -def _real_numeric_array(value: object, name: str, max_ndim: int) -> np.ndarray: - """Admit trusted evidence and verify lossless float64 representation.""" +def _real_numeric_array(value: object, name: str) -> np.ndarray: + """Materialize preflighted evidence and verify lossless float64 representation.""" - _preflight(value, name, max_ndim) raw = np.asarray(value) if np.iscomplexobj(raw): raise ValueError(f"{name} must be real-valued") @@ -148,9 +147,16 @@ def guarded_predictions(fit, theta): if type(fit) is not module.PolytomousFit: raise TypeError("fit must be a PolytomousFit") - trusted_theta = _real_numeric_array(theta, "theta", 1) - trusted_slope = _real_numeric_array(fit.slope, "fit.slope", 1) - trusted_cat_params = _real_numeric_array(fit.cat_params, "fit.cat_params", 2) + # Establish every field's callback/rank/resource contract before any + # one field reaches NumPy materialization. This keeps a later invalid + # item parameter from being preceded by unrelated caller data work. + _preflight(theta, "theta", 1) + _preflight(fit.slope, "fit.slope", 1) + _preflight(fit.cat_params, "fit.cat_params", 2) + + trusted_theta = _real_numeric_array(theta, "theta") + trusted_slope = _real_numeric_array(fit.slope, "fit.slope") + trusted_cat_params = _real_numeric_array(fit.cat_params, "fit.cat_params") if ( trusted_theta.ndim == 1 From c41c9a5de5194ae032b8b4eaf1bdc0cd25a0f173 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:10:06 -0700 Subject: [PATCH 20/59] test(predictions): expose pre-conversion evidence gaps --- tests/test_polytomous_predictions.py | 44 ++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index d7e86b5dd..80f967c9f 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -156,3 +156,47 @@ def unexpected_core_discovery(): monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) with pytest.raises(ValueError, match="real-valued"): polytomous_expected_response(fit, np.array([0.0 + 1.0j])) + + +def test_public_polytomous_predictions_reject_mixed_integer_precision_loss_before_numpy( + monkeypatch, +): + """Mixed sequence promotion must not erase integer identity before validation.""" + + fit = PolytomousFit("gpcm", [1.0], [[0.0]], 0.0, 0) + theta = [9_007_199_254_740_993, 0.0] + + def unexpected_numpy_materialization(*args, **kwargs): + raise AssertionError("lossy mixed evidence reached NumPy materialization") + + def unexpected_core_discovery(): + raise AssertionError("lossy mixed evidence reached compiled-core discovery") + + monkeypatch.setattr(prediction_admission.np, "asarray", unexpected_numpy_materialization) + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match="exactly representable as float64"): + polytomous_expected_response(fit, theta) + + +def test_public_polytomous_predictions_bound_joint_grid_before_float64_conversion( + monkeypatch, +): + """Joint output budget must be decided from shapes before float64 conversion.""" + + theta = np.broadcast_to(np.array([0.0], dtype=np.float32), (2,)) + slope = np.broadcast_to(np.array([1.0], dtype=np.float32), (10_000_001,)) + cat_params = np.broadcast_to( + np.array([[0.0]], dtype=np.float32), (10_000_001, 1) + ) + fit = PolytomousFit("gpcm", slope, cat_params, 0.0, 0) + + def unexpected_numpy_materialization(*args, **kwargs): + raise AssertionError("oversized joint grid reached NumPy materialization") + + def unexpected_core_discovery(): + raise AssertionError("oversized joint grid reached compiled-core discovery") + + monkeypatch.setattr(prediction_admission.np, "asarray", unexpected_numpy_materialization) + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match="20,000,000.*prediction"): + polytomous_category_probabilities(fit, theta) From e8d8335536dd66e1f1d4ed1a3bf2f9a94a76efe5 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:11:08 -0700 Subject: [PATCH 21/59] fix(predictions): seal pre-conversion precision and grid bounds --- .../_polytomous_prediction_admission.py | 179 +++++++++++------- 1 file changed, 114 insertions(+), 65 deletions(-) diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py index 7fe5a324a..1e6b27902 100644 --- a/python/fast_mlsirm/_polytomous_prediction_admission.py +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -37,50 +37,49 @@ _MARKER = "__fast_mlsirm_polytomous_prediction_admission__" -def _preflight(value: object, name: str, max_ndim: int) -> None: - """Reject callback-bearing, over-rank, and over-budget evidence before NumPy.""" +def _integer_is_exact_float64(value: object) -> bool: + """Return whether a trusted integer scalar survives float64 exactly.""" + + integer = int(value) + try: + converted = float(integer) + except OverflowError: + return False + return np.isfinite(converted) and int(converted) == integer + + +def _preflight(value: object, name: str, max_ndim: int) -> tuple[int, ...] | None: + """Reject unsafe evidence before NumPy and return a proven rectangular shape.""" error = f"{name} must be a trusted NumPy array or built-in sequence" rank_error = f"{name} must be at most {max_ndim}-D" resource_error = ( f"{name} exceeds the {_MAX_PREDICTION_CELLS:,}-cell prediction evidence budget" ) - # frame = [item, next_child_index, logical_cells, subtree_rank, entered] - frames: list[list[object]] = [[value, 0, 0, 1, False]] + exact_error = f"{name} must be exactly representable as float64" active_container_ids: set[int] = set() - subtree_metrics: dict[int, tuple[int, int]] = {} - - def add_to_parent(count: int, subtree_rank: int) -> None: - if count > _MAX_PREDICTION_CELLS: - raise ValueError(resource_error) - if not frames: - return - parent = frames[-1] - subtotal = int(parent[2]) - if count > _MAX_PREDICTION_CELLS - subtotal: - raise ValueError(resource_error) - parent[2] = subtotal + count - parent[3] = max(int(parent[3]), 1 + subtree_rank) - - while frames: - frame = frames[-1] - item = frame[0] + subtree_metrics: dict[int, tuple[int, int, tuple[int, ...] | None]] = {} + + def visit(item: object, depth: int) -> tuple[int, int, tuple[int, ...] | None]: item_type = type(item) - depth = len(frames) - 1 if item_type is np.ndarray: - if depth + int(item.ndim) > max_ndim: + item_rank = int(item.ndim) + if depth + item_rank > max_ndim: raise ValueError(rank_error) - frames.pop() - add_to_parent(int(item.size), int(item.ndim)) - continue + count = int(item.size) + if count > _MAX_PREDICTION_CELLS: + raise ValueError(resource_error) + return count, item_rank, tuple(int(size) for size in item.shape) if any(item_type is scalar_type for scalar_type in _TRUSTED_SCALAR_TYPES): if depth > max_ndim: raise ValueError(rank_error) - frames.pop() - add_to_parent(1, 0) - continue + if (item_type is int or item_type in _NUMPY_INTEGER_SCALAR_TYPES) and not ( + _integer_is_exact_float64(item) + ): + raise ValueError(exact_error) + return 1, 0, () if item_type is not list and item_type is not tuple: raise ValueError(error) @@ -88,32 +87,83 @@ def add_to_parent(count: int, subtree_rank: int) -> None: raise ValueError(rank_error) item_id = id(item) - if not bool(frame[4]): - cached = subtree_metrics.get(item_id) - if cached is not None: - count, subtree_rank = cached - if depth + subtree_rank > max_ndim: - raise ValueError(rank_error) - frames.pop() - add_to_parent(count, subtree_rank) - continue - if item_id in active_container_ids: - raise ValueError(error) - active_container_ids.add(item_id) - frame[4] = True - - child_index = int(frame[1]) - if child_index < len(item): - frame[1] = child_index + 1 - frames.append([item[child_index], 0, 0, 1, False]) - continue - - count = int(frame[2]) - subtree_rank = int(frame[3]) - active_container_ids.remove(item_id) - subtree_metrics[item_id] = (count, subtree_rank) - frames.pop() - add_to_parent(count, subtree_rank) + if item_id in active_container_ids: + raise ValueError(error) + cached = subtree_metrics.get(item_id) + if cached is not None: + count, subtree_rank, shape = cached + if depth + subtree_rank > max_ndim: + raise ValueError(rank_error) + return count, subtree_rank, shape + + active_container_ids.add(item_id) + total = 0 + max_child_rank = 0 + first_child_shape: tuple[int, ...] | None = None + rectangular = True + has_child = False + try: + for child in item: + count, child_rank, child_shape = visit(child, depth + 1) + if count > _MAX_PREDICTION_CELLS - total: + raise ValueError(resource_error) + total += count + max_child_rank = max(max_child_rank, child_rank) + if not has_child: + first_child_shape = child_shape + has_child = True + elif child_shape != first_child_shape: + rectangular = False + finally: + active_container_ids.remove(item_id) + + subtree_rank = 1 + max_child_rank + if not has_child: + shape: tuple[int, ...] | None = (0,) + elif rectangular and first_child_shape is not None: + shape = (len(item), *first_child_shape) + elif rectangular and first_child_shape == (): + shape = (len(item),) + else: + shape = None + result = (total, subtree_rank, shape) + subtree_metrics[item_id] = result + return result + + _, _, shape = visit(value, 0) + return shape + + +def _prediction_cells_from_shapes( + theta_shape: tuple[int, ...] | None, + slope_shape: tuple[int, ...] | None, + cat_params_shape: tuple[int, ...] | None, +) -> int | None: + """Return a provable joint output size without materializing evidence.""" + + if ( + theta_shape is None + or slope_shape is None + or cat_params_shape is None + or len(theta_shape) != 1 + or len(slope_shape) != 1 + or len(cat_params_shape) != 2 + or slope_shape[0] <= 0 + or cat_params_shape[0] != slope_shape[0] + or cat_params_shape[1] < 1 + ): + return None + return int(theta_shape[0]) * int(slope_shape[0]) * (int(cat_params_shape[1]) + 1) + + +def _raise_if_oversized_prediction_grid(cells: int | None) -> None: + """Fail closed when a provable prediction grid exceeds the public budget.""" + + if cells is not None and cells > _MAX_PREDICTION_CELLS: + raise ValueError( + f"prediction grid of {cells:,} cells exceeds the " + f"{_MAX_PREDICTION_CELLS:,} prediction-cell limit" + ) def _real_numeric_array(value: object, name: str) -> np.ndarray: @@ -148,11 +198,14 @@ def guarded_predictions(fit, theta): raise TypeError("fit must be a PolytomousFit") # Establish every field's callback/rank/resource contract before any - # one field reaches NumPy materialization. This keeps a later invalid - # item parameter from being preceded by unrelated caller data work. - _preflight(theta, "theta", 1) - _preflight(fit.slope, "fit.slope", 1) - _preflight(fit.cat_params, "fit.cat_params", 2) + # one field reaches NumPy materialization. Shape metadata also makes + # the joint output budget decidable before float64 copies are possible. + theta_shape = _preflight(theta, "theta", 1) + slope_shape = _preflight(fit.slope, "fit.slope", 1) + cat_params_shape = _preflight(fit.cat_params, "fit.cat_params", 2) + _raise_if_oversized_prediction_grid( + _prediction_cells_from_shapes(theta_shape, slope_shape, cat_params_shape) + ) trusted_theta = _real_numeric_array(theta, "theta") trusted_slope = _real_numeric_array(fit.slope, "fit.slope") @@ -168,11 +221,7 @@ def guarded_predictions(fit, theta): ): n_cat = int(trusted_cat_params.shape[1]) + 1 cells = int(trusted_theta.size) * int(trusted_slope.size) * n_cat - if cells > _MAX_PREDICTION_CELLS: - raise ValueError( - f"prediction grid of {cells:,} cells exceeds the " - f"{_MAX_PREDICTION_CELLS:,} prediction-cell limit" - ) + _raise_if_oversized_prediction_grid(cells) trusted_theta = np.ascontiguousarray(trusted_theta) trusted_slope = np.ascontiguousarray(trusted_slope) From a27c613f316514df9da857d155bf3ad341208814 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:18:37 -0700 Subject: [PATCH 22/59] test(predictions): preserve tied GRM boundary categories --- tests/test_polytomous_predictions.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index 80f967c9f..882e12b72 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -39,10 +39,24 @@ def test_public_polytomous_predictions_are_normalized_and_consistent(model, cat_ def test_public_polytomous_predictions_reject_invalid_grm_thresholds(): fit = PolytomousFit("grm", np.array([1.0]), np.array([[-1.0, 1.0]]), 0.0, 0) - with pytest.raises(ValueError, match="strictly decreasing"): + with pytest.raises(ValueError, match="non-increasing"): polytomous_category_probabilities(fit, np.array([0.0])) +def test_public_polytomous_predictions_accept_tied_grm_threshold_boundary(): + """A zero-width GRM category remains a valid fitted-boundary prediction.""" + + fit = PolytomousFit("grm", np.array([1.0]), np.array([[0.0, 0.0]]), 0.0, 0) + probabilities = polytomous_category_probabilities(fit, np.array([0.0])) + expected = polytomous_expected_response(fit, np.array([0.0])) + + assert probabilities.shape == (1, 1, 3) + assert np.all(np.isfinite(probabilities)) + assert np.allclose(probabilities.sum(axis=2), 1.0) + assert probabilities[0, 0, 1] == 0.0 + assert np.allclose(expected, probabilities @ np.arange(3, dtype=np.float64)) + + def test_public_polytomous_predictions_bound_output_before_native_dispatch(monkeypatch): """Reject oversized logical prediction grids before Rust/output allocation.""" From b6d6e99e931cb7dfede04b9bd7b3ca00b8ff50f5 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:22:52 -0700 Subject: [PATCH 23/59] test(predictions): defer GRM boundary-domain change --- tests/test_polytomous_predictions.py | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index 882e12b72..80f967c9f 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -39,24 +39,10 @@ def test_public_polytomous_predictions_are_normalized_and_consistent(model, cat_ def test_public_polytomous_predictions_reject_invalid_grm_thresholds(): fit = PolytomousFit("grm", np.array([1.0]), np.array([[-1.0, 1.0]]), 0.0, 0) - with pytest.raises(ValueError, match="non-increasing"): + with pytest.raises(ValueError, match="strictly decreasing"): polytomous_category_probabilities(fit, np.array([0.0])) -def test_public_polytomous_predictions_accept_tied_grm_threshold_boundary(): - """A zero-width GRM category remains a valid fitted-boundary prediction.""" - - fit = PolytomousFit("grm", np.array([1.0]), np.array([[0.0, 0.0]]), 0.0, 0) - probabilities = polytomous_category_probabilities(fit, np.array([0.0])) - expected = polytomous_expected_response(fit, np.array([0.0])) - - assert probabilities.shape == (1, 1, 3) - assert np.all(np.isfinite(probabilities)) - assert np.allclose(probabilities.sum(axis=2), 1.0) - assert probabilities[0, 0, 1] == 0.0 - assert np.allclose(expected, probabilities @ np.arange(3, dtype=np.float64)) - - def test_public_polytomous_predictions_bound_output_before_native_dispatch(monkeypatch): """Reject oversized logical prediction grids before Rust/output allocation.""" From a797fa61a53d2ce198f63de7189c0e6778b0babd Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:23:13 -0700 Subject: [PATCH 24/59] docs(changelog): record pre-conversion prediction bounds --- docs/changelog.d/1281-polytomous-prediction-admission.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/changelog.d/1281-polytomous-prediction-admission.md b/docs/changelog.d/1281-polytomous-prediction-admission.md index fc395aca0..fe2cd73b1 100644 --- a/docs/changelog.d/1281-polytomous-prediction-admission.md +++ b/docs/changelog.d/1281-polytomous-prediction-admission.md @@ -2,6 +2,6 @@ ## Fixed -- Reject callback-bearing, complex, non-numeric, and lossy GRM/GPCM prediction evidence before NumPy narrowing or compiled-core discovery. -- Preserve trusted exact NumPy and built-in sequence inputs while enforcing the 20,000,000-cell output ceiling before contiguous dense materialization. -- Keep category-probability and expected-score arithmetic in the Rust prediction kernel; Python performs validation, bounded materialization, and marshalling only. +- Reject callback-bearing, complex, non-numeric, and lossy GRM/GPCM prediction evidence before NumPy narrowing or compiled-core discovery, including mixed built-in sequences whose integer identity would be lost by NumPy promotion. +- Preflight theta, slope, and category-parameter rank/resource contracts together and enforce the 20,000,000-cell joint output ceiling from trusted shape metadata before any NumPy materialization or float64 copy. +- Preserve trusted exact NumPy and built-in sequence inputs while keeping category-probability and expected-score arithmetic in the Rust prediction kernel; Python performs validation, bounded materialization, and marshalling only. \ No newline at end of file From 454b50efdc7f77eefecfb8386e2986d25d197154 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:30:07 -0700 Subject: [PATCH 25/59] test(poly): bound prediction category domain --- .../tests/polytomous_predictions_resource.rs | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index 065f68bbd..483585705 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -14,6 +14,20 @@ fn native_predictions_reject_oversized_probability_grid_before_parameter_validat ); } +#[test] +fn native_predictions_reject_category_count_above_fitted_model_ceiling_before_parameter_validation() { + // Fitted polytomous models admit at most POLY_MAX_CAT=64 categories. Keep + // parameter slices empty so this regression also proves prediction replays + // the fitted-model category domain before parameter-dependent validation. + let error = polytomous_predictions(&[0.0], &[], &[], 1, 65, PolyModel::Gpcm) + .expect_err("prediction must reject category counts the fitter cannot produce"); + + assert!( + error.contains("n_cat") && error.contains("2..=64"), + "unexpected native prediction category-domain error: {error}" + ); +} + #[test] fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) From 208ddf78bcd9e96d3ceada13cf0a047e1e920c6a Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:34:10 -0700 Subject: [PATCH 26/59] test(poly): defer native category-domain bound --- .../tests/polytomous_predictions_resource.rs | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index 483585705..065f68bbd 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -14,20 +14,6 @@ fn native_predictions_reject_oversized_probability_grid_before_parameter_validat ); } -#[test] -fn native_predictions_reject_category_count_above_fitted_model_ceiling_before_parameter_validation() { - // Fitted polytomous models admit at most POLY_MAX_CAT=64 categories. Keep - // parameter slices empty so this regression also proves prediction replays - // the fitted-model category domain before parameter-dependent validation. - let error = polytomous_predictions(&[0.0], &[], &[], 1, 65, PolyModel::Gpcm) - .expect_err("prediction must reject category counts the fitter cannot produce"); - - assert!( - error.contains("n_cat") && error.contains("2..=64"), - "unexpected native prediction category-domain error: {error}" - ); -} - #[test] fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) From f69b11c0c982e158de0c7fbe0c816e197a828531 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:34:51 -0700 Subject: [PATCH 27/59] test(poly): reject prediction category domain drift --- tests/test_polytomous_predictions.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index 80f967c9f..439574547 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -43,6 +43,25 @@ def test_public_polytomous_predictions_reject_invalid_grm_thresholds(): polytomous_category_probabilities(fit, np.array([0.0])) +def test_public_polytomous_predictions_reject_category_count_above_fitter_ceiling_before_numpy( + monkeypatch, +): + """A manually constructed fit cannot exceed the fitter's 64-category domain.""" + + fit = PolytomousFit("gpcm", np.array([1.0]), np.zeros((1, 64)), 0.0, 0) + + def unexpected_numpy_materialization(*args, **kwargs): + raise AssertionError("out-of-domain category count reached NumPy materialization") + + def unexpected_core_discovery(): + raise AssertionError("out-of-domain category count reached compiled-core discovery") + + monkeypatch.setattr(prediction_admission.np, "asarray", unexpected_numpy_materialization) + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match=r"n_cat must be in 2\.\.=64"): + polytomous_category_probabilities(fit, np.array([0.0])) + + def test_public_polytomous_predictions_bound_output_before_native_dispatch(monkeypatch): """Reject oversized logical prediction grids before Rust/output allocation.""" From 2ee70059f198bc4186f57e90dfc8c5825410ae6e Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:35:36 -0700 Subject: [PATCH 28/59] fix(poly): replay fitted category ceiling in predictions --- .../_polytomous_prediction_admission.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py index 1e6b27902..7e1958c11 100644 --- a/python/fast_mlsirm/_polytomous_prediction_admission.py +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -6,6 +6,8 @@ import numpy as np +from .config import MAX_POLYTOMOUS_CATEGORIES + _NUMPY_INTEGER_SCALAR_TYPES = ( np.int8, @@ -156,6 +158,16 @@ def _prediction_cells_from_shapes( return int(theta_shape[0]) * int(slope_shape[0]) * (int(cat_params_shape[1]) + 1) +def _raise_if_invalid_category_shape(cat_params_shape: tuple[int, ...] | None) -> None: + """Replay the fitted-model category ceiling from trusted shape metadata.""" + + if cat_params_shape is None or len(cat_params_shape) != 2: + return + n_cat = int(cat_params_shape[1]) + 1 + if n_cat > MAX_POLYTOMOUS_CATEGORIES: + raise ValueError(f"n_cat must be in 2..={MAX_POLYTOMOUS_CATEGORIES}") + + def _raise_if_oversized_prediction_grid(cells: int | None) -> None: """Fail closed when a provable prediction grid exceeds the public budget.""" @@ -199,10 +211,12 @@ def guarded_predictions(fit, theta): # Establish every field's callback/rank/resource contract before any # one field reaches NumPy materialization. Shape metadata also makes - # the joint output budget decidable before float64 copies are possible. + # the fitted category domain and joint output budget decidable before + # float64 copies are possible. theta_shape = _preflight(theta, "theta", 1) slope_shape = _preflight(fit.slope, "fit.slope", 1) cat_params_shape = _preflight(fit.cat_params, "fit.cat_params", 2) + _raise_if_invalid_category_shape(cat_params_shape) _raise_if_oversized_prediction_grid( _prediction_cells_from_shapes(theta_shape, slope_shape, cat_params_shape) ) From fceaab5bb2462b5b0f4db3ab8f66a2bd4f523737 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 08:36:53 -0700 Subject: [PATCH 29/59] docs(poly): record public prediction category ceiling --- .../1292-polytomous-prediction-category-domain.md | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 docs/changelog.d/1292-polytomous-prediction-category-domain.md diff --git a/docs/changelog.d/1292-polytomous-prediction-category-domain.md b/docs/changelog.d/1292-polytomous-prediction-category-domain.md new file mode 100644 index 000000000..30c1da917 --- /dev/null +++ b/docs/changelog.d/1292-polytomous-prediction-category-domain.md @@ -0,0 +1,5 @@ +# Polytomous prediction category-domain admission + +## Fixed + +Public GRM/GPCM prediction admission now rejects manually constructed `PolytomousFit` evidence whose category-parameter shape implies more than the fitter-supported 64 categories before NumPy materialization or compiled-core discovery. Probability and expected-score arithmetic remain Rust-owned; direct native prediction-domain alignment remains tracked by issue #1292. From 6c22a7d5f32a4329b203b7bcaf6185159d450164 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 09:02:48 -0700 Subject: [PATCH 30/59] test(poly): require prediction category bound parity --- .../tests/polytomous_predictions_resource.rs | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index 065f68bbd..eb3d0d0d9 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -14,6 +14,32 @@ fn native_predictions_reject_oversized_probability_grid_before_parameter_validat ); } +#[test] +fn native_predictions_reject_category_count_above_fitter_limit_before_parameter_validation() { + let error = polytomous_predictions(&[0.0], &[], &[], 1, 65, PolyModel::Gpcm) + .expect_err("unsupported native category count must fail closed"); + + assert_eq!(error, "n_cat must be in 2..=64"); +} + +#[test] +fn native_prediction_category_guard_preserves_fitter_upper_boundary() { + let cat_params = vec![0.0; 63]; + let prediction = polytomous_predictions( + &[0.0], + &[1.0], + &cat_params, + 1, + 64, + PolyModel::Gpcm, + ) + .expect("fitter-supported 64-category prediction must remain accepted"); + + assert_eq!(prediction.probabilities.len(), 64); + assert_eq!(prediction.expected.len(), 1); + assert!((prediction.probabilities.iter().sum::() - 1.0).abs() < 1e-12); +} + #[test] fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) From 61f79075a15887dadd085cb5f842eb531a4849cd Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 09:05:08 -0700 Subject: [PATCH 31/59] test(poly): restore native category-bound probe --- .../tests/polytomous_predictions_resource.rs | 26 ------------------- 1 file changed, 26 deletions(-) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index eb3d0d0d9..065f68bbd 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -14,32 +14,6 @@ fn native_predictions_reject_oversized_probability_grid_before_parameter_validat ); } -#[test] -fn native_predictions_reject_category_count_above_fitter_limit_before_parameter_validation() { - let error = polytomous_predictions(&[0.0], &[], &[], 1, 65, PolyModel::Gpcm) - .expect_err("unsupported native category count must fail closed"); - - assert_eq!(error, "n_cat must be in 2..=64"); -} - -#[test] -fn native_prediction_category_guard_preserves_fitter_upper_boundary() { - let cat_params = vec![0.0; 63]; - let prediction = polytomous_predictions( - &[0.0], - &[1.0], - &cat_params, - 1, - 64, - PolyModel::Gpcm, - ) - .expect("fitter-supported 64-category prediction must remain accepted"); - - assert_eq!(prediction.probabilities.len(), 64); - assert_eq!(prediction.expected.len(), 1); - assert!((prediction.probabilities.iter().sum::() - 1.0).abs() < 1e-12); -} - #[test] fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) From 06e4893ae6543122f781407f4fa376ac468bd0dd Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:04:29 -0700 Subject: [PATCH 32/59] test(poly): require native prediction category bound parity --- .../tests/polytomous_predictions_resource.rs | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index 065f68bbd..eb3d0d0d9 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -14,6 +14,32 @@ fn native_predictions_reject_oversized_probability_grid_before_parameter_validat ); } +#[test] +fn native_predictions_reject_category_count_above_fitter_limit_before_parameter_validation() { + let error = polytomous_predictions(&[0.0], &[], &[], 1, 65, PolyModel::Gpcm) + .expect_err("unsupported native category count must fail closed"); + + assert_eq!(error, "n_cat must be in 2..=64"); +} + +#[test] +fn native_prediction_category_guard_preserves_fitter_upper_boundary() { + let cat_params = vec![0.0; 63]; + let prediction = polytomous_predictions( + &[0.0], + &[1.0], + &cat_params, + 1, + 64, + PolyModel::Gpcm, + ) + .expect("fitter-supported 64-category prediction must remain accepted"); + + assert_eq!(prediction.probabilities.len(), 64); + assert_eq!(prediction.expected.len(), 1); + assert!((prediction.probabilities.iter().sum::() - 1.0).abs() < 1e-12); +} + #[test] fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) From cba1fb5eb994b1c4237835648c0fcc66713ed0fa Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:06:16 -0700 Subject: [PATCH 33/59] test(poly): restore native category-bound probe --- .../tests/polytomous_predictions_resource.rs | 26 ------------------- 1 file changed, 26 deletions(-) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index eb3d0d0d9..065f68bbd 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -14,32 +14,6 @@ fn native_predictions_reject_oversized_probability_grid_before_parameter_validat ); } -#[test] -fn native_predictions_reject_category_count_above_fitter_limit_before_parameter_validation() { - let error = polytomous_predictions(&[0.0], &[], &[], 1, 65, PolyModel::Gpcm) - .expect_err("unsupported native category count must fail closed"); - - assert_eq!(error, "n_cat must be in 2..=64"); -} - -#[test] -fn native_prediction_category_guard_preserves_fitter_upper_boundary() { - let cat_params = vec![0.0; 63]; - let prediction = polytomous_predictions( - &[0.0], - &[1.0], - &cat_params, - 1, - 64, - PolyModel::Gpcm, - ) - .expect("fitter-supported 64-category prediction must remain accepted"); - - assert_eq!(prediction.probabilities.len(), 64); - assert_eq!(prediction.expected.len(), 1); - assert!((prediction.probabilities.iter().sum::() - 1.0).abs() < 1e-12); -} - #[test] fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) From fc198cf9579175c02bee2822a0dbf461c35e67d5 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:08:00 -0700 Subject: [PATCH 34/59] test(poly): require native prediction category bound parity --- .../tests/polytomous_predictions_resource.rs | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index 065f68bbd..eb3d0d0d9 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -14,6 +14,32 @@ fn native_predictions_reject_oversized_probability_grid_before_parameter_validat ); } +#[test] +fn native_predictions_reject_category_count_above_fitter_limit_before_parameter_validation() { + let error = polytomous_predictions(&[0.0], &[], &[], 1, 65, PolyModel::Gpcm) + .expect_err("unsupported native category count must fail closed"); + + assert_eq!(error, "n_cat must be in 2..=64"); +} + +#[test] +fn native_prediction_category_guard_preserves_fitter_upper_boundary() { + let cat_params = vec![0.0; 63]; + let prediction = polytomous_predictions( + &[0.0], + &[1.0], + &cat_params, + 1, + 64, + PolyModel::Gpcm, + ) + .expect("fitter-supported 64-category prediction must remain accepted"); + + assert_eq!(prediction.probabilities.len(), 64); + assert_eq!(prediction.expected.len(), 1); + assert!((prediction.probabilities.iter().sum::() - 1.0).abs() < 1e-12); +} + #[test] fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { let prediction = polytomous_predictions(&[0.0], &[1.0], &[0.0], 1, 2, PolyModel::Gpcm) From 0f4c78fe46a02b21e143b80cf4dfad858cfe9b7d Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:13:38 -0700 Subject: [PATCH 35/59] fix(poly): align native prediction category bound --- crates/mlsirm-core/src/poly.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/crates/mlsirm-core/src/poly.rs b/crates/mlsirm-core/src/poly.rs index 2ab9e5f4b..8a60d97b6 100644 --- a/crates/mlsirm-core/src/poly.rs +++ b/crates/mlsirm-core/src/poly.rs @@ -219,6 +219,9 @@ pub fn polytomous_predictions( if n_items == 0 || n_cat < 2 { return Err("n_items must be positive and n_cat must be at least 2".into()); } + if n_cat > POLY_MAX_CAT { + return Err(format!("n_cat must be in 2..={POLY_MAX_CAT}")); + } let cells = crate::checked_mul_usize( theta.len(), n_items, From bf2bbd4e5fe7f7a3e98fe0469fcda192c397b160 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:14:18 -0700 Subject: [PATCH 36/59] docs(changelog): record native prediction category parity --- docs/changelog.d/1292-polytomous-prediction-category-domain.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/changelog.d/1292-polytomous-prediction-category-domain.md b/docs/changelog.d/1292-polytomous-prediction-category-domain.md index 30c1da917..35c722aca 100644 --- a/docs/changelog.d/1292-polytomous-prediction-category-domain.md +++ b/docs/changelog.d/1292-polytomous-prediction-category-domain.md @@ -2,4 +2,4 @@ ## Fixed -Public GRM/GPCM prediction admission now rejects manually constructed `PolytomousFit` evidence whose category-parameter shape implies more than the fitter-supported 64 categories before NumPy materialization or compiled-core discovery. Probability and expected-score arithmetic remain Rust-owned; direct native prediction-domain alignment remains tracked by issue #1292. +GRM/GPCM prediction admission now enforces the fitter-supported `2..=64` category domain at both the public Python boundary and direct Rust `polytomous_predictions()` boundary. Manually constructed `PolytomousFit` evidence above 64 categories fails before NumPy/native work, while direct native requests above `POLY_MAX_CAT` fail before prediction-grid allocation or item-parameter validation. Probability and expected-score arithmetic remain Rust-owned and unchanged. From f1d417a0bba572da5a6cc3dff082589fea241b44 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:34:47 -0700 Subject: [PATCH 37/59] test(prediction): preserve admitted fit metadata --- tests/test_polytomous_predictions.py | 48 ++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/tests/test_polytomous_predictions.py b/tests/test_polytomous_predictions.py index 439574547..d4331b08e 100644 --- a/tests/test_polytomous_predictions.py +++ b/tests/test_polytomous_predictions.py @@ -1,5 +1,7 @@ """Public Rust-owned GRM/GPCM prediction contract.""" +from types import SimpleNamespace + import numpy as np import pytest @@ -37,6 +39,52 @@ def test_public_polytomous_predictions_are_normalized_and_consistent(model, cat_ assert np.allclose(expected, probabilities @ np.arange(3, dtype=np.float64)) +def test_prediction_admission_preserves_fit_metadata_when_normalizing_evidence(): + """Admission must not erase package-owned fit metadata while normalizing arrays.""" + + captured = {} + + def original(fit, theta): + captured["fit"] = fit + captured["theta"] = theta + return np.empty((0,), dtype=np.float64), np.empty((0,), dtype=np.float64) + + module = SimpleNamespace( + PolytomousFit=PolytomousFit, + _polytomous_predictions=original, + ) + prediction_admission.install(module) + + loglik_trace = np.array([-10.0, -8.5, -8.0], dtype=np.float64) + thresholds = np.array([[0.25]], dtype=np.float64) + fit = PolytomousFit( + "gpcm", + [np.float32(1.0)], + [[np.float32(0.0)]], + -8.0, + 7, + converged=True, + termination_reason="tolerance", + loglik_trace=loglik_trace, + final_delta=1e-7, + stopping_tolerance=1e-6, + thresholds=thresholds, + ) + + module._polytomous_predictions(fit, [np.float32(0.0)]) + trusted = captured["fit"] + + assert trusted.converged is True + assert trusted.termination_reason == "tolerance" + assert trusted.loglik_trace is loglik_trace + assert trusted.final_delta == pytest.approx(1e-7) + assert trusted.stopping_tolerance == pytest.approx(1e-6) + assert trusted.thresholds is thresholds + assert trusted.slope.dtype == np.float64 + assert trusted.cat_params.dtype == np.float64 + assert captured["theta"].dtype == np.float64 + + def test_public_polytomous_predictions_reject_invalid_grm_thresholds(): fit = PolytomousFit("grm", np.array([1.0]), np.array([[-1.0, 1.0]]), 0.0, 0) with pytest.raises(ValueError, match="strictly decreasing"): From 0664286a425620faba8c551c778127ece2856b81 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:35:23 -0700 Subject: [PATCH 38/59] fix(prediction): preserve admitted fit metadata --- .../fast_mlsirm/_polytomous_prediction_admission.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py index 7e1958c11..d49382e1c 100644 --- a/python/fast_mlsirm/_polytomous_prediction_admission.py +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -2,6 +2,7 @@ from __future__ import annotations +from dataclasses import replace from types import ModuleType import numpy as np @@ -240,12 +241,10 @@ def guarded_predictions(fit, theta): trusted_theta = np.ascontiguousarray(trusted_theta) trusted_slope = np.ascontiguousarray(trusted_slope) trusted_cat_params = np.ascontiguousarray(trusted_cat_params) - trusted_fit = module.PolytomousFit( - fit.model, - trusted_slope, - trusted_cat_params, - fit.loglik, - fit.n_iter, + trusted_fit = replace( + fit, + slope=trusted_slope, + cat_params=trusted_cat_params, ) return current(trusted_fit, trusted_theta) From cb9b4a66e33effa4588a8ef9e780c417cde9379e Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 11:35:41 -0700 Subject: [PATCH 39/59] docs(changelog): record prediction metadata preservation --- docs/changelog.d/1296-polytomous-prediction-fit-metadata.md | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 docs/changelog.d/1296-polytomous-prediction-fit-metadata.md diff --git a/docs/changelog.d/1296-polytomous-prediction-fit-metadata.md b/docs/changelog.d/1296-polytomous-prediction-fit-metadata.md new file mode 100644 index 000000000..1c4a27fb7 --- /dev/null +++ b/docs/changelog.d/1296-polytomous-prediction-fit-metadata.md @@ -0,0 +1,5 @@ +# Preserve polytomous prediction fit metadata + +## Fixed + +- Preserve all package-owned `PolytomousFit` convergence, trace, stopping, and threshold metadata when the public GRM/GPCM prediction admission boundary normalizes slope and category arrays before Rust dispatch. From 602c8dd7750d032214c6f3130941d5c12da13c6d Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 12:04:23 -0700 Subject: [PATCH 40/59] test(grm): reproduce collapsed initialization thresholds --- .../tests/polytomous_predictions_resource.rs | 45 ++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index eb3d0d0d9..3cc83a460 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -1,4 +1,4 @@ -use mlsirm_core::poly::{polytomous_predictions, PolyModel}; +use mlsirm_core::poly::{fit_poly_unidim, polytomous_predictions, PolyModel}; #[test] fn native_predictions_reject_oversized_probability_grid_before_parameter_validation() { @@ -50,3 +50,46 @@ fn native_prediction_resource_guard_preserves_small_gpcm_predictions() { assert!((prediction.probabilities.iter().sum::() - 1.0).abs() < 1e-12); assert!((prediction.expected[0] - 0.5).abs() < 1e-12); } + +#[test] +fn fitted_grm_with_empty_middle_category_stays_inside_prediction_domain() { + // Additive smoothing makes both cumulative probabilities finite, but the + // former independent 1e-4 clamp mapped them to the same upper boundary for + // this all-top-category item. A fitted GRM must never return that tied + // threshold vector because prediction/scoring correctly require strict + // threshold order. + let responses = vec![2usize; 32]; + let fit = fit_poly_unidim( + &responses, + None, + responses.len(), + 1, + 3, + PolyModel::Grm, + 7, + 2, + 1e-8, + ) + .expect("degenerate observed categories must not produce an invalid GRM fit"); + + let thresholds = &fit.cat_params[0]; + assert!( + thresholds.windows(2).all(|pair| pair[0] > pair[1]), + "fitted GRM thresholds must remain strictly decreasing: {thresholds:?}" + ); + + let prediction = polytomous_predictions( + &[-2.0, 0.0, 2.0], + &fit.slope, + thresholds, + 1, + 3, + PolyModel::Grm, + ) + .expect("a returned GRM fit must be directly accepted by prediction"); + + for probabilities in prediction.probabilities.chunks_exact(3) { + assert!(probabilities.iter().all(|value| value.is_finite() && *value >= 0.0)); + assert!((probabilities.iter().sum::() - 1.0).abs() < 1e-12); + } +} \ No newline at end of file From 2eeb3d42479779edcb58890b121ce73f4f321982 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 12:11:18 -0700 Subject: [PATCH 41/59] fix(grm): preserve ordered initialization thresholds --- crates/mlsirm-core/src/poly.rs | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/crates/mlsirm-core/src/poly.rs b/crates/mlsirm-core/src/poly.rs index 8a60d97b6..24a898925 100644 --- a/crates/mlsirm-core/src/poly.rs +++ b/crates/mlsirm-core/src/poly.rs @@ -551,8 +551,10 @@ pub fn fit_poly_unidim( let mut cum = 0.0_f64; for k in (1..n_cat).rev() { cum += freq[k]; - let c = cum.clamp(1e-4, 1.0 - 1e-4); - params[i][k] = (c / (1.0 - c)).ln(); + if !(0.0 < cum && cum < 1.0) { + return Err("smoothed GRM cumulative probability must be in (0, 1)".into()); + } + params[i][k] = (cum / (1.0 - cum)).ln(); } } } @@ -1426,8 +1428,10 @@ pub fn fit_poly_multigroup( let mut cum = 0.0_f64; for k in (1..n_cat).rev() { cum += freq[k]; - let c = cum.clamp(1e-4, 1.0 - 1e-4); - params[i][k] = (c / (1.0 - c)).ln(); + if !(0.0 < cum && cum < 1.0) { + return Err("smoothed GRM cumulative probability must be in (0, 1)".into()); + } + params[i][k] = (cum / (1.0 - cum)).ln(); } } } @@ -2362,9 +2366,9 @@ pub fn poly_s_x2( let lp = match model { PolyModel::Gpcm => { let scores: Vec = (0..n_cat).map(|c| c as f64).collect(); - let mut intercepts = vec![0.0_f64; n_cat]; - intercepts[1..].copy_from_slice(cp); - gpcm_logprobs(base, &scores, &intercepts) + let mut ic = vec![0.0_f64; n_cat]; + ic[1..].copy_from_slice(cp); + gpcm_logprobs(base, &scores, &ic) } PolyModel::Grm => grm_logprobs(base, cp), }; From 1587c501b7e2a79415d3b72cd6e0b0f8f7ea5505 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 12:15:38 -0700 Subject: [PATCH 42/59] chore(grm): restore unrelated S-X2 naming --- crates/mlsirm-core/src/poly.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/crates/mlsirm-core/src/poly.rs b/crates/mlsirm-core/src/poly.rs index 24a898925..c2b93e159 100644 --- a/crates/mlsirm-core/src/poly.rs +++ b/crates/mlsirm-core/src/poly.rs @@ -2366,9 +2366,9 @@ pub fn poly_s_x2( let lp = match model { PolyModel::Gpcm => { let scores: Vec = (0..n_cat).map(|c| c as f64).collect(); - let mut ic = vec![0.0_f64; n_cat]; - ic[1..].copy_from_slice(cp); - gpcm_logprobs(base, &scores, &ic) + let mut intercepts = vec![0.0_f64; n_cat]; + intercepts[1..].copy_from_slice(cp); + gpcm_logprobs(base, &scores, &intercepts) } PolyModel::Grm => grm_logprobs(base, cp), }; From 492176f648f40cb0b09aab5066f7875a45e722f7 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 12:17:07 -0700 Subject: [PATCH 43/59] docs(grm): record strict threshold initialization --- docs/changelog.d/1297-grm-threshold-initialization.md | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 docs/changelog.d/1297-grm-threshold-initialization.md diff --git a/docs/changelog.d/1297-grm-threshold-initialization.md b/docs/changelog.d/1297-grm-threshold-initialization.md new file mode 100644 index 000000000..343fffe6b --- /dev/null +++ b/docs/changelog.d/1297-grm-threshold-initialization.md @@ -0,0 +1,5 @@ +# Preserve strict GRM threshold order during initialization + +## Fixed + +- Preserve strictly decreasing finite GRM category thresholds for sparse or collapsed observed-category patterns by using the existing positive-pseudocount cumulative frequencies directly instead of independently clipping adjacent cumulative probabilities onto the same boundary. Returned GRM fits therefore remain inside the shared scoring and prediction parameter domain without changing Samejima category-probability arithmetic or GPCM behavior. From 1f93c471665fb9d7f87079fcb1593abe2acc9806 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 15:06:39 -0700 Subject: [PATCH 44/59] test(prediction): replay category ceiling in raw adapter --- .../test_polytomous_prediction_raw_domain.py | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 tests/test_polytomous_prediction_raw_domain.py diff --git a/tests/test_polytomous_prediction_raw_domain.py b/tests/test_polytomous_prediction_raw_domain.py new file mode 100644 index 000000000..5b58bea61 --- /dev/null +++ b/tests/test_polytomous_prediction_raw_domain.py @@ -0,0 +1,25 @@ +"""Defense-in-depth contract for the internal polytomous prediction adapter.""" + +import numpy as np +import pytest + +import fast_mlsirm.polytomous as polytomous_module + + +def test_raw_prediction_replays_fitter_category_ceiling_before_native(monkeypatch): + """The raw adapter must reject 65 categories before compiled-core discovery.""" + + fit = polytomous_module.PolytomousFit( + "gpcm", + np.array([1.0]), + np.zeros((1, 64)), + 0.0, + 0, + ) + + def unexpected_core_discovery(): + raise AssertionError("out-of-domain category count reached compiled-core discovery") + + monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + with pytest.raises(ValueError, match=r"n_cat must be in 2\.\.=64"): + polytomous_module._polytomous_predictions(fit, np.array([0.0])) From 6004e46e83e8d8c1b0ae7be5bf88f72a62ad29ec Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 15:08:02 -0700 Subject: [PATCH 45/59] test(prediction): exercise uninstalled raw category guard --- tests/test_polytomous_prediction_raw_domain.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/test_polytomous_prediction_raw_domain.py b/tests/test_polytomous_prediction_raw_domain.py index 5b58bea61..653304360 100644 --- a/tests/test_polytomous_prediction_raw_domain.py +++ b/tests/test_polytomous_prediction_raw_domain.py @@ -1,5 +1,7 @@ """Defense-in-depth contract for the internal polytomous prediction adapter.""" +from importlib import reload + import numpy as np import pytest @@ -7,9 +9,10 @@ def test_raw_prediction_replays_fitter_category_ceiling_before_native(monkeypatch): - """The raw adapter must reject 65 categories before compiled-core discovery.""" + """The uninstalled raw adapter must reject 65 categories before Rust discovery.""" - fit = polytomous_module.PolytomousFit( + raw_module = reload(polytomous_module) + fit = raw_module.PolytomousFit( "gpcm", np.array([1.0]), np.zeros((1, 64)), @@ -20,6 +23,6 @@ def test_raw_prediction_replays_fitter_category_ceiling_before_native(monkeypatc def unexpected_core_discovery(): raise AssertionError("out-of-domain category count reached compiled-core discovery") - monkeypatch.setattr(polytomous_module, "_core_module", unexpected_core_discovery) + monkeypatch.setattr(raw_module, "_core_module", unexpected_core_discovery) with pytest.raises(ValueError, match=r"n_cat must be in 2\.\.=64"): - polytomous_module._polytomous_predictions(fit, np.array([0.0])) + raw_module._polytomous_predictions(fit, np.array([0.0])) From 3de8dbcaf97a8aa48e0a5489b150034d64f3a819 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 15:09:30 -0700 Subject: [PATCH 46/59] test(prediction): defer raw adapter ceiling replay --- .../test_polytomous_prediction_raw_domain.py | 28 ------------------- 1 file changed, 28 deletions(-) delete mode 100644 tests/test_polytomous_prediction_raw_domain.py diff --git a/tests/test_polytomous_prediction_raw_domain.py b/tests/test_polytomous_prediction_raw_domain.py deleted file mode 100644 index 653304360..000000000 --- a/tests/test_polytomous_prediction_raw_domain.py +++ /dev/null @@ -1,28 +0,0 @@ -"""Defense-in-depth contract for the internal polytomous prediction adapter.""" - -from importlib import reload - -import numpy as np -import pytest - -import fast_mlsirm.polytomous as polytomous_module - - -def test_raw_prediction_replays_fitter_category_ceiling_before_native(monkeypatch): - """The uninstalled raw adapter must reject 65 categories before Rust discovery.""" - - raw_module = reload(polytomous_module) - fit = raw_module.PolytomousFit( - "gpcm", - np.array([1.0]), - np.zeros((1, 64)), - 0.0, - 0, - ) - - def unexpected_core_discovery(): - raise AssertionError("out-of-domain category count reached compiled-core discovery") - - monkeypatch.setattr(raw_module, "_core_module", unexpected_core_discovery) - with pytest.raises(ValueError, match=r"n_cat must be in 2\.\.=64"): - raw_module._polytomous_predictions(fit, np.array([0.0])) From a0f57302be0f929cbe2ffbd439fe0c90b721e702 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 15:09:44 -0700 Subject: [PATCH 47/59] test(prediction): reject nonfinite evidence before delegate --- .../test_polytomous_prediction_finiteness.py | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 tests/test_polytomous_prediction_finiteness.py diff --git a/tests/test_polytomous_prediction_finiteness.py b/tests/test_polytomous_prediction_finiteness.py new file mode 100644 index 000000000..e44a4270f --- /dev/null +++ b/tests/test_polytomous_prediction_finiteness.py @@ -0,0 +1,36 @@ +"""Finiteness ordering for public GRM/GPCM prediction admission.""" + +from types import SimpleNamespace + +import numpy as np +import pytest + +import fast_mlsirm._polytomous_prediction_admission as prediction_admission +from fast_mlsirm.polytomous import PolytomousFit + + +@pytest.mark.parametrize( + ("theta", "slope", "cat_params", "message"), + [ + ([np.inf], [1.0], [[0.0]], "theta must be a non-empty finite 1-D array"), + ([0.0], [np.inf], [[0.0]], "fit item parameters must be finite"), + ([0.0], [1.0], [[np.nan]], "fit item parameters must be finite"), + ], +) +def test_prediction_admission_rejects_nonfinite_evidence_before_delegate( + theta, slope, cat_params, message +): + """Non-finite trusted arrays must fail before the raw prediction delegate.""" + + def unexpected_delegate(*args, **kwargs): + raise AssertionError("non-finite evidence reached raw prediction delegate") + + module = SimpleNamespace( + PolytomousFit=PolytomousFit, + _polytomous_predictions=unexpected_delegate, + ) + prediction_admission.install(module) + fit = PolytomousFit("gpcm", slope, cat_params, 0.0, 0) + + with pytest.raises(ValueError, match=message): + module._polytomous_predictions(fit, theta) From b51a0720c9a6014656d8f6d07001dc05df4b722a Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 15:10:13 -0700 Subject: [PATCH 48/59] fix(prediction): reject nonfinite evidence before delegate --- python/fast_mlsirm/_polytomous_prediction_admission.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/python/fast_mlsirm/_polytomous_prediction_admission.py b/python/fast_mlsirm/_polytomous_prediction_admission.py index d49382e1c..ca21f0f1f 100644 --- a/python/fast_mlsirm/_polytomous_prediction_admission.py +++ b/python/fast_mlsirm/_polytomous_prediction_admission.py @@ -225,6 +225,12 @@ def guarded_predictions(fit, theta): trusted_theta = _real_numeric_array(theta, "theta") trusted_slope = _real_numeric_array(fit.slope, "fit.slope") trusted_cat_params = _real_numeric_array(fit.cat_params, "fit.cat_params") + if not np.all(np.isfinite(trusted_theta)): + raise ValueError("theta must be a non-empty finite 1-D array") + if not np.all(np.isfinite(trusted_slope)) or not np.all( + np.isfinite(trusted_cat_params) + ): + raise ValueError("fit item parameters must be finite") if ( trusted_theta.ndim == 1 From b60539fd1a9d72c48bdfddd2badac8116eabf07c Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 15:10:28 -0700 Subject: [PATCH 49/59] docs(changelog): record prediction finiteness admission --- docs/changelog.d/1281-polytomous-prediction-admission.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/changelog.d/1281-polytomous-prediction-admission.md b/docs/changelog.d/1281-polytomous-prediction-admission.md index fe2cd73b1..1b2061af0 100644 --- a/docs/changelog.d/1281-polytomous-prediction-admission.md +++ b/docs/changelog.d/1281-polytomous-prediction-admission.md @@ -2,6 +2,6 @@ ## Fixed -- Reject callback-bearing, complex, non-numeric, and lossy GRM/GPCM prediction evidence before NumPy narrowing or compiled-core discovery, including mixed built-in sequences whose integer identity would be lost by NumPy promotion. +- Reject callback-bearing, complex, non-numeric, lossy, and non-finite GRM/GPCM prediction evidence before the raw prediction delegate or compiled-core discovery, including mixed built-in sequences whose integer identity would be lost by NumPy promotion. - Preflight theta, slope, and category-parameter rank/resource contracts together and enforce the 20,000,000-cell joint output ceiling from trusted shape metadata before any NumPy materialization or float64 copy. - Preserve trusted exact NumPy and built-in sequence inputs while keeping category-probability and expected-score arithmetic in the Rust prediction kernel; Python performs validation, bounded materialization, and marshalling only. \ No newline at end of file From 9bf2a850a167950d96a5db288eb7a3054fbe233c Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 17:09:23 -0700 Subject: [PATCH 50/59] test(poly): prove GRM M-step preserves prediction domain --- .../tests/polytomous_predictions_resource.rs | 50 +++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs index 3cc83a460..357b8966c 100644 --- a/crates/mlsirm-core/tests/polytomous_predictions_resource.rs +++ b/crates/mlsirm-core/tests/polytomous_predictions_resource.rs @@ -92,4 +92,54 @@ fn fitted_grm_with_empty_middle_category_stays_inside_prediction_domain() { assert!(probabilities.iter().all(|value| value.is_finite() && *value >= 0.0)); assert!((probabilities.iter().sum::() - 1.0).abs() < 1e-12); } +} + +#[test] +fn fitted_grm_with_multiple_empty_middle_categories_remains_prediction_valid_after_em() { + // Stress more than one unobserved interior category and require the fitted + // parameters returned after at least one EM M-step to remain in the same + // strict threshold domain enforced by the public/native prediction path. + let responses: Vec = (0..64) + .map(|person| if person < 32 { 0 } else { 3 }) + .collect(); + let fit = fit_poly_unidim( + &responses, + None, + responses.len(), + 1, + 4, + PolyModel::Grm, + 7, + 2, + 1e-8, + ) + .expect("sparse four-category GRM evidence must retain a prediction-valid fit"); + + assert!( + fit.n_iter >= 1 && fit.loglik_trace.len() >= 2, + "regression must exercise at least one EM M-step: n_iter={}, trace_len={}", + fit.n_iter, + fit.loglik_trace.len() + ); + let thresholds = &fit.cat_params[0]; + assert_eq!(thresholds.len(), 3); + assert!( + thresholds.windows(2).all(|pair| pair[0] > pair[1]), + "post-M-step GRM thresholds must stay strictly decreasing: {thresholds:?}" + ); + + let prediction = polytomous_predictions( + &[-4.0, 0.0, 4.0], + &fit.slope, + thresholds, + 1, + 4, + PolyModel::Grm, + ) + .expect("a post-M-step GRM fit must be directly accepted by prediction"); + + for probabilities in prediction.probabilities.chunks_exact(4) { + assert!(probabilities.iter().all(|value| value.is_finite() && *value >= 0.0)); + assert!((probabilities.iter().sum::() - 1.0).abs() < 1e-12); + } } \ No newline at end of file From dc0d55b347f3cd208d02b9ced9c0099ce5d4c642 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 19:05:26 -0700 Subject: [PATCH 51/59] test(poly): require raw prediction category-domain replay --- .../test_polytomous_prediction_raw_domain.py | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 tests/test_polytomous_prediction_raw_domain.py diff --git a/tests/test_polytomous_prediction_raw_domain.py b/tests/test_polytomous_prediction_raw_domain.py new file mode 100644 index 000000000..ad72a5cd6 --- /dev/null +++ b/tests/test_polytomous_prediction_raw_domain.py @@ -0,0 +1,48 @@ +"""Defense-in-depth contract for the raw Python polytomous prediction helper.""" + +from __future__ import annotations + +import subprocess +import sys +import textwrap + + +def test_raw_polytomous_predictions_reject_category_count_above_fitter_domain() -> None: + """The raw helper must replay the 64-category fitter ceiling before core discovery.""" + + program = textwrap.dedent( + """ + import importlib + import numpy as np + import fast_mlsirm.polytomous as module + + module = importlib.reload(module) + + def unexpected_core_discovery(): + raise AssertionError("out-of-domain raw prediction reached compiled-core discovery") + + module._core_module = unexpected_core_discovery + fit = module.PolytomousFit( + "gpcm", + np.array([1.0]), + np.zeros((1, 64), dtype=np.float64), + 0.0, + 0, + ) + + try: + module._polytomous_predictions(fit, np.array([0.0])) + except ValueError as error: + assert str(error) == "n_cat must be in 2..=64" + else: + raise AssertionError("raw prediction accepted 65 categories") + """ + ) + completed = subprocess.run( + [sys.executable, "-c", program], + check=False, + capture_output=True, + text=True, + timeout=30, + ) + assert completed.returncode == 0, completed.stderr From af26250e30c774406aed31b9a42cdd577e49412c Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 19:06:05 -0700 Subject: [PATCH 52/59] test(poly): defer raw helper domain proof --- .../test_polytomous_prediction_raw_domain.py | 48 ------------------- 1 file changed, 48 deletions(-) delete mode 100644 tests/test_polytomous_prediction_raw_domain.py diff --git a/tests/test_polytomous_prediction_raw_domain.py b/tests/test_polytomous_prediction_raw_domain.py deleted file mode 100644 index ad72a5cd6..000000000 --- a/tests/test_polytomous_prediction_raw_domain.py +++ /dev/null @@ -1,48 +0,0 @@ -"""Defense-in-depth contract for the raw Python polytomous prediction helper.""" - -from __future__ import annotations - -import subprocess -import sys -import textwrap - - -def test_raw_polytomous_predictions_reject_category_count_above_fitter_domain() -> None: - """The raw helper must replay the 64-category fitter ceiling before core discovery.""" - - program = textwrap.dedent( - """ - import importlib - import numpy as np - import fast_mlsirm.polytomous as module - - module = importlib.reload(module) - - def unexpected_core_discovery(): - raise AssertionError("out-of-domain raw prediction reached compiled-core discovery") - - module._core_module = unexpected_core_discovery - fit = module.PolytomousFit( - "gpcm", - np.array([1.0]), - np.zeros((1, 64), dtype=np.float64), - 0.0, - 0, - ) - - try: - module._polytomous_predictions(fit, np.array([0.0])) - except ValueError as error: - assert str(error) == "n_cat must be in 2..=64" - else: - raise AssertionError("raw prediction accepted 65 categories") - """ - ) - completed = subprocess.run( - [sys.executable, "-c", program], - check=False, - capture_output=True, - text=True, - timeout=30, - ) - assert completed.returncode == 0, completed.stderr From cf6ac2b642fd6f4fb524b89d43c6547c150f28d1 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 21:07:43 -0700 Subject: [PATCH 53/59] test(poly): require raw helper category ceiling --- .../test_polytomous_prediction_raw_domain.py | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 tests/test_polytomous_prediction_raw_domain.py diff --git a/tests/test_polytomous_prediction_raw_domain.py b/tests/test_polytomous_prediction_raw_domain.py new file mode 100644 index 000000000..ad72a5cd6 --- /dev/null +++ b/tests/test_polytomous_prediction_raw_domain.py @@ -0,0 +1,48 @@ +"""Defense-in-depth contract for the raw Python polytomous prediction helper.""" + +from __future__ import annotations + +import subprocess +import sys +import textwrap + + +def test_raw_polytomous_predictions_reject_category_count_above_fitter_domain() -> None: + """The raw helper must replay the 64-category fitter ceiling before core discovery.""" + + program = textwrap.dedent( + """ + import importlib + import numpy as np + import fast_mlsirm.polytomous as module + + module = importlib.reload(module) + + def unexpected_core_discovery(): + raise AssertionError("out-of-domain raw prediction reached compiled-core discovery") + + module._core_module = unexpected_core_discovery + fit = module.PolytomousFit( + "gpcm", + np.array([1.0]), + np.zeros((1, 64), dtype=np.float64), + 0.0, + 0, + ) + + try: + module._polytomous_predictions(fit, np.array([0.0])) + except ValueError as error: + assert str(error) == "n_cat must be in 2..=64" + else: + raise AssertionError("raw prediction accepted 65 categories") + """ + ) + completed = subprocess.run( + [sys.executable, "-c", program], + check=False, + capture_output=True, + text=True, + timeout=30, + ) + assert completed.returncode == 0, completed.stderr From 71d2dba61a66687dcd473f2eb3e96f5c077b6526 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 21:09:01 -0700 Subject: [PATCH 54/59] test(poly): defer raw helper domain proof --- .../test_polytomous_prediction_raw_domain.py | 48 ------------------- 1 file changed, 48 deletions(-) delete mode 100644 tests/test_polytomous_prediction_raw_domain.py diff --git a/tests/test_polytomous_prediction_raw_domain.py b/tests/test_polytomous_prediction_raw_domain.py deleted file mode 100644 index ad72a5cd6..000000000 --- a/tests/test_polytomous_prediction_raw_domain.py +++ /dev/null @@ -1,48 +0,0 @@ -"""Defense-in-depth contract for the raw Python polytomous prediction helper.""" - -from __future__ import annotations - -import subprocess -import sys -import textwrap - - -def test_raw_polytomous_predictions_reject_category_count_above_fitter_domain() -> None: - """The raw helper must replay the 64-category fitter ceiling before core discovery.""" - - program = textwrap.dedent( - """ - import importlib - import numpy as np - import fast_mlsirm.polytomous as module - - module = importlib.reload(module) - - def unexpected_core_discovery(): - raise AssertionError("out-of-domain raw prediction reached compiled-core discovery") - - module._core_module = unexpected_core_discovery - fit = module.PolytomousFit( - "gpcm", - np.array([1.0]), - np.zeros((1, 64), dtype=np.float64), - 0.0, - 0, - ) - - try: - module._polytomous_predictions(fit, np.array([0.0])) - except ValueError as error: - assert str(error) == "n_cat must be in 2..=64" - else: - raise AssertionError("raw prediction accepted 65 categories") - """ - ) - completed = subprocess.run( - [sys.executable, "-c", program], - check=False, - capture_output=True, - text=True, - timeout=30, - ) - assert completed.returncode == 0, completed.stderr From fdce0d6fe45a6ff89f5ba6362435df7821d1516f Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 22:05:31 -0700 Subject: [PATCH 55/59] test(prediction): require raw category-domain replay --- ...lytomous_prediction_raw_category_domain.py | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 tests/test_polytomous_prediction_raw_category_domain.py diff --git a/tests/test_polytomous_prediction_raw_category_domain.py b/tests/test_polytomous_prediction_raw_category_domain.py new file mode 100644 index 000000000..412c73946 --- /dev/null +++ b/tests/test_polytomous_prediction_raw_category_domain.py @@ -0,0 +1,47 @@ +"""Regression for raw polytomous prediction category-domain replay.""" + +from __future__ import annotations + +import subprocess +import sys + + +def test_raw_prediction_replays_category_ceiling_before_core_discovery() -> None: + """The private raw helper must reject 65 categories before Rust discovery.""" + + script = r''' +import importlib +import numpy as np +import fast_mlsirm.polytomous as poly + +poly = importlib.reload(poly) + + +def core_forbidden(): + raise AssertionError("compiled core discovery must not occur") + + +poly._core_module = core_forbidden +fit = poly.PolytomousFit( + model="grm", + slope=np.ones(1, dtype=np.float64), + cat_params=np.zeros((1, 64), dtype=np.float64), + loglik=0.0, + n_iter=1, +) + +try: + poly._polytomous_predictions(fit, np.array([0.0], dtype=np.float64)) +except ValueError as exc: + assert str(exc) == "n_cat must be in 2..=64", str(exc) +else: + raise AssertionError("expected category-domain ValueError") +''' + completed = subprocess.run( + [sys.executable, "-c", script], + check=False, + capture_output=True, + text=True, + timeout=30, + ) + assert completed.returncode == 0, completed.stderr From b9ecf432bc257a9d99e2b3e1358331d72b971295 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 22:11:18 -0700 Subject: [PATCH 56/59] fix(prediction): replay raw category domain --- python/fast_mlsirm/polytomous.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/fast_mlsirm/polytomous.py b/python/fast_mlsirm/polytomous.py index 05e9e413d..cd4ce2354 100644 --- a/python/fast_mlsirm/polytomous.py +++ b/python/fast_mlsirm/polytomous.py @@ -105,12 +105,15 @@ def _polytomous_predictions( or cat_params.shape[1] < 1 ): raise ValueError("fit.cat_params must be n_items x (n_cat - 1)") + n_cat = int(cat_params.shape[1]) + 1 + if n_cat > MAX_POLYTOMOUS_CATEGORIES: + raise ValueError(f"n_cat must be in 2..={MAX_POLYTOMOUS_CATEGORIES}") if not np.all(np.isfinite(slope)) or not np.all(np.isfinite(cat_params)): raise ValueError("fit item parameters must be finite") model = fit.model.lower() if type(fit.model) is str else "" if model not in VALID_POLY_MODELS: raise ValueError(f"fit.model must be one of {sorted(VALID_POLY_MODELS)}") - prediction_cells = int(th.size) * int(slope.size) * int(cat_params.shape[1] + 1) + prediction_cells = int(th.size) * int(slope.size) * n_cat if prediction_cells > 20_000_000: raise ValueError( f"prediction grid of {prediction_cells:,} cells exceeds the " @@ -124,11 +127,11 @@ def _polytomous_predictions( slope, cat_params.reshape(-1), int(slope.size), - int(cat_params.shape[1] + 1), + n_cat, model, ) probabilities = np.asarray(result["probabilities"], dtype=np.float64).reshape( - th.size, slope.size, cat_params.shape[1] + 1 + th.size, slope.size, n_cat ) expected = np.asarray(result["expected"], dtype=np.float64).reshape( th.size, slope.size From 1dc38914439833b180983cc825e5340034c311b3 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 22:11:59 -0700 Subject: [PATCH 57/59] test(prediction): preserve raw 64-category boundary --- ...lytomous_prediction_raw_category_domain.py | 75 +++++++++++++++++-- 1 file changed, 67 insertions(+), 8 deletions(-) diff --git a/tests/test_polytomous_prediction_raw_category_domain.py b/tests/test_polytomous_prediction_raw_category_domain.py index 412c73946..05568e2da 100644 --- a/tests/test_polytomous_prediction_raw_category_domain.py +++ b/tests/test_polytomous_prediction_raw_category_domain.py @@ -6,10 +6,22 @@ import sys +def _run_isolated(script: str) -> None: + completed = subprocess.run( + [sys.executable, "-c", script], + check=False, + capture_output=True, + text=True, + timeout=30, + ) + assert completed.returncode == 0, completed.stderr + + def test_raw_prediction_replays_category_ceiling_before_core_discovery() -> None: """The private raw helper must reject 65 categories before Rust discovery.""" - script = r''' + _run_isolated( + r''' import importlib import numpy as np import fast_mlsirm.polytomous as poly @@ -37,11 +49,58 @@ def core_forbidden(): else: raise AssertionError("expected category-domain ValueError") ''' - completed = subprocess.run( - [sys.executable, "-c", script], - check=False, - capture_output=True, - text=True, - timeout=30, ) - assert completed.returncode == 0, completed.stderr + + +def test_raw_prediction_preserves_supported_64_category_boundary() -> None: + """The raw helper must still delegate the fitter-supported 64-category edge.""" + + _run_isolated( + r''' +import importlib +import numpy as np +import fast_mlsirm.polytomous as poly + +poly = importlib.reload(poly) + + +class Core: + def polytomous_predictions( + self, + theta, + slope, + cat_params, + n_items, + n_cat, + model, + ): + assert n_items == 1 + assert n_cat == 64 + assert model == "gpcm" + assert theta.shape == (1,) + assert slope.shape == (1,) + assert cat_params.shape == (63,) + return { + "probabilities": np.full(64, 1.0 / 64.0, dtype=np.float64), + "expected": np.array([31.5], dtype=np.float64), + } + + +poly._core_module = lambda: Core() +fit = poly.PolytomousFit( + model="gpcm", + slope=np.ones(1, dtype=np.float64), + cat_params=np.zeros((1, 63), dtype=np.float64), + loglik=0.0, + n_iter=1, +) +probabilities, expected = poly._polytomous_predictions( + fit, + np.array([0.0], dtype=np.float64), +) +assert probabilities.shape == (1, 1, 64) +assert np.allclose(probabilities.sum(axis=2), 1.0) +assert expected.shape == (1, 1) +assert expected[0, 0] == 31.5 +''' + ) From 84ffcb5bf407aaad79a3ebd89aa1100ebf6e7a54 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sun, 23 Aug 2026 22:12:26 -0700 Subject: [PATCH 58/59] docs(changelog): record raw prediction category replay --- docs/changelog.d/1307-polytomous-raw-category-domain.md | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 docs/changelog.d/1307-polytomous-raw-category-domain.md diff --git a/docs/changelog.d/1307-polytomous-raw-category-domain.md b/docs/changelog.d/1307-polytomous-raw-category-domain.md new file mode 100644 index 000000000..f415a8c52 --- /dev/null +++ b/docs/changelog.d/1307-polytomous-raw-category-domain.md @@ -0,0 +1,5 @@ +# Polytomous raw prediction category-domain replay + +## Fixed + +- Replayed the fitter-supported `2..=64` category domain in the package-private Python prediction helper before resource calculation or compiled-core discovery, while preserving the valid 64-category boundary and keeping GRM/GPCM probability arithmetic Rust-owned. From cbb0df16efcc194bcb47c75f96e8049a5ed47270 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Tue, 25 Aug 2026 04:22:08 -0700 Subject: [PATCH 59/59] feat(interaction-map): own residual factorization in Rust (#1387) * fix(release): resync CHANGELOG fragments after v0.8.0 cut PRs merged after the v0.8.0 release-cut PR (#966) added new docs/changelog.d fragments without re-running the renderer, and the already-released release-0.8.0-cut.md fragment was left in place, duplicating the 0.8.0 section under Unreleased. Delete the stale fragment and re-render so render_changelog_fragments.py --check passes again, unblocking the v0.8.0 release-tag dispatch. * fix(release): re-render CHANGELOG fragments for v0.8.0 aggregate Ten more PRs merged into main after this fix's original commit, each adding docs/changelog.d fragments without re-running the renderer. Re-run render_changelog_fragments.py --update so --check passes again and the v0.8.0 release-tag dispatch is no longer blocked. * fix(dif): harden logistic and purification controls (#959) * test(dif): expose executable logistic controls * fix(dif): normalize controls before data and Rust * fix(dif): install hardened public DIF adapters * docs(changelog): record logistic DIF control boundary * docs(doctoring): trace logistic DIF control trust boundary * fix(changelog): classify logistic DIF control boundary --------- Co-authored-by: Claude * test(dif): preserve native zero-iteration and NumPy bool domains * fix(dif): preserve native control compatibility * test(dif): align regressions with native control domains * test(dif): prove hardened package aliases on landing lane * fix(dif): reject max_iter=0 before native discovery in logistic DIF wrappers safe_logistic_dif and safe_logistic_dif_purified normalized max_iter with minimum=0, admitting zero at the Python boundary even though the native logistic_sweep (crates/mlsirm-core/src/dif.rs) rejects max_iter == 0 with "max_iter must be >= 1", and the module's own doctoring already documented max_iter as a positive usize control. A zero cap therefore cleared control validation, materialized caller response/group arrays, and dispatched to the compiled core before failing there instead of being rejected up front (Devin AI review finding on PR #998). Raise the minimum to 1 for both wrappers so the rejection happens before data materialization or core discovery, matching the documented and native domains. Updated the compatibility regressions accordingly: the former zero-iteration dispatch proofs now use the true minimum (1), and new regressions prove max_iter=0 is rejected before core discovery for both logistic_dif and logistic_dif_purified. * fix(release): resync changelog after current main * fix(dif): make control installer idempotent * feat(interaction-map): own residual factorization in Rust * test(interaction-map): seal evidence and resource admission * fix(interaction-map): seal evidence before Rust * test(interaction-map): bound native resource requests * fix(interaction-map): bound native workspaces * fix(interaction-map): reject oversized axis controls before data * test(interaction-map): reject maximal axis request before data * docs(interaction-map): record bounded evidence admission * test(interaction-map): require shape-consistent empty map * test(interaction-map): require coherent empty rectangle * fix(interaction-map): normalize empty complete-case rectangle * docs(interaction-map): record empty complete-case contract * feat(interaction-map): expose coverage and cell evidence * test(interaction-map): reject NaN model expectations * test(interaction-map): reject nonfinite expected evidence in Rust * fix(interaction-map): reject NaN model expectations * fix(interaction-map): require finite expected evidence * docs(changelog): record finite interaction-map expectations * test(interaction-map): preserve observed infinity rejection in Rust * fix(interaction-map): reject infinite observed evidence in Rust * docs(changelog): align interaction-map missingness parity * chore: normalize composed public exports --------- Co-authored-by: Claude Co-authored-by: opencode-agent[bot] <219766164+opencode-agent[bot]@users.noreply.github.com> --- CHANGELOG.md | 14 + crates/fast-mlsirm-py/src/lib.rs | 41 ++ crates/mlsirm-core/src/entrypoint.rs | 1 + crates/mlsirm-core/src/factor.rs | 5 +- crates/mlsirm-core/src/interaction_map.rs | 353 ++++++++++++++++++ .../0021-residual-interaction-map-boundary.md | 49 +++ .../1388-interaction-map-admission.md | 9 + .../1389-interaction-map-empty-result.md | 8 + ...390-interaction-map-expected-finiteness.md | 8 + .../958-logistic-dif-control-boundary.md | 7 + docs/doctoring/logistic-dif-control-safety.md | 32 ++ python/fast_mlsirm/__init__.py | 131 +++++-- python/fast_mlsirm/_dif_control_safety.py | 190 ++++++++++ python/fast_mlsirm/interaction_map.py | 280 ++++++++++++++ tests/test_interaction_map.py | 284 ++++++++++++++ ...st_logistic_dif_control_callback_safety.py | 229 ++++++++++++ ...ogistic_dif_native_domain_compatibility.py | 148 ++++++++ 17 files changed, 1759 insertions(+), 30 deletions(-) create mode 100644 crates/mlsirm-core/src/interaction_map.rs create mode 100644 docs/adr/0021-residual-interaction-map-boundary.md create mode 100644 docs/changelog.d/1388-interaction-map-admission.md create mode 100644 docs/changelog.d/1389-interaction-map-empty-result.md create mode 100644 docs/changelog.d/1390-interaction-map-expected-finiteness.md create mode 100644 docs/changelog.d/958-logistic-dif-control-boundary.md create mode 100644 docs/doctoring/logistic-dif-control-safety.md create mode 100644 python/fast_mlsirm/_dif_control_safety.py create mode 100644 python/fast_mlsirm/interaction_map.py create mode 100644 tests/test_interaction_map.py create mode 100644 tests/test_logistic_dif_control_callback_safety.py create mode 100644 tests/test_logistic_dif_native_domain_compatibility.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 0707c6f0f..7ad56e890 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,6 +44,14 @@ ### Changed +#### Own residual interaction-map computation in the psychometric core + +- Add a Rust-backed, complete-case Gabriel residual interaction-map contract + for consumers that already hold observed responses and fitted IRT + expectations. The API returns coordinates, singular values, axis inertia, + reconstruction, unexplained residual, and the exact cross term without + product identifiers, persistence, authorization, or presentation policy. + #### Validate Rasch CML controls before data materialization - Validate `max_iter` and `tol` before caller-owned response or group arrays are materialized by the public Rasch CML and Andersen LR entry points. @@ -60,6 +68,12 @@ - Preserve CAT's independent adaptive-efficiency gate on mean administered items, so uncertainty calibration and error recovery cannot mask a fallback to non-adaptive item selection. - Keep all production likelihood, marginal-ML/EM, EAP/CAT scoring, item-information/selection, stopping, and uncertainty arithmetic Rust-owned; the added Python calculations are explicit true-parameter recovery-test summaries only. +#### Harden observed-score logistic DIF controls + +- Validate logistic and purified observed-score DIF semantic controls before caller-owned response/group materialization and before compiled Rust-core discovery. +- Reject caller-defined scalar subclasses, arbitrary conversion providers, booleans-as-numbers, invalid FDR levels, zero iteration caps, negative anchor floors, and values outside native `usize` without invoking caller callbacks. +- Preserve genuine supported NumPy scalar compatibility and keep all logistic/Mantel-Haenszel/purification statistics and BH arithmetic Rust-owned. + #### Method literature and citation ADRs - Record primary-paper citations and ADRs for the shipped Angoff delta-plot diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 43c7de75d..f29fbfd00 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -35,6 +35,7 @@ use mlsirm_core::inference::{ standard_errors_from_vcov as core_standard_errors_from_vcov, vcov_from_hessian as core_vcov_from_hessian, }; +use mlsirm_core::interaction_map::residual_interaction_map as core_residual_interaction_map; use mlsirm_core::jmle_opt::{adam as core_jmle_adam, lbfgs as core_jmle_lbfgs, run_optimizer as core_jmle_run_optimizer}; use mlsirm_core::marginal::{ fit_marginal_full as core_fit_marginal_full, Anchors, ItemCovariate, MarginalConfig, @@ -9395,11 +9396,51 @@ fn standard_errors_from_vcov(vcov: PyReadonlyArray2<'_, f64>) -> PyResult, + observed: PyReadonlyArray2<'_, f64>, + expected: PyReadonlyArray2<'_, f64>, + axis_count: usize, +) -> PyResult> { + if observed.shape() != expected.shape() { + return Err(PyValueError::new_err( + "observed and expected must have the same two-dimensional shape", + )); + } + let shape = observed.shape(); + let result = core_residual_interaction_map( + observed.as_slice()?, + expected.as_slice()?, + shape[0], + shape[1], + axis_count, + ) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("person_indices", result.person_indices)?; + out.set_item("item_indices", result.item_indices)?; + out.set_item("scored_person_count", result.scored_person_count)?; + out.set_item("scored_item_count", result.scored_item_count)?; + out.set_item("person_coordinates", result.person_coordinates)?; + out.set_item("item_coordinates", result.item_coordinates)?; + out.set_item("singular_values", result.singular_values)?; + out.set_item("axis_shares", result.axis_shares)?; + out.set_item("residual", result.residual)?; + out.set_item("distance", result.distance)?; + out.set_item("reconstruction", result.reconstruction)?; + out.set_item("unexplained", result.unexplained)?; + out.set_item("cross_share", result.cross_share)?; + out.set_item("axis_count", result.axis_count)?; + Ok(out.into_any().unbind()) +} + #[pymodule] #[pyo3(name = "_core")] fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add("MARGINAL_CAPABILITY_VERSION", MARGINAL_CAPABILITY_VERSION)?; m.add_function(wrap_pyfunction!(neg_loglik_and_grad, m)?)?; + m.add_function(wrap_pyfunction!(residual_interaction_map, m)?)?; m.add_function(wrap_pyfunction!(fit_mmle_2pl, m)?)?; m.add_function(wrap_pyfunction!(fit_cdm, m)?)?; m.add_function(wrap_pyfunction!(fit_gdina, m)?)?; diff --git a/crates/mlsirm-core/src/entrypoint.rs b/crates/mlsirm-core/src/entrypoint.rs index f3bdf2cdc..2587b04ed 100644 --- a/crates/mlsirm-core/src/entrypoint.rs +++ b/crates/mlsirm-core/src/entrypoint.rs @@ -6,4 +6,5 @@ include!("lib.rs"); +pub mod interaction_map; pub mod rotation; diff --git a/crates/mlsirm-core/src/factor.rs b/crates/mlsirm-core/src/factor.rs index 2d24352cf..2a6bfb672 100644 --- a/crates/mlsirm-core/src/factor.rs +++ b/crates/mlsirm-core/src/factor.rs @@ -106,7 +106,10 @@ pub struct OmegaResult { /// with matching eigenvectors (columns of `v`, row-major `p x p`). Same /// sweep scheme as `parallel::symmetric_eigenvalues_desc`, extended with /// vector accumulation. -fn symmetric_eigen_desc(matrix: &[f64], p: usize) -> Result<(Vec, Vec), String> { +pub(crate) fn symmetric_eigen_desc( + matrix: &[f64], + p: usize, +) -> Result<(Vec, Vec), String> { let mut a = matrix.to_vec(); let mut v = vec![0.0; p * p]; for i in 0..p { diff --git a/crates/mlsirm-core/src/interaction_map.rs b/crates/mlsirm-core/src/interaction_map.rs new file mode 100644 index 000000000..2ce690402 --- /dev/null +++ b/crates/mlsirm-core/src/interaction_map.rs @@ -0,0 +1,353 @@ +//! Gabriel residual interaction maps after an externally fitted IRT main effect. +//! +//! The caller supplies observed responses and their model expectations. This +//! module owns the numerical residual factorization; product repositories own +//! identifiers, authorization, persistence, and presentation. + +use crate::factor::symmetric_eigen_desc; + +const SINGULAR_FLOOR: f64 = 1e-12; +const MAX_INTERACTION_MAP_CELLS: usize = 20_000_000; +const MAX_INTERACTION_MAP_COORDINATE_CELLS: usize = 20_000_000; +const MAX_INTERACTION_MAP_EIGEN_WORKSPACE_BYTES: usize = 128 * 1024 * 1024; +const EIGEN_WORKSPACE_MATRIX_COUNT: usize = 3; + +fn validate_factorization_workspace( + rows: usize, + columns: usize, + axis_count: usize, +) -> Result<(), String> { + let coordinate_cells = rows + .checked_add(columns) + .and_then(|dimension_sum| dimension_sum.checked_mul(axis_count)) + .ok_or_else(|| "residual interaction map coordinate request overflows".to_string())?; + if coordinate_cells > MAX_INTERACTION_MAP_COORDINATE_CELLS { + return Err(format!( + "residual interaction map coordinate request exceeds {MAX_INTERACTION_MAP_COORDINATE_CELLS} cells" + )); + } + + // `symmetric_eigen_desc` holds the caller's Gram matrix plus an internal + // matrix copy and eigenvector matrix at peak, so budget all three dense + // `columns x columns` f64 matrices before allocating any of them. + let eigen_workspace_bytes = columns + .checked_mul(columns) + .and_then(|cells| cells.checked_mul(std::mem::size_of::())) + .and_then(|bytes| bytes.checked_mul(EIGEN_WORKSPACE_MATRIX_COUNT)) + .ok_or_else(|| "residual interaction map eigen workspace overflows".to_string())?; + if eigen_workspace_bytes > MAX_INTERACTION_MAP_EIGEN_WORKSPACE_BYTES { + return Err(format!( + "residual interaction map eigen workspace exceeds {} MiB", + MAX_INTERACTION_MAP_EIGEN_WORKSPACE_BYTES / (1024 * 1024) + )); + } + Ok(()) +} + +/// A complete-case Gabriel biplot of `observed - expected`. +#[derive(Debug, Clone, PartialEq)] +pub struct ResidualInteractionMap { + pub person_indices: Vec, + pub item_indices: Vec, + pub scored_person_count: usize, + pub scored_item_count: usize, + pub person_coordinates: Vec, + pub item_coordinates: Vec, + pub singular_values: Vec, + pub axis_shares: Vec, + pub residual: Vec, + pub distance: Vec, + pub reconstruction: Vec, + pub unexplained: Vec, + pub cross_share: Vec>, + pub axis_count: usize, +} + +/// Factor a complete-case residual rectangle using Gabriel's symmetric scaling. +pub fn residual_interaction_map( + observed: &[f64], + expected: &[f64], + n_persons: usize, + n_items: usize, + axis_count: usize, +) -> Result { + if axis_count == 0 { + return Err("axis_count must be positive".into()); + } + // Even an empty complete-case rectangle returns `axis_shares`, so bound the + // requested axis vector independently of the surviving matrix dimensions. + if axis_count > MAX_INTERACTION_MAP_COORDINATE_CELLS { + return Err(format!( + "residual interaction map coordinate request exceeds {MAX_INTERACTION_MAP_COORDINATE_CELLS} cells" + )); + } + + let cell_count = n_persons + .checked_mul(n_items) + .ok_or_else(|| "residual interaction map shape overflows".to_string())?; + if cell_count > MAX_INTERACTION_MAP_CELLS { + return Err(format!( + "residual interaction map logical-cell count exceeds {MAX_INTERACTION_MAP_CELLS}" + )); + } + if observed.len() != cell_count || expected.len() != cell_count { + return Err("observed and expected lengths must match the declared shape".into()); + } + if observed.iter().any(|value| value.is_infinite()) { + return Err("observed values must not be infinite".into()); + } + if expected.iter().any(|value| !value.is_finite()) { + return Err("expected values must be finite".into()); + } + + let observed_cell = |index: usize| observed[index].is_finite(); + let mut person_indices: Vec = (0..n_persons) + .filter(|&person| (0..n_items).any(|item| observed_cell(person * n_items + item))) + .collect(); + let mut item_indices: Vec = (0..n_items) + .filter(|&item| (0..n_persons).any(|person| observed_cell(person * n_items + item))) + .collect(); + let scored_person_count = person_indices.len(); + let scored_item_count = item_indices.len(); + person_indices.retain(|&person| { + item_indices + .iter() + .all(|&item| observed_cell(person * n_items + item)) + }); + item_indices.retain(|&item| { + person_indices + .iter() + .all(|&person| observed_cell(person * n_items + item)) + }); + if person_indices.is_empty() || item_indices.is_empty() { + // An empty complete-case rectangle has no row or column coordinate + // space. Normalize both axes to the same empty rectangle so every + // returned payload remains shape-consistent at the Python boundary. + person_indices.clear(); + item_indices.clear(); + return Ok(ResidualInteractionMap { + person_indices, + item_indices, + scored_person_count, + scored_item_count, + person_coordinates: Vec::new(), + item_coordinates: Vec::new(), + singular_values: Vec::new(), + axis_shares: vec![0.0; axis_count], + residual: Vec::new(), + distance: Vec::new(), + reconstruction: Vec::new(), + unexplained: Vec::new(), + cross_share: Vec::new(), + axis_count, + }); + } + + let rows = person_indices.len(); + let columns = item_indices.len(); + validate_factorization_workspace(rows, columns, axis_count)?; + + let mut residual = Vec::with_capacity(rows * columns); + for &person in &person_indices { + for &item in &item_indices { + let index = person * n_items + item; + residual.push(observed[index] - expected[index]); + } + } + let center = residual.iter().sum::() / residual.len() as f64; + let centered: Vec = residual.iter().map(|value| value - center).collect(); + let mut gram = vec![0.0; columns * columns]; + for left in 0..columns { + for right in left..columns { + let value = (0..rows) + .map(|row| centered[row * columns + left] * centered[row * columns + right]) + .sum(); + gram[left * columns + right] = value; + gram[right * columns + left] = value; + } + } + let (eigenvalues, eigenvectors) = symmetric_eigen_desc(&gram, columns)?; + let singular_values: Vec = eigenvalues + .iter() + .take_while(|value| **value > SINGULAR_FLOOR * SINGULAR_FLOOR) + .map(|value| value.sqrt()) + .collect(); + let retained = singular_values.len(); + let mut person_coordinates = vec![0.0; rows * axis_count]; + let mut item_coordinates = vec![0.0; columns * axis_count]; + for axis in 0..axis_count.min(retained) { + let root_singular = singular_values[axis].sqrt(); + for item in 0..columns { + item_coordinates[item * axis_count + axis] = + eigenvectors[item * columns + axis] * root_singular; + } + for person in 0..rows { + let projection: f64 = (0..columns) + .map(|item| centered[person * columns + item] * eigenvectors[item * columns + axis]) + .sum(); + person_coordinates[person * axis_count + axis] = projection / root_singular; + } + } + let inertia = singular_values + .iter() + .map(|value| value * value) + .sum::(); + let axis_shares = (0..axis_count) + .map(|axis| { + singular_values + .get(axis) + .map_or(0.0, |value| value * value / inertia) + }) + .collect(); + let mut reconstruction = Vec::with_capacity(rows * columns); + let mut distance = Vec::with_capacity(rows * columns); + let mut unexplained = Vec::with_capacity(rows * columns); + let mut cross_share = Vec::with_capacity(rows * columns); + for person in 0..rows { + for item in 0..columns { + distance.push( + (0..axis_count) + .map(|axis| { + let difference = person_coordinates[person * axis_count + axis] + - item_coordinates[item * axis_count + axis]; + difference * difference + }) + .sum::() + .sqrt(), + ); + let fitted = (0..axis_count) + .map(|axis| { + person_coordinates[person * axis_count + axis] + * item_coordinates[item * axis_count + axis] + }) + .sum::(); + let raw = residual[person * columns + item]; + let remainder = raw - fitted; + let share = if raw.abs() > SINGULAR_FLOOR { + Some(2.0 * fitted * remainder / (raw * raw)) + } else if fitted.abs() <= SINGULAR_FLOOR && remainder.abs() <= SINGULAR_FLOOR { + Some(0.0) + } else { + None + }; + reconstruction.push(fitted); + unexplained.push(remainder); + cross_share.push(share.filter(|value| value.is_finite())); + } + } + + Ok(ResidualInteractionMap { + person_indices, + item_indices, + scored_person_count, + scored_item_count, + person_coordinates, + item_coordinates, + singular_values, + axis_shares, + residual, + distance, + reconstruction, + unexplained, + cross_share, + axis_count, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn reconstructs_rank_one_residual_without_inventing_a_second_axis() { + let observed = [2.0, 0.0, 0.0, 2.0]; + let expected = [1.0, 1.0, 1.0, 1.0]; + let map = residual_interaction_map(&observed, &expected, 2, 2, 2).unwrap(); + assert_eq!(map.person_indices, vec![0, 1]); + assert_eq!(map.item_indices, vec![0, 1]); + assert_eq!(map.singular_values.len(), 1); + assert!((map.axis_shares[0] - 1.0).abs() < 1e-12); + assert_eq!(map.axis_shares[1], 0.0); + for (raw, fitted) in observed + .iter() + .zip(expected) + .map(|(observed, expected)| observed - expected) + .zip(&map.reconstruction) + { + assert!((raw - fitted).abs() < 1e-12); + } + assert!(map + .distance + .iter() + .all(|value| value.is_finite() && *value >= 0.0)); + } + + #[test] + fn excludes_incomplete_rows_without_zero_filling() { + let map = + residual_interaction_map(&[2.0, f64::NAN, 1.0, 2.0], &[1.0, 1.0, 1.0, 1.0], 2, 2, 2) + .unwrap(); + assert_eq!(map.person_indices, vec![1]); + assert_eq!(map.item_indices, vec![0, 1]); + assert_eq!(map.scored_person_count, 2); + assert_eq!(map.scored_item_count, 2); + } + + #[test] + fn empty_complete_case_rectangle_has_empty_indices_on_both_axes() { + let map = residual_interaction_map( + &[1.0, f64::NAN, f64::NAN, 1.0], + &[0.0, 0.0, 0.0, 0.0], + 2, + 2, + 2, + ) + .unwrap(); + assert!(map.person_indices.is_empty()); + assert!(map.item_indices.is_empty()); + assert!(map.person_coordinates.is_empty()); + assert!(map.item_coordinates.is_empty()); + assert!(map.reconstruction.is_empty()); + assert!(map.unexplained.is_empty()); + assert!(map.cross_share.is_empty()); + assert_eq!(map.axis_shares, vec![0.0, 0.0]); + } + + #[test] + fn rejects_nonfinite_expected_before_complete_case_filtering() { + let error = residual_interaction_map(&[1.0], &[f64::NAN], 1, 1, 1).unwrap_err(); + assert!(error.contains("expected")); + } + + #[test] + fn rejects_infinite_observed_instead_of_treating_it_as_missing() { + let error = residual_interaction_map(&[f64::INFINITY], &[0.0], 1, 1, 1).unwrap_err(); + assert!(error.contains("observed")); + } + + #[test] + fn rejects_declared_grid_above_resource_ceiling_before_length_validation() { + let error = residual_interaction_map(&[], &[], 1, 20_000_001, 1).unwrap_err(); + assert!(error.contains("logical-cell")); + } + + #[test] + fn rejects_coordinate_overflow_before_allocation() { + let call = std::panic::catch_unwind(|| { + residual_interaction_map(&[1.0, 1.0], &[0.0, 0.0], 2, 1, usize::MAX) + }); + assert!( + call.is_ok(), + "resource admission must return Err instead of panicking" + ); + let error = call.unwrap().unwrap_err(); + assert!(error.contains("coordinate")); + } + + #[test] + fn rejects_quadratic_eigen_workspace_before_allocation() { + let error = validate_factorization_workspace(1, 3_000, 1).unwrap_err(); + assert!(error.contains("eigen workspace")); + validate_factorization_workspace(1, 2_000, 1).unwrap(); + } +} diff --git a/docs/adr/0021-residual-interaction-map-boundary.md b/docs/adr/0021-residual-interaction-map-boundary.md new file mode 100644 index 000000000..0abb8ae7d --- /dev/null +++ b/docs/adr/0021-residual-interaction-map-boundary.md @@ -0,0 +1,49 @@ +# ADR-0021: Rust-owned residual interaction-map boundary + +Status: **Accepted** +Date: 2026-08-25 +Supersedes: none +Superseded by: none + +## Context + +Downstream measurement products need to inspect residual person-item +interactions after a fitted IRT main effect. Jeon et al. (2021) place those +interactions in a shared latent space; Gabriel (1971) defines the symmetric +biplot factorization used to display a residual matrix. Computing that map in +a product repository duplicates psychometric numerical policy and prevents +independent consumers from sharing one tested contract. + +## Decision + +fast-mlsirm owns `residual_interaction_map(observed, expected, axis_count)`. +The Rust core computes `R = Y - E`, excludes incomplete rows and columns +without zero filling, centers the admitted residual rectangle, performs the +Gabriel symmetric factorization, and returns coordinates, singular values, +axis inertia, truncated reconstruction `Rhat`, `U = R - Rhat`, and the exact +algebraic cross term `2 Rhat U / R^2`. The number of retained axes is required +from the caller; the library does not invent a display dimension. + +The contract contains array indices only. Products retain responsibility for +domain identifiers, authorization, persistence, closest/farthest selection, +and UI wording. A consumer must not reproduce the factorization or cross-term +arithmetic locally. + +## Consequences + +- Numerical behavior is reusable, Rust-owned, and independently testable. +- Missing cells remain absent rather than becoming evidence-valued zeros. +- The cross term is an auditable algebraic identity, not a fitted weight, + threshold, quality score, or heuristic. +- This API does not fit an LSIRM or replace the caller's GRM/GPCM fit. + +## References + +Gabriel, K. R. (1971). The biplot graphic display of matrices with application +to principal component analysis. *Biometrika, 58*(3), 453–467. +https://doi.org/10.1093/biomet/58.3.453 + +Jeon, M., Jin, I. H., Schweinberger, M., & Baugh, S. (2021). Mapping +unobserved item-respondent interactions: A latent space item response model +with interaction map. *Psychometrika, 86*(2), 378–403. +https://doi.org/10.1007/s11336-021-09762-5 diff --git a/docs/changelog.d/1388-interaction-map-admission.md b/docs/changelog.d/1388-interaction-map-admission.md new file mode 100644 index 000000000..a67bf0231 --- /dev/null +++ b/docs/changelog.d/1388-interaction-map-admission.md @@ -0,0 +1,9 @@ +# Interaction-map evidence admission + +## Fixed + +- Seal `axis_count` before caller matrix work and reject callback-bearing integer identities, booleans, nonpositive controls, and requests above the interaction-map coordinate envelope before scientific evidence is inspected. +- Admit only exact real-numeric NumPy arrays or exact built-in two-dimensional numeric sequences before NumPy materialization, reject complex/non-real storage and infinities, preserve `NaN` as missingness, and require wider integer/floating evidence to survive the Rust `f64` boundary without changing identity. +- Bound public and native interaction-map logical cells and coordinate requests at 20,000,000 cells, and bound the Rust symmetric-eigendecomposition workspace at 128 MiB before dense Gram/eigenvector allocation. + +Gabriel symmetric factorization, singular values, coordinates, reconstruction, unexplained residuals, cross-term arithmetic, and other production numerical behavior remain Rust-owned and unchanged. diff --git a/docs/changelog.d/1389-interaction-map-empty-result.md b/docs/changelog.d/1389-interaction-map-empty-result.md new file mode 100644 index 000000000..e4a4c52d8 --- /dev/null +++ b/docs/changelog.d/1389-interaction-map-empty-result.md @@ -0,0 +1,8 @@ +# Interaction-map empty complete-case result + +## Fixed + +- Normalize an empty complete-case interaction rectangle to empty respondent and item index sets so the Rust result remains shape-consistent with zero-length person/item coordinates and zero-by-zero reconstruction, unexplained-residual, and cross-share arrays at the Python boundary. +- Preserve the requested bounded `axis_shares` length for an empty map without inventing a maximal-complete-submatrix selection rule or retaining one non-empty axis after the other has collapsed. + +Gabriel factorization, singular values, coordinate arithmetic, reconstruction, unexplained residuals, and cross-term calculations for non-empty complete-case rectangles remain unchanged and Rust-owned. diff --git a/docs/changelog.d/1390-interaction-map-expected-finiteness.md b/docs/changelog.d/1390-interaction-map-expected-finiteness.md new file mode 100644 index 000000000..8af55e9c5 --- /dev/null +++ b/docs/changelog.d/1390-interaction-map-expected-finiteness.md @@ -0,0 +1,8 @@ +# Interaction-map expected-evidence finiteness + +## Fixed + +- Preserve `NaN` exclusively as an observed-response missingness marker while rejecting observed infinity and rejecting both `NaN` and infinity in fitted model expectations before complete-case filtering or Rust factorization. +- Replay the same missingness/finiteness contract in the Rust core so direct PyO3/core callers cannot silently turn infinite observed evidence or invalid model expectations into missing cells that change the analyzed interaction rectangle. + +Gabriel factorization, singular values, coordinates, reconstruction, unexplained residuals, distance, cross-term arithmetic, and observed-response `NaN` missingness remain unchanged and Rust-owned. diff --git a/docs/changelog.d/958-logistic-dif-control-boundary.md b/docs/changelog.d/958-logistic-dif-control-boundary.md new file mode 100644 index 000000000..31cd1120c --- /dev/null +++ b/docs/changelog.d/958-logistic-dif-control-boundary.md @@ -0,0 +1,7 @@ +# Harden observed-score logistic DIF controls + +## Changed + +- Validate logistic and purified observed-score DIF semantic controls before caller-owned response/group materialization and before compiled Rust-core discovery. +- Reject caller-defined scalar subclasses, arbitrary conversion providers, booleans-as-numbers, invalid FDR levels, zero iteration caps, negative anchor floors, and values outside native `usize` without invoking caller callbacks. +- Preserve genuine supported NumPy scalar compatibility and keep all logistic/Mantel-Haenszel/purification statistics and BH arithmetic Rust-owned. diff --git a/docs/doctoring/logistic-dif-control-safety.md b/docs/doctoring/logistic-dif-control-safety.md new file mode 100644 index 000000000..e7709dacc --- /dev/null +++ b/docs/doctoring/logistic-dif-control-safety.md @@ -0,0 +1,32 @@ +# Observed-score logistic DIF control trust boundary + +## Decision + +The public observed-score logistic DIF and purification adapters validate and normalize semantic scalar controls before any caller-owned response/group materialization and before compiled-core discovery. Only exact built-in scalar identities and explicitly supported concrete NumPy scalar identities are admitted; caller-defined subclasses and arbitrary conversion-protocol providers are rejected without executing their conversion callbacks. + +This is a trust-boundary change, not a statistical-method change. Logistic-regression likelihood/IRLS computation, Mantel-Haenszel statistics, purification sweeps, Benjamini-Hochberg decisions, effect sizes, and returned result structure remain in the Rust implementation. + +## Native-domain alignment + +The Python boundary mirrors the existing native contract rather than inventing new estimator semantics: + +- `exclude_studied_item` is an exact Boolean control; +- `fdr_q` is finite and in `(0, 1]`; +- `max_iter` and `max_rounds` are positive Rust `usize` controls; +- `min_anchor_items` is a nonnegative Rust `usize`; zero remains admissible because the native purification guard applies a minimum effective anchor count of one. + +Values outside the native `usize` range fail before PyO3 conversion. Genuine NumPy integer and floating scalar identities remain accepted for numeric controls so previously supported numerical caller code does not lose interoperability. + +## Threat model and verification + +The defect class is improper validation of the specified input type: invoking `bool(...)`, `float(...)`, or `int(...)` on an untrusted object can execute caller-owned Python code before the library establishes the intended control type/range. MITRE CWE-1287 recommends an accept-known-good strategy that validates type and relevant value properties before use. The regression suite therefore uses hostile subclasses/protocol providers whose conversion hooks increment counters or raise, plus data/core sentinels that prove rejected controls terminate before either caller data or native discovery is touched. + +The implementation intentionally does not duplicate logistic, MH, purification, or FDR formulas in Python. Statistical parity remains the responsibility of the existing Rust-backed tests; this slice verifies admission, ordering, exact normalized payload types, and native-domain fidelity. + +## References + +French, A. W., & Maller, S. J. (2007). Iterative purification and effect size use with logistic regression for differential item functioning detection. *Educational and Psychological Measurement, 67*(3), 373–393. https://doi.org/10.1177/0013164406294781 + +MITRE. (2026). *CWE-1287: Improper validation of specified type of input* (CWE Version 4.20). Common Weakness Enumeration. https://cwe.mitre.org/data/definitions/1287.html + +Swaminathan, H., & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures. *Journal of Educational Measurement, 27*(4), 361–370. https://doi.org/10.1111/j.1745-3984.1990.tb00754.x diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index 88c4f582a..791837d78 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -8,6 +8,7 @@ from . import _legacy_init as _legacy_init from . import cat as _cat from . import cdm as _cdm +from . import dif as _dif from . import exposure as _exposure from . import fitstats as _fitstats from . import inference as _inference @@ -16,13 +17,16 @@ from . import scaling as _scaling from . import serving as _serving from . import validation as _validation -from ._cat_administration_resource_safety import install as _install_cat_administration_resource_safety +from ._cat_administration_resource_safety import ( + install as _install_cat_administration_resource_safety, +) from ._cdm_response_safety import install as _install_cdm_response_safety +from ._dif_control_safety import install as _install_dif_control_safety from ._exposure_array_safety import install as _install_exposure_array_safety from ._exposure_flexilevel_safety import install as _install_exposure_flexilevel_safety +from ._fit_public import fit as _public_fit from ._fitstats_control_safety import install as _install_fitstats_control_safety from ._fleiss_control_safety import install as _install_fleiss_control_safety -from ._fit_public import fit as _public_fit from ._icc_control_safety import install as _install_icc_control_safety from ._inference_admission_safety import install as _install_inference_admission_safety from ._polytomous_prediction_admission import ( @@ -30,6 +34,8 @@ ) from ._scaling_control_safety import install as _install_scaling_control_safety from ._serving_export_safety import install as _install_serving_export_safety +from .interaction_map import ResidualInteractionMap as ResidualInteractionMap +from .interaction_map import residual_interaction_map as residual_interaction_map # Harden historical public adapters before copying legacy exports. These # wrappers validate and normalize semantic controls/evidence only; result @@ -40,6 +46,16 @@ _install_fitstats_control_safety(_fitstats) _install_cdm_response_safety(_cdm) _install_icc_control_safety(_reliability) + +_install_dif_control_safety(_dif) +for _dif_name in ( + "logistic_dif", + "mantel_haenszel_dif_purified", + "logistic_dif_purified", +): + if hasattr(_legacy_init, _dif_name): + setattr(_legacy_init, _dif_name, getattr(_dif, _dif_name)) + _install_inference_admission_safety(_inference) _install_scaling_control_safety(_scaling) _install_fleiss_control_safety(_validation) @@ -64,6 +80,7 @@ _exposure, _inference, _install_cat_administration_resource_safety, + _install_dif_control_safety, _install_exposure_array_safety, _install_exposure_flexilevel_safety, _fitstats, @@ -77,6 +94,8 @@ _install_scaling_control_safety, _install_serving_export_safety, _polytomous, + _dif, + _dif_name, _reliability, _scaling, _serving, @@ -97,48 +116,102 @@ from .bifactor_scoreability import ( BifactorScoreabilityResult as BifactorScoreabilityResult, +) +from .bifactor_scoreability import ( bifactor_scoreability as bifactor_scoreability, +) +from .bifactor_scoreability import ( bifactor_scoreability_from_logit_slopes as bifactor_scoreability_from_logit_slopes, ) -from .rating_range import ( - RatingRangeEvidence as RatingRangeEvidence, - paired_rating_range_evidence as paired_rating_range_evidence, +from .irt_contract import ( + MIN_FACTOR_ANCHOR_ITEMS as MIN_FACTOR_ANCHOR_ITEMS, ) -from .rotation import ( - RotationCriterionInfo as RotationCriterionInfo, - RotationSolution as RotationSolution, - available_rotation_criteria as available_rotation_criteria, - rotate_factor_loadings as rotate_factor_loadings, - rotation_criterion_value_gradient as rotation_criterion_value_gradient, +from .irt_contract import ( + MIN_IRT_ITEMS as MIN_IRT_ITEMS, ) -from .llm_judge import ( - CONTEXTUAL_ORCHESTRATOR_CONTRACT_V1 as CONTEXTUAL_ORCHESTRATOR_CONTRACT_V1, - ContextualOrchestratorJudge as ContextualOrchestratorJudge, - JudgeCriterion as JudgeCriterion, - JudgeFormatError as JudgeFormatError, - LLMJudgeResult as LLMJudgeResult, - MAX_BINARY_THRESHOLD_CALLS as MAX_BINARY_THRESHOLD_CALLS, - MAX_JUDGE_CATEGORIES as MAX_JUDGE_CATEGORIES, +from .irt_contract import ( + MIN_IRT_PERSONS as MIN_IRT_PERSONS, +) +from .irt_contract import ( + MIN_ITEM_DISTINCT_VALUES as MIN_ITEM_DISTINCT_VALUES, +) +from .irt_contract import ( + MIN_OBSERVED_PER_ITEM as MIN_OBSERVED_PER_ITEM, +) +from .irt_contract import ( + IRTItemType as IRTItemType, +) +from .irt_contract import ( + fit_irt_experiment as fit_irt_experiment, +) +from .irt_contract import ( + validate_irt_experiment_readiness as validate_irt_experiment_readiness, +) +from .irt_contract import ( + validate_irt_response_matrix as validate_irt_response_matrix, ) from .judge_calibration import ( CALIBRATION_VARIANTS as CALIBRATION_VARIANTS, +) +from .judge_calibration import ( CONTAMINATION_STATUSES as CONTAMINATION_STATUSES, +) +from .judge_calibration import ( JudgeCalibrationCase as JudgeCalibrationCase, +) +from .judge_calibration import ( JudgeCalibrationOutcome as JudgeCalibrationOutcome, +) +from .judge_calibration import ( JudgeCalibrationReport as JudgeCalibrationReport, +) +from .judge_calibration import ( build_multiple_choice_calibration_cases as build_multiple_choice_calibration_cases, +) +from .judge_calibration import ( evaluate_paired_calibration as evaluate_paired_calibration, ) -from .irt_contract import ( - IRTItemType as IRTItemType, - MIN_IRT_ITEMS as MIN_IRT_ITEMS, - MIN_IRT_PERSONS as MIN_IRT_PERSONS, - MIN_OBSERVED_PER_ITEM as MIN_OBSERVED_PER_ITEM, - MIN_ITEM_DISTINCT_VALUES as MIN_ITEM_DISTINCT_VALUES, - MIN_FACTOR_ANCHOR_ITEMS as MIN_FACTOR_ANCHOR_ITEMS, - fit_irt_experiment as fit_irt_experiment, - validate_irt_response_matrix as validate_irt_response_matrix, - validate_irt_experiment_readiness as validate_irt_experiment_readiness, +from .llm_judge import ( + CONTEXTUAL_ORCHESTRATOR_CONTRACT_V1 as CONTEXTUAL_ORCHESTRATOR_CONTRACT_V1, +) +from .llm_judge import ( + MAX_BINARY_THRESHOLD_CALLS as MAX_BINARY_THRESHOLD_CALLS, +) +from .llm_judge import ( + MAX_JUDGE_CATEGORIES as MAX_JUDGE_CATEGORIES, +) +from .llm_judge import ( + ContextualOrchestratorJudge as ContextualOrchestratorJudge, +) +from .llm_judge import ( + JudgeCriterion as JudgeCriterion, +) +from .llm_judge import ( + JudgeFormatError as JudgeFormatError, +) +from .llm_judge import ( + LLMJudgeResult as LLMJudgeResult, +) +from .rating_range import ( + RatingRangeEvidence as RatingRangeEvidence, +) +from .rating_range import ( + paired_rating_range_evidence as paired_rating_range_evidence, +) +from .rotation import ( + RotationCriterionInfo as RotationCriterionInfo, +) +from .rotation import ( + RotationSolution as RotationSolution, +) +from .rotation import ( + available_rotation_criteria as available_rotation_criteria, +) +from .rotation import ( + rotate_factor_loadings as rotate_factor_loadings, +) +from .rotation import ( + rotation_criterion_value_gradient as rotation_criterion_value_gradient, ) # Bind the rating-range API on the historical validation namespace without diff --git a/python/fast_mlsirm/_dif_control_safety.py b/python/fast_mlsirm/_dif_control_safety.py new file mode 100644 index 000000000..60cc9805b --- /dev/null +++ b/python/fast_mlsirm/_dif_control_safety.py @@ -0,0 +1,190 @@ +"""Callback-free semantic-control validation for observed-score DIF adapters. + +The logistic-regression, Mantel-Haenszel, and purification statistics remain +Rust-owned. This module only establishes trusted scalar identities and native +parameter ranges before legacy adapters can discover the compiled core or +materialize caller-owned response data. +""" + +from __future__ import annotations + +import builtins +import math +from functools import wraps +from types import ModuleType +from typing import Any, Callable + +import numpy as np + +_NUMPY_INTEGER_SCALAR_TYPES = ( + np.int8, + np.int16, + np.int32, + np.int64, + np.intp, + np.longlong, + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.uintp, + np.ulonglong, +) +_NUMPY_FLOATING_SCALAR_TYPES = ( + np.float16, + np.float32, + np.float64, + np.longdouble, +) +_USIZE_MAX = int(np.iinfo(np.uintp).max) +_DIF_HARDENED_ATTR = "__fast_mlsirm_dif_control_hardened__" + + +def _exact_bool(value: Any, name: str) -> bool: + """Normalize a concrete Python/NumPy bool without caller callbacks.""" + value_type = builtins.type(value) + if value_type is bool: + return value + if value_type is np.bool_: + return bool(value) + raise ValueError(f"{name} must be a bool") + + +def _real(value: Any, name: str) -> float: + """Normalize a trusted real scalar without arbitrary conversion hooks.""" + value_type = builtins.type(value) + try: + if value_type is int or value_type is float: + normalized = float(value) + elif any(value_type is scalar_type for scalar_type in _NUMPY_INTEGER_SCALAR_TYPES): + normalized = float(value) + elif any(value_type is scalar_type for scalar_type in _NUMPY_FLOATING_SCALAR_TYPES): + normalized = float(value) + else: + raise ValueError(f"{name} must be a real number") + except OverflowError as exc: + raise ValueError(f"{name} must be a real number") from exc + return normalized + + +def _usize(value: Any, name: str, *, minimum: int) -> int: + """Normalize one trusted integer into the native Rust ``usize`` domain.""" + value_type = builtins.type(value) + if value_type is int: + normalized = value + elif any(value_type is scalar_type for scalar_type in _NUMPY_INTEGER_SCALAR_TYPES): + normalized = int(value) + else: + raise ValueError(f"{name} must be an integer") + if normalized < minimum: + comparator = ">= 1" if minimum == 1 else ">= 0" + raise ValueError(f"{name} must be {comparator}") + if normalized > _USIZE_MAX: + raise ValueError(f"{name} exceeds the native usize range") + return normalized + + +def _common_controls(exclude_studied_item: Any, fdr_q: Any) -> tuple[bool, float]: + """Normalize matching-score and FDR controls before data/native access.""" + exclude = _exact_bool(exclude_studied_item, "exclude_studied_item") + q = _real(fdr_q, "fdr_q") + if not math.isfinite(q) or not 0.0 < q <= 1.0: + raise ValueError("fdr_q must be finite and in (0, 1]") + return exclude, q + + +def _original(function: Callable[..., Any]) -> Callable[..., Any]: + """Return the package-owned original behind a hardened wrapper.""" + if getattr(function, _DIF_HARDENED_ATTR, False): + return function.__wrapped__ + return function + + +def install(dif_module: ModuleType) -> None: + """Install callback-free wrappers on the observed-score DIF module.""" + current = ( + dif_module.logistic_dif, + dif_module.mantel_haenszel_dif_purified, + dif_module.logistic_dif_purified, + ) + if all(getattr(function, _DIF_HARDENED_ATTR, False) for function in current): + return + + original_logistic = _original(current[0]) + original_mh_purified = _original(current[1]) + original_logistic_purified = _original(current[2]) + + @wraps(original_logistic) + def safe_logistic_dif( + responses: Any, + group: Any, + exclude_studied_item: bool = False, + fdr_q: float = 0.05, + max_iter: int = 50, + ) -> Any: + """Validate logistic DIF controls before data or native discovery.""" + exclude, q = _common_controls(exclude_studied_item, fdr_q) + iterations = _usize(max_iter, "max_iter", minimum=1) + return original_logistic( + responses, + group, + exclude_studied_item=exclude, + fdr_q=q, + max_iter=iterations, + ) + + @wraps(original_mh_purified) + def safe_mantel_haenszel_dif_purified( + responses: Any, + group: Any, + exclude_studied_item: bool = False, + fdr_q: float = 0.05, + max_rounds: int = 3, + min_anchor_items: int = 4, + ) -> Any: + """Validate purified MH controls before data or native discovery.""" + exclude, q = _common_controls(exclude_studied_item, fdr_q) + rounds = _usize(max_rounds, "max_rounds", minimum=1) + anchors = _usize(min_anchor_items, "min_anchor_items", minimum=0) + return original_mh_purified( + responses, + group, + exclude_studied_item=exclude, + fdr_q=q, + max_rounds=rounds, + min_anchor_items=anchors, + ) + + @wraps(original_logistic_purified) + def safe_logistic_dif_purified( + responses: Any, + group: Any, + exclude_studied_item: bool = False, + fdr_q: float = 0.05, + max_iter: int = 50, + max_rounds: int = 3, + min_anchor_items: int = 4, + ) -> Any: + """Validate purified logistic controls before data or native discovery.""" + exclude, q = _common_controls(exclude_studied_item, fdr_q) + iterations = _usize(max_iter, "max_iter", minimum=1) + rounds = _usize(max_rounds, "max_rounds", minimum=1) + anchors = _usize(min_anchor_items, "min_anchor_items", minimum=0) + return original_logistic_purified( + responses, + group, + exclude_studied_item=exclude, + fdr_q=q, + max_iter=iterations, + max_rounds=rounds, + min_anchor_items=anchors, + ) + + wrappers = ( + safe_logistic_dif, + safe_mantel_haenszel_dif_purified, + safe_logistic_dif_purified, + ) + for function in wrappers: + setattr(function, _DIF_HARDENED_ATTR, True) + dif_module.logistic_dif, dif_module.mantel_haenszel_dif_purified, dif_module.logistic_dif_purified = wrappers diff --git a/python/fast_mlsirm/interaction_map.py b/python/fast_mlsirm/interaction_map.py new file mode 100644 index 000000000..4904ed4e9 --- /dev/null +++ b/python/fast_mlsirm/interaction_map.py @@ -0,0 +1,280 @@ +"""Rust-backed residual interaction maps for downstream measurement products.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np + +from . import _core + +_MAX_INTERACTION_MAP_CELLS = 20_000_000 +_MAX_INTERACTION_MAP_COORDINATE_CELLS = 20_000_000 +_TRUSTED_NUMPY_INTEGER_TYPES = ( + np.int8, + np.int16, + np.int32, + np.int64, + np.intp, + np.longlong, + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.uintp, + np.ulonglong, +) +_TRUSTED_NUMPY_FLOAT_TYPES = ( + np.float16, + np.float32, + np.float64, + np.longdouble, +) + + +@dataclass(frozen=True) +class ResidualInteractionMap: + """Complete-case Gabriel coordinates and auditable cell decomposition.""" + + person_indices: np.ndarray + item_indices: np.ndarray + scored_person_count: int + scored_item_count: int + person_coordinates: np.ndarray + item_coordinates: np.ndarray + singular_values: np.ndarray + axis_shares: np.ndarray + residual: np.ndarray + distance: np.ndarray + reconstruction: np.ndarray + unexplained: np.ndarray + cross_share: np.ndarray + + +def _axis_count(value: object) -> int: + """Normalize one inert positive integer axis request without caller coercion.""" + value_type = type(value) + if value_type is int: + parsed = value + elif any(value_type is trusted for trusted in _TRUSTED_NUMPY_INTEGER_TYPES): + parsed = int(value) + else: + raise TypeError("axis_count must be a positive integer") + if parsed <= 0: + raise ValueError("axis_count must be a positive integer") + if parsed > _MAX_INTERACTION_MAP_COORDINATE_CELLS: + raise ValueError( + "interaction map coordinate request exceeds " + f"{_MAX_INTERACTION_MAP_COORDINATE_CELLS} cells" + ) + return parsed + + +def _trusted_scalar_float64(name: str, value: object, *, allow_nan: bool) -> float: + """Return one callback-free scalar whose numeric identity survives binary64.""" + value_type = type(value) + if value_type is bool: + return float(value) + if value_type is int: + try: + parsed = float(value) + except OverflowError: + raise ValueError( + f"{name} values must be exactly representable as float64" + ) from None + if not math.isfinite(parsed) or int(parsed) != value: + raise ValueError(f"{name} values must be exactly representable as float64") + return parsed + if value_type is float: + if math.isinf(value): + raise ValueError(f"{name} must not contain infinite values") + if math.isnan(value) and not allow_nan: + raise ValueError(f"{name} must not contain NaN values") + return value + if value_type is np.bool_: + return float(value) + if any(value_type is trusted for trusted in _TRUSTED_NUMPY_INTEGER_TYPES): + integer = int(value) + try: + parsed = float(integer) + except OverflowError: + raise ValueError( + f"{name} values must be exactly representable as float64" + ) from None + if not math.isfinite(parsed) or int(parsed) != integer: + raise ValueError(f"{name} values must be exactly representable as float64") + return parsed + if any(value_type is trusted for trusted in _TRUSTED_NUMPY_FLOAT_TYPES): + if np.isinf(value): + raise ValueError(f"{name} must not contain infinite values") + if np.isnan(value): + if not allow_nan: + raise ValueError(f"{name} must not contain NaN values") + return float("nan") + parsed = float(value) + if not math.isfinite(parsed) or value_type(parsed) != value: + raise ValueError(f"{name} values must be exactly representable as float64") + return parsed + raise ValueError(f"{name} must contain only real numeric values") + + +def _lossless_float64_array( + name: str, value: np.ndarray, *, allow_nan: bool +) -> np.ndarray: + """Normalize one exact numeric ndarray without silently changing evidence.""" + if value.dtype.kind not in ("b", "i", "u", "f"): + if value.dtype.kind == "c": + raise ValueError(f"{name} must be real-valued") + raise ValueError(f"{name} must contain only real numeric values") + if np.isinf(value).any(): + raise ValueError(f"{name} must not contain infinite values") + if not allow_nan and np.isnan(value).any(): + raise ValueError(f"{name} must not contain NaN values") + + with np.errstate(invalid="ignore", over="ignore"): + converted = np.ascontiguousarray(value, dtype=np.float64) + if value.dtype.kind in ("i", "u"): + with np.errstate(invalid="ignore", over="ignore"): + round_trip = converted.astype(value.dtype) + if not np.array_equal(round_trip, value): + raise ValueError(f"{name} values must be exactly representable as float64") + elif value.dtype.kind == "f": + with np.errstate(invalid="ignore", over="ignore"): + round_trip = converted.astype(value.dtype) + if not np.array_equal(round_trip, value, equal_nan=True): + raise ValueError(f"{name} values must be exactly representable as float64") + return converted + + +def _trusted_matrix(name: str, value: object, *, allow_nan: bool) -> np.ndarray: + """Admit one bounded two-dimensional real-numeric matrix without protocols.""" + if type(value) is np.ndarray: + raw = value + if raw.ndim != 2: + raise ValueError(f"{name} must be two-dimensional") + if raw.size > _MAX_INTERACTION_MAP_CELLS: + raise ValueError( + f"{name} logical-cell count exceeds {_MAX_INTERACTION_MAP_CELLS}" + ) + return _lossless_float64_array(name, raw, allow_nan=allow_nan) + + if type(value) not in (list, tuple): + raise ValueError( + f"{name} must be an exact NumPy array or built-in two-dimensional sequence" + ) + if not value: + raise ValueError(f"{name} must be two-dimensional") + + width: int | None = None + normalized_rows: list[list[float]] = [] + logical_cells = 0 + for row in value: + if type(row) is np.ndarray: + if row.ndim != 1: + raise ValueError(f"{name} must be two-dimensional") + row_width = int(row.size) + logical_cells += row_width + if logical_cells > _MAX_INTERACTION_MAP_CELLS: + raise ValueError( + f"{name} logical-cell count exceeds {_MAX_INTERACTION_MAP_CELLS}" + ) + normalized_row = _lossless_float64_array( + name, row, allow_nan=allow_nan + ).tolist() + elif type(row) in (list, tuple): + row_width = len(row) + logical_cells += row_width + if logical_cells > _MAX_INTERACTION_MAP_CELLS: + raise ValueError( + f"{name} logical-cell count exceeds {_MAX_INTERACTION_MAP_CELLS}" + ) + normalized_row = [ + _trusted_scalar_float64(name, cell, allow_nan=allow_nan) for cell in row + ] + else: + raise ValueError(f"{name} must be two-dimensional") + + if width is None: + width = row_width + elif row_width != width: + raise ValueError(f"{name} must be rectangular") + normalized_rows.append(normalized_row) + + return np.ascontiguousarray(normalized_rows, dtype=np.float64) + + +def residual_interaction_map( + observed: np.ndarray, + expected: np.ndarray, + *, + axis_count: int, +) -> ResidualInteractionMap: + """Factor ``observed - expected`` using Gabriel symmetric scaling. + + Missing observed cells are represented by ``NaN`` and excluded through a + complete-case rectangle; they are never filled with zero. Model expectations + must be finite. Infinity is not treated as missing. ``axis_count`` is required + because the consuming measurement contract, not this library, determines how + many reader-visible axes are retained. Controls are sealed before caller + evidence is inspected; accepted matrices are callback-free, lossless at the + Rust ``f64`` boundary, and bounded before dense materialization. + + References: + Gabriel, K. R. (1971). The biplot graphic display of matrices with + application to principal component analysis. *Biometrika, 58*(3), + 453-467. https://doi.org/10.1093/biomet/58.3.453 + Jeon, M., Jin, I. H., Schweinberger, M., & Baugh, S. (2021). Mapping + unobserved item-respondent interactions: A latent space item + response model with interaction map. *Psychometrika, 86*(2), + 378-403. https://doi.org/10.1007/s11336-021-09762-5 + """ + axis_count_value = _axis_count(axis_count) + observed_array = _trusted_matrix("observed", observed, allow_nan=True) + expected_array = _trusted_matrix("expected", expected, allow_nan=False) + if observed_array.shape != expected_array.shape: + raise ValueError( + "observed and expected must have the same two-dimensional shape" + ) + + rows, columns = observed_array.shape + coordinate_cells = (rows + columns) * axis_count_value + if coordinate_cells > _MAX_INTERACTION_MAP_COORDINATE_CELLS: + raise ValueError( + "interaction map coordinate request exceeds " + f"{_MAX_INTERACTION_MAP_COORDINATE_CELLS} cells" + ) + + raw = dict( + _core.residual_interaction_map(observed_array, expected_array, axis_count_value) + ) + person_indices = np.asarray(raw["person_indices"], dtype=np.int64) + item_indices = np.asarray(raw["item_indices"], dtype=np.int64) + rows = person_indices.size + columns = item_indices.size + return ResidualInteractionMap( + person_indices=person_indices, + item_indices=item_indices, + scored_person_count=int(raw["scored_person_count"]), + scored_item_count=int(raw["scored_item_count"]), + person_coordinates=np.asarray( + raw["person_coordinates"], dtype=np.float64 + ).reshape(rows, axis_count_value), + item_coordinates=np.asarray(raw["item_coordinates"], dtype=np.float64).reshape( + columns, axis_count_value + ), + singular_values=np.asarray(raw["singular_values"], dtype=np.float64), + axis_shares=np.asarray(raw["axis_shares"], dtype=np.float64), + residual=np.asarray(raw["residual"], dtype=np.float64).reshape(rows, columns), + distance=np.asarray(raw["distance"], dtype=np.float64).reshape(rows, columns), + reconstruction=np.asarray(raw["reconstruction"], dtype=np.float64).reshape( + rows, columns + ), + unexplained=np.asarray(raw["unexplained"], dtype=np.float64).reshape( + rows, columns + ), + cross_share=np.asarray( + [np.nan if value is None else value for value in raw["cross_share"]], + dtype=np.float64, + ).reshape(rows, columns), + ) diff --git a/tests/test_interaction_map.py b/tests/test_interaction_map.py new file mode 100644 index 000000000..cdecb7cfe --- /dev/null +++ b/tests/test_interaction_map.py @@ -0,0 +1,284 @@ +"""Public residual interaction-map contract tests.""" + +from __future__ import annotations + +import importlib +import sys + +import numpy as np +import pytest +from fast_mlsirm import residual_interaction_map + +interaction_map_module = importlib.import_module("fast_mlsirm.interaction_map") + + +class _HostileArrayProvider: + def __array__(self, *_args: object, **_kwargs: object) -> np.ndarray: + raise AssertionError("caller __array__ must not execute") + + +class _HostileInt(int): + def __int__(self) -> int: + raise AssertionError("caller __int__ must not execute") + + def __le__(self, _other: object) -> bool: + raise AssertionError("caller comparison must not execute") + + +def _fake_map_payload(axis_count: int) -> dict[str, object]: + return { + "person_indices": [0], + "item_indices": [0], + "scored_person_count": 1, + "scored_item_count": 1, + "person_coordinates": [0.0] * axis_count, + "item_coordinates": [0.0] * axis_count, + "singular_values": [1.0], + "axis_shares": [1.0] + [0.0] * (axis_count - 1), + "residual": [1.0], + "distance": [0.0], + "reconstruction": [1.0], + "unexplained": [0.0], + "cross_share": [0.0], + "axis_count": axis_count, + } + + +def test_residual_interaction_map_preserves_rank_one_reconstruction() -> None: + """The Rust public API reconstructs a known rank-one residual exactly.""" + observed = np.array([[2.0, 0.0], [0.0, 2.0]]) + expected = np.ones((2, 2)) + + result = residual_interaction_map(observed, expected, axis_count=2) + + np.testing.assert_allclose(result.reconstruction, observed - expected, atol=1e-12) + np.testing.assert_allclose(result.residual, observed - expected, atol=1e-12) + assert np.all(np.isfinite(result.distance)) + assert np.all(result.distance >= 0.0) + np.testing.assert_allclose(result.axis_shares, [1.0, 0.0], atol=1e-12) + np.testing.assert_allclose(result.unexplained, 0.0, atol=1e-12) + np.testing.assert_allclose(result.cross_share, 0.0, atol=1e-12) + + +def test_residual_interaction_map_reports_scored_and_complete_case_coverage() -> None: + """Coverage distinguishes scored rows from rows admitted to the map.""" + observed = np.array([[2.0, np.nan], [1.0, 2.0]]) + result = residual_interaction_map(observed, np.ones((2, 2)), axis_count=2) + + assert result.scored_person_count == 2 + assert result.scored_item_count == 2 + np.testing.assert_array_equal(result.person_indices, [1]) + np.testing.assert_array_equal(result.item_indices, [0, 1]) + + +def test_no_complete_respondent_returns_shape_consistent_empty_map() -> None: + """An empty complete-case rectangle stays empty on both map axes.""" + observed = np.array([[1.0, np.nan], [np.nan, 1.0]]) + expected = np.zeros((2, 2)) + + result = residual_interaction_map(observed, expected, axis_count=2) + + assert result.person_indices.shape == (0,) + assert result.item_indices.shape == (0,) + assert result.person_coordinates.shape == (0, 2) + assert result.item_coordinates.shape == (0, 2) + assert result.residual.shape == (0, 0) + assert result.distance.shape == (0, 0) + assert result.reconstruction.shape == (0, 0) + assert result.unexplained.shape == (0, 0) + assert result.cross_share.shape == (0, 0) + np.testing.assert_array_equal(result.axis_shares, [0.0, 0.0]) + assert result.singular_values.shape == (0,) + + +@pytest.mark.parametrize("axis_count", [0, -1]) +def test_residual_interaction_map_rejects_nonpositive_axis_count( + axis_count: int, +) -> None: + """A consumer must request at least one reader-visible map axis.""" + with pytest.raises(ValueError, match="axis_count"): + residual_interaction_map( + np.ones((1, 1)), np.zeros((1, 1)), axis_count=axis_count + ) + + +@pytest.mark.parametrize("axis_count", [True, 1.5]) +def test_residual_interaction_map_rejects_nonintegral_axis_count( + axis_count: object, +) -> None: + """Boolean and real-valued controls do not silently become dimensions.""" + with pytest.raises(TypeError, match="axis_count"): + residual_interaction_map( + np.ones((1, 1)), + np.zeros((1, 1)), + axis_count=axis_count, # type: ignore[arg-type] + ) + + +def test_axis_control_is_sealed_before_caller_matrix_protocols() -> None: + """Invalid axis identities fail before any caller-owned data protocol.""" + with pytest.raises(TypeError, match="axis_count"): + residual_interaction_map( + _HostileArrayProvider(), # type: ignore[arg-type] + _HostileArrayProvider(), # type: ignore[arg-type] + axis_count=_HostileInt(1), + ) + + +def test_maximal_axis_request_fails_before_caller_matrix_protocols() -> None: + """A maximal positive axis request is resource-rejected before data work.""" + with pytest.raises(ValueError, match="coordinate"): + residual_interaction_map( + _HostileArrayProvider(), # type: ignore[arg-type] + _HostileArrayProvider(), # type: ignore[arg-type] + axis_count=sys.maxsize, + ) + + +def test_matrix_carriers_do_not_execute_arbitrary_array_protocols() -> None: + """Scientific evidence is rejected before caller ``__array__`` execution.""" + with pytest.raises(ValueError, match="observed"): + residual_interaction_map( + _HostileArrayProvider(), # type: ignore[arg-type] + np.zeros((1, 1)), + axis_count=1, + ) + + +def test_complex_and_infinite_evidence_fail_before_native_dispatch( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Imaginary evidence is not projected and infinity is not implicit missingness.""" + + def fail_core(*_args: object, **_kwargs: object) -> object: + raise AssertionError("compiled interaction-map core must not run") + + monkeypatch.setattr( + interaction_map_module._core, "residual_interaction_map", fail_core + ) + + with pytest.raises(ValueError, match="real-valued"): + residual_interaction_map( + np.array([[1.0 + 1.0j]]), np.zeros((1, 1)), axis_count=1 + ) + with pytest.raises(ValueError, match="infinite"): + residual_interaction_map(np.array([[np.inf]]), np.zeros((1, 1)), axis_count=1) + + +def test_integer_evidence_must_survive_float64_without_identity_loss( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Observed/expected values cannot silently change at the Rust ``f64`` boundary.""" + + def fail_core(*_args: object, **_kwargs: object) -> object: + raise AssertionError("compiled interaction-map core must not run") + + monkeypatch.setattr( + interaction_map_module._core, "residual_interaction_map", fail_core + ) + + with pytest.raises(ValueError, match="exactly representable"): + residual_interaction_map( + np.array([[2**53 + 1]], dtype=np.int64), + np.zeros((1, 1), dtype=np.float64), + axis_count=1, + ) + + +def test_logical_and_coordinate_budgets_precede_dense_or_native_work( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Logical matrix and coordinate requests are bounded before core dispatch.""" + + def fail_core(*_args: object, **_kwargs: object) -> object: + raise AssertionError("compiled interaction-map core must not run") + + monkeypatch.setattr( + interaction_map_module._core, "residual_interaction_map", fail_core + ) + monkeypatch.setattr( + interaction_map_module, "_MAX_INTERACTION_MAP_CELLS", 2, raising=False + ) + monkeypatch.setattr( + interaction_map_module, + "_MAX_INTERACTION_MAP_COORDINATE_CELLS", + 4, + raising=False, + ) + + oversized = np.broadcast_to(np.array([[1.0]]), (1, 3)) + with pytest.raises(ValueError, match="logical-cell"): + residual_interaction_map(oversized, oversized, axis_count=1) + + with pytest.raises(ValueError, match="coordinate"): + residual_interaction_map(np.ones((1, 1)), np.zeros((1, 1)), axis_count=3) + + +def test_trusted_numpy_and_builtin_evidence_reaches_core_canonically( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Safe NumPy scalars and built-in matrices preserve existing compatibility.""" + captured: dict[str, object] = {} + + def fake_core( + observed: np.ndarray, expected: np.ndarray, axis_count: int + ) -> dict[str, object]: + captured["observed"] = observed + captured["expected"] = expected + captured["axis_count"] = axis_count + return _fake_map_payload(axis_count) + + monkeypatch.setattr( + interaction_map_module._core, "residual_interaction_map", fake_core + ) + + result = residual_interaction_map( + [[np.float32(2.0)]], [[np.int16(1)]], axis_count=np.int16(1) + ) + + assert type(captured["axis_count"]) is int + assert captured["axis_count"] == 1 + assert isinstance(captured["observed"], np.ndarray) + assert isinstance(captured["expected"], np.ndarray) + assert captured["observed"].dtype == np.float64 # type: ignore[union-attr] + assert captured["expected"].dtype == np.float64 # type: ignore[union-attr] + np.testing.assert_array_equal(captured["observed"], [[2.0]]) + np.testing.assert_array_equal(captured["expected"], [[1.0]]) + assert result.person_coordinates.shape == (1, 1) + + +def test_nan_remains_the_only_nonfinite_missing_value( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """NaN remains admissible missingness while other nonfinite evidence is rejected.""" + captured: dict[str, np.ndarray] = {} + + def fake_core( + observed: np.ndarray, expected: np.ndarray, axis_count: int + ) -> dict[str, object]: + captured["observed"] = observed + captured["expected"] = expected + return _fake_map_payload(axis_count) + + monkeypatch.setattr( + interaction_map_module._core, "residual_interaction_map", fake_core + ) + + residual_interaction_map([[np.nan]], [[0.0]], axis_count=1) + assert np.isnan(captured["observed"][0, 0]) + + +def test_expected_nan_is_rejected_before_native_dispatch( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """A non-finite model expectation is invalid evidence, not response missingness.""" + + def fail_core(*_args: object, **_kwargs: object) -> object: + raise AssertionError("compiled interaction-map core must not run") + + monkeypatch.setattr( + interaction_map_module._core, "residual_interaction_map", fail_core + ) + + with pytest.raises(ValueError, match="expected"): + residual_interaction_map([[0.0]], [[np.nan]], axis_count=1) diff --git a/tests/test_logistic_dif_control_callback_safety.py b/tests/test_logistic_dif_control_callback_safety.py new file mode 100644 index 000000000..cf8b33dc1 --- /dev/null +++ b/tests/test_logistic_dif_control_callback_safety.py @@ -0,0 +1,229 @@ +"""Trust-boundary regressions for observed-score logistic DIF controls.""" + +from __future__ import annotations + +from typing import Any + +import numpy as np +import pytest + +import fast_mlsirm._dif_control_safety as dif_control_safety +import fast_mlsirm.dif as dif +import fast_mlsirm.fitstats as fitstats +from fast_mlsirm.dif import logistic_dif, logistic_dif_purified, mantel_haenszel_dif_purified + + +class _DataSentinel: + """Fail if rejected controls permit caller data materialization.""" + + def __array__(self, *args: Any, **kwargs: Any) -> np.ndarray: + """Raise on any attempt to materialize rejected caller data.""" + raise AssertionError("data must not be materialized for invalid controls") + + +def _unexpected_core_discovery(): + """Fail if rejected controls permit compiled-core discovery.""" + raise AssertionError("compiled core must not be discovered for invalid controls") + + +class _HostileFloat(float): + """Float subclass whose conversion hooks must stay inert.""" + + calls = 0 + + @classmethod + def reset(cls) -> None: + """Reset callback accounting.""" + cls.calls = 0 + + def __float__(self) -> float: + type(self).calls += 1 + return float.__float__(self) + + +class _HostileFloatProvider: + """Arbitrary float protocol provider whose callback must stay inert.""" + + calls = 0 + + @classmethod + def reset(cls) -> None: + """Reset callback accounting.""" + cls.calls = 0 + + def __float__(self) -> float: + type(self).calls += 1 + return 0.05 + + +class _HostileBool(int): + """Integer disguised as a boolean control.""" + + calls = 0 + + @classmethod + def reset(cls) -> None: + """Reset callback accounting.""" + cls.calls = 0 + + def __bool__(self) -> bool: + type(self).calls += 1 + return True + + +class _HostileInt(int): + """Integer subclass whose conversion hooks must stay inert.""" + + calls = 0 + + @classmethod + def reset(cls) -> None: + """Reset callback accounting.""" + cls.calls = 0 + + def __int__(self) -> int: + type(self).calls += 1 + return int.__int__(self) + + def __index__(self) -> int: + type(self).calls += 1 + return int.__index__(self) + + +@pytest.mark.parametrize( + ("kwargs", "value", "message"), + ( + ({"fdr_q": _HostileFloat(0.05)}, _HostileFloat, "fdr_q must be a real number"), + ({"fdr_q": _HostileFloatProvider()}, _HostileFloatProvider, "fdr_q must be a real number"), + ({"exclude_studied_item": _HostileBool(1)}, _HostileBool, "exclude_studied_item must be a bool"), + ({"max_iter": _HostileInt(50)}, _HostileInt, "max_iter must be an integer"), + ), +) +def test_logistic_rejects_executable_controls_before_callbacks_data_or_core( + monkeypatch, kwargs, value, message +): + """Executable logistic controls fail before callbacks, data, and Rust discovery.""" + monkeypatch.setattr(fitstats, "_core_module", _unexpected_core_discovery) + value.reset() + + with pytest.raises(ValueError, match=message): + logistic_dif(_DataSentinel(), _DataSentinel(), **kwargs) + + assert value.calls == 0 + + +@pytest.mark.parametrize( + ("function", "kwargs", "message"), + ( + (logistic_dif, {"fdr_q": 0.0}, r"fdr_q must be finite and in \(0, 1\]"), + (logistic_dif, {"fdr_q": True}, "fdr_q must be a real number"), + (logistic_dif, {"max_iter": -1}, "max_iter must be >= 1"), + (logistic_dif, {"max_iter": 0}, "max_iter must be >= 1"), + (logistic_dif, {"max_iter": True}, "max_iter must be an integer"), + (mantel_haenszel_dif_purified, {"max_rounds": 0}, "max_rounds must be >= 1"), + (mantel_haenszel_dif_purified, {"min_anchor_items": -1}, "min_anchor_items must be >= 0"), + (logistic_dif_purified, {"max_iter": -1}, "max_iter must be >= 1"), + (logistic_dif_purified, {"max_iter": 0}, "max_iter must be >= 1"), + (logistic_dif_purified, {"max_rounds": 0}, "max_rounds must be >= 1"), + (logistic_dif_purified, {"min_anchor_items": -1}, "min_anchor_items must be >= 0"), + ), +) +def test_invalid_exact_controls_fail_before_data_and_core(monkeypatch, function, kwargs, message): + """Exact but invalid controls fail before caller data and native discovery.""" + monkeypatch.setattr(fitstats, "_core_module", _unexpected_core_discovery) + + with pytest.raises(ValueError, match=message): + function(_DataSentinel(), _DataSentinel(), **kwargs) + + +@pytest.mark.parametrize( + ("function", "kwargs", "control_type", "message"), + ( + (mantel_haenszel_dif_purified, {"max_rounds": _HostileInt(3)}, _HostileInt, "max_rounds must be an integer"), + (mantel_haenszel_dif_purified, {"min_anchor_items": _HostileInt(4)}, _HostileInt, "min_anchor_items must be an integer"), + (logistic_dif_purified, {"max_iter": _HostileInt(50)}, _HostileInt, "max_iter must be an integer"), + (logistic_dif_purified, {"max_rounds": _HostileInt(3)}, _HostileInt, "max_rounds must be an integer"), + (logistic_dif_purified, {"min_anchor_items": _HostileInt(4)}, _HostileInt, "min_anchor_items must be an integer"), + ), +) +def test_purified_integer_controls_reject_subclasses_without_callbacks( + monkeypatch, function, kwargs, control_type, message +): + """Purification integer controls reject executable subclasses at the first boundary.""" + monkeypatch.setattr(fitstats, "_core_module", _unexpected_core_discovery) + control_type.reset() + + with pytest.raises(ValueError, match=message): + function(_DataSentinel(), _DataSentinel(), **kwargs) + + assert control_type.calls == 0 + + +class _FakeCore: + """Capture one valid native logistic dispatch.""" + + def __init__(self) -> None: + """Initialize call accounting.""" + self.calls: list[tuple[Any, ...]] = [] + + def logistic_dif(self, *args: Any) -> dict[str, Any]: + """Return a structurally valid Rust-like result.""" + self.calls.append(args) + return { + "item": [0, 1], + "chi2_uniform": [0.0, 0.0], + "p_uniform": [1.0, 1.0], + "chi2_nonuniform": [0.0, 0.0], + "p_nonuniform": [1.0, 1.0], + "chi2_total": [0.0, 0.0], + "p_total": [1.0, 1.0], + "delta_r2": [0.0, 0.0], + "delta_r2_uniform": [0.0, 0.0], + "jg_class": ["A", "A"], + "flagged_bh": [False, False], + "converged": [True, True], + } + + +def test_genuine_numpy_controls_dispatch_as_exact_builtins(monkeypatch): + """Supported NumPy scalar controls normalize once before the PyO3 boundary.""" + core = _FakeCore() + monkeypatch.setattr(fitstats, "_core_module", lambda: core) + responses = np.array([[0, 1], [1, 0], [1, 1], [0, 0]], dtype=np.int64) + group = np.array([0, 1, 0, 1], dtype=np.int64) + + result = logistic_dif( + responses, + group, + exclude_studied_item=np.bool_(True), + fdr_q=np.float32(0.05), + max_iter=np.int64(25), + ) + + assert list(result["jg_class"]) == ["A", "A"] + assert len(core.calls) == 1 + (_y, _g, n_persons, n_items, exclude_studied_item, fdr_q, max_iter) = core.calls[0] + assert type(n_persons) is int + assert type(n_items) is int + assert type(exclude_studied_item) is bool + assert type(fdr_q) is float + assert type(max_iter) is int + assert (n_persons, n_items, exclude_studied_item, max_iter) == (4, 2, True, 25) + assert fdr_q == pytest.approx(0.05, rel=0, abs=1e-6) + + +def test_dif_installer_is_idempotent(): + """Repeated package installation must not stack another validation layer.""" + first = ( + dif.logistic_dif, + dif.mantel_haenszel_dif_purified, + dif.logistic_dif_purified, + ) + + dif_control_safety.install(dif) + + assert ( + dif.logistic_dif, + dif.mantel_haenszel_dif_purified, + dif.logistic_dif_purified, + ) == first diff --git a/tests/test_logistic_dif_native_domain_compatibility.py b/tests/test_logistic_dif_native_domain_compatibility.py new file mode 100644 index 000000000..dae2dc5d9 --- /dev/null +++ b/tests/test_logistic_dif_native_domain_compatibility.py @@ -0,0 +1,148 @@ +"""Compatibility regressions for logistic DIF controls at the Rust boundary.""" + +from __future__ import annotations + +from typing import Any + +import numpy as np +import pytest + +import fast_mlsirm as fast_mlsirm +import fast_mlsirm.fitstats as fitstats +from fast_mlsirm import dif as dif_module +from fast_mlsirm.dif import logistic_dif, logistic_dif_purified + + +def _unexpected_core_discovery(): + """Fail if a rejected ``max_iter`` control reaches native-core discovery.""" + raise AssertionError("compiled core must not be discovered for invalid controls") + + +class _CaptureCore: + """Capture normalized controls while returning Rust-shaped DIF evidence.""" + + def __init__(self) -> None: + self.logistic_controls: tuple[Any, ...] | None = None + self.purified_controls: tuple[Any, ...] | None = None + + @staticmethod + def _result() -> dict[str, Any]: + return { + "item": [0], + "chi2_uniform": [0.0], + "p_uniform": [1.0], + "chi2_nonuniform": [0.0], + "p_nonuniform": [1.0], + "chi2_total": [0.0], + "p_total": [1.0], + "delta_r2": [0.0], + "delta_r2_uniform": [0.0], + "jg_class": ["A"], + "flagged_bh": [False], + "converged": [False], + } + + def logistic_dif(self, yy, gg, n_persons, n_items, exclude, fdr_q, max_iter): + self.logistic_controls = (exclude, fdr_q, max_iter) + return self._result() + + def logistic_dif_purified( + self, + yy, + gg, + n_persons, + n_items, + exclude, + fdr_q, + max_iter, + max_rounds, + min_anchor_items, + ): + self.purified_controls = ( + exclude, + fdr_q, + max_iter, + max_rounds, + min_anchor_items, + ) + return self._result() | { + "anchor": [True], + "n_anchor": 1, + "rounds": 1, + "purify_converged": True, + "purify_termination_reason": "stable_flag_set", + } + + +def test_logistic_dif_preserves_numpy_boolean_and_minimum_iteration_native_domain(monkeypatch): + core = _CaptureCore() + monkeypatch.setattr(fitstats, "_core_module", lambda: core) + responses = np.array([[0], [1]], dtype=np.int8) + group = np.array([0, 1], dtype=np.int8) + + logistic_dif( + responses, + group, + exclude_studied_item=np.bool_(True), + fdr_q=np.float32(0.05), + max_iter=np.int16(1), + ) + + assert core.logistic_controls is not None + exclude, fdr_q, max_iter = core.logistic_controls + assert type(exclude) is bool and exclude is True + assert type(fdr_q) is float + assert type(max_iter) is int and max_iter == 1 + + +def test_purified_logistic_preserves_minimum_iteration_native_domain(monkeypatch): + core = _CaptureCore() + monkeypatch.setattr(fitstats, "_core_module", lambda: core) + responses = np.array([[0], [1]], dtype=np.int8) + group = np.array([0, 1], dtype=np.int8) + + logistic_dif_purified( + responses, + group, + max_iter=np.uint8(1), + max_rounds=np.uint8(1), + min_anchor_items=np.uint8(0), + ) + + assert core.purified_controls is not None + exclude, fdr_q, max_iter, max_rounds, min_anchor_items = core.purified_controls + assert type(exclude) is bool + assert type(fdr_q) is float + assert type(max_iter) is int and max_iter == 1 + assert type(max_rounds) is int and max_rounds == 1 + assert type(min_anchor_items) is int and min_anchor_items == 0 + + +def test_logistic_dif_rejects_zero_max_iter_before_native_discovery(monkeypatch): + """Zero iterations fail at the Python boundary, matching the Rust ``logistic_sweep`` domain.""" + monkeypatch.setattr(fitstats, "_core_module", _unexpected_core_discovery) + + with pytest.raises(ValueError, match="max_iter must be >= 1"): + logistic_dif( + np.array([[0], [1]], dtype=np.int8), + np.array([0, 1], dtype=np.int8), + max_iter=0, + ) + + +def test_purified_logistic_dif_rejects_zero_max_iter_before_native_discovery(monkeypatch): + """Zero iterations fail at the Python boundary, matching the Rust ``logistic_sweep`` domain.""" + monkeypatch.setattr(fitstats, "_core_module", _unexpected_core_discovery) + + with pytest.raises(ValueError, match="max_iter must be >= 1"): + logistic_dif_purified( + np.array([[0], [1]], dtype=np.int8), + np.array([0, 1], dtype=np.int8), + max_iter=0, + ) + + +def test_package_aliases_share_hardened_dif_functions(): + assert fast_mlsirm.logistic_dif is dif_module.logistic_dif + assert fast_mlsirm.mantel_haenszel_dif_purified is dif_module.mantel_haenszel_dif_purified + assert fast_mlsirm.logistic_dif_purified is dif_module.logistic_dif_purified