Add joint MAP hierarchical CT-AR Rasch slice (stacked on #976)

ARCHITECTURE.md

@@ -206,10 +206,14 @@ hierarchical, multiply affiliated, or longitudinal:
   temporal validity rules.
 
 The `fast_mlsirm.multilevel` contracts are content-addressed and fail closed.
-Nested estimators that consume multiple-membership and longitudinal contracts
-remain explicitly paper-scoped until their Rust implementation and recovery
-evidence are complete; the presence of a contract is not a claim that the
-estimator is already production-ready.
+The ADR-0018 state layer fits independent OLS trends and caller-supplied
+discrete AR predictions. ADR-0019 adds a separate joint MAP hierarchical
+continuous-time AR(1) Rasch slice with estimated `(mu, tau, lambda)`, elapsed-day
+transitions, and Wald observed-information intervals. That slice excludes
+estimated multiple-membership `u_h` and does not claim GPU parity. Remaining
+nested/crossed estimators stay paper-scoped until their own Rust
+implementation and recovery evidence are complete; the presence of a contract
+is not a claim that every estimator is production-ready.
 
 ## 5. Numerical and scientific architecture
 

CHANGELOG.md

@@ -3,8 +3,53 @@
 ## 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.
+
 ### Changed
 
+#### Govern non-psychometric item-bank suspension concerns
+
+- Governed item-bank suspension and reactivation can now bind exact non-psychometric concern evidence for evidence/content validity and security/privacy findings, alongside existing DIF, drift, exposure, and linking evidence, without fabricating psychometric drift evidence.
+- Suspended records bind the exact newly asserted concern classes into their content-addressed identity, and reactivation requires fresh evidence for those same classes so unrelated evidence cannot clear a quarantine.
+- Reactivation rejects a historical approval or concern fingerprint even when it is presented under a replacement evidence identifier; every required reactivation artifact must bind new evidence content.
+
+#### Harden remaining equating controls before native discovery
+
+- Validate circle-arc method/point/scalar controls, nominal-weights score ceilings and synthetic-population weight, and the composite-linking exponent before compiled-core discovery.
+- Reject caller-defined scalar/container subclasses and arbitrary conversion providers without executing their conversion, comparison, representation, hashing, or iteration callbacks.
+- Preserve exact built-in and genuine NumPy scalar compatibility while keeping circle-arc geometry, nominal-weights moments, composite-linking weight arithmetic, and all result-affecting equating mathematics in Rust.
+
 #### Release cut 0.8.0
 
 - Project version is bumped to 0.8.0 in `pyproject.toml`, `crates/mlsirm-core`,
@@ -26,6 +71,49 @@
   Python 3.14 CI).
 - Released authoritative fragments are removed from `docs/changelog.d`; the
   directory again holds only genuinely unreleased notes.
