diff --git a/ARCHITECTURE.md b/ARCHITECTURE.md index b17fcb3e4..be6cff26e 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. +ADR-0018 proposes a state layer for independent OLS trends and caller-supplied +discrete AR predictions. ADR-0019 proposes 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/crates/fast-mlsirm-py/src/multilevel_bindings.rs b/crates/fast-mlsirm-py/src/multilevel_bindings.rs index f08eb8021..c709bc4ea 100644 --- a/crates/fast-mlsirm-py/src/multilevel_bindings.rs +++ b/crates/fast-mlsirm-py/src/multilevel_bindings.rs @@ -6,10 +6,17 @@ //! counts, and its input validation are owned by //! `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; use pyo3::types::PyModule; use pyo3::wrap_pyfunction; @@ -17,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 @@ -87,9 +96,249 @@ fn py_weighted_contextual_effect<'py>( Ok(result.to_pyarray(py)) } +/// Fit the Rust-owned repeated-measurement state layer. +/// +/// Parameters +/// ---------- +/// row_offsets : numpy.ndarray[uint64] +/// CSR-style respondent pointer, length ``n_respondents + 1``. +/// sequence_indices : numpy.ndarray[uint64] +/// Discrete occasion indices aligned with ``values``. +/// time_offsets_milliseconds : numpy.ndarray[int64] +/// Exact millisecond offsets aligned with ``values``. +/// values : numpy.ndarray[float64] +/// Observed states; ``NaN`` marks a missing occasion. +/// state_kind : str +/// Compatibility wire label for the requested state predictor. +/// ar_coefficient : float or None +/// Caller-supplied discrete AR coefficient, or ``None`` for OLS trends. +/// worker_count : int +/// Number of deterministic worker threads (``>= 1``). +/// +/// Returns +/// ------- +/// dict +/// Predicted states, respondent intercepts/slopes, RMSE, and counts. +#[pyfunction(name = "fit_longitudinal_state")] +fn py_fit_longitudinal_state<'py>( + py: Python<'py>, + row_offsets: PyReadonlyArray1<'_, u64>, + sequence_indices: PyReadonlyArray1<'_, u64>, + time_offsets_milliseconds: PyReadonlyArray1<'_, i64>, + values: PyReadonlyArray1<'_, f64>, + state_kind: &str, + ar_coefficient: Option, + worker_count: usize, +) -> PyResult> { + let row_offsets = checked_usize_values(row_offsets.as_slice()?, "row_offsets")?; + let sequence_indices = checked_usize_values(sequence_indices.as_slice()?, "sequence_indices")?; + let time_offsets = time_offsets_milliseconds.as_slice()?.to_vec(); + let values = values.as_slice()?.to_vec(); + let state_kind = state_kind.to_owned(); + let fit = py + .detach(move || { + core_fit_longitudinal_state( + &row_offsets, + &sequence_indices, + &time_offsets, + &values, + &state_kind, + ar_coefficient, + worker_count, + ) + }) + .map_err(PyValueError::new_err)?; + let result = PyDict::new(py); + result.set_item("state", fit.state.to_pyarray(py))?; + result.set_item("intercepts", fit.intercepts.to_pyarray(py))?; + result.set_item("slopes", fit.slopes.to_pyarray(py))?; + result.set_item("ar_coefficient", fit.ar_coefficient)?; + result.set_item("rmse", fit.rmse)?; + result.set_item("observed_count", fit.observed_count)?; + result.set_item("transition_count", fit.transition_count)?; + result.set_item("engine", "rust_cpu_multithreaded")?; + 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_view = item_intercepts.as_slice()?; + if item_intercepts_view.len() > MAX_HIERARCHICAL_ITEMS { + return Err(PyValueError::new_err(format!( + "item_intercepts exceeds maximum supported length of {MAX_HIERARCHICAL_ITEMS}" + ))); + } + let item_intercepts = item_intercepts_view.to_vec(); + let n_items = item_intercepts.len(); + let row_offsets = checked_usize_values(row_offsets.as_slice()?, "row_offsets")?; + let time_offsets_view = time_offsets_milliseconds.as_slice()?; + if time_offsets_view.len() > MAX_HIERARCHICAL_OCCASIONS { + return Err(PyValueError::new_err(format!( + "time offsets exceed maximum supported length of {MAX_HIERARCHICAL_OCCASIONS}" + ))); + } + let time_offsets = time_offsets_view.to_vec(); + 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 f2719a800..6a3b64a52 100644 --- a/crates/mlsirm-core/src/lib.rs +++ b/crates/mlsirm-core/src/lib.rs @@ -18,6 +18,8 @@ pub mod gtheory; 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.rs b/crates/mlsirm-core/src/longitudinal.rs new file mode 100644 index 000000000..dd72ebb84 --- /dev/null +++ b/crates/mlsirm-core/src/longitudinal.rs @@ -0,0 +1,806 @@ +//! Rust-owned longitudinal state estimation for repeated psychometric observations. +//! +//! The estimator is intentionally small and explicit. The compatibility wire +//! label `random_intercept_slope` denotes an independent per-respondent +//! ordinary-least-squares trend; it does not estimate a population random- +//! effects distribution or apply shrinkage. A stationary AR(1) state uses the +//! caller-supplied discrete-occasion coefficient and produces one-step latent +//! predictions. Respondents are independent, so the CPU path shards them +//! across scoped threads and reduces diagnostics in respondent order for +//! deterministic results. The item/factor likelihood remains in the existing +//! CPU/GPU Rust kernels; this module owns only the repeated-measurement state +//! layer described by the multilevel RFC. +//! +//! Missing observations are represented by `NaN` and are excluded from fitting +//! while retaining a predicted state at every declared occasion. Time offsets +//! are exact milliseconds at the boundary and are converted to days only for +//! the OLS design matrix. The AR coefficient is a discrete occasion +//! parameter, not a continuous-time decay parameter. + +use std::thread; + +const MILLIS_PER_DAY: f64 = 86_400_000.0; +const MAX_ABS_TIME_DAYS: f64 = 10_000_000.0; +const MAX_AR_SEQUENCE_GAP: usize = i32::MAX as usize; + +/// Result of a validated independent-OLS-trend or stationary-AR state fit. +#[derive(Clone, Debug, PartialEq)] +pub struct LongitudinalStateFit { + /// Predicted latent state aligned with the flattened occasion input. + pub state: Vec, + /// Respondent-level intercept estimates, in respondent input order. + pub intercepts: Vec, + /// Respondent-level slope estimates in outcome units per day. + pub slopes: Vec, + /// The fixed or validated AR coefficient. Zero for intercept/slope fits. + pub ar_coefficient: f64, + /// RMSE over observed values or one-step AR predictions. + pub rmse: f64, + /// Number of finite outcome observations used in fitting. + pub observed_count: usize, + /// Number of observed respondent transitions used by the AR diagnostic. + pub transition_count: usize, +} + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum StateKind { + IndependentRespondentOlsTrend, + StationaryAutoregressive, +} + +#[derive(Clone, Debug)] +struct RespondentFit { + state: Vec, + intercept: f64, + slope: f64, + squared_error: f64, + observed_count: usize, + transition_count: usize, +} + +fn parse_state_kind(value: &str) -> Result { + match value { + "random_intercept_slope" => Ok(StateKind::IndependentRespondentOlsTrend), + "stationary_autoregressive" => Ok(StateKind::StationaryAutoregressive), + _ => Err( + "state_kind must be random_intercept_slope or stationary_autoregressive".to_string(), + ), + } +} + +fn validate_offsets(row_offsets: &[usize], n_values: 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_values) { + return Err("row_offsets must end at the value count".to_string()); + } + Ok(row_offsets.len() - 1) +} + +fn checked_ar_gap(gap: usize) -> Result { + i32::try_from(gap) + .map_err(|_| "accumulated AR occasion gap exceeds the supported range".to_string()) +} + +fn respondent_sequence_span(sequences: &[usize]) -> Result { + let Some(first) = sequences.first().copied() else { + return Ok(0); + }; + sequences[sequences.len() - 1] + .checked_sub(first) + .ok_or_else(|| "sequence indices must increase within the supported AR gap".to_string()) +} + +fn validate_inputs( + row_offsets: &[usize], + sequence_indices: &[usize], + time_offsets_milliseconds: &[i64], + values: &[f64], + state_kind: StateKind, + ar_coefficient: Option, +) -> Result<(usize, f64), String> { + if sequence_indices.len() != values.len() { + return Err("sequence_indices and values must have equal length".to_string()); + } + if time_offsets_milliseconds.len() != values.len() { + return Err("time_offsets_milliseconds and values must have equal length".to_string()); + } + let respondents = validate_offsets(row_offsets, values.len())?; + for row in 0..respondents { + let start = row_offsets[row]; + let end = row_offsets[row + 1]; + let sequences = &sequence_indices[start..end]; + if sequences + .windows(2) + .any(|window| window[1] <= window[0] || window[1] - window[0] > MAX_AR_SEQUENCE_GAP) + { + return Err("sequence indices must increase within the supported AR gap".to_string()); + } + if !sequences.is_empty() { + checked_ar_gap(respondent_sequence_span(sequences)?)?; + } + 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()); + } + } + let phi = match state_kind { + StateKind::IndependentRespondentOlsTrend => { + if ar_coefficient.is_some() { + return Err("random_intercept_slope does not accept an AR coefficient".to_string()); + } + 0.0 + } + StateKind::StationaryAutoregressive => { + let value = ar_coefficient.ok_or_else(|| { + "stationary_autoregressive requires an AR coefficient".to_string() + })?; + if !value.is_finite() || !(-1.0 < value && value < 1.0) { + return Err( + "AR coefficient must be finite and strictly between -1 and 1".to_string(), + ); + } + value + } + }; + for &value in values { + if !value.is_finite() && !value.is_nan() { + return Err("values must be finite or NaN for missing observations".to_string()); + } + } + for &offset in time_offsets_milliseconds { + 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((respondents, phi)) +} + +fn slope_is_identified(denominator: f64, max_abs_deviation: f64) -> bool { + denominator > f64::EPSILON * max_abs_deviation * max_abs_deviation +} + +fn fit_intercept_slope(times: &[i64], values: &[f64]) -> Result { + let first_time = times.first().copied().unwrap_or(0); + let x: Vec = times + .iter() + .map(|value| ((*value - first_time) as f64) / MILLIS_PER_DAY) + .collect(); + let observed: Vec<(f64, f64)> = x + .iter() + .zip(values) + .filter_map(|(time, value)| value.is_finite().then_some((*time, *value))) + .collect(); + if observed.is_empty() { + return Ok(RespondentFit { + state: vec![0.0; values.len()], + intercept: 0.0, + slope: 0.0, + squared_error: 0.0, + observed_count: 0, + transition_count: 0, + }); + } + let count = observed.len() as f64; + let mean_x = observed.iter().map(|(time, _)| time).sum::() / count; + let mean_y = observed.iter().map(|(_, value)| value).sum::() / count; + let denominator = observed + .iter() + .map(|(time, _)| (time - mean_x).powi(2)) + .sum::(); + let max_abs_deviation = observed + .iter() + .map(|(time, _)| (time - mean_x).abs()) + .fold(0.0_f64, f64::max); + let numerator = observed + .iter() + .map(|(time, value)| (time - mean_x) * (value - mean_y)) + .sum::(); + // Time offsets are strictly increasing before this private fitter runs. + // Degeneracy is therefore either an identified intercept-only case + // (fewer than two finite observations) or a genuine scale-relative + // collapse that must fail closed rather than invent a zero slope. + let slope = if observed.len() < 2 { + 0.0 + } else if slope_is_identified(denominator, max_abs_deviation) { + numerator / denominator + } else { + return Err("OLS slope is degenerate relative to the time scale".to_string()); + }; + let intercept = mean_y - slope * mean_x; + let state: Vec = x.iter().map(|time| intercept + slope * time).collect(); + let squared_error = observed + .iter() + .map(|(time, value)| (value - (intercept + slope * time)).powi(2)) + .sum::(); + Ok(RespondentFit { + state, + intercept, + slope, + squared_error, + observed_count: observed.len(), + transition_count: 0, + }) +} + +fn fit_ar(sequence_indices: &[usize], values: &[f64], phi: f64) -> Result { + let mut state = vec![0.0; values.len()]; + let observed: Vec = values + .iter() + .copied() + .filter(|value| value.is_finite()) + .collect(); + let Some(first) = values.iter().position(|value| value.is_finite()) else { + return Ok(RespondentFit { + state, + intercept: 0.0, + slope: 0.0, + squared_error: 0.0, + observed_count: 0, + transition_count: 0, + }); + }; + // Psychometric latent scores use the identified zero-centered origin. The + // state specification has no free intercept, so estimating a respondent + // mean here would reintroduce an unanchored location parameter. + let mean = 0.0; + let mut previous_index = first; + let mut previous_observed = values[first]; + state[first] = previous_observed; + let mut squared_error = 0.0; + let mut transition_count = 0; + for index in (first + 1)..values.len() { + let gap = sequence_indices[index] + .checked_sub(sequence_indices[previous_index]) + .ok_or_else(|| "AR sequence indices must increase".to_string())?; + let exponent = checked_ar_gap(gap)?; + let prediction = mean + phi.powi(exponent) * (previous_observed - mean); + state[index] = prediction; + if values[index].is_finite() { + squared_error += (values[index] - prediction).powi(2); + transition_count += 1; + previous_index = index; + previous_observed = values[index]; + } + } + Ok(RespondentFit { + state, + intercept: mean, + slope: 0.0, + squared_error, + observed_count: observed.len(), + transition_count, + }) +} + +fn map_worker_join(joined: thread::Result) -> Result { + joined.map_err(|_| "longitudinal worker failed".to_string()) +} + +fn require_respondent_fit(fit: Option) -> Result { + fit.ok_or_else(|| "a respondent state fit is missing".to_string()) +} + +/// Fit respondent-level repeated-measurement states on the Rust CPU path. +/// +/// `row_offsets` partitions flattened occasions by respondent. The caller +/// supplies a validated state specification from the Python contract, while +/// this function independently rechecks all array and numeric invariants at +/// the trust boundary. `worker_count` controls deterministic respondent +/// sharding; a value larger than the respondent count is capped. +pub fn fit_longitudinal_state( + row_offsets: &[usize], + sequence_indices: &[usize], + time_offsets_milliseconds: &[i64], + values: &[f64], + state_kind: &str, + ar_coefficient: Option, + worker_count: usize, +) -> Result { + if worker_count == 0 { + return Err("worker_count must be at least one".to_string()); + } + let kind = parse_state_kind(state_kind)?; + let (respondent_count, phi) = validate_inputs( + row_offsets, + sequence_indices, + time_offsets_milliseconds, + values, + kind, + ar_coefficient, + )?; + let mut fits: Vec> = (0..respondent_count).map(|_| None).collect(); + if respondent_count > 0 { + let workers = worker_count.min(respondent_count); + let chunk = respondent_count.div_ceil(workers); + 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(respondent_count); + if start >= end { + continue; + } + handles.push(scope.spawn(move || { + (start..end) + .map(|row| { + let value_start = row_offsets[row]; + let value_end = row_offsets[row + 1]; + let sequences = &sequence_indices[value_start..value_end]; + let times = &time_offsets_milliseconds[value_start..value_end]; + let row_values = &values[value_start..value_end]; + let fit = match kind { + StateKind::IndependentRespondentOlsTrend => { + fit_intercept_slope(times, row_values) + } + StateKind::StationaryAutoregressive => { + fit_ar(sequences, row_values, phi) + } + }?; + Ok((row, fit)) + }) + .collect::, String>>() + })); + } + for handle in handles { + let rows = map_worker_join(handle.join())??; + for (row, fit) in rows { + fits[row] = Some(fit); + } + } + Ok(()) + }); + joined?; + } + let mut state = Vec::with_capacity(values.len()); + let mut intercepts = Vec::with_capacity(respondent_count); + let mut slopes = Vec::with_capacity(respondent_count); + let mut squared_error = 0.0; + let mut observed_count = 0; + let mut transition_count = 0; + for fit in fits { + let fit = require_respondent_fit(fit)?; + state.extend(fit.state); + intercepts.push(fit.intercept); + slopes.push(fit.slope); + squared_error += fit.squared_error; + observed_count += fit.observed_count; + transition_count += fit.transition_count; + } + let denominator = if kind == StateKind::StationaryAutoregressive { + transition_count + } else { + observed_count + }; + Ok(LongitudinalStateFit { + state, + intercepts, + slopes, + ar_coefficient: phi, + rmse: if denominator == 0 { + 0.0 + } else { + (squared_error / denominator as f64).sqrt() + }, + observed_count, + transition_count, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn recovers_two_respondent_intercept_slopes_and_is_worker_deterministic() { + let offsets = [0, 3, 6]; + let sequences = [0, 1, 2, 0, 1, 2]; + let times = [0, 86_400_000, 172_800_000, 0, 86_400_000, 172_800_000]; + let values = [2.0, 3.5, 5.0, -1.0, -3.0, -5.0]; + let one = fit_longitudinal_state( + &offsets, + &sequences, + ×, + &values, + "random_intercept_slope", + None, + 1, + ) + .unwrap(); + let many = fit_longitudinal_state( + &offsets, + &sequences, + ×, + &values, + "random_intercept_slope", + None, + 8, + ) + .unwrap(); + assert_eq!(one, many); + assert_eq!(one.intercepts, vec![2.0, -1.0]); + assert_eq!(one.slopes, vec![1.5, -2.0]); + assert!(one.rmse < 1e-12); + } + + #[test] + fn recovers_noisy_ols_trends_with_bounded_rmse() { + let respondents = 12_usize; + let occasions = 6_usize; + let mut offsets = vec![0]; + let mut sequences = Vec::new(); + let mut times = Vec::new(); + let mut values = Vec::new(); + let mut true_intercepts = Vec::new(); + let mut true_slopes = Vec::new(); + for respondent in 0..respondents { + let intercept = (respondent as f64) * 0.25 - 1.0; + let slope = 0.5 - (respondent as f64) * 0.05; + true_intercepts.push(intercept); + true_slopes.push(slope); + for occasion in 0..occasions { + let days = occasion as f64; + sequences.push(occasion); + times.push((days * MILLIS_PER_DAY) as i64); + let noise = if occasion % 2 == 0 { 0.01 } else { -0.01 }; + values.push(intercept + slope * days + noise); + } + offsets.push(values.len()); + } + let fit = fit_longitudinal_state( + &offsets, + &sequences, + ×, + &values, + "random_intercept_slope", + None, + 4, + ) + .unwrap(); + let intercept_rmse = true_intercepts + .iter() + .zip(&fit.intercepts) + .map(|(truth, estimate)| (truth - estimate).powi(2)) + .sum::() + .sqrt() + / (respondents as f64).sqrt(); + let slope_rmse = true_slopes + .iter() + .zip(&fit.slopes) + .map(|(truth, estimate)| (truth - estimate).powi(2)) + .sum::() + .sqrt() + / (respondents as f64).sqrt(); + assert!(intercept_rmse < 0.02, "{intercept_rmse}"); + assert!(slope_rmse < 0.02, "{slope_rmse}"); + assert!(fit.rmse < 0.02, "{}", fit.rmse); + assert_eq!(fit.observed_count, respondents * occasions); + } + + #[test] + fn ar_state_starts_from_first_finite_observation() { + let fit = fit_longitudinal_state( + &[0, 3], + &[0, 1, 2], + &[0, 86_400_000, 172_800_000], + &[f64::NAN, 1.0, 0.4], + "stationary_autoregressive", + Some(0.4), + 1, + ) + .unwrap(); + assert_eq!(fit.state[0], 0.0); + assert_eq!(fit.state[1], 1.0); + assert!((fit.state[2] - 0.4).abs() < 1e-12); + assert_eq!(fit.observed_count, 2); + assert_eq!(fit.transition_count, 1); + assert!(fit.rmse < 1e-12); + } + + #[test] + fn ar_state_preserves_missing_occasion_and_recovers_prediction_error() { + let offsets = [0, 4]; + let sequences = [0, 1, 2, 3]; + let times = [0, 86_400_000, 172_800_000, 259_200_000]; + let values = [1.0, f64::NAN, 0.25, 0.125]; + let fit = fit_longitudinal_state( + &offsets, + &sequences, + ×, + &values, + "stationary_autoregressive", + Some(0.5), + 2, + ) + .unwrap(); + assert_eq!(fit.ar_coefficient, 0.5); + assert_eq!(fit.observed_count, 3); + assert_eq!(fit.transition_count, 2); + assert_eq!(fit.state[1], 0.5); + assert_eq!(fit.state[2], 0.25); + assert_eq!(fit.state[3], 0.125); + assert!(fit.rmse < 1e-12); + } + + #[test] + fn recovers_caller_supplied_ar_series_with_bounded_rmse() { + let phi = 0.6; + let start = 1.25; + let mut values = vec![start]; + for _ in 0..7 { + let next = phi * values[values.len() - 1]; + values.push(next); + } + let sequences: Vec = (0..values.len()).collect(); + let times: Vec = sequences + .iter() + .map(|step| (*step as i64) * 86_400_000) + .collect(); + let fit = fit_longitudinal_state( + &[0, values.len()], + &sequences, + ×, + &values, + "stationary_autoregressive", + Some(phi), + 1, + ) + .unwrap(); + assert_eq!(fit.ar_coefficient, phi); + assert_eq!(fit.transition_count, values.len() - 1); + assert!(fit.rmse < 1e-12, "{}", fit.rmse); + assert!((fit.state[0] - start).abs() < 1e-12); + } + + #[test] + fn intercept_only_and_all_missing_rows_are_identified() { + let empty = + fit_longitudinal_state(&[0], &[], &[], &[], "random_intercept_slope", None, 3).unwrap(); + assert!(empty.state.is_empty()); + assert_eq!(empty.observed_count, 0); + assert_eq!(empty.rmse, 0.0); + + let intercept_only = fit_longitudinal_state( + &[0, 1, 1], + &[0], + &[0], + &[4.0], + "random_intercept_slope", + None, + 2, + ) + .unwrap(); + assert_eq!(intercept_only.intercepts, vec![4.0, 0.0]); + assert_eq!(intercept_only.slopes, vec![0.0, 0.0]); + assert_eq!(intercept_only.observed_count, 1); + assert_eq!(intercept_only.rmse, 0.0); + + let missing_ar = fit_longitudinal_state( + &[0, 2], + &[0, 1], + &[0, 1], + &[f64::NAN, f64::NAN], + "stationary_autoregressive", + Some(0.3), + 1, + ) + .unwrap(); + assert_eq!(missing_ar.observed_count, 0); + assert_eq!(missing_ar.transition_count, 0); + assert_eq!(missing_ar.rmse, 0.0); + assert_eq!(missing_ar.state, vec![0.0, 0.0]); + } + + #[test] + fn rejects_invalid_contracts_without_panicking() { + let values = [1.0, 2.0]; + let sequences = [0, 1]; + let times = [0, 1]; + for (offsets, kind, phi, message) in [ + (&[1, 2][..], "random_intercept_slope", None, "row_offsets"), + (&[0, 2][..], "stationary_autoregressive", None, "requires"), + ( + &[0, 2][..], + "stationary_autoregressive", + Some(1.0), + "strictly", + ), + ] { + let error = fit_longitudinal_state(offsets, &sequences, ×, &values, kind, phi, 1) + .unwrap_err(); + assert!(error.contains(message), "{error}"); + } + assert!( + fit_longitudinal_state(&[0, 2], &sequences, ×, &values, "unknown", None, 1) + .unwrap_err() + .contains("state_kind") + ); + assert!(fit_longitudinal_state( + &[0, 2], + &sequences, + ×, + &values, + "random_intercept_slope", + None, + 0 + ) + .unwrap_err() + .contains("worker_count")); + assert!(fit_longitudinal_state( + &[0, 2], + &[0, MAX_AR_SEQUENCE_GAP + 1], + ×, + &values, + "stationary_autoregressive", + Some(0.5), + 1 + ) + .unwrap_err() + .contains("supported AR gap")); + assert!(fit_longitudinal_state( + &[0, 2], + &[0], + ×, + &values, + "random_intercept_slope", + None, + 1 + ) + .unwrap_err() + .contains("equal length")); + assert!(fit_longitudinal_state( + &[0, 2], + &sequences, + &[0], + &values, + "random_intercept_slope", + None, + 1 + ) + .unwrap_err() + .contains("time_offsets_milliseconds")); + assert!(validate_offsets(&[0, 1, 0], 0) + .unwrap_err() + .contains("non-decreasing")); + assert!(validate_offsets(&[0, 1], 2) + .unwrap_err() + .contains("end at the value count")); + assert!(fit_longitudinal_state( + &[0, 2], + &[1, 0], + ×, + &values, + "random_intercept_slope", + None, + 1 + ) + .unwrap_err() + .contains("sequence indices")); + assert!(fit_longitudinal_state( + &[0, 2], + &sequences, + &[2, 1], + &values, + "random_intercept_slope", + None, + 1 + ) + .unwrap_err() + .contains("time offsets")); + assert!(fit_longitudinal_state( + &[0, 2], + &sequences, + ×, + &values, + "random_intercept_slope", + Some(0.1), + 1 + ) + .unwrap_err() + .contains("does not accept")); + assert!(fit_longitudinal_state( + &[0, 2], + &sequences, + ×, + &values, + "stationary_autoregressive", + Some(f64::NAN), + 1 + ) + .unwrap_err() + .contains("strictly")); + assert!(fit_longitudinal_state( + &[0, 2], + &sequences, + ×, + &[1.0, f64::INFINITY], + "random_intercept_slope", + None, + 1 + ) + .unwrap_err() + .contains("finite or NaN")); + assert!(fit_longitudinal_state( + &[0, 2], + &sequences, + &[0, 86_400_000_i64.saturating_mul(10_000_001 * 2)], + &values, + "random_intercept_slope", + None, + 1 + ) + .unwrap_err() + .contains("supported finite range")); + } + + #[test] + fn unused_worker_chunk_is_skipped_without_changing_estimates() { + // Four respondents and three workers yield chunk=2, so the third + // worker starts at index 4 and is skipped. The remaining shards still + // produce the independent intercept-only estimates. + let offsets = [0, 1, 2, 3, 4]; + let sequences = [0, 0, 0, 0]; + let times = [0, 0, 0, 0]; + let values = [1.0, 2.0, 3.0, 4.0]; + let fit = fit_longitudinal_state( + &offsets, + &sequences, + ×, + &values, + "random_intercept_slope", + None, + 3, + ) + .unwrap(); + assert_eq!(fit.intercepts, vec![1.0, 2.0, 3.0, 4.0]); + assert_eq!(fit.slopes, vec![0.0, 0.0, 0.0, 0.0]); + assert_eq!(fit.observed_count, 4); + assert_eq!(fit.rmse, 0.0); + } + + #[test] + fn package_owned_join_and_missing_fit_errors_are_stable() { + let join_error = map_worker_join::<()>(Err(Box::new("boom"))); + assert_eq!(join_error.unwrap_err(), "longitudinal worker failed"); + assert_eq!( + require_respondent_fit(None).unwrap_err(), + "a respondent state fit is missing" + ); + assert!(checked_ar_gap(MAX_AR_SEQUENCE_GAP + 1) + .unwrap_err() + .contains("accumulated AR occasion gap")); + assert_eq!(respondent_sequence_span(&[]).unwrap(), 0); + assert_eq!(respondent_sequence_span(&[5]).unwrap(), 0); + assert_eq!(respondent_sequence_span(&[0, 3]).unwrap(), 3); + assert!(respondent_sequence_span(&[5, 2]) + .unwrap_err() + .contains("sequence indices must increase")); + assert!(!slope_is_identified(0.0, 0.0)); + assert!(slope_is_identified(1e-16, 1e-8)); + } + + #[test] + fn fit_ar_rejects_decreasing_sequence_after_validation_bypass() { + let error = fit_ar(&[2, 1], &[1.0, 0.5], 0.4).unwrap_err(); + assert!( + error.contains("AR sequence indices must increase"), + "{error}" + ); + } + + #[test] + fn intercept_slope_fails_closed_on_scale_relative_degeneracy() { + let error = fit_intercept_slope(&[0, 0], &[1.0, 2.0]).unwrap_err(); + assert!(error.contains("degenerate"), "{error}"); + let empty = fit_intercept_slope(&[], &[]).unwrap(); + assert_eq!(empty.observed_count, 0); + assert_eq!(empty.slope, 0.0); + } +} diff --git a/crates/mlsirm-core/src/longitudinal_irt.rs b/crates/mlsirm-core/src/longitudinal_irt.rs new file mode 100644 index 000000000..6f3c1c73f --- /dev/null +++ b/crates/mlsirm-core/src/longitudinal_irt.rs @@ -0,0 +1,1677 @@ +//! 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] = if params[1] < MIN_LOG_SD || params[1] > MAX_LOG_SD { + 0.0 + } else { + reduced.log_sd_grad + }; + grad[2] = if params[2] < MIN_LOG_DECAY || params[2] > MAX_LOG_DECAY { + 0.0 + } else { + 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; + if params[1] - step <= MIN_LOG_SD + || params[1] + step >= MAX_LOG_SD + || params[2] - step <= MIN_LOG_DECAY + || params[2] + step >= MAX_LOG_DECAY + { + return Err( + "hyperparameter Hessian is not identified at the supported log-scale boundary" + .to_string(), + ); + } + 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 clamped_hyperparameters_have_flat_raw_gradients_and_unidentified_boundary_hessian() { + let offsets = [0, 3]; + let times = [0, 86_400_000, 172_800_000]; + let responses = [1.0, 0.0, 0.0, 1.0, 1.0, 0.0]; + let base = 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, 0.1], + }; + let params = pack(&base); + for (index, first_raw, second_raw) in [ + (1, MAX_LOG_SD + 1.0, MAX_LOG_SD + 2.0), + (1, MIN_LOG_SD - 1.0, MIN_LOG_SD - 2.0), + (2, MAX_LOG_DECAY + 1.0, MAX_LOG_DECAY + 2.0), + (2, MIN_LOG_DECAY - 1.0, MIN_LOG_DECAY - 2.0), + ] { + let mut first = params.clone(); + first[index] = first_raw; + let (first_nll, first_grad, _, _, _) = + joint_objective(&offsets, ×, &responses, 2, &first, 1).unwrap(); + let mut second = params.clone(); + second[index] = second_raw; + let (second_nll, _, _, _, _) = + joint_objective(&offsets, ×, &responses, 2, &second, 1).unwrap(); + assert!((first_nll - second_nll).abs() < 1e-12); + assert_eq!(first_grad[index], 0.0); + } + + for (index, boundary_near) in [ + (1, MAX_LOG_SD - 0.5e-3), + (2, MIN_LOG_DECAY + 0.5e-3), + ] { + let mut near_boundary = params.clone(); + near_boundary[index] = boundary_near; + let err = hyperparameter_hessian( + &offsets, + ×, + &responses, + 2, + &near_boundary, + 1, + 1e-3, + ) + .unwrap_err(); + assert!(err.contains("supported log-scale boundary")); + let unpacked = unpack(&near_boundary, 2, 3).unwrap(); + let (se, lower, upper, identified) = hyperparameter_interval_estimates( + &offsets, + ×, + &responses, + 2, + &near_boundary, + &unpacked, + 1, + 1e-3, + ) + .unwrap(); + assert!(!identified); + assert!(se.iter().all(|value| value.is_nan())); + assert!(lower.iter().all(|value| value.is_nan())); + assert!(upper.iter().all(|value| value.is_nan())); + } + } + + #[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/crates/mlsirm-core/tests/longitudinal_fail_closed.rs b/crates/mlsirm-core/tests/longitudinal_fail_closed.rs new file mode 100644 index 000000000..1a506e3e0 --- /dev/null +++ b/crates/mlsirm-core/tests/longitudinal_fail_closed.rs @@ -0,0 +1,17 @@ +use mlsirm_core::longitudinal::fit_longitudinal_state; + +#[test] +fn rejects_accumulated_ar_gap_beyond_i32_exponent_range() { + let max_gap = i32::MAX as usize; + let error = fit_longitudinal_state( + &[0, 3], + &[0, max_gap, max_gap + 1], + &[0, 1, 2], + &[1.0, f64::NAN, 0.25], + "stationary_autoregressive", + Some(0.5), + 1, + ) + .unwrap_err(); + assert!(error.contains("accumulated AR occasion gap"), "{error}"); +} diff --git a/crates/mlsirm-core/tests/longitudinal_short_interval.rs b/crates/mlsirm-core/tests/longitudinal_short_interval.rs new file mode 100644 index 000000000..e2074d99b --- /dev/null +++ b/crates/mlsirm-core/tests/longitudinal_short_interval.rs @@ -0,0 +1,21 @@ +use mlsirm_core::longitudinal::fit_longitudinal_state; + +#[test] +fn recovers_one_millisecond_respondent_trend() { + let fit = fit_longitudinal_state( + &[0, 2], + &[0, 1], + &[0, 1], + &[2.0, 3.0], + "random_intercept_slope", + None, + 1, + ) + .unwrap(); + + assert!((fit.intercepts[0] - 2.0).abs() < 1e-12); + assert!((fit.slopes[0] - 86_400_000.0).abs() < 1e-3); + assert!((fit.state[0] - 2.0).abs() < 1e-12); + assert!((fit.state[1] - 3.0).abs() < 1e-12); + assert!(fit.rmse < 1e-12); +} 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 1b65624d4..63d525645 100644 --- a/docs/adr/0007-multilevel-multiple-membership-temporal.md +++ b/docs/adr/0007-multilevel-multiple-membership-temporal.md @@ -7,7 +7,12 @@ 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. -A current open PR contains reusable contract work for nested, cross-classified, multiple-membership and longitudinal designs, but it is not yet protected-integrated. Numerical estimators for the full structures are not accepted production behavior. Therefore this ADR remains Proposed. +Reusable contracts plus a Rust-owned respondent state layer define a supported +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 @@ -23,6 +28,19 @@ The architecture treats the following as distinct, explicit structures: - future continuous-time state transitions; - rater/model/prompt drift. +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. 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`. @@ -61,6 +79,8 @@ This architecture avoids forcing product-specific tenant/org structures into the - governed contracts merged to protected main; - architecture/serialization tests pass; - at least one Rust estimator or clear handoff contract exists for a supported multilevel/temporal inference use case; +- the state-layer recovery fixture reports slope recovery, missing-occasion + behavior, AR transition RMSE, and equality across worker counts; - recovery evidence meets the numerical release rule. ## References diff --git a/docs/adr/0018-rust-longitudinal-state-engine.md b/docs/adr/0018-rust-longitudinal-state-engine.md new file mode 100644 index 000000000..69e5582dd --- /dev/null +++ b/docs/adr/0018-rust-longitudinal-state-engine.md @@ -0,0 +1,145 @@ +# ADR-0018: Rust-owned longitudinal state layer + +Status: **Proposed** +Date: 2026-08-17 +Supersedes: none +Superseded by: none + +## Context + +The sealed `fast_mlsirm.multilevel` contracts preserve respondent identity, +sequence, exact time offsets, revision provenance, weighted contexts, and +missing occasions. A contract without a numerical consumer still leaves a +product integration gap: a report or evaluation workflow cannot distinguish a +stable respondent level from a time-varying state. At the same time, a full +multilevel IRT estimator would introduce a much larger identification, +uncertainty, GPU, and recovery surface. + +Live protected main already occupies ADR-0015 (multi-item IRT fit boundary). +PR #948 independently records Angoff delta-plot as ADR-0016 and Bradley–Terry +Hunter MM as ADR-0017. This decision therefore uses ADR-0018 so the +longitudinal state layer does not collide with those citation records. + +## Ownership and dependency direction + +`ContextualWisdomLab/fast-mlsirm` owns this reusable measurement state layer. +Time-flow / longitudinal state arithmetic belongs here, not in `kaefa`, +`aFIPC`, or Psychometrics Commons. Downstream products consume the sealed +design and returned diagnostics; they do not reimplement the estimator. + +## Decision + +Add a narrow Rust-owned state layer behind the existing Python contract and +PyO3 boundary: + +1. `random_intercept_slope` is retained as the state-specification wire label, + but the fitted estimand is an **independent per-respondent OLS trend**. Each + respondent's finite observations are fitted by ordinary least squares on + exact millisecond offsets converted to days. With + `x_pi = (t_pi - t_p1) / 86,400,000`, the state is + \(\hat\eta_{pi}=\hat\alpha_p+\hat\beta_p x_{pi}\). A respondent with fewer + than two distinct observed times receives a zero slope rather than an + invented trend. A scale-relative degeneracy with two or more observations + fails closed. This path does not estimate a population random-effects + distribution and applies no multilevel shrinkage. +2. `stationary_autoregressive` is likewise a state-specification wire label for + a **discrete AR state predictor**. It uses the sealed caller-supplied + `-1 < phi < 1` value and actual `sequence_index` gaps. Missing observations + remain in the output state and do not reset the last observed state. The + coefficient is a discrete-occasion parameter; this layer does not estimate + \(\phi\) or its uncertainty. Cumulative gaps use checked `i32` conversion + before exponentiation. Exact calendar offsets are retained for audit and + are not silently converted into continuous-time decay. +3. Respondents are independently sharded over scoped Rust CPU threads. Results + are reduced in respondent order, making worker-count changes deterministic. + Worker join failures are converted to the stable package-owned Rust error + `longitudinal worker failed` rather than unwinding across the PyO3 boundary. +4. PyO3 performs only bounded array marshalling, detaching the numeric fit + from the Python GIL. Python exposes the returned state, intercepts, slopes, + RMSE, counts, engine identity, and aligned respondent/occasion identifiers. + +The Python result metadata is normative for interpretation. The OLS path emits +`estimand_scope="independent_respondent_ols_trend"`, +`population_random_effects_estimated=False`, and +`ar_coefficient_source="not_applicable"`. The AR path emits +`estimand_scope="discrete_ar_state_prediction"`, +`population_random_effects_estimated=False`, `ar_coefficient_estimated=False`, +and `ar_coefficient_source="caller_supplied"`. Consumers must use these fields +rather than infer an estimand from the compatibility wire label alone. + +The weighted multiple-membership contextual predictor remains a separate Rust +kernel. This ADR does not claim Bayesian random-effect integration, joint item +and context likelihood estimation, interval uncertainty, continuous-time +transitions, or GPU parity for the recurrent state recurrence. Those require +their own equations, identification rules, and recovery evidence. + +## Invariants / acceptance evidence + +1. Independent OLS intercepts and slopes recover known generating parameters + with bounded RMSE; worker counts do not change the result. +2. Caller-supplied discrete AR predictions recover a known series with bounded + one-step RMSE and preserve missing occasions. +3. Invalid worker counts, foreign designs, hostile observation mappings, + non-real values, and unsupported AR gaps fail with package-owned errors. +4. Closely spaced valid occasions keep a finite slope; genuine scale-relative + degeneracy fails closed. + +## Non-goals and claims not made + +- No population random-effects distribution or shrinkage. +- No AR coefficient or uncertainty estimation. +- No continuous-time or interval-adjusted transitions. +- No joint multilevel IRT likelihood or GPU recurrent-state parity. +- No hosted product session, consent, or persistence behavior. + +## Consequences + +LineageWeave and other consumers can call one provider-neutral state estimator +without copying arithmetic into Python or treating missing occasions as absent +rows. A single-period report group can be represented honestly as a state +observation, but it cannot be reported as evidence of temporal trend. Consumers +must persist the state specification, design fingerprint, engine version, and +RMSE/count diagnostics before interpreting a score. + +The current state layer is intentionally not a replacement for fast-mlsirm's +full psychometric calibration path. Production promotion remains gated on +protected integration, true-parameter recovery across nested/crossed/weighted +membership and unbalanced temporal designs, interval coverage, and GPU/CPU +parity where a GPU state implementation is justified. + +## Verification + +- Rust unit tests recover two respondent slopes exactly and compare one and + many worker counts. +- Rust recovery fixtures report intercept/slope RMSE and AR prediction RMSE + against known generating parameters. +- Rust/Python tests preserve a missing occasion and reject malformed offsets, + state kinds, coefficients, and worker counts. +- Rust regression coverage preserves valid slopes at millisecond-scale time + intervals rather than treating a small centered sum of squares as a missing + trend. +- Python tests use a non-contiguous sequence gap and verify that AR prediction + uses the declared discrete gap rather than array position or calendar days. +- The PyO3 extension is built from the root Maturin project and exposes + `fit_longitudinal_state` through the existing multilevel loader. + +## Research basis (APA 7) + +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 + +Embretson, S. E. (1991). A multidimensional latent trait model for measuring +learning and change. *Psychometrika, 56*, 495–515. +https://doi.org/10.1007/BF02294487 + +Fox, J.-P., & Glas, C. A. W. (2001). Bayesian estimation of a multilevel IRT +model using Gibbs sampling. *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 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 8575be25d..27a2ce988 100644 --- a/docs/adr/README.md +++ b/docs/adr/README.md @@ -32,6 +32,8 @@ A conversation, issue, PR body, design note, or paper summary is not an Accepted | [0013](0013-continuous-execution-and-documentation-governance.md) | Proposed | Keep autonomous work work-conserving and enforce one canonical cross-cutting documentation writer with explicit maturity states. | | [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/848-rust-longitudinal-state-engine.md b/docs/changelog.d/848-rust-longitudinal-state-engine.md new file mode 100644 index 000000000..253bf34ee --- /dev/null +++ b/docs/changelog.d/848-rust-longitudinal-state-engine.md @@ -0,0 +1,15 @@ +# Rust longitudinal state layer + +## Added + +- 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. 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 7499558d7..30cbb68cc 100644 --- a/docs/doctoring/multilevel_longitudinal_measurement.md +++ b/docs/doctoring/multilevel_longitudinal_measurement.md @@ -5,13 +5,23 @@ `fast_mlsirm.multilevel` introduces immutable design contracts for nested, 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. - -This first slice performs no statistical estimation. Python owns validation, -content identity, bounded collection handling, replay protection, sparse design -marshalling, and serialization only. Likelihood, integration, gradients, -optimization, uncertainty, CPU multithreading, and any justified GPU batching -remain in Rust. +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. 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, 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 @@ -63,7 +73,14 @@ The implementation: - canonicalizes input order without changing assignments; - retains exact membership and occasion revision fingerprints; - requires strict respondent-level sequence and time ordering; -- distinguishes a random-intercept/slope state from stationary AR(1); +- 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; - bounds all collections before aggregate allocation; - rejects Boolean-as-number coercion, non-finite values, duplicate cells, and @@ -73,17 +90,17 @@ 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 +discrete AR path must not be described as continuous-time or interval-adjusted; +its coefficient is tied to sequence gaps. The OLS path uses exact day-scaled +offsets and does not estimate a continuous-time transition. ## Identification and interpretation limits @@ -111,18 +128,23 @@ must use explicit schema migration rather than mutating the accepted contract. ## Verification boundary -The current evidence is limited to contract validation, deterministic identity, +The current evidence includes contract validation, deterministic identity, child replay, resource bounds, dimension-scoped assignment, strict temporal -ordering, and source-text-free serialization. It is not evidence of: - -- estimator correctness; -- variance-component identification; -- true-parameter recovery; -- interval coverage; +ordering, Rust state-layer slope recovery, missing-occasion preservation, +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 performance or parity; or +- GPU recurrent-state performance or parity; or - high-stakes deployment readiness. Those claims require separate Rust implementations and same-head recovery, @@ -153,6 +175,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 61be2ca2d..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 | PARTIAL | governed contracts/design primitives exist; continuous-time or richer estimator claims require separate recovery evidence | +| 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/multilevel_multiple_membership_longitudinal_rfc.md b/docs/multilevel_multiple_membership_longitudinal_rfc.md index 23acff62e..87d4cc334 100644 --- a/docs/multilevel_multiple_membership_longitudinal_rfc.md +++ b/docs/multilevel_multiple_membership_longitudinal_rfc.md @@ -2,9 +2,12 @@ ## Status -This RFC defines the provider-neutral contract boundary delivered by issue #565. -It does not yet add a numerical estimator. All future psychometric arithmetic -remains a Rust-core responsibility. +This RFC defines the provider-neutral contract boundary delivered by issue #565 +and the first Rust-owned repeated-measurement state layer. The state layer is a +small, identified handoff for independent per-respondent OLS trends and +discrete-step AR(1) state predictions; it is not a claim that the full +multilevel IRT likelihood, Bayesian random-effect estimation, uncertainty, or +GPU state recurrence is complete. ## Product outcome @@ -83,13 +86,21 @@ and revision identity. Within a respondent: - time offsets increase strictly with sequence order; - irregular spacing is preserved as provenance. -The initial state specifications are: +The initial state-specification wire labels are: - `random_intercept_slope`; - `stationary_autoregressive` with \(-1<\phi<1\). -The current `autoregressive_coefficient` is a **discrete occasion-step AR(1) -parameter**. The millisecond offsets do not transform \(\phi\), and the contract +These labels are compatibility identifiers, not claims about the fitted +estimand. A `random_intercept_slope` result reports +`estimand_scope="independent_respondent_ols_trend"` and +`population_random_effects_estimated=False`; there is no population variance +component or shrinkage in this state predictor. A `stationary_autoregressive` +result reports `estimand_scope="discrete_ar_state_prediction"`, +`ar_coefficient_estimated=False`, and `ar_coefficient_source="caller_supplied"`. +The current `autoregressive_coefficient` is therefore a **discrete occasion-step +AR(1) parameter supplied by the caller**, not a coefficient estimated by this +state layer. The millisecond offsets do not transform \(\phi\), and the contract does not claim continuous-time or interval-adjusted transitions. A later Rust estimator may use irregular gaps only after a separate, explicit continuous-time or elapsed-gap parameterization and recovery contract is introduced. @@ -116,16 +127,17 @@ use `NVIDIA_NIM_API_KEY` rather than `COPILOT_GITHUB_TOKEN`. ## Numerical boundary Python performs validation, canonicalization, hashing, sparse design marshalling, -and serialization. Future Rust PRs own: - -- multilevel and cross-classified predictors; -- random-effect integration; -- multiple-membership weighting; -- longitudinal state transitions; -- likelihood and gradients; -- optimization and uncertainty; -- CPU multithreading and justified GPU batching; -- true-parameter recovery. +and serialization. Rust owns: + +- the multilevel and cross-classified weighted predictor; +- independent per-respondent OLS state estimates; +- discrete-step stationary AR(1) state prediction using caller-supplied \(\phi\); +- deterministic CPU respondent sharding and diagnostics for those state paths. + +The following remain explicit future boundaries: full multilevel IRT random-effect +integration and estimation, uncertainty/intervals, joint item/context/state +likelihood and gradients, GPU batching for the recurrent state path, continuous +time transitions, and true-parameter recovery for the full joint model. A Python fallback estimator is explicitly out of scope. 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 b15ee0540..bcaa3b866 100644 --- a/docs/traceability/research-basis.md +++ b/docs/traceability/research-basis.md @@ -134,6 +134,20 @@ Primary basis: - Fox, J.-P., & Glas, C. A. W. (2001). Bayesian estimation of a multilevel IRT model. *Psychometrika, 66*, 271–288. - Uto, M. (2022). A Bayesian many-facet Rasch model with Markov modeling for rater severity drift. *Behavior Research Methods, 55*, 3910–3928. +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 Architecture effect: 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 214bb1104..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 weighted_contextual_effect +from .estimation import ( + fit_hierarchical_longitudinal_irt, + fit_longitudinal_state, + simulate_hierarchical_longitudinal_irt, + weighted_contextual_effect, +) __all__ = [ "ContextMembership", @@ -29,5 +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/contracts.py b/python/fast_mlsirm/multilevel/contracts.py index 84bea83f4..8dc471367 100644 --- a/python/fast_mlsirm/multilevel/contracts.py +++ b/python/fast_mlsirm/multilevel/contracts.py @@ -35,7 +35,12 @@ class LongitudinalStateKind(str, Enum): - """Supported initial latent-state structures for repeated measurement.""" + """Supported repeated-measurement states and their compatibility labels. + + ``RANDOM_INTERCEPT_SLOPE`` is retained as a wire-compatible input label, + while the current implementation reports its estimand as an independent + respondent OLS trend rather than a population random-effects fit. + """ RANDOM_INTERCEPT_SLOPE = "random_intercept_slope" STATIONARY_AUTOREGRESSIVE = "stationary_autoregressive" diff --git a/python/fast_mlsirm/multilevel/estimation.py b/python/fast_mlsirm/multilevel/estimation.py index e0e6b3cba..a988e7c8c 100644 --- a/python/fast_mlsirm/multilevel/estimation.py +++ b/python/fast_mlsirm/multilevel/estimation.py @@ -1,26 +1,77 @@ -"""Typed Python access to the Rust-native contextual-effects predictor. +"""Typed Python access to the Rust-native contextual and longitudinal predictors. This module performs marshalling only: converting a validated ``ContextMembershipDesign`` (see ``fast_mlsirm.multilevel.contracts``) into the flat CSR arrays ``mlsirm_core::multilevel::weighted_contextual_effect`` -expects, and converting the caller's per-context random-effect values into -the matching flat vector. The additive sum, its determinism across worker -counts, and its numerical input validation are owned by the Rust core; see -that module's docstring for the full linear-predictor context and the -Browne, Goldstein, and Rasbash (2001) citation. +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 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 from .._multilevel_core_loader import multilevel_core -from .contracts import ContextMembershipDesign +from .contracts import ContextMembershipDesign, LongitudinalDesign, LongitudinalStateKind ContextKey = tuple[str, str] +_NUMPY_INTEGER_SCALAR_TYPES = ( + np.int8, + np.int16, + np.int32, + np.int64, + np.intp, + np.longlong, + np.uint8, + np.uint16, + np.uint32, + np.uint64, + np.uintp, + np.ulonglong, +) +_NUMPY_REAL_SCALAR_TYPES = _NUMPY_INTEGER_SCALAR_TYPES + ( + np.float16, + np.float32, + np.float64, + np.longdouble, +) + + +def _trusted_positive_integer(value: object, name: str) -> int: + """Return one positive execution integer without caller coercion callbacks.""" + value_type = type(value) + if value_type is int: + normalized = value + elif any(value_type is trusted_type for trusted_type in _NUMPY_INTEGER_SCALAR_TYPES): + normalized = int(value) + else: + raise ValueError(f"{name} must be an integer") + if normalized < 1: + raise ValueError(f"{name} must be at least one") + return normalized + + +def _trusted_positive_real(value: object, name: str) -> float: + """Return one positive finite execution real without caller coercion callbacks.""" + value_type = type(value) + if value_type is int or value_type is float: + normalized = float(value) + elif any(value_type is trusted_type for trusted_type in _NUMPY_REAL_SCALAR_TYPES): + normalized = float(value) + else: + raise ValueError(f"{name} must be a real number") + if not np.isfinite(normalized) or normalized <= 0.0: + raise ValueError(f"{name} must be finite and strictly positive") + return normalized + def _snapshot_context_effects( context_keys: tuple[ContextKey, ...], @@ -127,4 +178,393 @@ def weighted_contextual_effect( ) -__all__ = ["ContextKey", "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: + raw = values.get(occasion_id, np.nan) + except Exception: + raise ValueError("values must be a plain read-only mapping") from None + if isinstance(raw, (bool, np.bool_)) or not isinstance( + raw, (int, float, np.integer, np.floating) + ): + raise ValueError(f"values[{occasion_id!r}] must be a real number") + return float(raw) + + +def fit_longitudinal_state( + design: LongitudinalDesign, + values: Mapping[str, float], + *, + worker_count: int = 1, +) -> dict[str, object]: + """Fit the Rust-owned respondent state predictor for a sealed design. + + ``values`` maps exact occasion identifiers to observed factor scores. A + missing identifier is represented as ``NaN`` so the design remains intact + while the Rust fitter excludes that observation from estimation. The + returned state is aligned with ``design.occasions`` sorted by respondent + and sequence, and includes normative estimand metadata so callers cannot + mistake independent respondent OLS trends for population random effects or + a caller-supplied AR coefficient for an estimated parameter. + + Parameters + ---------- + design: + A package-built ``LongitudinalDesign``. Integrity is verified before + marshalling so a tampered or hand-constructed design raises here. + values: + Mapping from occasion identifier to a real observed state. Absent + keys become ``NaN``. Boolean, complex, and non-numeric values are + rejected before native dispatch. + worker_count: + Number of deterministic worker threads (``>= 1``). The numerical + result does not depend on this value. + + Returns + ------- + dict + Predicted states, respondent intercepts and slopes, RMSE, counts, + engine identity, fingerprints, and normative estimand metadata. + + Raises + ------ + ValueError + If ``worker_count < 1``, the design is not an exact sealed + ``LongitudinalDesign``, a caller observation cannot be read or + converted safely, or the Rust-side state contract is invalid. + """ + if worker_count < 1: + raise ValueError("worker_count must be at least one") + 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) + 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: + ar_coefficient = None + estimand_scope = "independent_respondent_ols_trend" + ar_coefficient_source = "not_applicable" + else: + estimand_scope = "discrete_ar_state_prediction" + ar_coefficient_source = "caller_supplied" + core = multilevel_core() + result = core.fit_longitudinal_state( + np.asarray(row_offsets, dtype=np.uint64), + np.asarray(sequence_indices, dtype=np.uint64), + np.asarray(time_offsets, dtype=np.int64), + np.asarray(observations, dtype=np.float64), + state_kind.value, + ar_coefficient, + worker_count, + ) + return { + "state_kind": state_kind.value, + "estimand_scope": estimand_scope, + "population_random_effects_estimated": False, + "ar_coefficient_estimated": False, + "ar_coefficient_source": ar_coefficient_source, + "state_spec_fingerprint": design.state_spec.state_spec_fingerprint, + "design_fingerprint": design.design_fingerprint, + "state": np.asarray(result["state"], dtype=np.float64), + "intercepts": np.asarray(result["intercepts"], dtype=np.float64), + "slopes": np.asarray(result["slopes"], dtype=np.float64), + "ar_coefficient": float(result["ar_coefficient"]), + "rmse": float(result["rmse"]), + "observed_count": int(result["observed_count"]), + "transition_count": int(result["transition_count"]), + "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 _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. + """ + worker_count = _trusted_positive_integer(worker_count, "worker_count") + max_iter = _trusted_positive_integer(max_iter, "max_iter") + tolerance = _trusted_positive_real(tolerance, "tolerance") + hessian_step = _trusted_positive_real(hessian_step, "hessian_step") + 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, + tolerance, + 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", +] \ No newline at end of file diff --git a/tests/test_hierarchical_longitudinal_control_boundary.py b/tests/test_hierarchical_longitudinal_control_boundary.py new file mode 100644 index 000000000..28b29bf80 --- /dev/null +++ b/tests/test_hierarchical_longitudinal_control_boundary.py @@ -0,0 +1,75 @@ +"""Trust-boundary regressions for hierarchical longitudinal execution controls.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from fast_mlsirm.multilevel import fit_hierarchical_longitudinal_irt + + +class _HostileOrder: + """Control object whose comparison callback must never execute.""" + + def __init__(self) -> None: + self.calls = 0 + + def __lt__(self, other: object) -> bool: + self.calls += 1 + raise AssertionError("caller comparison callback executed") + + +class _HostileFloat: + """Control object whose float conversion callback must never execute.""" + + def __init__(self) -> None: + self.calls = 0 + + def __float__(self) -> float: + self.calls += 1 + raise AssertionError("caller float callback executed") + + +@pytest.mark.parametrize("name", ["worker_count", "max_iter"]) +def test_integer_execution_controls_fail_closed_before_callbacks(name: str) -> None: + """Integer controls reject alien ordering protocols with package errors.""" + value = _HostileOrder() + with pytest.raises(ValueError): + fit_hierarchical_longitudinal_irt( + object(), + np.zeros((1, 2), dtype=np.float64), + **{name: value}, + ) + assert value.calls == 0 + + +@pytest.mark.parametrize("name", ["tolerance", "hessian_step"]) +def test_real_execution_controls_fail_closed_before_callbacks(name: str) -> None: + """Real controls reject alien conversion protocols with package errors.""" + value = _HostileFloat() + with pytest.raises(ValueError): + fit_hierarchical_longitudinal_irt( + object(), + np.zeros((1, 2), dtype=np.float64), + **{name: value}, + ) + assert value.calls == 0 + + +@pytest.mark.parametrize( + ("name", "value"), + [ + ("worker_count", "1"), + ("max_iter", "250"), + ("tolerance", "1e-5"), + ("hessian_step", "1e-3"), + ], +) +def test_nonnumeric_execution_controls_raise_value_error(name: str, value: str) -> None: + """Documented execution-control failures normalize to ValueError.""" + with pytest.raises(ValueError): + fit_hierarchical_longitudinal_irt( + object(), + np.zeros((1, 2), dtype=np.float64), + **{name: value}, + ) diff --git a/tests/test_hierarchical_longitudinal_irt.py b/tests/test_hierarchical_longitudinal_irt.py new file mode 100644 index 000000000..e73adbf7c --- /dev/null +++ b/tests/test_hierarchical_longitudinal_irt.py @@ -0,0 +1,297 @@ +"""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) + with pytest.raises( + ValueError, match="at least one respondent must have two or more occasions" + ): + fit_hierarchical_longitudinal_irt( + _design(1, 1), + np.array([[1.0, 0.0]], dtype=np.float64), + ) + + 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] + ) + + +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_longitudinal_state_estimation.py b/tests/test_longitudinal_state_estimation.py new file mode 100644 index 000000000..325b30fb3 --- /dev/null +++ b/tests/test_longitudinal_state_estimation.py @@ -0,0 +1,247 @@ +"""End-to-end tests for the Rust longitudinal state estimator.""" + +from __future__ import annotations + +from collections.abc import Mapping + +import numpy as np +import pytest + +from fast_mlsirm.multilevel import ( + LongitudinalStateKind, + build_longitudinal_design, + build_longitudinal_state_spec, + build_temporal_occasion, + fit_longitudinal_state, +) + + +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 _single_occasion_design(): + """Return the smallest sealed design for public-boundary validation tests.""" + return build_longitudinal_design( + occasions=[_occasion("respondent_a", "guard", 0, 0)], + state_spec=build_longitudinal_state_spec( + state_kind=LongitudinalStateKind.RANDOM_INTERCEPT_SLOPE, + ), + ) + + +def test_rust_state_fit_recovers_slopes_and_missing_values() -> None: + """A two-person synthetic recovery fixture exercises multithreading.""" + occasions = [ + _occasion("respondent_a", "a0", 0, 0), + _occasion("respondent_a", "a1", 1, 86_400_000), + _occasion("respondent_a", "a2", 2, 172_800_000), + _occasion("respondent_b", "b0", 0, 0), + _occasion("respondent_b", "b1", 1, 86_400_000), + _occasion("respondent_b", "b2", 2, 172_800_000), + ] + design = build_longitudinal_design( + occasions=occasions, + state_spec=build_longitudinal_state_spec( + state_kind=LongitudinalStateKind.RANDOM_INTERCEPT_SLOPE, + ), + ) + values = { + "occasion_a0": 2.0, + "occasion_a1": 3.5, + "occasion_a2": 5.0, + "occasion_b0": -1.0, + "occasion_b1": -3.0, + } + result = fit_longitudinal_state(design, values, worker_count=4) + single = fit_longitudinal_state(design, values, worker_count=1) + np.testing.assert_allclose(result["intercepts"], [2.0, -1.0], atol=1e-12) + np.testing.assert_allclose(result["slopes"], [1.5, -2.0], atol=1e-12) + np.testing.assert_array_equal(single["state"], result["state"]) + assert single["rmse"] == result["rmse"] + assert result["observed_count"] == 5 + assert result["engine"] == "rust_cpu_multithreaded" + assert result["state_kind"] == "random_intercept_slope" + assert result["estimand_scope"] == "independent_respondent_ols_trend" + assert result["population_random_effects_estimated"] is False + assert result["ar_coefficient_estimated"] is False + assert result["ar_coefficient_source"] == "not_applicable" + assert len(result["design_fingerprint"]) == 64 + assert result["occasion_records"][0]["sequence_index"] == 0 + + +def test_rust_state_fit_recovers_true_ols_parameters_with_rmse() -> None: + """Known intercepts and slopes are recovered with bounded RMSE.""" + true_intercepts = (1.25, -0.5, 0.0) + true_slopes = (0.75, -1.0, 0.25) + occasions = [] + values: dict[str, float] = {} + for respondent_index, respondent in enumerate( + ("respondent_r0", "respondent_r1", "respondent_r2") + ): + intercept = true_intercepts[respondent_index] + slope = true_slopes[respondent_index] + for occasion_index in range(5): + occasion_id = f"{respondent}{occasion_index}" + occasions.append( + _occasion( + respondent, + occasion_id, + occasion_index, + occasion_index * 86_400_000, + ) + ) + values[f"occasion_{occasion_id}"] = intercept + slope * occasion_index + design = build_longitudinal_design( + occasions=occasions, + state_spec=build_longitudinal_state_spec( + state_kind=LongitudinalStateKind.RANDOM_INTERCEPT_SLOPE, + ), + ) + result = fit_longitudinal_state(design, values, worker_count=3) + intercept_rmse = float( + np.sqrt(np.mean((np.asarray(result["intercepts"]) - true_intercepts) ** 2)) + ) + slope_rmse = float(np.sqrt(np.mean((np.asarray(result["slopes"]) - true_slopes) ** 2))) + assert intercept_rmse < 1e-12 + assert slope_rmse < 1e-12 + assert result["rmse"] < 1e-12 + assert result["observed_count"] == 15 + + +def test_all_missing_respondent_has_zero_rmse() -> None: + """An all-missing design uses the intercept-only empty-observation branch.""" + design = build_longitudinal_design( + occasions=[ + _occasion("respondent_a", "a0", 0, 0), + _occasion("respondent_a", "a1", 1, 86_400_000), + ], + state_spec=build_longitudinal_state_spec( + state_kind=LongitudinalStateKind.RANDOM_INTERCEPT_SLOPE, + ), + ) + empty = fit_longitudinal_state(design, {}, worker_count=3) + assert empty["observed_count"] == 0 + assert empty["rmse"] == 0.0 + np.testing.assert_allclose(empty["intercepts"], [0.0]) + np.testing.assert_allclose(empty["slopes"], [0.0]) + + +def test_rust_state_fit_preserves_discrete_ar_and_irregular_time() -> None: + """An AR fixture verifies the explicit discrete-step contract.""" + design = build_longitudinal_design( + occasions=[ + _occasion("respondent_a", "a0", 0, 0), + _occasion("respondent_a", "a1", 1, 86_400_000), + _occasion("respondent_a", "a2", 4, 259_200_000), + ], + state_spec=build_longitudinal_state_spec( + state_kind=LongitudinalStateKind.STATIONARY_AUTOREGRESSIVE, + autoregressive_coefficient=0.5, + ), + ) + result = fit_longitudinal_state( + design, + {"occasion_a0": 1.0, "occasion_a1": 0.5, "occasion_a2": 0.125}, + ) + assert result["ar_coefficient"] == 0.5 + assert result["transition_count"] == 2 + assert result["estimand_scope"] == "discrete_ar_state_prediction" + assert result["population_random_effects_estimated"] is False + assert result["ar_coefficient_estimated"] is False + assert result["ar_coefficient_source"] == "caller_supplied" + # The third declared occasion is sequence three steps after the second; + # the discrete-step AR(1) therefore predicts 0.5**3 * 0.5 = 0.0625. + np.testing.assert_allclose(result["state"], [1.0, 0.5, 0.0625], atol=1e-12) + + +def test_rust_state_fit_recovers_true_ar_predictions_with_rmse() -> None: + """A caller-supplied AR series is recovered with near-zero prediction RMSE.""" + phi = 0.4 + start = 1.6 + values = {"occasion_a0": start} + occasions = [_occasion("respondent_a", "a0", 0, 0)] + current = start + for step in range(1, 6): + current *= phi + occasion_id = f"a{step}" + occasions.append(_occasion("respondent_a", occasion_id, step, step * 86_400_000)) + values[f"occasion_{occasion_id}"] = current + design = build_longitudinal_design( + occasions=occasions, + state_spec=build_longitudinal_state_spec( + state_kind=LongitudinalStateKind.STATIONARY_AUTOREGRESSIVE, + autoregressive_coefficient=phi, + ), + ) + result = fit_longitudinal_state(design, values) + assert result["rmse"] < 1e-12 + assert result["transition_count"] == 5 + np.testing.assert_allclose(result["state"][0], start, atol=1e-12) + + +def test_state_fit_rejects_invalid_worker_count_and_foreign_design() -> None: + """The public estimator validates execution controls before Rust dispatch.""" + design = _single_occasion_design() + with pytest.raises(ValueError, match="worker_count"): + fit_longitudinal_state(design, {}, worker_count=0) + + class ForeignDesign: + """Represent an object that was not produced by the package factory.""" + + with pytest.raises(ValueError, match="LongitudinalDesign"): + fit_longitudinal_state(ForeignDesign(), {}) # type: ignore[arg-type] + + +@pytest.mark.parametrize("value", [True, np.bool_(True), 1 + 2j, "not-a-number"]) +def test_state_fit_rejects_non_real_observation_values(value: object) -> None: + """Caller-controlled observations fail with the package-owned exception.""" + with pytest.raises(ValueError, match="must be a real number"): + fit_longitudinal_state( + _single_occasion_design(), + {"occasion_guard": value}, # type: ignore[dict-item] + ) + + +def test_state_fit_translates_hostile_mapping_reads_to_value_error() -> None: + """A hostile mapping cannot leak its implementation exception.""" + + class ExplodingMapping(Mapping[str, float]): + """Raise from every read path to model an adversarial mapping.""" + + def __getitem__(self, key: str) -> float: + raise RuntimeError(f"blocked read: {key}") + + def __iter__(self): + return iter(()) + + def __len__(self) -> int: + return 0 + + def get(self, key: str, default=None): + raise RuntimeError(f"blocked get: {key}") + + with pytest.raises(ValueError, match="plain read-only mapping"): + fit_longitudinal_state(_single_occasion_design(), ExplodingMapping()) + + +@pytest.mark.parametrize("value", [2, np.int64(2), 2.5, np.float64(2.5)]) +def test_state_fit_accepts_real_numeric_scalars(value: object) -> None: + """Supported Python and NumPy real scalars convert through the helper.""" + result = fit_longitudinal_state( + _single_occasion_design(), + {"occasion_guard": value}, # type: ignore[dict-item] + ) + np.testing.assert_allclose(result["intercepts"], [float(value)]) + assert result["observed_count"] == 1 diff --git a/tests/test_multilevel_core_loader.py b/tests/test_multilevel_core_loader.py index 5327463b0..1dc81c1cd 100644 --- a/tests/test_multilevel_core_loader.py +++ b/tests/test_multilevel_core_loader.py @@ -29,6 +29,9 @@ 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") + assert hasattr(first, "fit_longitudinal_state") second = loader.multilevel_core() assert second is first