feat(irt): expose Rust polytomous predictions

CHANGELOG.md

@@ -3,8 +3,55 @@
 ## Unreleased
 
 <!-- BEGIN AUTHORITATIVE CHANGELOG FRAGMENTS -->
+### Added
+
+#### Rust longitudinal state layer
+
+- Added a Rust-owned independent per-respondent OLS trend and discrete-sequence
+  AR(1) state predictor behind the sealed `fast_mlsirm.multilevel` contract.
+- Preserved exact sequence gaps, missing-occasion output state, deterministic
+  respondent sharding, RMSE/count diagnostics, and PyO3/Python marshalling.
+- Documented the compatibility wire label `random_intercept_slope` as independent
+  OLS with no population random-effects distribution or shrinkage, and the AR
+  path as caller-supplied `phi` without coefficient estimation.
+- Added slope-recovery, missingness, irregular-calendar/non-contiguous-sequence,
+  worker-determinism, and fail-closed contract tests with APA 7 doctoring.
+- This fragment does not claim full multilevel IRT random-effect integration,
+  uncertainty, continuous-time transitions, or GPU recurrent-state parity.
+
+#### Joint MAP hierarchical continuous-time AR(1) Rasch
+
+- Added a Rust-owned joint MAP hierarchical continuous-time AR(1) Rasch
+  estimator behind `fit_hierarchical_longitudinal_irt`, stacked on the
+  `#976` longitudinal design handoff.
+- Estimated shared population hyperparameters `(mu, tau, lambda)` and person-
+  occasion states from exact millisecond elapsed-day gaps. State intervals
+  are Wald intervals from measurement observed information; short series
+  leave `lambda` weakly identified under joint MAP.
+- Documented the estimand as joint MAP, not independent OLS, not caller-
+  supplied discrete AR, not Fox and Glas Gibbs, and not estimated
+  multiple-membership `u_h`. GPU parity is reported false because the
+  existing wgpu path owns a different MLSIRM objective.
+- Added multi-seed true-parameter recovery, irregular-time, missing-response,
+  worker-determinism, and fail-closed marshalling tests with APA 7 ADR and
+  doctoring.
+- Normalized oversized integer and non-finite real execution controls before
+  native dispatch so Python callers receive package-owned validation errors.
+- Enforced finite, identified sum-zero Rasch item intercepts in the simulator
+  before generating recovery data.
+- Labeled hyperparameter intervals as conditional on fixed item/state nuisance
+  blocks and aligned CT-AR gradients with the active variance branch.
+
 ### 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.
@@ -21,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
@@ -56,6 +109,12 @@
 
 ### Fixed
 
+#### Harden S-X² scalar control admission
+
+- Reject caller-defined integer and floating subclasses at the public S-X² control boundary before numeric conversion or compiled-core dispatch, while preserving exact built-in and concrete NumPy scalar compatibility and leaving all S-X²/G², quadrature, and BH/FDR arithmetic Rust-owned.
+- Reject built-in or concrete NumPy integer-valued real controls when float64 normalization would change the integer identity, so `min_expected`, `fdr_q`, and `min_effect` cannot be silently rounded before domain validation.
+- Reject finite extended-precision NumPy floating controls when conversion to the Rust `f64` boundary would change their value, while preserving exact float16/float32/float64 and lossless long-double controls.
+
 #### Validate G-theory controls and score evidence before Rust discovery
 
 - G-theory D-study sizes and `Phi(lambda)` scalar controls now fail closed before caller-owned score-array materialization and before compiled Rust capability discovery when invalid, while preserving the existing callback-free Python/NumPy scalar contract.
@@ -163,6 +222,40 @@
 
 - Reject Rating Scale Model response evidence above 20,000,000 logical cells during the callback-free source preflight, including oversized exact NumPy arrays and exact NumPy rows nested in trusted built-in sequences, before NumPy stacking or contiguous `float64` allocation while preserving existing response semantics and Rust-owned Andrich arithmetic.
 
+#### Bound polytomous prediction resources
+
+- Reject public GRM/GPCM prediction grids above 20,000,000 dense probability cells before compiled-core discovery or output allocation.
+- 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.
+
+#### Polytomous prediction evidence admission
+
+- 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.
+
+#### Polytomous prediction evidence rank admission
+
+- 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.
+
+#### Polytomous prediction category-domain admission
+
+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.
+
+#### Preserve polytomous prediction fit metadata
+
+- 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.
+
+#### Preserve strict GRM threshold order during initialization
+
+- 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.
+
+#### Polytomous raw prediction category-domain replay
+
+- 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.
+
 #### Linking evidence admission
 
 - Seal fixed-anchor and common-item IRT linking numeric evidence before NumPy materialization so arbitrary array/container/numeric protocols cannot execute while scientific inputs are being admitted.

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,
@@ -6303,6 +6304,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<Py<pyo3::types::PyDict>> {
+    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]
@@ -9368,11 +9396,51 @@ fn standard_errors_from_vcov(vcov: PyReadonlyArray2<'_, f64>) -> PyResult<Vec<f6
 /// Version of the Python-to-Rust marginal-MMLE call contract.
 const MARGINAL_CAPABILITY_VERSION: u32 = 1;
 
+#[pyfunction]
+fn residual_interaction_map(
+    py: Python<'_>,
+    observed: PyReadonlyArray2<'_, f64>,
+    expected: PyReadonlyArray2<'_, f64>,
+    axis_count: usize,
+) -> PyResult<Py<PyAny>> {
+    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)?)?;
@@ -9561,6 +9629,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)?)?;

crates/mlsirm-core/src/entrypoint.rs

@@ -6,4 +6,5 @@
 
 include!("lib.rs");
 
+pub mod interaction_map;
 pub mod rotation;

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<f64>, Vec<f64>), String> {
+pub(crate) fn symmetric_eigen_desc(
+    matrix: &[f64],
+    p: usize,
+) -> Result<(Vec<f64>, Vec<f64>), String> {
     let mut a = matrix.to_vec();
     let mut v = vec![0.0; p * p];
     for i in 0..p {