+
+### Fixed
+
+#### Fail closed on unsafe multilevel contextual effects
+
+- Multilevel contextual-effect evaluation now fails closed when any referenced context random-effect value is NaN or infinite and when finite inputs overflow the weighted sum, preventing non-finite predictor results from escaping the Rust boundary while leaving unreferenced table capacity outside sparse validation work.
+- Python context-effect marshalling snapshots each required mapping value once without caller-defined membership probes and normalizes hostile lookup callbacks to non-reflective package errors before native dispatch.
+
+#### Workflow-registry audit transport retry hardening
+
+- Expanded the read-only Actions-registry audit transport's bounded retry classifier to cover transient HTTP 403, 404, 429, and all 5xx responses, while preserving fail-closed exhaustion and immediate failure for non-transient authentication errors such as HTTP 401.
+- Added direct transport regression coverage so incident audits do not misclassify one transient GitHub control-plane response as a completed inventory failure.
+
+#### Harden RAG metadata callback safety
+
+- Validate caller-provided RAG metadata keys exactly once before reading any values, then freeze only the captured allowlisted values. Hostile membership, key/value, duplicate-key, and key-reiteration callbacks now fail through non-reflective package errors without granting new metadata authority.
+
+#### Exposure-control scalar callback safety
+
+- Validate CAT/exposure integer controls from exact built-in Python and genuine NumPy scalar types before caller-dispatchable coercion or Rust-core discovery, preserving integral built-in/NumPy floating controls, package-owned bounds/errors, and Rust-owned exposure, routing, scoring, posterior, recovery, and simulation arithmetic.
+
+#### Harden scoring-policy integer callback boundaries
+
+- Reject caller-defined integer coercion at scoring-policy positive-integer boundaries before any `__index__` callback can run, while preserving exact built-in and genuine NumPy integer scalar compatibility and existing bounded `AssessmentSpecError` semantics.
+
+#### ATA integer callback safety
+
+- Automated test assembly now admits only exact built-in integers and explicitly supported genuine NumPy integer scalar identities for public length, seed, exposure, content-count, and exclusion controls before normalization.
+- Caller-defined Python and NumPy integer subclasses fail closed before conversion callbacks or item-information work, while existing finite-domain validation and genuine NumPy scalar compatibility are preserved.
+- Added focused public-boundary regressions for hostile scalar and container controls without changing ATA information, selection, or scoring arithmetic.
+
+#### Fleiss kappa control trust boundary
+
+- Hardened the public Fleiss/Conger kappa control boundary so explicit category counts and exact-mode selection are validated without executing caller-defined integer, index, or truthiness callbacks before ratings materialization or compiled-core discovery.
+- Preserved genuine Python/NumPy scalar compatibility, capped explicit and inferred category counts at the Rust contract maximum of 10,000, and kept all agreement arithmetic Rust-owned.
+
+### Security
+
+#### Parallel-analysis control trust hardening
+
+- Validate `n_iterations`, `centile`, and `seed` before native-core discovery, accepting only exact built-in integers and supported concrete NumPy integer scalars while rejecting booleans, `np.bool_`, caller-defined subclasses, and conversion providers without executing their callbacks. Workspace and `u64` seed limits fail at the same pre-discovery boundary.
+- Normalize nonnumeric `data` conversion failures to a package-owned `ValueError` before native-core discovery while preserving dimensionality and workspace validation for successfully converted arrays.
+- Preserve the existing positive-iteration, centile `0..99`, Rust `u64` seed, and 128 MiB random-benchmark workspace limits without changing Rust-owned Horn/Glorfeld factor-retention arithmetic.
 <!-- END AUTHORITATIVE CHANGELOG FRAGMENTS -->
 ## [0.8.0] - 2026-08-17
 

crates/fast-mlsirm-py/src/multilevel_bindings.rs

@@ -7,8 +7,13 @@
 //! `mlsirm_core::multilevel::weighted_contextual_effect`.
 
 use mlsirm_core::longitudinal::fit_longitudinal_state as core_fit_longitudinal_state;
+use mlsirm_core::longitudinal_irt::{
+    fit_hierarchical_ctar_rasch as core_fit_hierarchical_ctar_rasch,
+    simulate_hierarchical_ctar_rasch as core_simulate_hierarchical_ctar_rasch,
+    HierarchicalCtarRaschConfig,
+};
 use mlsirm_core::multilevel::weighted_contextual_effect as core_weighted_contextual_effect;
-use numpy::{PyArray1, PyReadonlyArray1, ToPyArray};
+use numpy::{PyArray1, PyReadonlyArray1, PyReadonlyArray2, PyUntypedArrayMethods, ToPyArray};
 use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
 use pyo3::types::PyDict;
@@ -19,6 +24,8 @@ use pyo3::wrap_pyfunction;
 // contract. The regression tests import the Python constant so drift fails CI.
 const MAX_CONTEXT_MEMBERSHIPS: usize = 100_000;
 const MAX_ROW_OFFSETS: usize = MAX_CONTEXT_MEMBERSHIPS + 1;
+const MAX_HIERARCHICAL_OCCASIONS: usize = 100_000;
+const MAX_HIERARCHICAL_ITEMS: usize = 4_096;
 
 fn checked_usize_values(values: &[u64], name: &str) -> PyResult<Vec<usize>> {
     values
@@ -153,10 +160,183 @@ fn py_fit_longitudinal_state<'py>(
     Ok(result)
 }
 
