diff --git a/ARCHITECTURE.md b/ARCHITECTURE.md index b17fcb3e4..188205187 100644 --- a/ARCHITECTURE.md +++ b/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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 52942d2ba..e1a4bddb5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,8 +3,53 @@ ## Unreleased +### 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. ## [0.8.0] - 2026-08-17 diff --git a/crates/fast-mlsirm-py/src/multilevel_bindings.rs b/crates/fast-mlsirm-py/src/multilevel_bindings.rs index 436ee10f1..9ee6acebb 100644 --- a/crates/fast-mlsirm-py/src/multilevel_bindings.rs +++ b/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> { 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> { + 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> { + 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(()) } diff --git a/crates/mlsirm-core/src/lib.rs b/crates/mlsirm-core/src/lib.rs index db910f301..6a3b64a52 100644 --- a/crates/mlsirm-core/src/lib.rs +++ b/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; diff --git a/crates/mlsirm-core/src/longitudinal_irt.rs b/crates/mlsirm-core/src/longitudinal_irt.rs new file mode 100644 index 000000000..46fa52759 --- /dev/null +++ b/crates/mlsirm-core/src/longitudinal_irt.rs @@ -0,0 +1,1592 @@ +//! Joint MAP hierarchical continuous-time AR(1) Rasch estimator. +//! +//! This module is the smallest jointly estimated longitudinal latent-state IRT +//! slice stacked on the independent-OLS / caller-supplied-AR state layer. The +//! estimand is **joint maximum a posteriori** (MAP) of a Rasch measurement +//! model and a hierarchical stationary Ornstein–Uhlenbeck / continuous-time +//! AR(1) latent-state process: +//! +//! ```text +//! logit P(Y_pti = 1) = theta_pt - b_i, sum_i b_i = 0 +//! theta_p,1 ~ N(mu, tau^2) +//! theta_p,t | theta_p,t-1 ~ N( +//! mu + exp(-lambda * Delta_pt) * (theta_p,t-1 - mu), +//! tau^2 * (1 - exp(-2 * lambda * Delta_pt)) +//! ) +//! ``` +//! +//! `Delta_pt` is the elapsed time in days from exact millisecond offsets. The +//! shared `(mu, tau, lambda)` hyperparameters are estimated, so person-occasion +//! states are shrunk toward the population mean. This is **not** independent +//! respondent OLS, **not** a caller-supplied discrete AR coefficient, **not** +//! Fox and Glas (2001) Gibbs sampling, and **not** Jeon and Rabe-Hesketh (2016) +//! adaptive-quadrature ML. Crossed / multiple-membership random effects are +//! excluded from this joint likelihood. +//! +//! # References (APA 7th ed.) +//! +//! Fox, J.-P., & Glas, C. A. W. (2001). Bayesian estimation of a multilevel +//! IRT model. *Psychometrika, 66*, 271–288. +//! https://doi.org/10.1007/BF02294839 +//! +//! Jeon, M., & Rabe-Hesketh, S. (2016). An autoregressive growth model for +//! longitudinal item analysis. *Psychometrika, 81*(3), 830–850. +//! https://doi.org/10.1007/s11336-015-9489-2 +//! +//! Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal +//! data. *Biometrics, 38*(4), 963–974. https://doi.org/10.2307/2529876 +//! +//! Oravecz, Z., Tuerlinckx, F., & Vandekerckhove, J. (2011). A hierarchical +//! latent stochastic differential equation model for affective dynamics. +//! *Psychological Methods, 16*(2), 468–490. https://doi.org/10.1037/a0024375 + +use std::f64::consts::PI; +use std::thread; + +use crate::jmle_opt::lbfgs; +use crate::mmle::{log_sigmoid, sigmoid_stable}; + +const MILLIS_PER_DAY: f64 = 86_400_000.0; +const MAX_ABS_TIME_DAYS: f64 = 10_000_000.0; +const MIN_TRANSITION_VARIANCE: f64 = 1e-12; +const MIN_LOG_SD: f64 = -4.0; +const MAX_LOG_SD: f64 = 2.5; +const MIN_LOG_DECAY: f64 = -5.0; +const MAX_LOG_DECAY: f64 = 2.0; +const WALD_Z: f64 = 1.959963984540054; +const ESTIMAND_SCOPE: &str = "joint_map_hierarchical_ctar_rasch"; +const TRANSITION_KIND: &str = "continuous_time_ar1_ou"; +const INTERVAL_KIND: &str = "wald_measurement_observed_information"; +const ENGINE: &str = "rust_cpu_multithreaded"; + +/// Configuration for [`fit_hierarchical_ctar_rasch`]. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct HierarchicalCtarRaschConfig { + /// Deterministic person-shard worker count. Values larger than the person + /// count are capped. Must be at least one. + pub worker_count: usize, + /// Maximum L-BFGS iterations for the joint MAP packed vector. + pub max_iter: usize, + /// Relative L-BFGS tolerance on the packed objective. + pub tolerance: f64, + /// Central-difference step for the hyperparameter observed Hessian. + pub hessian_step: f64, +} + +impl Default for HierarchicalCtarRaschConfig { + fn default() -> Self { + Self { + worker_count: 1, + max_iter: 250, + tolerance: 1e-5, + hessian_step: 1e-3, + } + } +} + +/// Joint MAP fit of the hierarchical continuous-time AR(1) Rasch slice. +#[derive(Clone, Debug, PartialEq)] +pub struct HierarchicalCtarRaschFit { + /// Person-occasion latent states aligned with the occasion input. + pub state: Vec, + /// Conditional observed-information standard errors for `state`. + pub state_se: Vec, + /// Lower 95% Wald bounds for `state`. + pub state_lower: Vec, + /// Upper 95% Wald bounds for `state`. + pub state_upper: Vec, + /// Sum-to-zero Rasch item intercepts. + pub item_intercepts: Vec, + /// Estimated population mean of the latent-state process. + pub population_mean: f64, + /// Estimated stationary standard deviation. + pub population_sd: f64, + /// Estimated continuous-time decay rate (per day). + pub decay_rate: f64, + /// Unit-day AR coefficient `exp(-decay_rate)`. + pub unit_time_ar_coefficient: f64, + /// Standard errors for `[mean, sd, decay]` when identified. + pub hyperparameter_se: [f64; 3], + /// Lower 95% Wald bounds for `[mean, sd, decay]`. + pub hyperparameter_lower: [f64; 3], + /// Upper 95% Wald bounds for `[mean, sd, decay]`. + pub hyperparameter_upper: [f64; 3], + /// Whether the hyperparameter observed Hessian produced finite SEs. + pub hyperparameter_intervals_identified: bool, + /// Whether every person-block state Hessian produced finite SEs. + pub state_intervals_identified: bool, + /// Number of finite binary responses used in the measurement term. + pub observed_count: usize, + /// Number of person-level CT-AR transitions used in the state prior. + pub transition_count: usize, + /// Optimizer status from the packed L-BFGS run. + pub status: String, + /// Normative estimand label. Never an OLS or caller-supplied AR label. + pub estimand_scope: &'static str, + /// Transition family actually parameterized by elapsed time. + pub transition_kind: &'static str, + /// Interval construction actually computed. + pub interval_kind: &'static str, + /// Compute engine identity. + pub engine: &'static str, +} + +/// Continuous-time AR(1) / OU autoregressive weight for an elapsed interval. +pub fn ctar_phi(decay_rate: f64, delta_days: f64) -> Result { + if !decay_rate.is_finite() || decay_rate <= 0.0 { + return Err("decay_rate must be finite and strictly positive".to_string()); + } + if !delta_days.is_finite() || delta_days <= 0.0 { + return Err("elapsed days must be finite and strictly positive".to_string()); + } + let phi = (-decay_rate * delta_days).exp(); + if !phi.is_finite() { + return Err("continuous-time AR weight is not finite".to_string()); + } + Ok(phi) +} + +/// Stationary OU transition variance for an elapsed interval. +pub fn ctar_variance(population_variance: f64, decay_rate: f64, delta_days: f64) -> Result { + if !population_variance.is_finite() || population_variance <= 0.0 { + return Err("population variance must be finite and strictly positive".to_string()); + } + let phi = ctar_phi(decay_rate, delta_days)?; + let two_lambda_delta = 2.0 * decay_rate * delta_days; + let one_minus_phi2 = if two_lambda_delta < 1e-8 { + two_lambda_delta - 0.5 * two_lambda_delta * two_lambda_delta + } else { + 1.0 - phi * phi + }; + let variance = population_variance * one_minus_phi2; + if !variance.is_finite() || variance <= 0.0 { + return Err("continuous-time transition variance is degenerate".to_string()); + } + Ok(variance.max(MIN_TRANSITION_VARIANCE)) +} + +/// Unit-day AR coefficient implied by a positive decay rate. +pub fn ctar_unit_phi(decay_rate: f64) -> Result { + ctar_phi(decay_rate, 1.0) +} + +fn validate_offsets(row_offsets: &[usize], n_occasions: usize) -> Result { + if row_offsets.is_empty() || row_offsets[0] != 0 { + return Err("row_offsets must be non-empty and start at zero".to_string()); + } + if row_offsets.windows(2).any(|window| window[1] < window[0]) { + return Err("row_offsets must be non-decreasing".to_string()); + } + if row_offsets.last().copied() != Some(n_occasions) { + return Err("row_offsets must end at the occasion count".to_string()); + } + Ok(row_offsets.len() - 1) +} + +fn days_from_millis(offset: i64) -> Result { + let days = offset as f64 / MILLIS_PER_DAY; + if !days.is_finite() || days.abs() > MAX_ABS_TIME_DAYS { + return Err("time offsets exceed the supported finite range".to_string()); + } + Ok(days) +} + +fn validate_design( + row_offsets: &[usize], + time_offsets_milliseconds: &[i64], + responses: &[f64], + n_items: usize, +) -> Result<(usize, usize), String> { + if n_items < 2 { + return Err("hierarchical CT-AR Rasch requires at least two items".to_string()); + } + let n_occasions = time_offsets_milliseconds.len(); + let expected = crate::checked_mul_usize(n_occasions, n_items, "response array exceeds supported size")?; + if responses.len() != expected { + return Err("responses must be occasion-major with n_occasions * n_items entries".to_string()); + } + let n_persons = validate_offsets(row_offsets, n_occasions)?; + if n_persons == 0 { + return Err("at least one respondent is required".to_string()); + } + let mut has_transition = false; + let mut item_observed = vec![false; n_items]; + for person in 0..n_persons { + let start = row_offsets[person]; + let end = row_offsets[person + 1]; + if start >= end { + return Err("each respondent must have at least one occasion".to_string()); + } + if end - start >= 2 { + has_transition = true; + } + if time_offsets_milliseconds[start..end] + .windows(2) + .any(|window| window[1] <= window[0]) + { + return Err("time offsets must increase strictly within each respondent".to_string()); + } + for &offset in &time_offsets_milliseconds[start..end] { + days_from_millis(offset)?; + } + let mut person_observed = false; + for occasion in start..end { + for item in 0..n_items { + let value = responses[occasion * n_items + item]; + if value.is_nan() { + continue; + } + if value != 0.0 && value != 1.0 { + return Err("responses must be 0, 1, or NaN".to_string()); + } + item_observed[item] = true; + person_observed = true; + } + } + if !person_observed { + return Err("each respondent must have at least one observed response".to_string()); + } + } + if !has_transition { + return Err("at least one respondent must have two or more occasions".to_string()); + } + if item_observed.iter().any(|seen| !seen) { + return Err("each item must have at least one observed response".to_string()); + } + Ok((n_persons, n_occasions)) +} + +fn validate_config(config: HierarchicalCtarRaschConfig) -> Result { + if config.worker_count == 0 { + return Err("worker_count must be at least one".to_string()); + } + if config.max_iter == 0 { + return Err("max_iter must be at least one".to_string()); + } + if !config.tolerance.is_finite() || config.tolerance <= 0.0 { + return Err("tolerance must be finite and strictly positive".to_string()); + } + if !config.hessian_step.is_finite() || config.hessian_step <= 0.0 { + return Err("hessian_step must be finite and strictly positive".to_string()); + } + Ok(config) +} + +fn map_worker_join(joined: thread::Result) -> Result { + joined.map_err(|_| "hierarchical longitudinal worker failed".to_string()) +} + +#[derive(Clone, Debug)] +struct Unpacked { + mean: f64, + log_sd: f64, + log_decay: f64, + items: Vec, + state: Vec, +} + +fn n_hyper(n_items: usize) -> usize { + 3 + n_items +} + +fn unpack(params: &[f64], n_items: usize, n_occasions: usize) -> Result { + let expected = n_hyper(n_items) + n_occasions; + if params.len() != expected { + return Err("packed parameter length does not match the design".to_string()); + } + if params.iter().any(|value| !value.is_finite()) { + return Err("packed parameters must be finite".to_string()); + } + let mut items = params[3..3 + n_items].to_vec(); + let item_mean = items.iter().sum::() / n_items as f64; + for item in &mut items { + *item -= item_mean; + } + Ok(Unpacked { + mean: params[0], + log_sd: params[1].clamp(MIN_LOG_SD, MAX_LOG_SD), + log_decay: params[2].clamp(MIN_LOG_DECAY, MAX_LOG_DECAY), + items, + state: params[3 + n_items..].to_vec(), + }) +} + +fn pack(unpacked: &Unpacked) -> Vec { + let mut params = Vec::with_capacity(3 + unpacked.items.len() + unpacked.state.len()); + params.push(unpacked.mean); + params.push(unpacked.log_sd); + params.push(unpacked.log_decay); + params.extend_from_slice(&unpacked.items); + params.extend_from_slice(&unpacked.state); + params +} + +fn sd_from_log(log_sd: f64) -> Result { + let sd = log_sd.exp(); + if !sd.is_finite() || sd <= 0.0 { + return Err("population sd is not a finite positive value".to_string()); + } + Ok(sd) +} + +fn decay_from_log(log_decay: f64) -> Result { + let decay = log_decay.exp(); + if !decay.is_finite() || decay <= 0.0 { + return Err("decay rate is not a finite positive value".to_string()); + } + Ok(decay) +} + +#[derive(Clone, Debug)] +struct PersonNll { + nll: f64, + observed_count: usize, + transition_count: usize, + mean_grad: f64, + log_sd_grad: f64, + log_decay_grad: f64, + item_grad: Vec, + state_grad: Vec, +} + +fn person_objective( + times: &[i64], + responses: &[f64], + n_items: usize, + unpacked: &Unpacked, + state: &[f64], +) -> Result { + let sd = sd_from_log(unpacked.log_sd)?; + let variance = sd * sd; + let decay = decay_from_log(unpacked.log_decay)?; + let mut nll = 0.0; + let mut observed_count = 0; + let mut item_grad = vec![0.0; n_items]; + let mut state_grad = vec![0.0; state.len()]; + for (occasion, &theta) in state.iter().enumerate() { + for item in 0..n_items { + let value = responses[occasion * n_items + item]; + if value.is_nan() { + continue; + } + let eta = theta - unpacked.items[item]; + nll += -value * log_sigmoid(eta) - (1.0 - value) * log_sigmoid(-eta); + let residual = sigmoid_stable(eta) - value; + state_grad[occasion] += residual; + item_grad[item] -= residual; + observed_count += 1; + } + } + let first_resid = state[0] - unpacked.mean; + nll += 0.5 * first_resid * first_resid / variance + unpacked.log_sd; + state_grad[0] += first_resid / variance; + let mut mean_grad = -first_resid / variance; + let mut log_sd_grad = 1.0 - (first_resid * first_resid) / variance; + let mut log_decay_grad = 0.0; + let mut transition_count = 0; + for occasion in 1..state.len() { + let delta = days_from_millis(times[occasion])? - days_from_millis(times[occasion - 1])?; + let phi = ctar_phi(decay, delta)?; + let transition_variance = ctar_variance(variance, decay, delta)?; + let mean = unpacked.mean + phi * (state[occasion - 1] - unpacked.mean); + let resid = state[occasion] - mean; + nll += 0.5 * resid * resid / transition_variance + 0.5 * transition_variance.ln(); + let d_nll_d_theta = resid / transition_variance; + let d_nll_d_mean = -d_nll_d_theta; + let d_nll_d_var = -0.5 * resid * resid / (transition_variance * transition_variance) + + 0.5 / transition_variance; + state_grad[occasion] += d_nll_d_theta; + state_grad[occasion - 1] += d_nll_d_mean * phi; + mean_grad += d_nll_d_mean * (1.0 - phi); + let d_mean_d_phi = state[occasion - 1] - unpacked.mean; + let d_var_d_phi = -2.0 * phi * variance; + let d_nll_d_phi = d_nll_d_mean * d_mean_d_phi + d_nll_d_var * d_var_d_phi; + log_decay_grad += d_nll_d_phi * (-delta * phi * decay); + log_sd_grad += d_nll_d_var * (1.0 - phi * phi) * 2.0 * variance; + transition_count += 1; + } + if !nll.is_finite() { + return Err("hierarchical CT-AR objective is not finite".to_string()); + } + Ok(PersonNll { + nll, + observed_count, + transition_count, + mean_grad, + log_sd_grad, + log_decay_grad, + item_grad, + state_grad, + }) +} + +fn reduce_person_nll(parts: Vec, n_items: usize, n_occasions: usize) -> PersonNll { + let mut total = PersonNll { + nll: 0.0, + observed_count: 0, + transition_count: 0, + mean_grad: 0.0, + log_sd_grad: 0.0, + log_decay_grad: 0.0, + item_grad: vec![0.0; n_items], + state_grad: vec![0.0; n_occasions], + }; + let mut cursor = 0; + for part in parts { + total.nll += part.nll; + total.observed_count += part.observed_count; + total.transition_count += part.transition_count; + total.mean_grad += part.mean_grad; + total.log_sd_grad += part.log_sd_grad; + total.log_decay_grad += part.log_decay_grad; + for (dst, src) in total.item_grad.iter_mut().zip(&part.item_grad) { + *dst += src; + } + total.state_grad[cursor..cursor + part.state_grad.len()].copy_from_slice(&part.state_grad); + cursor += part.state_grad.len(); + } + let item_mean = total.item_grad.iter().sum::() / n_items as f64; + for value in &mut total.item_grad { + *value -= item_mean; + } + total +} + +fn joint_objective( + row_offsets: &[usize], + times: &[i64], + responses: &[f64], + n_items: usize, + params: &[f64], + worker_count: usize, +) -> Result<(f64, Vec, f64, usize, usize), String> { + let n_occasions = times.len(); + let unpacked = unpack(params, n_items, n_occasions)?; + let n_persons = row_offsets.len() - 1; + let workers = worker_count.min(n_persons).max(1); + let chunk = n_persons.div_ceil(workers); + let mut parts: Vec> = (0..n_persons).map(|_| None).collect(); + let joined: Result<(), String> = thread::scope(|scope| { + let mut handles = Vec::with_capacity(workers); + for worker in 0..workers { + let start = worker * chunk; + let end = (start + chunk).min(n_persons); + if start >= end { + continue; + } + let unpacked = &unpacked; + handles.push(scope.spawn(move || { + (start..end) + .map(|person| { + let occ_start = row_offsets[person]; + let occ_end = row_offsets[person + 1]; + let fit = person_objective( + ×[occ_start..occ_end], + &responses[occ_start * n_items..occ_end * n_items], + n_items, + unpacked, + &unpacked.state[occ_start..occ_end], + )?; + Ok((person, fit)) + }) + .collect::, String>>() + })); + } + for handle in handles { + let rows = map_worker_join(handle.join())??; + for (person, fit) in rows { + parts[person] = Some(fit); + } + } + Ok(()) + }); + joined?; + let ordered = collect_person_fits(parts)?; + let reduced = reduce_person_nll(ordered, n_items, n_occasions); + let mut grad = vec![0.0; params.len()]; + grad[0] = reduced.mean_grad; + grad[1] = reduced.log_sd_grad; + grad[2] = reduced.log_decay_grad; + grad[3..3 + n_items].copy_from_slice(&reduced.item_grad); + grad[3 + n_items..].copy_from_slice(&reduced.state_grad); + Ok(( + reduced.nll, + grad, + -reduced.nll, + reduced.observed_count, + reduced.transition_count, + )) +} + +fn initialize_params( + _row_offsets: &[usize], + responses: &[f64], + n_items: usize, + n_occasions: usize, +) -> Vec { + let mut item_success = vec![0.0; n_items]; + let mut item_count = vec![0.0; n_items]; + let mut state = vec![0.0; n_occasions]; + for occasion in 0..n_occasions { + let mut success = 0.0; + let mut count = 0.0; + for item in 0..n_items { + let value = responses[occasion * n_items + item]; + if value.is_nan() { + continue; + } + success += value; + count += 1.0; + item_success[item] += value; + item_count[item] += 1.0; + } + let proportion = if count == 0.0 { + 0.5 + } else { + (success + 0.5) / (count + 1.0) + }; + state[occasion] = (proportion / (1.0 - proportion)).ln(); + } + let mut items = vec![0.0; n_items]; + for item in 0..n_items { + let proportion = if item_count[item] == 0.0 { + 0.5 + } else { + (item_success[item] + 0.5) / (item_count[item] + 1.0) + }; + items[item] = -((proportion / (1.0 - proportion)).ln()); + } + let item_mean = items.iter().sum::() / n_items as f64; + for item in &mut items { + *item -= item_mean; + } + let mean = state.iter().sum::() / n_occasions.max(1) as f64; + let var = state + .iter() + .map(|value| (value - mean).powi(2)) + .sum::() + / n_occasions.max(1) as f64; + let sd = var.sqrt().clamp(0.25, 2.0); + pack(&Unpacked { + mean, + log_sd: sd.ln(), + log_decay: (0.5_f64).ln(), + items, + state, + }) +} + +#[cfg(test)] +fn empirical_state_sd(state: &[f64]) -> f64 { + if state.is_empty() { + return 0.0; + } + let mean = state.iter().sum::() / state.len() as f64; + let var = state + .iter() + .map(|value| (value - mean).powi(2)) + .sum::() + / state.len() as f64; + var.sqrt() +} + +#[cfg(test)] +fn interval_population_sd(unpacked: &Unpacked) -> Result { + let fitted = sd_from_log(unpacked.log_sd)?; + Ok(fitted.max(empirical_state_sd(&unpacked.state)).max(0.25)) +} + +fn person_state_hessian( + times: &[i64], + responses: &[f64], + n_items: usize, + unpacked: &Unpacked, + state: &[f64], +) -> Result<(Vec, Vec), String> { + let n = state.len(); + let mut diag = vec![0.0; n]; + let off = vec![0.0; n.saturating_sub(1)]; + for (occasion, &theta) in state.iter().enumerate() { + for item in 0..n_items { + let value = responses[occasion * n_items + item]; + if value.is_nan() { + continue; + } + let pi = sigmoid_stable(theta - unpacked.items[item]); + diag[occasion] += pi * (1.0 - pi); + } + } + let _times = times; + if diag.iter().all(|value| *value <= 0.0) { + return Err("measurement observed information is empty".to_string()); + } + Ok((diag, off)) +} + +fn tridiagonal_inverse_diagonal(diag: &[f64], off: &[f64]) -> Result, String> { + let n = diag.len(); + if n == 0 { + return Err("state Hessian must be non-empty".to_string()); + } + if off.len() != n.saturating_sub(1) { + return Err("state Hessian off-diagonal length is inconsistent".to_string()); + } + if diag.iter().any(|value| !value.is_finite()) || off.iter().any(|value| !value.is_finite()) { + return Err("state Hessian entries must be finite".to_string()); + } + let mut variances = vec![0.0; n]; + for column in 0..n { + let mut rhs = vec![0.0; n]; + rhs[column] = 1.0; + let solved = solve_tridiagonal(diag, off, &rhs)?; + variances[column] = solved[column]; + if !variances[column].is_finite() || variances[column] <= 0.0 { + return Err("conditional state information is not positive definite".to_string()); + } + } + Ok(variances) +} + +fn solve_tridiagonal(diag: &[f64], off: &[f64], rhs: &[f64]) -> Result, String> { + let n = diag.len(); + let mut c_prime = vec![0.0; n]; + let mut d_prime = vec![0.0; n]; + let mut denom = diag[0]; + if !denom.is_finite() || denom.abs() < 1e-12 { + return Err("conditional state information is singular".to_string()); + } + c_prime[0] = if n > 1 { off[0] / denom } else { 0.0 }; + d_prime[0] = rhs[0] / denom; + for i in 1..n { + denom = diag[i] - off[i - 1] * c_prime[i - 1]; + if !denom.is_finite() || denom.abs() < 1e-12 { + return Err("conditional state information is singular".to_string()); + } + c_prime[i] = if i + 1 < n { off[i] / denom } else { 0.0 }; + d_prime[i] = (rhs[i] - off[i - 1] * d_prime[i - 1]) / denom; + } + let mut x = vec![0.0; n]; + x[n - 1] = d_prime[n - 1]; + for i in (0..n - 1).rev() { + x[i] = d_prime[i] - c_prime[i] * x[i + 1]; + } + Ok(x) +} + +fn wald_interval(estimate: f64, se: f64) -> (f64, f64) { + (estimate - WALD_Z * se, estimate + WALD_Z * se) +} + +fn hyperparameter_hessian( + row_offsets: &[usize], + times: &[i64], + responses: &[f64], + n_items: usize, + params: &[f64], + worker_count: usize, + step: f64, +) -> Result, String> { + let n = 3; + let base = joint_objective(row_offsets, times, responses, n_items, params, worker_count)?.0; + let mut hessian = vec![0.0; n * n]; + for i in 0..n { + let mut plus = params.to_vec(); + let mut minus = params.to_vec(); + plus[i] += step; + minus[i] -= step; + let f_plus = joint_objective(row_offsets, times, responses, n_items, &plus, worker_count)?.0; + let f_minus = joint_objective(row_offsets, times, responses, n_items, &minus, worker_count)?.0; + hessian[i * n + i] = (f_plus - 2.0 * base + f_minus) / (step * step); + for j in (i + 1)..n { + let mut pp = params.to_vec(); + let mut pm = params.to_vec(); + let mut mp = params.to_vec(); + let mut mm = params.to_vec(); + pp[i] += step; + pp[j] += step; + pm[i] += step; + pm[j] -= step; + mp[i] -= step; + mp[j] += step; + mm[i] -= step; + mm[j] -= step; + let f_pp = joint_objective(row_offsets, times, responses, n_items, &pp, worker_count)?.0; + let f_pm = joint_objective(row_offsets, times, responses, n_items, &pm, worker_count)?.0; + let f_mp = joint_objective(row_offsets, times, responses, n_items, &mp, worker_count)?.0; + let f_mm = joint_objective(row_offsets, times, responses, n_items, &mm, worker_count)?.0; + let value = (f_pp - f_pm - f_mp + f_mm) / (4.0 * step * step); + hessian[i * n + j] = value; + hessian[j * n + i] = value; + } + } + Ok(hessian) +} + +fn invert_three(hessian: &[f64]) -> Result, String> { + crate::inference::vcov_from_hessian(hessian, 3, 1e-10) +} + +fn collect_person_fits(parts: Vec>) -> Result, String> { + parts + .into_iter() + .map(|part| part.ok_or_else(|| "a respondent hierarchical fit is missing".to_string())) + .collect() +} + +fn state_interval_estimates( + row_offsets: &[usize], + times: &[i64], + responses: &[f64], + n_items: usize, + unpacked: &Unpacked, +) -> (Vec, Vec, Vec, bool) { + let n_occasions = unpacked.state.len(); + let mut state_se = vec![f64::NAN; n_occasions]; + let mut identified = true; + for person in 0..(row_offsets.len() - 1) { + let start = row_offsets[person]; + let end = row_offsets[person + 1]; + match person_state_hessian( + ×[start..end], + &responses[start * n_items..end * n_items], + n_items, + unpacked, + &unpacked.state[start..end], + ) + .and_then(|(diag, off)| tridiagonal_inverse_diagonal(&diag, &off)) + { + Ok(variances) => { + for (offset, variance) in variances.into_iter().enumerate() { + state_se[start + offset] = variance.sqrt(); + } + } + Err(_) => { + identified = false; + } + } + } + let mut state_lower = vec![f64::NAN; n_occasions]; + let mut state_upper = vec![f64::NAN; n_occasions]; + if identified { + for occasion in 0..n_occasions { + let (lower, upper) = wald_interval(unpacked.state[occasion], state_se[occasion]); + state_lower[occasion] = lower; + state_upper[occasion] = upper; + } + } + (state_se, state_lower, state_upper, identified) +} + +fn hyperparameter_interval_estimates( + row_offsets: &[usize], + times: &[i64], + responses: &[f64], + n_items: usize, + params: &[f64], + unpacked: &Unpacked, + worker_count: usize, + step: f64, +) -> Result<([f64; 3], [f64; 3], [f64; 3], bool), String> { + let mut hyper_se = [f64::NAN; 3]; + let mut hyper_lower = [f64::NAN; 3]; + let mut hyper_upper = [f64::NAN; 3]; + let mut identified = false; + if let Ok(hessian) = hyperparameter_hessian( + row_offsets, + times, + responses, + n_items, + params, + worker_count, + step, + ) { + if let Ok(vcov) = invert_three(&hessian) { + let mean_se = vcov[0].sqrt(); + let log_sd_se = vcov[4].sqrt(); + let log_decay_se = vcov[8].sqrt(); + if mean_se.is_finite() && log_sd_se.is_finite() && log_decay_se.is_finite() { + let sd = sd_from_log(unpacked.log_sd)?; + let decay = decay_from_log(unpacked.log_decay)?; + hyper_se = [ + mean_se, + delta_method_sd_se(unpacked.log_sd, log_sd_se), + delta_method_decay_se(unpacked.log_decay, log_decay_se), + ]; + let mean_int = wald_interval(unpacked.mean, hyper_se[0]); + let sd_int = wald_interval(sd, hyper_se[1]); + let decay_int = wald_interval(decay, hyper_se[2]); + hyper_lower = [mean_int.0, sd_int.0, decay_int.0]; + hyper_upper = [mean_int.1, sd_int.1, decay_int.1]; + identified = true; + } + } + } + Ok((hyper_se, hyper_lower, hyper_upper, identified)) +} + +fn delta_method_sd_se(log_sd: f64, log_sd_se: f64) -> f64 { + log_sd.exp() * log_sd_se +} + +fn delta_method_decay_se(log_decay: f64, log_decay_se: f64) -> f64 { + log_decay.exp() * log_decay_se +} + +/// Fit the joint MAP hierarchical continuous-time AR(1) Rasch slice. +pub fn fit_hierarchical_ctar_rasch( + row_offsets: &[usize], + time_offsets_milliseconds: &[i64], + responses: &[f64], + n_items: usize, + config: HierarchicalCtarRaschConfig, +) -> Result { + let config = validate_config(config)?; + let (_n_persons, n_occasions) = validate_design( + row_offsets, + time_offsets_milliseconds, + responses, + n_items, + )?; + let mut params = initialize_params(row_offsets, responses, n_items, n_occasions); + let worker_count = config.worker_count; + let (fitted, _trace, _loglik, status) = lbfgs( + ¶ms, + &mut |candidate| { + let (nll, grad, loglik, _, _) = joint_objective( + row_offsets, + time_offsets_milliseconds, + responses, + n_items, + candidate, + worker_count, + )?; + Ok((nll, grad, loglik)) + }, + config.max_iter, + config.tolerance, + 12, + )?; + params = fitted; + let unpacked = unpack(¶ms, n_items, n_occasions)?; + let (_, _, _, observed_count, transition_count) = joint_objective( + row_offsets, + time_offsets_milliseconds, + responses, + n_items, + ¶ms, + worker_count, + )?; + let (state_se, state_lower, state_upper, state_identified) = state_interval_estimates( + row_offsets, + time_offsets_milliseconds, + responses, + n_items, + &unpacked, + ); + let (hyper_se, hyper_lower, hyper_upper, hyper_identified) = hyperparameter_interval_estimates( + row_offsets, + time_offsets_milliseconds, + responses, + n_items, + ¶ms, + &unpacked, + worker_count, + config.hessian_step, + )?; + let decay = decay_from_log(unpacked.log_decay)?; + Ok(HierarchicalCtarRaschFit { + state: unpacked.state, + state_se, + state_lower, + state_upper, + item_intercepts: unpacked.items, + population_mean: unpacked.mean, + population_sd: sd_from_log(unpacked.log_sd)?, + decay_rate: decay, + unit_time_ar_coefficient: ctar_unit_phi(decay)?, + hyperparameter_se: hyper_se, + hyperparameter_lower: hyper_lower, + hyperparameter_upper: hyper_upper, + hyperparameter_intervals_identified: hyper_identified, + state_intervals_identified: state_identified, + observed_count, + transition_count, + status, + estimand_scope: ESTIMAND_SCOPE, + transition_kind: TRANSITION_KIND, + interval_kind: INTERVAL_KIND, + engine: ENGINE, + }) +} + +/// Deterministic simulator for hierarchical CT-AR Rasch recovery fixtures. +pub fn simulate_hierarchical_ctar_rasch( + row_offsets: &[usize], + time_offsets_milliseconds: &[i64], + n_items: usize, + population_mean: f64, + population_sd: f64, + decay_rate: f64, + item_intercepts: &[f64], + seed: u64, +) -> Result<(Vec, Vec), String> { + if item_intercepts.len() != n_items { + return Err("item_intercepts length must equal n_items".to_string()); + } + if !population_mean.is_finite() || !population_sd.is_finite() || population_sd <= 0.0 { + return Err("simulator population parameters must be finite with positive sd".to_string()); + } + let n_occasions = time_offsets_milliseconds.len(); + let n_persons = validate_offsets(row_offsets, n_occasions)?; + let mut rng = Lcg(seed | 1); + let mut state = vec![0.0; n_occasions]; + let mut responses = vec![f64::NAN; n_occasions * n_items]; + for person in 0..n_persons { + let start = row_offsets[person]; + let end = row_offsets[person + 1]; + if start >= end { + return Err("each respondent must have at least one occasion".to_string()); + } + state[start] = population_mean + population_sd * rng.standard_normal(); + for occasion in (start + 1)..end { + let delta = days_from_millis(time_offsets_milliseconds[occasion])? + - days_from_millis(time_offsets_milliseconds[occasion - 1])?; + let phi = ctar_phi(decay_rate, delta)?; + let variance = ctar_variance(population_sd * population_sd, decay_rate, delta)?; + let mean = population_mean + phi * (state[occasion - 1] - population_mean); + state[occasion] = mean + variance.sqrt() * rng.standard_normal(); + } + for occasion in start..end { + for item in 0..n_items { + let eta = state[occasion] - item_intercepts[item]; + let draw = if rng.next_f64() < sigmoid_stable(eta) { + 1.0 + } else { + 0.0 + }; + responses[occasion * n_items + item] = draw; + } + } + } + Ok((state, responses)) +} + +struct Lcg(u64); + +impl Lcg { + fn next_f64(&mut self) -> f64 { + self.0 = self + .0 + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + ((self.0 >> 11) as f64) / ((1u64 << 53) as f64) + } + + fn standard_normal(&mut self) -> f64 { + let u1 = self.next_f64().max(1e-12); + let u2 = self.next_f64(); + (-2.0 * u1.ln()).sqrt() * (2.0 * PI * u2).cos() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn regular_offsets(n_persons: usize, n_occasions_each: usize) -> (Vec, Vec) { + let mut offsets = vec![0]; + let mut times = Vec::new(); + for _ in 0..n_persons { + for occasion in 0..n_occasions_each { + let days = occasion as f64 + if occasion == 2 { 0.5 } else { 0.0 }; + times.push((days * MILLIS_PER_DAY) as i64); + } + offsets.push(times.len()); + } + (offsets, times) + } + + #[test] + fn ctar_helpers_match_the_ou_transition_and_reject_invalid_inputs() { + let phi = ctar_phi(0.5, 2.0).unwrap(); + assert!((phi - (-1.0_f64).exp()).abs() < 1e-12); + let variance = ctar_variance(1.0, 0.5, 2.0).unwrap(); + assert!((variance - (1.0 - phi * phi)).abs() < 1e-12); + let tiny = ctar_variance(1.0, 1e-6, 1e-6).unwrap(); + assert!(tiny > 0.0 && tiny.is_finite()); + assert!((ctar_unit_phi(0.4).unwrap() - (-0.4_f64).exp()).abs() < 1e-12); + assert!(ctar_phi(0.0, 1.0).unwrap_err().contains("strictly positive")); + assert!(ctar_phi(0.5, 0.0).unwrap_err().contains("elapsed days")); + assert!(ctar_variance(0.0, 0.5, 1.0) + .unwrap_err() + .contains("population variance")); + assert_eq!(ctar_phi(1e9, 1e9).unwrap(), 0.0); + let floored = ctar_variance(1.0, 1e-16, 1e-16).unwrap(); + assert!((floored - MIN_TRANSITION_VARIANCE).abs() < 1e-18); + assert!(ctar_variance(f64::INFINITY, 0.5, 1.0) + .unwrap_err() + .contains("population variance")); + assert!(ctar_variance(f64::NAN, 0.5, 1.0) + .unwrap_err() + .contains("population variance")); + } + + #[test] + fn measurement_and_transition_gradients_match_finite_differences() { + let times = [0, 86_400_000, 172_800_000]; + let responses = [ + 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, + ]; + let unpacked = Unpacked { + mean: 0.1, + log_sd: 0.0, + log_decay: (0.4_f64).ln(), + items: vec![-0.2, 0.1, 0.1], + state: vec![0.3, -0.1, 0.2], + }; + let analytic = person_objective(×, &responses, 3, &unpacked, &unpacked.state).unwrap(); + let step = 1e-5; + let mut plus = unpacked.clone(); + plus.state[1] += step; + let mut minus = unpacked.clone(); + minus.state[1] -= step; + let f_plus = person_objective(×, &responses, 3, &plus, &plus.state) + .unwrap() + .nll; + let f_minus = person_objective(×, &responses, 3, &minus, &minus.state) + .unwrap() + .nll; + let numeric = (f_plus - f_minus) / (2.0 * step); + assert!((analytic.state_grad[1] - numeric).abs() < 1e-6, "{numeric}"); + let mut plus_mean = unpacked.clone(); + plus_mean.mean += step; + let mut minus_mean = unpacked.clone(); + minus_mean.mean -= step; + let mean_numeric = (person_objective(×, &responses, 3, &plus_mean, &plus_mean.state) + .unwrap() + .nll + - person_objective(×, &responses, 3, &minus_mean, &minus_mean.state) + .unwrap() + .nll) + / (2.0 * step); + assert!((analytic.mean_grad - mean_numeric).abs() < 1e-6); + let mut plus_sd = unpacked.clone(); + plus_sd.log_sd += step; + let mut minus_sd = unpacked.clone(); + minus_sd.log_sd -= step; + let sd_numeric = (person_objective(×, &responses, 3, &plus_sd, &plus_sd.state) + .unwrap() + .nll + - person_objective(×, &responses, 3, &minus_sd, &minus_sd.state) + .unwrap() + .nll) + / (2.0 * step); + assert!((analytic.log_sd_grad - sd_numeric).abs() < 1e-5); + let mut plus_decay = unpacked.clone(); + plus_decay.log_decay += step; + let mut minus_decay = unpacked.clone(); + minus_decay.log_decay -= step; + let decay_numeric = ( + person_objective(×, &responses, 3, &plus_decay, &plus_decay.state) + .unwrap() + .nll + - person_objective(×, &responses, 3, &minus_decay, &minus_decay.state) + .unwrap() + .nll + ) / (2.0 * step); + assert!((analytic.log_decay_grad - decay_numeric).abs() < 1e-5); + } + + #[test] + fn packed_centering_and_unpack_errors_are_stable() { + let unpacked = Unpacked { + mean: 0.0, + log_sd: 0.0, + log_decay: 0.0, + items: vec![1.0, -1.0], + state: vec![0.2, -0.1], + }; + let packed = pack(&unpacked); + let replayed = unpack(&packed, 2, 2).unwrap(); + assert_eq!(replayed.items, vec![1.0, -1.0]); + assert!(unpack(&[0.0; 3], 2, 2) + .unwrap_err() + .contains("packed parameter length")); + assert!(unpack(&[0.0, 0.0, 0.0, f64::NAN, 0.0, 0.0, 0.0], 2, 2) + .unwrap_err() + .contains("finite")); + assert_eq!(n_hyper(4), 7); + assert!(sd_from_log(1e9).unwrap_err().contains("population sd")); + assert!(decay_from_log(1e9).unwrap_err().contains("decay rate")); + } + + #[test] + fn rejects_invalid_designs_and_config_without_panicking() { + let offsets = [0, 2]; + let times = [0, 86_400_000]; + let responses = [1.0, 0.0, 0.0, 1.0]; + assert!(fit_hierarchical_ctar_rasch( + &offsets, + ×, + &responses, + 1, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("at least two items")); + assert!(fit_hierarchical_ctar_rasch( + &[1, 2], + ×, + &responses, + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("start at zero")); + assert!(validate_offsets(&[0, 1, 0], 0) + .unwrap_err() + .contains("non-decreasing")); + assert!(validate_offsets(&[0, 1], 2) + .unwrap_err() + .contains("end at the occasion count")); + assert!(fit_hierarchical_ctar_rasch( + &[0], + &[], + &[], + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("at least one respondent")); + assert!(fit_hierarchical_ctar_rasch( + &[0, 0], + &[], + &[], + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("at least one occasion")); + assert!(fit_hierarchical_ctar_rasch( + &[0, 1], + &[0], + &[1.0, 0.0], + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("two or more occasions")); + assert!(fit_hierarchical_ctar_rasch( + &offsets, + &[2, 1], + &responses, + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("increase strictly")); + assert!(fit_hierarchical_ctar_rasch( + &offsets, + ×, + &[1.0, 0.0, 2.0, 0.0], + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("0, 1, or NaN")); + assert!(fit_hierarchical_ctar_rasch( + &offsets, + ×, + &[1.0, f64::NAN, 0.0, f64::NAN], + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("each item")); + assert!(fit_hierarchical_ctar_rasch( + &[0, 2, 4], + &[0, 1, 0, 1], + &[1.0, 0.0, 0.0, 1.0, f64::NAN, f64::NAN, f64::NAN, f64::NAN], + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("at least one observed")); + assert!(fit_hierarchical_ctar_rasch( + &offsets, + ×, + &[1.0, 0.0], + 2, + HierarchicalCtarRaschConfig::default() + ) + .unwrap_err() + .contains("occasion-major")); + let mut bad = HierarchicalCtarRaschConfig::default(); + bad.worker_count = 0; + assert!(validate_config(bad).unwrap_err().contains("worker_count")); + bad = HierarchicalCtarRaschConfig::default(); + bad.max_iter = 0; + assert!(validate_config(bad).unwrap_err().contains("max_iter")); + bad = HierarchicalCtarRaschConfig::default(); + bad.tolerance = 0.0; + assert!(validate_config(bad).unwrap_err().contains("tolerance")); + bad = HierarchicalCtarRaschConfig::default(); + bad.hessian_step = -1.0; + assert!(validate_config(bad).unwrap_err().contains("hessian_step")); + assert!(days_from_millis(i64::MAX).unwrap_err().contains("supported finite range")); + } + + #[test] + fn worker_count_does_not_change_the_joint_map() { + let (offsets, times) = regular_offsets(4, 3); + let items = [-0.6, -0.2, 0.2, 0.6]; + let (_state, responses) = simulate_hierarchical_ctar_rasch( + &offsets, + ×, + 4, + 0.0, + 0.6, + 0.4, + &items, + 17, + ) + .unwrap(); + let mut config = HierarchicalCtarRaschConfig { + worker_count: 1, + max_iter: 80, + tolerance: 1e-4, + hessian_step: 1e-3, + }; + let one = fit_hierarchical_ctar_rasch(&offsets, ×, &responses, 4, config).unwrap(); + config.worker_count = 3; + let many = fit_hierarchical_ctar_rasch(&offsets, ×, &responses, 4, config).unwrap(); + for (left, right) in one.state.iter().zip(&many.state) { + assert!((left - right).abs() < 1e-8, "{left} vs {right}"); + } + assert_eq!(one.estimand_scope, ESTIMAND_SCOPE); + assert_eq!(one.engine, ENGINE); + assert_eq!(one.transition_kind, TRANSITION_KIND); + assert_eq!(one.interval_kind, INTERVAL_KIND); + } + + #[test] + fn unused_worker_shard_is_skipped_and_join_errors_are_package_owned() { + let (offsets, times) = regular_offsets(2, 2); + let items = [-0.4, 0.4]; + let (_state, responses) = simulate_hierarchical_ctar_rasch( + &offsets, + ×, + 2, + 0.0, + 0.5, + 0.5, + &items, + 3, + ) + .unwrap(); + let config = HierarchicalCtarRaschConfig { + worker_count: 8, + max_iter: 40, + tolerance: 1e-4, + hessian_step: 1e-3, + }; + let fit = fit_hierarchical_ctar_rasch(&offsets, ×, &responses, 2, config).unwrap(); + assert_eq!(fit.state.len(), 4); + assert_eq!( + map_worker_join::<()>(Err(Box::new("boom"))).unwrap_err(), + "hierarchical longitudinal worker failed" + ); + } + + #[test] + fn recovers_true_states_and_transition_parameters_across_seeds() { + let n_persons = 16; + let n_occasions = 8; + let n_items = 8; + let (offsets, times) = regular_offsets(n_persons, n_occasions); + let items: Vec = (0..n_items) + .map(|item| (item as f64 - 3.5) * 0.25) + .collect(); + let true_mean = 0.0; + let true_sd = 0.7; + let true_decay = 0.35; + let seeds = [11_u64, 23, 41, 59, 73]; + let mut state_sse = 0.0; + let mut state_count = 0.0; + let mut covered = 0.0; + let mut mean_err = 0.0; + let mut sd_err = 0.0; + let mut decay_values = Vec::new(); + for seed in seeds { + let (true_state, responses) = simulate_hierarchical_ctar_rasch( + &offsets, + ×, + n_items, + true_mean, + true_sd, + true_decay, + &items, + seed, + ) + .unwrap(); + let config = HierarchicalCtarRaschConfig { + worker_count: 4, + max_iter: 200, + tolerance: 1e-5, + hessian_step: 1e-3, + }; + let fit = fit_hierarchical_ctar_rasch(&offsets, ×, &responses, n_items, config) + .unwrap(); + assert_eq!(fit.observed_count, n_persons * n_occasions * n_items); + assert_eq!(fit.transition_count, n_persons * (n_occasions - 1)); + assert!(fit.state_intervals_identified); + for (estimate, truth, lower, upper) in fit + .state + .iter() + .zip(&true_state) + .zip(&fit.state_lower) + .zip(&fit.state_upper) + .map(|(((estimate, truth), lower), upper)| (estimate, truth, lower, upper)) + { + state_sse += (estimate - truth).powi(2); + state_count += 1.0; + if *lower <= *truth && *truth <= *upper { + covered += 1.0; + } + } + mean_err += (fit.population_mean - true_mean).powi(2); + sd_err += (fit.population_sd - true_sd).powi(2); + assert!(fit.decay_rate.is_finite() && fit.decay_rate > 0.0); + assert!( + fit.unit_time_ar_coefficient > 0.0 && fit.unit_time_ar_coefficient < 1.0 + ); + decay_values.push(fit.decay_rate); + } + let state_rmse = (state_sse / state_count).sqrt(); + let coverage = covered / state_count; + let mean_rmse = (mean_err / seeds.len() as f64).sqrt(); + let sd_rmse = (sd_err / seeds.len() as f64).sqrt(); + let decay_rmse = (decay_values + .iter() + .map(|value| (value - true_decay).powi(2)) + .sum::() + / seeds.len() as f64) + .sqrt(); + assert!(state_rmse < 0.85, "state RMSE {state_rmse}"); + assert!(coverage > 0.80, "state coverage {coverage}"); + assert!(mean_rmse < 0.35, "mean RMSE {mean_rmse}"); + // Joint MAP shrinks tau; this bound is not an unbiased-ML claim. + assert!(sd_rmse < 0.70, "sd RMSE {sd_rmse}"); + // Short irregular series leave lambda weakly identified under joint MAP. + // The recovery claim is a finite positive decay and a unit-day phi in + // (0, 1), not a tight RMSE for lambda itself. + assert!(decay_rmse.is_finite(), "decay RMSE {decay_rmse}"); + } + + #[test] + fn missing_responses_are_excluded_and_irregular_gaps_change_phi() { + let offsets = [0, 3]; + let times = [0, 86_400_000, 3 * 86_400_000]; + let responses = [ + 1.0, 0.0, 1.0, 0.0, f64::NAN, f64::NAN, 0.0, 1.0, 0.0, 1.0, f64::NAN, f64::NAN, + ]; + let config = HierarchicalCtarRaschConfig { + worker_count: 1, + max_iter: 80, + tolerance: 1e-4, + hessian_step: 1e-3, + }; + let fit = fit_hierarchical_ctar_rasch(&offsets, ×, &responses, 4, config).unwrap(); + assert_eq!(fit.observed_count, 8); + assert_eq!(fit.transition_count, 2); + let one_day = ctar_phi(fit.decay_rate, 1.0).unwrap(); + let two_day = ctar_phi(fit.decay_rate, 2.0).unwrap(); + assert!(two_day < one_day); + assert_eq!(fit.unit_time_ar_coefficient, one_day); + } + + #[test] + fn tridiagonal_and_hessian_helpers_cover_failure_branches() { + assert!(tridiagonal_inverse_diagonal(&[], &[]) + .unwrap_err() + .contains("non-empty")); + assert!(tridiagonal_inverse_diagonal(&[1.0, 1.0], &[1.0, 1.0]) + .unwrap_err() + .contains("off-diagonal")); + assert!(tridiagonal_inverse_diagonal(&[f64::NAN], &[]) + .unwrap_err() + .contains("finite")); + assert!(solve_tridiagonal(&[0.0], &[], &[1.0]) + .unwrap_err() + .contains("singular")); + assert!(solve_tridiagonal(&[1.0, 0.0], &[0.0], &[1.0, 1.0]) + .unwrap_err() + .contains("singular")); + let solved = solve_tridiagonal(&[2.0, 2.0], &[-1.0], &[1.0, 0.0]).unwrap(); + assert!((solved[0] - 2.0 / 3.0).abs() < 1e-12); + let (lower, upper) = wald_interval(0.0, 1.0); + assert!((upper - lower - 2.0 * WALD_Z).abs() < 1e-12); + assert!( + simulate_hierarchical_ctar_rasch(&[0, 1], &[0], 2, 0.0, 0.5, 0.4, &[0.0], 1) + .unwrap_err() + .contains("item_intercepts") + ); + assert!( + simulate_hierarchical_ctar_rasch(&[0, 1], &[0], 2, 0.0, 0.0, 0.4, &[0.0, 0.0], 1) + .unwrap_err() + .contains("population parameters") + ); + assert!( + simulate_hierarchical_ctar_rasch(&[0, 0], &[], 2, 0.0, 0.5, 0.4, &[0.0, 0.0], 1) + .unwrap_err() + .contains("at least one occasion") + ); + assert_eq!(empirical_state_sd(&[]), 0.0); + assert!(empirical_state_sd(&[1.0, -1.0]) > 0.0); + let scale = interval_population_sd(&Unpacked { + mean: 0.0, + log_sd: (0.2_f64).ln(), + log_decay: 0.0, + items: vec![0.0, 0.0], + state: vec![1.0, -1.0], + }) + .unwrap(); + assert!(scale >= 0.25); + assert!(person_state_hessian( + &[0, 86_400_000], + &[f64::NAN, f64::NAN, f64::NAN, f64::NAN], + 2, + &Unpacked { + mean: 0.0, + log_sd: 0.0, + log_decay: 0.0, + items: vec![0.0, 0.0], + state: vec![0.0, 0.0], + }, + &[0.0, 0.0], + ) + .unwrap_err() + .contains("measurement observed information")); + assert!(delta_method_sd_se(0.0, 0.2) > 0.0); + assert!(delta_method_decay_se(0.0, 0.2) > 0.0); + let identity = invert_three(&[1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]).unwrap(); + assert!((identity[0] - 1.0).abs() < 1e-12); + } + + #[test] + fn default_config_and_metadata_are_normative() { + let config = HierarchicalCtarRaschConfig::default(); + assert_eq!(config.worker_count, 1); + assert_eq!(config.max_iter, 250); + assert_eq!(config.tolerance, 1e-5); + assert_eq!(config.hessian_step, 1e-3); + assert_eq!(ESTIMAND_SCOPE, "joint_map_hierarchical_ctar_rasch"); + assert_ne!(ESTIMAND_SCOPE, "independent_respondent_ols_trend"); + assert_ne!(ESTIMAND_SCOPE, "discrete_ar_state_prediction"); + } + + #[test] + fn interval_helpers_and_defensive_branches_are_covered() { + let offsets = [0, 2]; + let times = [0, 86_400_000]; + let responses = [1.0, 0.0, 0.0, 1.0]; + let good = Unpacked { + mean: 0.0, + log_sd: 0.0, + log_decay: (0.4_f64).ln(), + items: vec![-0.1, 0.1], + state: vec![0.2, -0.1], + }; + let (se, lower, upper, identified) = + state_interval_estimates(&offsets, ×, &responses, 2, &good); + assert!(identified); + assert!(se.iter().all(|value| value.is_finite() && *value > 0.0)); + assert!(lower[0] < upper[0]); + let missing = [f64::NAN, f64::NAN, f64::NAN, f64::NAN]; + let (_se, lower, upper, identified) = + state_interval_estimates(&offsets, ×, &missing, 2, &good); + assert!(!identified); + assert!(lower.iter().all(|value| value.is_nan())); + assert!(upper.iter().all(|value| value.is_nan())); + let packed = pack(&good); + let (hyper_se, _, _, hyper_ok) = hyperparameter_interval_estimates( + &offsets, + ×, + &responses, + 2, + &packed, + &good, + 1, + 1e-3, + ) + .unwrap(); + if hyper_ok { + assert!(hyper_se.iter().all(|value| value.is_finite() && *value > 0.0)); + } + let (hyper_se, _, _, hyper_ok) = hyperparameter_interval_estimates( + &offsets, + ×, + &responses, + 2, + &packed, + &good, + 1, + 50.0, + ) + .unwrap(); + assert!(!hyper_ok); + assert!(hyper_se.iter().all(|value| value.is_nan())); + assert!(collect_person_fits(vec![None]) + .unwrap_err() + .contains("hierarchical fit is missing")); + assert!(joint_objective(&offsets, ×, &responses, 2, &[0.0; 3], 1) + .unwrap_err() + .contains("packed parameter length")); + let mut exploding = good.clone(); + exploding.state = vec![1e200, -1e200]; + assert!(person_objective(×, &responses, 2, &exploding, &exploding.state) + .unwrap_err() + .contains("not finite")); + assert!(invert_three(&[f64::NAN; 9]).is_err()); + assert!(validate_design( + &[0, 2], + &[0, 1], + &[], + (usize::MAX / 2) + 1 + ) + .unwrap_err() + .contains("exceeds supported size")); + let init = initialize_params( + &[0, 2], + &[f64::NAN, f64::NAN, 1.0, 0.0], + 2, + 2, + ); + assert!(init.iter().all(|value| value.is_finite())); + let empty_item_init = initialize_params(&[0, 1], &[f64::NAN, f64::NAN], 2, 1); + assert!(empty_item_init.iter().all(|value| value.is_finite())); + let all_zero = initialize_params(&[0, 1], &[0.0, 0.0], 2, 1); + assert!(all_zero.iter().all(|value| value.is_finite())); + assert!(hyperparameter_hessian( + &offsets, + ×, + &responses, + 2, + &[0.0; 3], + 1, + 1e-3 + ) + .unwrap_err() + .contains("packed parameter length")); + let mut negative_diag = good.clone(); + negative_diag.log_sd = 0.0; + let (diag, off) = + person_state_hessian(×, &responses, 2, &negative_diag, &negative_diag.state) + .unwrap(); + let mut broken = diag; + broken[0] = -1.0; + assert!(tridiagonal_inverse_diagonal(&broken, &off) + .unwrap_err() + .contains("positive definite")); + } +} diff --git a/docs/PRD.md b/docs/PRD.md index cb62e08b5..834d56b88 100644 --- a/docs/PRD.md +++ b/docs/PRD.md @@ -162,7 +162,7 @@ Open PRs and issues may contain additional capabilities. They are **not** consid **PRD-FR-060** Reusable contracts shall represent explicit context dimensions, context identities, membership weights, repeated occasions, and temporal state specifications without inferring random-effect families from labels. -**PRD-FR-061** Multiple-membership weights and temporal ordering shall be provenance-bound. Elapsed-time effects shall not be claimed unless the fitted model actually parameterizes elapsed-time transitions. +**PRD-FR-061** Multiple-membership weights and temporal ordering shall be provenance-bound. Elapsed-time effects shall not be claimed unless the fitted model actually parameterizes elapsed-time transitions. The joint MAP hierarchical CT-AR Rasch slice parameterizes elapsed days; the independent OLS and caller-supplied discrete AR layer does not. **PRD-FR-062** Future Rust estimators for these contracts shall establish identification and true-parameter recovery before release as production estimators. diff --git a/docs/TRD.md b/docs/TRD.md index 40f7acdb0..5949ab414 100644 --- a/docs/TRD.md +++ b/docs/TRD.md @@ -182,9 +182,9 @@ docs/ PRD/TRD, method docs, ADRs, doctoring, diagrams, **TRD-MLT-003** Cross-classified designs shall maintain dimension-qualified identities. -**TRD-MLT-004** Temporal occasions shall retain ordering/time provenance and explicitly separate discrete occasion-step AR effects from continuous-time parameterizations. +**TRD-MLT-004** Temporal occasions shall retain ordering/time provenance and explicitly separate discrete occasion-step AR effects from continuous-time parameterizations. The hierarchical CT-AR Rasch slice uses elapsed-day Ornstein–Uhlenbeck transitions; the OLS/AR state layer does not. -**TRD-MLT-005** Numerical multilevel/longitudinal estimators shall remain proposed until Rust implementations pass identification and true-parameter recovery studies. +**TRD-MLT-005** Numerical multilevel/longitudinal estimators shall remain proposed until Rust implementations pass identification and true-parameter recovery studies. The hierarchical CT-AR Rasch slice reports multi-seed state RMSE/coverage as joint MAP evidence, not as Fox–Glas Gibbs or estimated MMMC recovery. ### 4.12 Testing and scientific evidence diff --git a/docs/adr/0007-multilevel-multiple-membership-temporal.md b/docs/adr/0007-multilevel-multiple-membership-temporal.md index 6446900a2..63d525645 100644 --- a/docs/adr/0007-multilevel-multiple-membership-temporal.md +++ b/docs/adr/0007-multilevel-multiple-membership-temporal.md @@ -8,10 +8,11 @@ Date: 2026-08-09 Psychometric and AI-evaluation observations commonly sit inside schools, teams, organizations, prompts, testlets, documents, clients, time periods or other overlapping contexts. Repeated observations also evolve over time. Flattening those structures into independent rows can produce atomistic fallacy, understate uncertainty, confound stable traits with context effects and drift, and misinterpret temporal dependence. Reusable contracts plus a Rust-owned respondent state layer define a supported -longitudinal handoff. The ADR remains Proposed until that boundary is -protected-integrated. Full multilevel IRT random-effect estimation, uncertainty, -joint likelihood recovery, and GPU recurrent-state parity are not accepted -production behavior. +longitudinal handoff. ADR-0019 adds a separate joint MAP hierarchical +continuous-time AR(1) Rasch slice on that handoff. This ADR remains Proposed +until those boundaries are protected-integrated. Fox and Glas Gibbs sampling, +Jeon and Rabe-Hesketh adaptive-quadrature ML, estimated multiple-membership +`u_h`, and GPU recurrent-state parity are not accepted production behavior. ## Decision @@ -27,14 +28,19 @@ The architecture treats the following as distinct, explicit structures: - future continuous-time state transitions; - rater/model/prompt drift. -The accepted state-layer boundary is deliberately narrower: independent +The ADR-0018 state-layer boundary remains deliberately narrower: independent per-respondent OLS trends are fitted by Rust on exact day-scaled offsets, and stationary AR(1) states produce discrete-sequence predictions from a caller-supplied coefficient. The latter uses sequence gaps, not elapsed milliseconds, so irregular calendar spacing cannot be silently treated as a -continuous-time decay. These predictors do not estimate population +continuous-time decay. Those predictors do not estimate population random-effects distributions or AR-coefficient uncertainty. +ADR-0019 is a separate joint MAP slice. It estimates shared +`(mu, tau, lambda)`, shrinks person-occasion states toward `mu`, and uses +elapsed days in an Ornstein–Uhlenbeck / continuous-time AR(1) transition. +It does not estimate crossed or multiple-membership `u_h`. + ### Contract rules 1. Context membership names an explicit `context_dimension_id` and dimension-scoped `context_id`. diff --git a/docs/adr/0019-joint-hierarchical-ctar-rasch.md b/docs/adr/0019-joint-hierarchical-ctar-rasch.md new file mode 100644 index 000000000..ddd5178cc --- /dev/null +++ b/docs/adr/0019-joint-hierarchical-ctar-rasch.md @@ -0,0 +1,212 @@ +# ADR-0019: Joint MAP hierarchical continuous-time AR(1) Rasch slice + +Status: **Proposed** +Date: 2026-08-17 +Supersedes: none +Superseded by: none + +## Context + +ADR-0018 on the live longitudinal state-engine branch (`#976`) owns an +independent per-respondent OLS trend and a caller-supplied discrete AR +predictor. Those kernels are scientifically valid for their stated estimands +and must not be relabeled as population random effects, estimated +autoregression, interval coverage, or continuous-time transitions. + +The next product gap is the smallest jointly estimated longitudinal +latent-state IRT slice that *does* estimate shared population hyperparameters, +uses elapsed time in the transition, and reports uncertainty. Open PRs `#978`, +`#979`, and `#981` do not implement that gap. Closed `#848` is the same OLS/AR +tree as `#976`. Combining estimated multiple-membership context effects +`u_h` with person-occasion states is a larger identification problem and is +not this slice. + +Journal PDFs for the governing methods are copyrighted. This record cites, +links, and summarizes them rather than attaching the PDFs. + +## Decision drivers + +- Honest estimand labels: do not call OLS/AR random effects if they are not. +- Hierarchical pooling requires shared population parameters and shrinkage. +- Irregular occasion spacing requires a continuous-time or explicitly + time-scaled transition, not a discrete `phi` reused across unequal gaps. +- Uncertainty must be computed and interval coverage tested against known + states, without inventing 95% hyperparameter coverage from five seeds. +- Multiple-membership / crossed grouping is compatible only if the joint + likelihood actually estimates those effects; otherwise exclude it in + contract language. +- GPU parity is honest only when an existing GPU abstraction implements the + same estimand. + +## Ownership and dependency direction + +`ContextualWisdomLab/fast-mlsirm` owns this reusable measurement kernel. +Psychometrics Commons and other hosted products consume the sealed design and +returned diagnostics. This slice depends on the `#976` longitudinal design +handoff and must not recreate hosted session, consent, persistence, or HTTP +runtime. `kaefa` and `aFIPC` are not oracles. + +## Decision + +Add a separate Rust-owned joint MAP estimator behind a new Python entry point +`fit_hierarchical_longitudinal_irt`: + +```text +logit P(Y_pti = 1) = theta_pt - b_i, sum_i b_i = 0 +theta_p,1 ~ N(mu, tau^2) +theta_p,t | theta_p,t-1 ~ N( + mu + exp(-lambda * Delta_pt) * (theta_p,t-1 - mu), + tau^2 * (1 - exp(-2 * lambda * Delta_pt)) +) +``` + +`Delta_pt` is elapsed time in days from exact millisecond offsets. Packed +parameters are `[mu, log tau, log lambda, b, theta]`. Optimization is the +existing Rust L-BFGS path. Person shards run on scoped CPU threads and reduce +in respondent order. + +Normative metadata: + +- `estimand_scope = "joint_map_hierarchical_ctar_rasch"` +- `transition_kind = "continuous_time_ar1_ou"` +- `interval_kind = "wald_measurement_observed_information"` +- `engine = "rust_cpu_multithreaded"` +- `population_random_effects_estimated = True` +- `ar_coefficient_estimated = True` +- `ar_coefficient_source = "joint_map"` +- `multiple_membership_estimated = False` +- `gpu_parity = False` + +State standard errors use the person-block measurement observed information +only. The hierarchical prior regularizes the joint MAP point estimates; it is +not treated as known truth when forming Wald state intervals. +Hyperparameter standard errors use a 3×3 finite-difference Hessian on +`(mu, log tau, log lambda)` plus the delta method for `tau` and `lambda`. +Intervals are Wald 95% intervals. Identification flags are returned rather +than invented when the Hessian is not usable. + +The `#976` OLS/AR entry point and its metadata remain unchanged. + +## Invariants / acceptance evidence + +1. Worker counts do not change the joint MAP state vector. +2. Multi-seed recovery against known generating states and transition + parameters reports RMSE/bias and state-interval coverage inside + documented MAP bounds. Those bounds are not a claim of unbiased ML or + exact 95% hyperparameter coverage. +3. Missing responses are excluded; irregular gaps change `exp(-lambda Delta)`. +4. Invalid designs, non-binary responses, and invalid optimizer controls fail + closed with package-owned errors. +5. Metadata never uses OLS or caller-supplied-AR estimand labels. + +## Non-goals and claims not made + +- Not independent respondent OLS and not caller-supplied discrete AR. +- Not Fox and Glas (2001) Gibbs sampling. +- Not Jeon and Rabe-Hesketh (2016) adaptive-quadrature ML. +- Not estimated multiple-membership or crossed `u_h` in this joint likelihood. +- Not GPU parity. The existing wgpu path owns MLSIRM distance/likelihood + kernels in f32, a different estimand. +- Not a claim that five recovery seeds establish frequentist 95% coverage + for `(mu, tau, lambda)`. + +## Consequences and trade-offs + +### Benefits + +- Smallest scientifically valid joint longitudinal IRT slice on top of `#976`. +- Elapsed-time transitions and hierarchical shrinkage are named honestly. +- Python remains marshalling-only. + +### Costs / risks + +- Joint MAP shrinks `tau` and person states toward `mu`. +- Wald observed-information intervals are local and can be unidentified. +- Crossed / multiple-membership structure still requires a later joint model. + +## Alternatives considered + +### Relabel the ADR-0018 OLS/AR kernels as random effects + +Rejected. Those kernels do not estimate a population distribution or `phi`. + +### Discrete time-scaled AR with caller `phi` + +Rejected as the next slice. It would still leave the coefficient +caller-supplied. + +### Full MMMC + longitudinal joint likelihood + +Deferred. Estimating context-level `u_h` together with CT-AR person states +needs a separate identification and recovery design (Fox & Glas; Browne, +Goldstein, & Rasbash). + +### GPU kernel for this objective + +Rejected. The current GPU abstraction does not implement this estimand. + +## Failure, degraded, and recovery behavior + +Invalid designs, non-finite packed parameters, degenerate OU transitions, and +worker join failures return package-owned errors. Unidentified Hessians set +interval flags to false and leave the corresponding bounds as NaN. Rollback +is removal of the new module and entry point; the `#976` OLS/AR layer remains. + +## Security and privacy implications + +No new credentials, network access, or raw source retention. Response +matrices are bounded before native allocation. Hostile Python mappings and +Boolean-as-number coercion fail closed. + +## Compatibility, migration, and rollback + +New public names only. Existing `fit_longitudinal_state` metadata is +unchanged. The sealed `LongitudinalDesign` is reused for occasion identity +and is not reinterpreted as this IRT estimand. Stacked on `#976`; do not +merge to `main` without that dependency. + +## Verification and release evidence + +- Rust unit tests: OU helpers, analytic vs finite-difference gradients, + fail-closed validation, worker determinism, unused-shard join errors, + missing/irregular time, tridiagonal/Hessian failures, multi-seed recovery. +- Python tests: public marshalling, recovery, metadata honesty, binding + bounds, and simulator controls. +- Cargo, pytest, package, existing fuzz, and existing GPU-smoke evidence. + This slice adds no GPU kernel, so GPU-smoke remains the MLSIRM path. + +## Research and standards basis + +Fox, J.-P., & Glas, C. A. W. (2001). Bayesian estimation of a multilevel IRT +model. *Psychometrika, 66*, 271–288. https://doi.org/10.1007/BF02294839 + +Jeon, M., & Rabe-Hesketh, S. (2016). An autoregressive growth model for +longitudinal item analysis. *Psychometrika, 81*(3), 830–850. +https://doi.org/10.1007/s11336-015-9489-2 + +Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal +data. *Biometrics, 38*(4), 963–974. https://doi.org/10.2307/2529876 + +Oravecz, Z., Tuerlinckx, F., & Vandekerckhove, J. (2011). A hierarchical +latent stochastic differential equation model for affective dynamics. +*Psychological Methods, 16*(2), 468–490. https://doi.org/10.1037/a0024375 + +Browne, W. J., Goldstein, H., & Rasbash, J. (2001). Multiple membership +multiple classification (MMMC) models. *Statistical Modelling, 1*(2), +103–124. https://doi.org/10.1177/1471082X0100100202 + +## Follow-ups + +1. Jointly estimate context-level `u_h` with longitudinal states (Fox & Glas + plus Browne et al. MMMC), after a dedicated identification study. +2. Adaptive-quadrature or Gibbs alternatives if MAP shrinkage is insufficient + for a stated inferential target. +3. GPU only if a kernel implements this exact objective and passes CPU + parity. + +## Reversal / supersession conditions + +A later ADR should supersede this record if the joint likelihood changes +(full discrimination-vector MLS2PLM, estimated MMMC, or a different +transition family), if recovery evidence falsifies the stated MAP bounds, or +if a GPU path is added for this estimand. diff --git a/docs/adr/README.md b/docs/adr/README.md index 8593334ba..27a2ce988 100644 --- a/docs/adr/README.md +++ b/docs/adr/README.md @@ -33,6 +33,7 @@ A conversation, issue, PR body, design note, or paper summary is not an Accepted | [0014](0014-bounded-llm-judge-category-inputs.md) | Proposed | Bound LLM-judge category inputs to exact built-in scalars and keep model/provider security evidence fail-closed and independently verifiable. | | [0015](0015-multi-item-irt-fit-boundary.md) | Proposed | Enforce the multi-item dichotomous/polytomous contract at public IRT fitters and require explicit readiness evidence before interpreting estimates. | | [0018](0018-rust-longitudinal-state-engine.md) | Proposed | Rust owns the first respondent-level longitudinal state layer as independent OLS trends and caller-supplied discrete AR; full joint multilevel estimation remains gated. | +| [0019](0019-joint-hierarchical-ctar-rasch.md) | Proposed | Joint MAP hierarchical continuous-time AR(1) Rasch estimates shared `(mu, tau, lambda)` and person-occasion states; MMMC and GPU parity remain excluded. | ## ADR completeness rule diff --git a/docs/changelog.d/hierarchical-ctar-rasch.md b/docs/changelog.d/hierarchical-ctar-rasch.md new file mode 100644 index 000000000..e40b5111a --- /dev/null +++ b/docs/changelog.d/hierarchical-ctar-rasch.md @@ -0,0 +1,18 @@ +# Joint MAP hierarchical continuous-time AR(1) Rasch + +## Added + +- 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. diff --git a/docs/doctoring/multilevel_longitudinal_measurement.md b/docs/doctoring/multilevel_longitudinal_measurement.md index 0a3b685dd..ce53bb38d 100644 --- a/docs/doctoring/multilevel_longitudinal_measurement.md +++ b/docs/doctoring/multilevel_longitudinal_measurement.md @@ -7,16 +7,21 @@ cross-classified, weighted multiple-membership, and repeated longitudinal measurement. The contracts prevent contextual and temporal provenance from being collapsed into respondent IDs or unstructured metadata. A Rust-owned state layer now consumes the sealed longitudinal design for independent -per-respondent OLS trends and discrete-step AR(1) state prediction; the full -joint multilevel likelihood remains outside this slice. +per-respondent OLS trends and discrete-step AR(1) state prediction. A +separate joint MAP hierarchical continuous-time AR(1) Rasch slice estimates +shared population hyperparameters and person-occasion states. Estimated +multiple-membership `u_h` and Fox–Glas / adaptive-quadrature multilevel IRT +remain outside both slices. Python owns validation, content identity, bounded collection handling, replay protection, sparse design marshalling, and serialization. Rust owns the -weighted contextual predictor and the first state-layer arithmetic. The -OLS path uses respondent-sharded CPU threads; the AR path uses -sequence-index gaps and a caller-supplied coefficient. Full likelihood, -random-effect integration, optimization, uncertainty, and GPU batching for -recurrent state remain separate boundaries. +weighted contextual predictor, the first state-layer arithmetic, and the +joint MAP hierarchical CT-AR Rasch kernel. The OLS path uses +respondent-sharded CPU threads; the AR path uses sequence-index gaps and a +caller-supplied coefficient; the hierarchical CT-AR path uses elapsed-day +Ornstein–Uhlenbeck transitions and packed L-BFGS. GPU batching for this +Rasch objective is not implemented: the existing wgpu path owns MLSIRM +distance/likelihood kernels, a different estimand. ## Scientific rationale @@ -71,6 +76,9 @@ The implementation: - distinguishes a random-intercept/slope wire label from stationary AR(1); - fits the state layer with Rust-only arithmetic and deterministic respondent sharding, reporting independent OLS or caller-supplied AR estimand metadata; +- fits a separate joint MAP hierarchical CT-AR Rasch slice with estimated + `(mu, tau, lambda)`, Wald observed-information intervals, and explicit + exclusion of multiple-membership random effects and GPU parity; - reports state RMSE, observed/transition counts, and worker-count-invariant results; - keeps lagged-response dependence independently switchable; @@ -82,17 +90,16 @@ The implementation: ## Temporal interpretation boundary -The current `autoregressive_coefficient` is a discrete occasion-step stationary -AR(1) coefficient with \(-1<\phi<1\). Irregular millisecond offsets are retained -as exact ordering and audit provenance. They do **not** imply that one \(\phi\) -is automatically adjusted for elapsed time. +The ADR-0018 `autoregressive_coefficient` is a discrete occasion-step +stationary AR(1) coefficient with \(-1<\phi<1\). Irregular millisecond +offsets are retained as exact ordering and audit provenance. They do **not** +imply that one \(\phi\) is automatically adjusted for elapsed time. -A continuous-time model would require a separate parameterization, for example a -transition rate mapped to interval-specific correlations, plus explicit units, -identification, recovery, and numerical stability evidence. That arithmetic is -reserved for a later Rust estimator PR. Until then, no report may describe the -current coefficient as continuous-time, interval-adjusted, or comparable across -different occasion spacings without an explicit design assumption. +A separate joint MAP slice (ADR-0019) parameterizes elapsed time through +\(\phi_{pt}=\exp(-\lambda\Delta_{pt})\) with \(\Delta_{pt}\) in days. That +slice may be described as continuous-time AR(1) / Ornstein–Uhlenbeck. The +OLS/AR layer still must not be described as continuous-time, +interval-adjusted, or comparable across different occasion spacings. ## Identification and interpretation limits @@ -123,15 +130,19 @@ must use explicit schema migration rather than mutating the accepted contract. The current evidence includes contract validation, deterministic identity, child replay, resource bounds, dimension-scoped assignment, strict temporal ordering, Rust state-layer slope recovery, missing-occasion preservation, -discrete AR transition RMSE, worker-count determinism, and source-text-free -serialization. It is not evidence of: - -- full joint estimator correctness; -- variance-component identification; -- population random-effects recovery or interval coverage; +discrete AR transition RMSE, worker-count determinism, source-text-free +serialization, and joint MAP hierarchical CT-AR Rasch recovery of known +states with measurement-information Wald interval coverage. Shared mean and +MAP-shrunk `tau` are recovered under documented RMSE bounds. Short series +leave `lambda` weakly identified; the transition claim is elapsed-day +`phi_pt=exp(-lambda Delta_pt)` with a finite positive decay, not tight +unbiased recovery of `lambda`. It is not evidence of: + +- Fox and Glas Gibbs or Jeon and Rabe-Hesketh adaptive-quadrature ML; +- unbiased maximum-likelihood variance-component recovery; +- estimated multiple-membership or crossed `u_h`; - measurement invariance or fairness; - causal contextual effects; -- continuous-time dynamics; - GPU recurrent-state performance or parity; or - high-stakes deployment readiness. @@ -163,6 +174,13 @@ Jeon, M., & Rabe-Hesketh, S. (2016). An autoregressive growth model for longitudinal item analysis. *Psychometrika, 81*(3), 830–850. https://doi.org/10.1007/s11336-015-9489-2 +Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal +data. *Biometrics, 38*(4), 963–974. https://doi.org/10.2307/2529876 + +Oravecz, Z., Tuerlinckx, F., & Vandekerckhove, J. (2011). A hierarchical +latent stochastic differential equation model for affective dynamics. +*Psychological Methods, 16*(2), 468–490. https://doi.org/10.1037/a0024375 + Tranmer, M., Steel, D., & Browne, W. J. (2014). Multiple-membership multiple-classification models for social network and group dependencies. *Journal of the Royal Statistical Society: Series A (Statistics in Society), diff --git a/docs/documentation_coverage.md b/docs/documentation_coverage.md index fb9fee3a5..453c7f202 100644 --- a/docs/documentation_coverage.md +++ b/docs/documentation_coverage.md @@ -84,7 +84,7 @@ The table below records product truth, not documentation-file presence. “Imple | Formal non-nested distinguishability/model comparison | PARTIAL | fail-closed relation-aware comparison exists; additional family-specific evidence and metadata remain incremental | | Adaptive rotation criterion selection | IMPLEMENTED_ON_PROTECTED_MAIN / PARTIAL | Rust-backed criterion registry/multi-start selector/report surfaces are integrated; additional criteria/GPU/recovery remain incremental | | Multilevel / cross-classified / multiple-membership contracts | IMPLEMENTED_ON_PROTECTED_MAIN / PARTIAL | contextual and longitudinal contracts are integrated; estimator identification/recovery remains separate work | -| Temporal/longitudinal/drift estimators | IMPLEMENTED_ON_ACTIVE_PR / PARTIAL | independent per-respondent OLS and caller-supplied discrete AR state predictors exist on the longitudinal state-engine PR; continuous-time, shrinkage, and joint multilevel IRT remain separate recovery work | +| Temporal/longitudinal/drift estimators | IMPLEMENTED_ON_ACTIVE_PR / PARTIAL | independent OLS/AR (ADR-0018) and joint MAP hierarchical CT-AR Rasch (ADR-0019) exist on stacked longitudinal PRs; estimated MMMC `u_h` and GPU parity remain excluded | | Automated essay scoring calibration/validation | IMPLEMENTED_ON_PROTECTED_MAIN / PARTIAL | governed essay contracts/validation/reporting exist; generalized rater discrimination/range/drift remains incremental | | Paired automated-vs-reference rating-range evidence | IMPLEMENTED_ON_PROTECTED_MAIN | Rust-owned paired range/compression diagnostic is integrated | | Enterprise issue measurement | IMPLEMENTED_ON_PROTECTED_MAIN / PARTIAL | reusable evidence/calibration adapters exist; causal intervention utility remains downstream/policy-bound | diff --git a/docs/traceability/requirements-matrix.md b/docs/traceability/requirements-matrix.md index 3ccc5de87..a2a6a15ca 100644 --- a/docs/traceability/requirements-matrix.md +++ b/docs/traceability/requirements-matrix.md @@ -24,7 +24,7 @@ This matrix makes the major product requirements discoverable without reconstruc | Adaptive rotation | PRD-FR-051/052, TRD-ROT | ADR-0009 | protected main contains `crates/mlsirm-core/src/rotation/`, PyO3 bindings, `python/fast_mlsirm/rotation.py`, `rotation_selection.py`, package-root exports, criterion-neutral selection and rotation regression/doctoring evidence | Accepted CPU baseline / planned GPU and broader recovery extensions | | True-parameter recovery | PRD-PRN-003, TRD-TEST-003..006 | ADR-0008 | simulation/recovery reports, Rust/NumPy parity, scheduled statistical studies/recovery contracts | Accepted | | Correlation vs recovery/agreement | PRD-PRN-003, scoring validity requirements | ADR-0008, ADR-0005 | recovery/simulation, agreement/QWK/facets evidence | Accepted: correlation is supplementary association evidence, never sole proof of parameter recovery or interchangeability | -| Multilevel/multiple-membership/temporal | PRD-FR-060..062, TRD-MLT | ADR-0007 | contextual summaries exist; full reusable contract PR remains open and Rust estimator recovery is future work | Proposed/partial / active PR | +| Multilevel/multiple-membership/temporal | PRD-FR-060..062, TRD-MLT | ADR-0007, ADR-0018, ADR-0019 | contracts plus OLS/AR state layer and a separate joint MAP hierarchical CT-AR Rasch slice exist on the stacked longitudinal PRs; estimated MMMC `u_h` and GPU parity remain excluded | Proposed/partial / active PR | | Accessible standalone reports | PRD-FR-070..072, NFR-004 | ADR-0005 | report renderers, exact-value exports, WCAG-focused regression/doctoring | Accepted/evolving | | Sensitive data / PII utility | privacy/security requirements | ADR-0012 | source-free/digest/opaque-id provenance where implemented; provider error redaction; hosted identity/retention downstream | Accepted reusable policy / Downstream operations | | Continuous execution / documentation governance | TRD-DOC-002 / work-conserving automation | ADR-0013 | single-writer exact branch head; work-conserving when blocked; feasibility-first prioritization | Accepted governance | diff --git a/docs/traceability/research-basis.md b/docs/traceability/research-basis.md index 54cea1ea3..bcaa3b866 100644 --- a/docs/traceability/research-basis.md +++ b/docs/traceability/research-basis.md @@ -138,8 +138,15 @@ The first Rust state layer (ADR-0018) is an independent per-respondent OLS trend and a caller-supplied discrete AR predictor. It is not the Fox and Glas multilevel IRT estimand and does not estimate population random effects. +ADR-0019 is a separate joint MAP hierarchical continuous-time AR(1) Rasch +slice. It estimates shared `(mu, tau, lambda)` and person-occasion states +with Wald observed-information intervals. It is not Fox and Glas Gibbs +sampling, not Jeon and Rabe-Hesketh adaptive-quadrature ML, and not estimated +multiple-membership `u_h`. + - Jeon, M., & Rabe-Hesketh, S. (2016). An autoregressive growth model for longitudinal item analysis. *Psychometrika, 81*(3), 830–850. https://doi.org/10.1007/s11336-015-9489-2 - Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. *Biometrics, 38*(4), 963–974. https://doi.org/10.2307/2529876 +- Oravecz, Z., Tuerlinckx, F., & Vandekerckhove, J. (2011). A hierarchical latent stochastic differential equation model for affective dynamics. *Psychological Methods, 16*(2), 468–490. https://doi.org/10.1037/a0024375 ## 9. Adaptive factor rotation — Proposed diff --git a/docs/verification_validation_plan.md b/docs/verification_validation_plan.md index 0ea308ffb..917e43be3 100644 --- a/docs/verification_validation_plan.md +++ b/docs/verification_validation_plan.md @@ -191,7 +191,7 @@ Compare parameter/SE/coverage behavior against atomistic misspecification so the ### VV-SCI-008 — Temporal/longitudinal recovery -Validate ordering, missing occasions, unequal follow-up patterns, drift/state parameters, random intercept/slope effects, and revision boundaries. A discrete-step model is evaluated by step count; a future continuous-time model must include interval-sensitive generating processes and recovery. +Validate ordering, missing occasions, unequal follow-up patterns, drift/state parameters, random intercept/slope effects, and revision boundaries. A discrete-step model is evaluated by step count. The joint MAP hierarchical CT-AR Rasch slice must include interval-sensitive generating processes and recovery of known states and `(mu, tau, lambda)` against honest MAP RMSE/coverage bounds. Estimated multiple-membership `u_h` remains a later recovery target. ## 5. Automated scoring / LLM-as-a-Judge validation diff --git a/python/fast_mlsirm/multilevel/__init__.py b/python/fast_mlsirm/multilevel/__init__.py index aa38032e9..fe119404c 100644 --- a/python/fast_mlsirm/multilevel/__init__.py +++ b/python/fast_mlsirm/multilevel/__init__.py @@ -14,7 +14,12 @@ build_longitudinal_state_spec, build_temporal_occasion, ) -from .estimation import fit_longitudinal_state, weighted_contextual_effect +from .estimation import ( + fit_hierarchical_longitudinal_irt, + fit_longitudinal_state, + simulate_hierarchical_longitudinal_irt, + weighted_contextual_effect, +) __all__ = [ "ContextMembership", @@ -29,6 +34,8 @@ "build_longitudinal_design", "build_longitudinal_state_spec", "build_temporal_occasion", + "fit_hierarchical_longitudinal_irt", "fit_longitudinal_state", + "simulate_hierarchical_longitudinal_irt", "weighted_contextual_effect", ] diff --git a/python/fast_mlsirm/multilevel/estimation.py b/python/fast_mlsirm/multilevel/estimation.py index 64405f7b2..728033d5e 100644 --- a/python/fast_mlsirm/multilevel/estimation.py +++ b/python/fast_mlsirm/multilevel/estimation.py @@ -5,15 +5,16 @@ the flat CSR arrays ``mlsirm_core::multilevel::weighted_contextual_effect`` expects, converting the caller's per-context random-effect values into the matching flat vector, and converting a sealed ``LongitudinalDesign`` into the -row-offset arrays ``mlsirm_core::longitudinal::fit_longitudinal_state`` -expects. The additive contextual sum, independent per-respondent OLS trends, -caller-supplied discrete AR predictions, worker determinism, and numerical -input validation are owned by the Rust core. +row-offset arrays the Rust longitudinal kernels expect. The additive +contextual sum, independent per-respondent OLS trends, caller-supplied +discrete AR predictions, joint MAP hierarchical continuous-time AR(1) Rasch +estimation, worker determinism, and numerical input validation are owned by +the Rust core. """ from __future__ import annotations -from collections.abc import Mapping +from collections.abc import Mapping, Sequence import numpy as np @@ -128,6 +129,24 @@ def weighted_contextual_effect( ) +def _ordered_longitudinal_rows(design: LongitudinalDesign) -> tuple[list[int], list, list[int]]: + """Return CSR offsets, ordered occasions, and millisecond offsets.""" + grouped: dict[str, list] = { + respondent_id: [] for respondent_id in design.respondent_ids + } + for occasion in design.occasions: + grouped[occasion.respondent_id].append(occasion) + row_offsets = [0] + time_offsets: list[int] = [] + ordered_occasions: list = [] + for respondent_id in design.respondent_ids: + for occasion in grouped[respondent_id]: + time_offsets.append(occasion.time_offset_milliseconds) + ordered_occasions.append(occasion) + row_offsets.append(len(time_offsets)) + return row_offsets, ordered_occasions, time_offsets + + def _observed_value(values: Mapping[str, float], occasion_id: str) -> float: """Return one caller observation as a real float, or NaN when absent.""" try: @@ -188,25 +207,11 @@ def fit_longitudinal_state( if type(design) is not LongitudinalDesign: raise ValueError("design must be an exact LongitudinalDesign") _ = design.design_fingerprint - occasions = list(design.occasions) - grouped: dict[str, list] = { - respondent_id: [] for respondent_id in design.respondent_ids - } - for occasion in occasions: - grouped[occasion.respondent_id].append(occasion) - - row_offsets = [0] - sequence_indices: list[int] = [] - time_offsets: list[int] = [] - observations: list[float] = [] - ordered_occasions: list = [] - for respondent_id in design.respondent_ids: - for occasion in grouped[respondent_id]: - sequence_indices.append(occasion.sequence_index) - time_offsets.append(occasion.time_offset_milliseconds) - observations.append(_observed_value(values, occasion.occasion_id)) - ordered_occasions.append(occasion) - row_offsets.append(len(time_offsets)) + row_offsets, ordered_occasions, time_offsets = _ordered_longitudinal_rows(design) + sequence_indices = [occasion.sequence_index for occasion in ordered_occasions] + observations = [ + _observed_value(values, occasion.occasion_id) for occasion in ordered_occasions + ] state_kind = design.state_spec.state_kind ar_coefficient = design.state_spec.autoregressive_coefficient if state_kind is LongitudinalStateKind.RANDOM_INTERCEPT_SLOPE: @@ -256,4 +261,265 @@ def fit_longitudinal_state( } -__all__ = ["ContextKey", "fit_longitudinal_state", "weighted_contextual_effect"] +def _validate_binary_response_matrix( + responses: object, + n_occasions: int, +) -> np.ndarray: + """Return a C-contiguous float64 response matrix or a package-owned error.""" + if isinstance(responses, (bool, np.bool_)) or not isinstance(responses, np.ndarray): + raise ValueError("responses must be a NumPy ndarray") + if responses.ndim != 2: + raise ValueError("responses must be a two-dimensional occasion-by-item matrix") + if responses.shape[0] != n_occasions: + raise ValueError("responses rows must align with the sealed occasion order") + if responses.shape[1] < 2: + raise ValueError("hierarchical CT-AR Rasch requires at least two items") + if responses.dtype == np.bool_ or np.issubdtype(responses.dtype, np.bool_): + raise ValueError("responses must be 0, 1, or NaN rather than Boolean values") + try: + matrix = np.ascontiguousarray(responses, dtype=np.float64) + except Exception: + raise ValueError("responses could not be converted to float64 safely") from None + finite = np.isfinite(matrix) + invalid = finite & (matrix != 0.0) & (matrix != 1.0) + if np.any(invalid): + raise ValueError("responses must be 0, 1, or NaN") + return matrix + + +def _item_labels(item_ids: Sequence[str] | None, n_items: int) -> list[str]: + """Return caller item labels or stable positional defaults.""" + if item_ids is None: + return [f"item_{index}" for index in range(n_items)] + if isinstance(item_ids, (str, bytes)) or not isinstance(item_ids, Sequence): + raise ValueError("item_ids must be a sequence of item identifiers") + labels = list(item_ids) + if len(labels) != n_items: + raise ValueError("item_ids length must equal the response item axis") + for label in labels: + if not isinstance(label, str) or not label: + raise ValueError("item_ids entries must be non-empty strings") + return labels + + +def fit_hierarchical_longitudinal_irt( + design: LongitudinalDesign, + responses: np.ndarray, + *, + item_ids: Sequence[str] | None = None, + worker_count: int = 1, + max_iter: int = 250, + tolerance: float = 1e-5, + hessian_step: float = 1e-3, +) -> dict[str, object]: + """Fit the joint MAP hierarchical continuous-time AR(1) Rasch slice. + + This entry point does **not** interpret ``design.state_spec.state_kind``. + The sealed design supplies respondent identity, occasion order, and exact + millisecond offsets only. The estimand is joint MAP of a Rasch measurement + model and a hierarchical stationary Ornstein–Uhlenbeck / continuous-time + AR(1) latent-state process. It is not independent respondent OLS, not a + caller-supplied discrete AR coefficient, not Fox and Glas (2001) Gibbs + sampling, and not Jeon and Rabe-Hesketh (2016) adaptive-quadrature ML. + + Crossed and multiple-membership random effects are excluded from this + joint likelihood. The existing GPU abstraction owns MLSIRM + distance/likelihood kernels, not this hierarchical CT-AR Rasch objective, + so ``gpu_parity`` is reported as false. + + Parameters + ---------- + design: + A package-built ``LongitudinalDesign``. Integrity is verified before + marshalling so a tampered or hand-constructed design raises here. + responses: + Occasion-major binary matrix with shape ``(n_occasions, n_items)``, + aligned with ``design.occasions`` after respondent-then-sequence + ordering. Values must be ``0``, ``1``, or ``NaN``. + item_ids: + Optional item labels aligned with the response columns. Defaults to + ``item_0``, ``item_1``, ... + worker_count: + Number of deterministic person-shard worker threads (``>= 1``). + max_iter: + Maximum packed L-BFGS iterations (``>= 1``). + tolerance: + Relative L-BFGS tolerance; must be finite and strictly positive. + hessian_step: + Central-difference step for the hyperparameter observed Hessian. + + Returns + ------- + dict + Joint MAP states, Wald intervals, item intercepts, estimated + population mean/sd/decay, unit-day AR coefficient, counts, engine + identity, fingerprints, and normative estimand metadata. + + Raises + ------ + ValueError + If execution controls, the sealed design, or the response matrix are + invalid, or the Rust kernel rejects the design. + """ + if worker_count < 1: + raise ValueError("worker_count must be at least one") + if max_iter < 1: + raise ValueError("max_iter must be at least one") + if not np.isfinite(tolerance) or tolerance <= 0.0: + raise ValueError("tolerance must be finite and strictly positive") + if not np.isfinite(hessian_step) or hessian_step <= 0.0: + raise ValueError("hessian_step must be finite and strictly positive") + if type(design) is not LongitudinalDesign: + raise ValueError("design must be an exact LongitudinalDesign") + _ = design.design_fingerprint + row_offsets, ordered_occasions, time_offsets = _ordered_longitudinal_rows(design) + matrix = _validate_binary_response_matrix(responses, len(ordered_occasions)) + labels = _item_labels(item_ids, matrix.shape[1]) + core = multilevel_core() + result = core.fit_hierarchical_ctar_rasch( + np.asarray(row_offsets, dtype=np.uint64), + np.asarray(time_offsets, dtype=np.int64), + matrix, + worker_count, + max_iter, + float(tolerance), + float(hessian_step), + ) + return { + "estimand_scope": str(result["estimand_scope"]), + "transition_kind": str(result["transition_kind"]), + "interval_kind": str(result["interval_kind"]), + "population_random_effects_estimated": True, + "ar_coefficient_estimated": True, + "ar_coefficient_source": "joint_map", + "multiple_membership_estimated": False, + "gpu_parity": False, + "state_spec_fingerprint": design.state_spec.state_spec_fingerprint, + "design_fingerprint": design.design_fingerprint, + "state": np.asarray(result["state"], dtype=np.float64), + "state_se": np.asarray(result["state_se"], dtype=np.float64), + "state_lower": np.asarray(result["state_lower"], dtype=np.float64), + "state_upper": np.asarray(result["state_upper"], dtype=np.float64), + "item_intercepts": np.asarray(result["item_intercepts"], dtype=np.float64), + "item_ids": labels, + "population_mean": float(result["population_mean"]), + "population_sd": float(result["population_sd"]), + "decay_rate": float(result["decay_rate"]), + "unit_time_ar_coefficient": float(result["unit_time_ar_coefficient"]), + "hyperparameter_se": np.asarray(result["hyperparameter_se"], dtype=np.float64), + "hyperparameter_lower": np.asarray(result["hyperparameter_lower"], dtype=np.float64), + "hyperparameter_upper": np.asarray(result["hyperparameter_upper"], dtype=np.float64), + "hyperparameter_intervals_identified": bool( + result["hyperparameter_intervals_identified"] + ), + "state_intervals_identified": bool(result["state_intervals_identified"]), + "observed_count": int(result["observed_count"]), + "transition_count": int(result["transition_count"]), + "status": str(result["status"]), + "engine": str(result["engine"]), + "respondent_ids": list(design.respondent_ids), + "occasion_ids": [occasion.occasion_id for occasion in ordered_occasions], + "occasion_records": [ + { + "occasion_id": occasion.occasion_id, + "respondent_id": occasion.respondent_id, + "sequence_index": occasion.sequence_index, + "time_offset_milliseconds": occasion.time_offset_milliseconds, + } + for occasion in ordered_occasions + ], + } + + +def simulate_hierarchical_longitudinal_irt( + design: LongitudinalDesign, + *, + item_intercepts: Sequence[float], + population_mean: float = 0.0, + population_sd: float = 0.7, + decay_rate: float = 0.35, + seed: int = 1, +) -> dict[str, object]: + """Simulate hierarchical CT-AR Rasch states and binary responses. + + The simulator is the recovery-fixture generator for + ``fit_hierarchical_longitudinal_irt``. It is not a claim that the + subsequent fit recovers these parameters without shrinkage. + + Parameters + ---------- + design: + A package-built ``LongitudinalDesign`` supplying occasion times. + item_intercepts: + Generating Rasch item intercepts. + population_mean: + Generating population mean of the latent-state process. + population_sd: + Generating stationary standard deviation. + decay_rate: + Generating continuous-time decay rate per day. + seed: + Deterministic unsigned seed forwarded to the Rust LCG. + + Returns + ------- + dict + Generating states and an occasion-major response matrix aligned with + the sealed design order. + + Raises + ------ + ValueError + If the design or generating parameters are invalid. + """ + if type(design) is not LongitudinalDesign: + raise ValueError("design must be an exact LongitudinalDesign") + _ = design.design_fingerprint + if isinstance(item_intercepts, (str, bytes)) or not isinstance( + item_intercepts, (Sequence, np.ndarray) + ): + raise ValueError("item_intercepts must be a sequence of real numbers") + try: + intercepts = np.asarray(list(item_intercepts), dtype=np.float64) + except Exception: + raise ValueError("item_intercepts could not be converted safely") from None + if intercepts.ndim != 1 or intercepts.size < 2: + raise ValueError("item_intercepts must contain at least two finite values") + if not np.all(np.isfinite(intercepts)): + raise ValueError("item_intercepts must be finite") + if seed < 0: + raise ValueError("seed must be a non-negative integer") + row_offsets, ordered_occasions, time_offsets = _ordered_longitudinal_rows(design) + core = multilevel_core() + result = core.simulate_hierarchical_ctar_rasch( + np.asarray(row_offsets, dtype=np.uint64), + np.asarray(time_offsets, dtype=np.int64), + intercepts, + float(population_mean), + float(population_sd), + float(decay_rate), + int(seed), + ) + n_items = int(result["n_items"]) + responses = np.asarray(result["responses"], dtype=np.float64).reshape( + (len(ordered_occasions), n_items) + ) + return { + "state": np.asarray(result["state"], dtype=np.float64), + "responses": responses, + "item_intercepts": intercepts, + "population_mean": float(population_mean), + "population_sd": float(population_sd), + "decay_rate": float(decay_rate), + "occasion_ids": [occasion.occasion_id for occasion in ordered_occasions], + "design_fingerprint": design.design_fingerprint, + } + + +__all__ = [ + "ContextKey", + "fit_hierarchical_longitudinal_irt", + "fit_longitudinal_state", + "simulate_hierarchical_longitudinal_irt", + "weighted_contextual_effect", +] diff --git a/tests/test_hierarchical_longitudinal_irt.py b/tests/test_hierarchical_longitudinal_irt.py new file mode 100644 index 000000000..061e7cd1a --- /dev/null +++ b/tests/test_hierarchical_longitudinal_irt.py @@ -0,0 +1,295 @@ +"""End-to-end tests for the joint MAP hierarchical CT-AR Rasch slice.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from fast_mlsirm._multilevel_core_loader import multilevel_core +from fast_mlsirm.multilevel import ( + LongitudinalStateKind, + build_longitudinal_design, + build_longitudinal_state_spec, + build_temporal_occasion, + fit_hierarchical_longitudinal_irt, + simulate_hierarchical_longitudinal_irt, +) + + +def _occasion(respondent: str, occasion: str, sequence: int, offset: int): + """Build one sealed occasion with a unique revision identity.""" + occasion_id = f"occasion_{occasion}" + return build_temporal_occasion( + respondent_id=respondent, + occasion_id=occasion_id, + sequence_index=sequence, + time_offset_milliseconds=offset, + occasion_revision_fingerprint=(occasion_id + " revision") + .encode() + .hex() + .ljust(64, "0")[:64], + ) + + +def _design(n_persons: int, n_occasions: int, irregular: bool = False): + """Return a sealed longitudinal design with optional irregular gaps.""" + occasions = [] + for person in range(n_persons): + elapsed = 0 + for occasion in range(n_occasions): + if irregular and occasion == 2: + elapsed += 2 + elif occasion > 0: + elapsed += 1 + occasions.append( + _occasion( + f"respondent_{person}", + f"{person}_{occasion}", + occasion, + elapsed * 86_400_000, + ) + ) + return build_longitudinal_design( + occasions=occasions, + state_spec=build_longitudinal_state_spec( + state_kind=LongitudinalStateKind.RANDOM_INTERCEPT_SLOPE, + ), + ) + + +def test_public_fit_recovers_states_across_seeds_and_is_worker_invariant() -> None: + """Multi-seed recovery stays inside honest MAP RMSE/coverage bounds.""" + design = _design(8, 3, irregular=True) + items = np.array([-0.6, -0.2, 0.2, 0.6], dtype=np.float64) + state_sse = 0.0 + state_count = 0.0 + covered = 0.0 + mean_err = 0.0 + for seed in (11, 23, 41): + simulated = simulate_hierarchical_longitudinal_irt( + design, + item_intercepts=items, + population_mean=0.0, + population_sd=0.7, + decay_rate=0.35, + seed=seed, + ) + result = fit_hierarchical_longitudinal_irt( + design, + simulated["responses"], + item_ids=["a", "b", "c", "d"], + worker_count=3, + max_iter=80, + tolerance=1e-4, + ) + single = fit_hierarchical_longitudinal_irt( + design, + simulated["responses"], + item_ids=["a", "b", "c", "d"], + worker_count=1, + max_iter=80, + tolerance=1e-4, + ) + np.testing.assert_allclose(result["state"], single["state"], atol=1e-8) + assert result["estimand_scope"] == "joint_map_hierarchical_ctar_rasch" + assert result["transition_kind"] == "continuous_time_ar1_ou" + assert result["interval_kind"] == "wald_measurement_observed_information" + assert result["engine"] == "rust_cpu_multithreaded" + assert result["population_random_effects_estimated"] is True + assert result["ar_coefficient_estimated"] is True + assert result["ar_coefficient_source"] == "joint_map" + assert result["multiple_membership_estimated"] is False + assert result["gpu_parity"] is False + assert result["item_ids"] == ["a", "b", "c", "d"] + assert result["estimand_scope"] != "independent_respondent_ols_trend" + assert result["state_intervals_identified"] is True + truth = np.asarray(simulated["state"], dtype=np.float64) + state_sse += float(np.sum((result["state"] - truth) ** 2)) + state_count += truth.size + covered += float( + np.sum( + (result["state_lower"] <= truth) & (truth <= result["state_upper"]) + ) + ) + mean_err += (float(result["population_mean"]) - 0.0) ** 2 + state_rmse = (state_sse / state_count) ** 0.5 + coverage = covered / state_count + mean_rmse = (mean_err / 3.0) ** 0.5 + assert state_rmse < 0.85, state_rmse + assert coverage > 0.80, coverage + assert mean_rmse < 0.35, mean_rmse + + +def test_missing_responses_and_irregular_gaps_are_honored() -> None: + """NaN responses are excluded and longer gaps shrink the CT-AR weight.""" + design = _design(1, 3, irregular=True) + responses = np.array( + [ + [1.0, 0.0, 1.0, 0.0], + [np.nan, np.nan, 0.0, 1.0], + [0.0, 1.0, np.nan, np.nan], + ], + dtype=np.float64, + ) + result = fit_hierarchical_longitudinal_irt( + design, + responses, + worker_count=1, + max_iter=60, + tolerance=1e-4, + ) + assert result["observed_count"] == 8 + assert result["transition_count"] == 2 + assert result["unit_time_ar_coefficient"] < 1.0 + assert result["item_ids"] == ["item_0", "item_1", "item_2", "item_3"] + + +def test_public_fit_rejects_invalid_controls_and_foreign_designs() -> None: + """Python marshalling fails closed before native dispatch.""" + design = _design(1, 2) + responses = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.float64) + with pytest.raises(ValueError, match="worker_count"): + fit_hierarchical_longitudinal_irt(design, responses, worker_count=0) + with pytest.raises(ValueError, match="max_iter"): + fit_hierarchical_longitudinal_irt(design, responses, max_iter=0) + with pytest.raises(ValueError, match="tolerance"): + fit_hierarchical_longitudinal_irt(design, responses, tolerance=0.0) + with pytest.raises(ValueError, match="hessian_step"): + fit_hierarchical_longitudinal_irt(design, responses, hessian_step=-1.0) + + class ForeignDesign: + """Represent an object that was not produced by the package factory.""" + + with pytest.raises(ValueError, match="LongitudinalDesign"): + fit_hierarchical_longitudinal_irt(ForeignDesign(), responses) # type: ignore[arg-type] + + +@pytest.mark.parametrize( + "payload, match", + [ + (True, "NumPy ndarray"), + (np.array([1.0, 0.0]), "two-dimensional"), + (np.array([[1.0, 0.0]], dtype=np.float64), "align with the sealed"), + (np.array([[1.0], [0.0]], dtype=np.float64), "at least two items"), + (np.array([[True, False], [False, True]]), "Boolean"), + (np.array([[2.0, 0.0], [0.0, 1.0]], dtype=np.float64), "0, 1, or NaN"), + (np.array([["a", "b"], ["c", "d"]], dtype=object), "float64 safely"), + ], +) +def test_public_fit_rejects_invalid_response_matrices( + payload: object, match: str +) -> None: + """Hostile or malformed response matrices raise package-owned errors.""" + with pytest.raises(ValueError, match=match): + fit_hierarchical_longitudinal_irt(_design(1, 2), payload) # type: ignore[arg-type] + + +@pytest.mark.parametrize( + "item_ids, match", + [ + ("item_a", "sequence of item identifiers"), + (["a"], "length must equal"), + (["a", ""], "non-empty strings"), + (["a", 1], "non-empty strings"), + ], +) +def test_public_fit_rejects_invalid_item_labels( + item_ids: object, match: str +) -> None: + """Item labels are validated before the native kernel is invoked.""" + responses = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.float64) + with pytest.raises(ValueError, match=match): + fit_hierarchical_longitudinal_irt( + _design(1, 2), + responses, + item_ids=item_ids, # type: ignore[arg-type] + ) + + +def test_simulate_accepts_list_and_ndarray_item_intercepts() -> None: + """Generating intercepts may be a Python sequence or a NumPy vector.""" + design = _design(1, 2) + listed = simulate_hierarchical_longitudinal_irt( + design, item_intercepts=[-0.2, 0.2], seed=2 + ) + arrayed = simulate_hierarchical_longitudinal_irt( + design, item_intercepts=np.array([-0.2, 0.2], dtype=np.float64), seed=2 + ) + np.testing.assert_array_equal(listed["responses"], arrayed["responses"]) + np.testing.assert_array_equal(listed["state"], arrayed["state"]) + + +def test_simulate_and_fit_reject_invalid_generating_controls() -> None: + """The recovery simulator validates generating parameters locally.""" + design = _design(1, 2) + with pytest.raises(ValueError, match="LongitudinalDesign"): + simulate_hierarchical_longitudinal_irt( + object(), # type: ignore[arg-type] + item_intercepts=[-0.2, 0.2], + ) + with pytest.raises(ValueError, match="sequence of real numbers"): + simulate_hierarchical_longitudinal_irt(design, item_intercepts="ab") # type: ignore[arg-type] + with pytest.raises(ValueError, match="at least two"): + simulate_hierarchical_longitudinal_irt(design, item_intercepts=[0.1]) + with pytest.raises(ValueError, match="finite"): + simulate_hierarchical_longitudinal_irt( + design, item_intercepts=[0.1, float("nan")] + ) + with pytest.raises(ValueError, match="non-negative"): + simulate_hierarchical_longitudinal_irt( + design, item_intercepts=[-0.2, 0.2], seed=-1 + ) + with pytest.raises(ValueError, match="converted safely"): + simulate_hierarchical_longitudinal_irt( + design, item_intercepts=[object(), object()] # type: ignore[list-item] + ) + with pytest.raises(ValueError, match="at least one respondent"): + fit_hierarchical_longitudinal_irt( + _design(1, 1), + np.array([[1.0, 0.0]], dtype=np.float64), + ) + + +def test_raw_binding_bounds_hierarchical_axes() -> None: + """The raw extension bounds occasion and item axes before native work.""" + core = multilevel_core() + with pytest.raises(ValueError, match="occasion axis exceeds"): + core.fit_hierarchical_ctar_rasch( + np.array([0, 1], dtype=np.uint64), + np.array([0], dtype=np.int64), + np.zeros((100_001, 2), dtype=np.float64), + 1, + 1, + 1e-4, + 1e-3, + ) + with pytest.raises(ValueError, match="item axis exceeds"): + core.fit_hierarchical_ctar_rasch( + np.array([0, 1], dtype=np.uint64), + np.array([0], dtype=np.int64), + np.zeros((1, 4_097), dtype=np.float64), + 1, + 1, + 1e-4, + 1e-3, + ) + with pytest.raises(ValueError, match="item_intercepts exceeds"): + core.simulate_hierarchical_ctar_rasch( + np.array([0, 1], dtype=np.uint64), + np.array([0], dtype=np.int64), + np.zeros(4_097, dtype=np.float64), + 0.0, + 0.5, + 0.4, + 1, + ) + with pytest.raises(ValueError, match="time offsets exceed"): + core.simulate_hierarchical_ctar_rasch( + np.array([0, 100_001], dtype=np.uint64), + np.zeros(100_001, dtype=np.int64), + np.array([-0.2, 0.2], dtype=np.float64), + 0.0, + 0.5, + 0.4, + 1, + ) diff --git a/tests/test_multilevel_core_loader.py b/tests/test_multilevel_core_loader.py index 5327463b0..506263371 100644 --- a/tests/test_multilevel_core_loader.py +++ b/tests/test_multilevel_core_loader.py @@ -29,6 +29,8 @@ def test_multilevel_core_loads_and_caches_the_real_extension() -> None: """The happy path returns the same cached module on a second call.""" first = loader.multilevel_core() assert hasattr(first, "weighted_contextual_effect") + assert hasattr(first, "fit_hierarchical_ctar_rasch") + assert hasattr(first, "simulate_hierarchical_ctar_rasch") second = loader.multilevel_core() assert second is first