+/// Fit the joint MAP hierarchical continuous-time AR(1) Rasch slice.
+///
+/// Parameters
+/// ----------
+/// row_offsets : numpy.ndarray[uint64]
+///     CSR-style respondent pointer, length ``n_respondents + 1``.
+/// time_offsets_milliseconds : numpy.ndarray[int64]
+///     Exact millisecond offsets aligned with the occasion axis of
+///     ``responses``.
+/// responses : numpy.ndarray[float64]
+///     Occasion-major binary matrix with shape
+///     ``(n_occasions, n_items)``. ``NaN`` marks a missing response.
+/// worker_count : int
+///     Number of deterministic person-shard worker threads (``>= 1``).
+/// max_iter : int
+///     Maximum packed L-BFGS iterations (``>= 1``).
+/// tolerance : float
+///     Relative L-BFGS tolerance; must be finite and strictly positive.
+/// hessian_step : float
+///     Central-difference step for the hyperparameter Hessian.
+///
+/// Returns
+/// -------
+/// dict
+///     Joint MAP states, Wald intervals, item intercepts, estimated
+///     ``(mu, tau, lambda)``, and normative estimand metadata.
+#[pyfunction(name = "fit_hierarchical_ctar_rasch")]
+fn py_fit_hierarchical_ctar_rasch<'py>(
+    py: Python<'py>,
+    row_offsets: PyReadonlyArray1<'_, u64>,
+    time_offsets_milliseconds: PyReadonlyArray1<'_, i64>,
+    responses: PyReadonlyArray2<'_, f64>,
+    worker_count: usize,
+    max_iter: usize,
+    tolerance: f64,
+    hessian_step: f64,
+) -> PyResult<Bound<'py, PyDict>> {
+    let shape = responses.shape();
+    if shape[0] > MAX_HIERARCHICAL_OCCASIONS {
+        return Err(PyValueError::new_err(format!(
+            "responses occasion axis exceeds maximum supported length of {MAX_HIERARCHICAL_OCCASIONS}"
+        )));
+    }
+    if shape[1] > MAX_HIERARCHICAL_ITEMS {
+        return Err(PyValueError::new_err(format!(
+            "responses item axis exceeds maximum supported length of {MAX_HIERARCHICAL_ITEMS}"
+        )));
+    }
+    let row_offsets = checked_usize_values(row_offsets.as_slice()?, "row_offsets")?;
+    let time_offsets = time_offsets_milliseconds.as_slice()?.to_vec();
+    let responses = responses.as_slice()?.to_vec();
+    let n_items = shape[1];
+    let config = HierarchicalCtarRaschConfig {
+        worker_count,
+        max_iter,
+        tolerance,
+        hessian_step,
+    };
+    let fit = py
+        .detach(move || {
+            core_fit_hierarchical_ctar_rasch(
+                &row_offsets,
+                &time_offsets,
+                &responses,
+                n_items,
+                config,
+            )
+        })
+        .map_err(PyValueError::new_err)?;
+    let result = PyDict::new(py);
+    result.set_item("state", fit.state.to_pyarray(py))?;
+    result.set_item("state_se", fit.state_se.to_pyarray(py))?;
+    result.set_item("state_lower", fit.state_lower.to_pyarray(py))?;
+    result.set_item("state_upper", fit.state_upper.to_pyarray(py))?;
+    result.set_item("item_intercepts", fit.item_intercepts.to_pyarray(py))?;
+    result.set_item("population_mean", fit.population_mean)?;
+    result.set_item("population_sd", fit.population_sd)?;
+    result.set_item("decay_rate", fit.decay_rate)?;
+    result.set_item("unit_time_ar_coefficient", fit.unit_time_ar_coefficient)?;
+    result.set_item("hyperparameter_se", fit.hyperparameter_se.to_vec())?;
+    result.set_item("hyperparameter_lower", fit.hyperparameter_lower.to_vec())?;
+    result.set_item("hyperparameter_upper", fit.hyperparameter_upper.to_vec())?;
+    result.set_item(
+        "hyperparameter_intervals_identified",
+        fit.hyperparameter_intervals_identified,
+    )?;
+    result.set_item("state_intervals_identified", fit.state_intervals_identified)?;
+    result.set_item("observed_count", fit.observed_count)?;
+    result.set_item("transition_count", fit.transition_count)?;
+    result.set_item("status", fit.status)?;
+    result.set_item("estimand_scope", fit.estimand_scope)?;
+    result.set_item("transition_kind", fit.transition_kind)?;
+    result.set_item("interval_kind", fit.interval_kind)?;
+    result.set_item("engine", fit.engine)?;
+    result.set_item("population_random_effects_estimated", true)?;
+    result.set_item("ar_coefficient_estimated", true)?;
+    result.set_item("ar_coefficient_source", "joint_map")?;
+    result.set_item("multiple_membership_estimated", false)?;
+    result.set_item("gpu_parity", false)?;
+    Ok(result)
+}
+
+/// Simulate hierarchical continuous-time AR(1) Rasch responses.
+///
+/// Parameters
+/// ----------
+/// row_offsets : numpy.ndarray[uint64]
+///     CSR-style respondent pointer, length ``n_respondents + 1``.
+/// time_offsets_milliseconds : numpy.ndarray[int64]
+///     Exact millisecond offsets aligned with the generated occasions.
+/// item_intercepts : numpy.ndarray[float64]
+///     Sum-to-zero Rasch item intercepts used as generating values.
+/// population_mean : float
+///     Generating population mean.
+/// population_sd : float
+///     Generating stationary standard deviation.
+/// decay_rate : float
+///     Generating continuous-time decay rate per day.
+/// seed : int
+///     Deterministic LCG seed.
+///
+/// Returns
+/// -------
+/// dict
+///     Generating latent states and occasion-major binary responses.
+#[pyfunction(name = "simulate_hierarchical_ctar_rasch")]
+fn py_simulate_hierarchical_ctar_rasch<'py>(
+    py: Python<'py>,
+    row_offsets: PyReadonlyArray1<'_, u64>,
+    time_offsets_milliseconds: PyReadonlyArray1<'_, i64>,
+    item_intercepts: PyReadonlyArray1<'_, f64>,
+    population_mean: f64,
+    population_sd: f64,
+    decay_rate: f64,
+    seed: u64,
+) -> PyResult<Bound<'py, PyDict>> {
+    let item_intercepts = item_intercepts.as_slice()?.to_vec();
+    if item_intercepts.len() > MAX_HIERARCHICAL_ITEMS {
+        return Err(PyValueError::new_err(format!(
+            "item_intercepts exceeds maximum supported length of {MAX_HIERARCHICAL_ITEMS}"
+        )));
+    }
+    let n_items = item_intercepts.len();
+    let row_offsets = checked_usize_values(row_offsets.as_slice()?, "row_offsets")?;
+    let time_offsets = time_offsets_milliseconds.as_slice()?.to_vec();
+    if time_offsets.len() > MAX_HIERARCHICAL_OCCASIONS {
+        return Err(PyValueError::new_err(format!(
+            "time offsets exceed maximum supported length of {MAX_HIERARCHICAL_OCCASIONS}"
+        )));
+    }
+    let (state, responses) = py
+        .detach(move || {
+            core_simulate_hierarchical_ctar_rasch(
+                &row_offsets,
+                &time_offsets,
+                n_items,
+                population_mean,
+                population_sd,
+                decay_rate,
+                &item_intercepts,
+                seed,
+            )
+        })
+        .map_err(PyValueError::new_err)?;
+    let result = PyDict::new(py);
+    result.set_item("state", state.to_pyarray(py))?;
+    result.set_item("responses", responses.to_pyarray(py))?;
+    result.set_item("n_items", n_items)?;
+    Ok(result)
+}
+
 #[pymodule]
 #[pyo3(name = "_multilevel_core")]
 fn fast_mlsirm_multilevel_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(py_weighted_contextual_effect, m)?)?;
     m.add_function(wrap_pyfunction!(py_fit_longitudinal_state, m)?)?;
+    m.add_function(wrap_pyfunction!(py_fit_hierarchical_ctar_rasch, m)?)?;
+    m.add_function(wrap_pyfunction!(py_simulate_hierarchical_ctar_rasch, m)?)?;
     Ok(())
 }

crates/mlsirm-core/src/lib.rs

@@ -19,6 +19,7 @@ pub mod ksirt;
 pub mod linking;
 pub mod lltm;
 pub mod longitudinal;
+pub mod longitudinal_irt;
 pub mod marginal;
 pub mod mhrm;
 pub mod mixed;