diff --git a/CHANGELOG.md b/CHANGELOG.md index bc9855d5a..acc6f377a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -119,6 +119,9 @@ ### Added +- `fleiss_kappa` Fleiss' multi-rater kappa for nominal agreement with the exact (Conger) chance-agreement variant (CRAN irr 0.85 `kappam.fleiss()` `R/kappam.fleiss.R` — READ and normative; Fleiss 1971 and Conger 1980 NOT READ, cited as model origins only): classification-table agreement `agreeP = (1/m)Σᵢ(Σⱼtᵢⱼ² − nr)/(nr(nr−1))`, classic chance `Σⱼpⱼ²` vs exact `Σⱼpⱼ² − (1/nr)Σⱼs²ⱼ` (sample variance over per-rater category proportions; algebra verified against R's `apply(rtab,2,var)` form), Fleiss' large-sample z test and category-wise kappas (classic mode; NaN for empty categories, matching R's 0/0), and listwise row drop for missing ratings. API deviations documented: index codes 0..k-1 with explicit/inferred k, negative-or-NaN = missing, error on degenerate `1 − chanceP = 0` (R returns NaN). Exact-Fraction oracle anchors FK1–FK5 (classic κ = 139/399, exact κ = 37/102, category κ = [1/21, 31/91, 43/63]), a 6-mutant EXECUTED kill map (row-vs-column chance sums, missing-as-category, variance-sign, pjk-centering), and an MC-500 subject/rater permutation-invariance test. +- `bratt_mm` Bradley-Terry model with ties fitted by MM (VGAM 1.1-14 `bratt()` family, `R/family.categorical.R` — READ and normative; Bradley & Terry 1952 NOT READ, cited as model origin): `P(i>j) = αᵢ/(αᵢ+αⱼ+α₀)`, `P(tie) = α₀/(αᵢ+αⱼ+α₀)` with a hand-derived supporting-hyperplane MM ascent (same pattern as the crate's `bradley_terry_mm`) and a joint reference rescale of α AND α₀ (likelihood-preserving; verified identity `Σ wins + T = Σ n_ij`). This is the additive-α₀ ties model, NOT Rao-Kupper/Davidson (neither read; disambiguation only). Contract: fractional weighted counts accepted, symmetric ties matrix required, tie-free data rejected (use `bradley_terry_mm`; an API contract, not VGAM behavior), zero-win contestants rejected, n capped at 10000 (O(n²) guard). Exact-Fraction oracle anchors B1–B4 (iteration-1 pins `[1, 27/40, 3/4]`, α₀ = 9/14), a 5-mutant EXECUTED kill map (incl. a tol-separated convergence anchor killing an α₀-blind convergence check), and an MC-500 log-likelihood dominance test. +- `predict_rating` / `predict_rating_multi` game-outcome prediction from fitted ratings (CRAN PlayerRatings 1.1-0 `predict.rating` `R/ratings.R` lines 1056–1133 — READ and normative; no journal paper exists for this dispatch, CRAN package provenance only): Elo logistic branch, deviation-shrunk Glicko/Glicko-2/Stephenson branch (`qip3 = 3(ln10/400π)²`, joint shrink over BOTH players' squared deviations), and multi-player EloM branch (`(rating − rowmean)/40`, na.rm row means, optional min-tie placing ranks with NaN kept). R semantics preserved: strict `games < tng` unrated cutoff, `trat` replacement of ALL missing extracted values (unmatched, low-games, stored-NA), `pred >= thresh` binarization with NaN propagation. REDUCED-SCOPE vs R: index-based (−1 = unmatched; caller does name matching), per-game/scalar gamma only. Exact-oracle fixtures P1–P9 and a 7-mutant EXECUTED kill map (incl. both branches' tng comparisons). - `fide_rating` FIDE-style Elo ratings (CRAN PlayerRatings 1.1-0 `fide()` `R/ratings.R` lines 125–272 + `kfide()` lines 959–972 — READ and normative; no journal paper exists for this variant, CRAN package provenance only): per-period batch Elo with the kfide K-factor schedule (K = kv[0] elite / kv[1] ≥30 games / kv[2] novice, evaluated from PERIOD-START state), sticky elite flag set from POST-update ratings ≥ 2400, and per-player running mean of POST-update opponent ratings. REDUCED-SCOPE vs R: no status/history frames, kfide-only K schedule, self-play rejected, thresholds 30/2400 hard-coded. kv=(k,k,k) reduces bitwise to `elo_rating(kfac=k)` (MC-500 anchor); exact-oracle fixtures F1–F5 and a 5-mutant EXECUTED kill map. - `metrics_rating` prediction-quality metrics for binary-outcome forecasts (binomial deviance on capped predictions, RMSE/MAE on raw predictions, optional 0.5-baseline scaling), a Rust reimplementation of CRAN PlayerRatings 1.1-0 `metrics()` with its cap quirk and elementwise NaN semantics preserved; exact-oracle anchor tests, mutation-kill map, and MC-500 invariants. - **Multiplayer Elo rating (CRAN PlayerRatings 1.1-0 `elom()` diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 7d2ddd264..efd3f0c6e 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -5741,6 +5741,142 @@ fn fide_rating( Ok(d.into()) } +/// Predicted game outcomes from fitted ratings, PlayerRatings +/// `predict.rating` two-player branches (see +/// `mlsirm_core::scaling::predict_rating_two`). `white`/`black` are player +/// indices with -1 = unmatched. Pass `deviations` for the +/// Glicko/Glicko-2/Stephenson deviation-shrunk branch. +#[pyfunction] +#[pyo3(signature = (ratings, deviations, games, white, black, gamma, tng, trat_rating, trat_deviation, thresh))] +#[allow(clippy::too_many_arguments)] +fn predict_rating_two( + py: Python<'_>, + ratings: PyReadonlyArray1<'_, f64>, + deviations: Option>, + games: PyReadonlyArray1<'_, u64>, + white: PyReadonlyArray1<'_, i64>, + black: PyReadonlyArray1<'_, i64>, + gamma: PyReadonlyArray1<'_, f64>, + tng: u64, + trat_rating: Option, + trat_deviation: Option, + thresh: Option, +) -> PyResult>> { + let dev_slice = deviations.as_ref().map(|d| d.as_slice()).transpose()?; + let trat = trat_rating.map(|t1| (t1, trat_deviation.unwrap_or(f64::NAN))); + let out = mlsirm_core::scaling::predict_rating_two( + ratings.as_slice()?, + dev_slice, + games.as_slice()?, + white.as_slice()?, + black.as_slice()?, + gamma.as_slice()?, + tng, + trat, + thresh, + ) + .map_err(PyValueError::new_err)?; + Ok(PyArray1::from_slice(py, &out).into()) +} + +/// Predicted expected scores for multi-player (EloM) events, PlayerRatings +/// `predict.rating` EloM branch (see +/// `mlsirm_core::scaling::predict_rating_multi`). `players` is flattened +/// row-major nr x np with -1 = empty seat; `placing` returns min-tie ranks. +#[pyfunction] +#[pyo3(signature = (ratings, games, players, nr, np, tng, trat, placing))] +#[allow(clippy::too_many_arguments)] +fn predict_rating_multi( + py: Python<'_>, + ratings: PyReadonlyArray1<'_, f64>, + games: PyReadonlyArray1<'_, u64>, + players: PyReadonlyArray1<'_, i64>, + nr: usize, + np: usize, + tng: u64, + trat: Option, + placing: bool, +) -> PyResult>> { + let out = mlsirm_core::scaling::predict_rating_multi( + ratings.as_slice()?, + games.as_slice()?, + players.as_slice()?, + nr, + np, + tng, + trat, + placing, + ) + .map_err(PyValueError::new_err)?; + Ok(PyArray1::from_slice(py, &out).into()) +} + + +/// Bradley-Terry model with ties (additive alpha0, VGAM `bratt`) fitted +/// by MM (see `mlsirm_core::scaling::bratt_mm`). `wins` and `ties` are +/// flat row-major n*n matrices (ties symmetric); returns dict with alpha +/// (worths, alpha[ref_index] == ref_value), alpha0 (tie parameter), +/// iterations, log_likelihood. +#[pyfunction] +#[pyo3(signature = (wins, ties, n, ref_index=0, ref_value=1.0, max_iter=10000, tol=1e-10))] +#[allow(clippy::too_many_arguments)] +fn bratt_mm( + py: Python<'_>, + wins: PyReadonlyArray1<'_, f64>, + ties: PyReadonlyArray1<'_, f64>, + n: usize, + ref_index: usize, + ref_value: f64, + max_iter: usize, + tol: f64, +) -> PyResult> { + let res = mlsirm_core::scaling::bratt_mm( + wins.as_slice()?, + ties.as_slice()?, + n, + ref_index, + ref_value, + max_iter, + tol, + ) + .map_err(PyValueError::new_err)?; + let d = pyo3::types::PyDict::new(py); + d.set_item("alpha", PyArray1::from_slice(py, &res.alpha))?; + d.set_item("alpha0", res.alpha0)?; + d.set_item("iterations", res.iterations as u64)?; + d.set_item("log_likelihood", res.log_likelihood)?; + Ok(d.into()) +} + +/// Fleiss' kappa for nominal agreement among nr raters over ns subjects, +/// with the exact (Conger) variant (irr 0.85 `kappam.fleiss`; see +/// `mlsirm_core::agreement::fleiss_kappa`). `ratings` is flat row-major +/// ns*nr of codes 0..k-1; negative = missing (listwise row drop). Returns +/// dict with kappa, subjects_used, z, p_value, category_kappa/z/p +/// (empty arrays and NaN z/p in exact mode). +#[pyfunction] +#[pyo3(signature = (ratings, ns, nr, k, exact=false))] +fn fleiss_kappa( + py: Python<'_>, + ratings: PyReadonlyArray1<'_, i64>, + ns: usize, + nr: usize, + k: usize, + exact: bool, +) -> PyResult> { + let res = mlsirm_core::agreement::fleiss_kappa(ratings.as_slice()?, ns, nr, k, exact) + .map_err(PyValueError::new_err)?; + let d = pyo3::types::PyDict::new(py); + d.set_item("kappa", res.kappa)?; + d.set_item("subjects_used", res.subjects_used as u64)?; + d.set_item("z", res.z)?; + d.set_item("p_value", res.p_value)?; + d.set_item("category_kappa", PyArray1::from_slice(py, &res.category_kappa))?; + d.set_item("category_z", PyArray1::from_slice(py, &res.category_z))?; + d.set_item("category_p", PyArray1::from_slice(py, &res.category_p))?; + Ok(d.into()) +} + /// GPCM/nominal softmax cell log-probabilities at one node (parity surface for /// the NumPy `category_logprobs` reference). #[pyfunction] @@ -8095,6 +8231,10 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(elom_rating, m)?)?; m.add_function(wrap_pyfunction!(metrics_rating, m)?)?; m.add_function(wrap_pyfunction!(fide_rating, m)?)?; + m.add_function(wrap_pyfunction!(predict_rating_two, m)?)?; + m.add_function(wrap_pyfunction!(predict_rating_multi, m)?)?; + m.add_function(wrap_pyfunction!(bratt_mm, m)?)?; + m.add_function(wrap_pyfunction!(fleiss_kappa, m)?)?; m.add_function(wrap_pyfunction!(circle_arc_middle_anchor, m)?)?; m.add_function(wrap_pyfunction!(loglinear_smooth, m)?)?; m.add_function(wrap_pyfunction!(person_fit_stat, m)?)?; diff --git a/crates/mlsirm-core/src/agreement.rs b/crates/mlsirm-core/src/agreement.rs index 3b88a1de0..c075bc58f 100644 --- a/crates/mlsirm-core/src/agreement.rs +++ b/crates/mlsirm-core/src/agreement.rs @@ -241,6 +241,198 @@ pub fn validate_scoring( }) } +/// Result of Fleiss' multi-rater kappa (`fleiss_kappa`). +/// +/// In exact (Conger) mode `z`/`p_value` are NaN and the category vectors +/// are empty, mirroring irr's `kappam.fleiss(exact=TRUE)` which returns +/// neither a test statistic nor category detail. +#[derive(Clone, Debug)] +pub struct FleissKappaResult { + pub kappa: f64, + /// Subjects remaining after listwise deletion of rows with missing codes. + pub subjects_used: usize, + pub z: f64, + pub p_value: f64, + /// Category-wise kappas (classic mode only); NaN for empty categories. + pub category_kappa: Vec, + pub category_z: Vec, + pub category_p: Vec, +} + +/// Fleiss' kappa for nominal agreement among `nr` raters over `ns` subjects, +/// with the exact (Conger) chance-agreement variant. +/// +/// Reimplements `kappam.fleiss()` from CRAN irr 0.85 (`R/kappam.fleiss.R`, +/// READ in full; algorithm source of truth). The model originates in +/// Fleiss, J. L. (1971), "Measuring nominal scale agreement among many +/// raters," Psychological Bulletin, 76(5), 378-382, and the exact variant in +/// Conger, A. J. (1980), "Integration and generalization of kappas for +/// multiple raters," Psychological Bulletin, 88(2), 322-328 — both cited as +/// origins only (NOT READ); every formula below was verified against the irr +/// R source. +/// +/// `ratings` is row-major `ns x nr` with category codes `0..k-1`; a negative +/// code marks a missing rating and drops the whole subject row (listwise, as +/// in R: `ratings[apply(is.na(ratings),1,sum)==0,]`). With `m` used subjects +/// and `ttab[i][j]` = raters assigning subject `i` to category `j`: +/// +/// - `agreeP = (1/m) sum_i (sum_j ttab_ij^2 - nr) / (nr(nr-1))` +/// - classic `chanceP = sum_j p_j^2` with `p_j = C_j/(m nr)`, `C_j` column sums +/// - exact `chanceP = sum_j p_j^2 - (1/nr) sum_j s2_j`, `s2_j` the sample +/// variance (divisor `nr-1`) over raters of per-rater category proportions +/// (algebraically equal to R's `sum(apply(rtab,2,var)*(nr-1)/nr)/(nr-1)`) +/// - `kappa = (agreeP - chanceP)/(1 - chanceP)` +/// +/// Classic mode adds Fleiss' large-sample test +/// `var = 2[(sum p_j q_j)^2 - sum p_j q_j (q_j - p_j)] / +/// [(sum p_j q_j)^2 m nr (nr-1)]`, `z = kappa/sqrt(var)`, +/// `p = 2(1 - Phi(|z|))`, and category-wise kappas +/// `pjk_j = (sum_i ttab_ij^2 - m nr p_j)/(m nr (nr-1) p_j)`, +/// `kappa_j = (pjk_j - p_j)/(1 - p_j)`, `var_j = 2/(m nr (nr-1))` +/// (computed unconditionally; identical to irr's `detail=TRUE`). Empty +/// categories yield NaN, matching R's 0/0. +/// +/// API deviations from R (documented contract, not transcription): codes are +/// index-based `0..k-1` with explicit `k` (R derives factor levels; negative +/// numeric labels must be remapped by the caller since negative = missing +/// here), degenerate `1 - chanceP == 0` is an error (R returns NaN), and the +/// size caps below are safety bounds. +pub fn fleiss_kappa( + ratings: &[i64], + ns: usize, + nr: usize, + k: usize, + exact: bool, +) -> Result { + if ns == 0 { + return Err("need at least one subject".into()); + } + if nr < 2 { + return Err("need at least 2 raters".into()); + } + if k < 2 { + return Err("need at least 2 categories".into()); + } + if ns > 1_000_000 || nr > 10_000 || k > 10_000 { + return Err("size caps: ns <= 1e6, nr <= 1e4, k <= 1e4".into()); + } + if ratings.len() != ns * nr { + return Err(format!( + "ratings length {} != ns*nr = {}", + ratings.len(), + ns * nr + )); + } + for &c in ratings { + if c >= k as i64 { + return Err(format!("category code {c} out of range 0..{k}")); + } + } + // Listwise drop of rows containing any negative (missing) code, then + // classification table ttab[i][j] and per-rater counts. Counts are exact + // in f64: entries <= nr <= 1e4, sums of squares <= m*nr^2 <= 1e14 < 2^53. + let mut ttab: Vec> = Vec::new(); + let mut rater_counts = vec![vec![0.0_f64; k]; nr]; + for i in 0..ns { + let row = &ratings[i * nr..(i + 1) * nr]; + if row.iter().any(|&c| c < 0) { + continue; + } + let mut t = vec![0.0_f64; k]; + for (r, &c) in row.iter().enumerate() { + t[c as usize] += 1.0; + rater_counts[r][c as usize] += 1.0; + } + ttab.push(t); + } + let m = ttab.len(); + if m == 0 { + return Err("all subject rows dropped for missing ratings".into()); + } + let mf = m as f64; + let nrf = nr as f64; + + let agree_p: f64 = ttab + .iter() + .map(|t| (t.iter().map(|&v| v * v).sum::() - nrf) / (nrf * (nrf - 1.0))) + .sum::() + / mf; + + let col: Vec = (0..k) + .map(|j| ttab.iter().map(|t| t[j]).sum::()) + .collect(); + let p: Vec = col.iter().map(|&c| c / (mf * nrf)).collect(); + let mut chance_p: f64 = p.iter().map(|&v| v * v).sum(); + if exact { + // Sample variance over raters of rtab[r][j] = rater_counts[r][j]/m. + let mut s2_sum = 0.0_f64; + for j in 0..k { + let props: Vec = (0..nr).map(|r| rater_counts[r][j] / mf).collect(); + let mean = props.iter().sum::() / nrf; + let s2 = props.iter().map(|&v| (v - mean) * (v - mean)).sum::() / (nrf - 1.0); + s2_sum += s2; + } + chance_p -= s2_sum / nrf; + } + let denom = 1.0 - chance_p; + if denom.abs() < 1e-12 { + return Err("degenerate marginals: no chance-corrected agreement is defined".into()); + } + let kappa = (agree_p - chance_p) / denom; + + if exact { + return Ok(FleissKappaResult { + kappa, + subjects_used: m, + z: f64::NAN, + p_value: f64::NAN, + category_kappa: Vec::new(), + category_z: Vec::new(), + category_p: Vec::new(), + }); + } + + let sqrt2 = std::f64::consts::SQRT_2; + let pq: f64 = p.iter().map(|&v| v * (1.0 - v)).sum(); + let var = 2.0 + * (pq * pq + - p.iter() + .map(|&v| v * (1.0 - v) * (1.0 - 2.0 * v)) + .sum::()) + / (pq * pq * mf * nrf * (nrf - 1.0)); + let z = kappa / var.sqrt(); + let p_value = crate::fitstats::erfc(z.abs() / sqrt2); + + let var_k = 2.0 / (mf * nrf * (nrf - 1.0)); + let mut category_kappa = Vec::with_capacity(k); + let mut category_z = Vec::with_capacity(k); + let mut category_p = Vec::with_capacity(k); + for j in 0..k { + let sum_sq: f64 = ttab.iter().map(|t| t[j] * t[j]).sum(); + // Empty category: p_j = 0 gives R's 0/0 = NaN, preserved here. + let pjk = (sum_sq - mf * nrf * p[j]) / (mf * nrf * (nrf - 1.0) * p[j]); + let kj = (pjk - p[j]) / (1.0 - p[j]); + let zj = kj / var_k.sqrt(); + category_kappa.push(kj); + category_z.push(zj); + category_p.push(if zj.is_finite() { + crate::fitstats::erfc(zj.abs() / sqrt2) + } else { + f64::NAN + }); + } + + Ok(FleissKappaResult { + kappa, + subjects_used: m, + z, + p_value, + category_kappa, + category_z, + category_p, + }) +} + #[cfg(test)] #[path = "../../../tests/unit/agreement_tests.rs"] mod tests; diff --git a/crates/mlsirm-core/src/scaling.rs b/crates/mlsirm-core/src/scaling.rs index 8a88c80e4..b82104b6a 100644 --- a/crates/mlsirm-core/src/scaling.rs +++ b/crates/mlsirm-core/src/scaling.rs @@ -3701,6 +3701,552 @@ pub fn fide_rating( }) } +/// Game-outcome prediction from fitted ratings — two-player branches of +/// `predict.rating` from CRAN PlayerRatings 1.1-0 (`R/ratings.R` lines +/// 1056-1133; source READ). No journal paper exists for this dispatch +/// function; the CRAN R source is the normative reference (provenance as +/// for `stephenson_rating`). Two branches are implemented here: +/// +/// * Elo branch (fitted `type == "Elo"`, produced by `elo()`/`fide()`): +/// `pred = 1 / (1 + 10^((brat - wrat - gamma)/400))` (R line 1118). +/// * Deviation branch (`Glicko`/`Glicko-2`/`Stephenson`), engaged when +/// `deviations` is supplied: with `qv = ln(10)/400` and +/// `qip3 = 3 (qv/pi)^2`, +/// `vec = 1/sqrt(1 + qip3 (wdev^2 + bdev^2))` and +/// `pred = 1 / (1 + 10^(vec (brat - wrat - gamma)/400))` +/// (R lines 1071, 1119-1121; note the joint shrink sums BOTH players' +/// squared deviations). +/// +/// Pre-processing (verified against R 1069-1114 and the executed oracle): +/// players with `games < tng` (strict `<`) have their stored rating (and +/// deviation) set to NA before extraction; when `trat` is supplied it +/// then replaces ALL extracted-NA values — unmatched players (index +/// sentinel `-1`), low-games players, and matched players whose stored +/// rating/deviation is already NaN. Without `trat`, NA propagates to a +/// NaN prediction. `thresh` (R 1127-1128) maps `pred >= thresh` to 1 else +/// 0, leaving NaN predictions NaN. +/// +/// REDUCED SCOPE relative to R: index-based (caller does name matching; +/// `-1` = unmatched), per-game `gamma` (scalar recycling is a wrapper +/// concern), no data-frame interface or `object$type` dispatch. +/// +/// All formulas were verified against an executed exact/float oracle +/// transcribed from the R source; claims are limited to that source. +pub fn predict_rating_two( + ratings: &[f64], + deviations: Option<&[f64]>, + games: &[u64], + white: &[i64], + black: &[i64], + gamma: &[f64], + tng: u64, + trat: Option<(f64, f64)>, + thresh: Option, +) -> Result, String> { + let n = ratings.len(); + if n < 2 || n > 10_000 { + return Err(format!( + "predict_rating_two: number of players must be in 2..=10000, got {n}" + )); + } + if games.len() != n { + return Err(format!( + "predict_rating_two: games length {} != number of players {n}", + games.len() + )); + } + if let Some(dev) = deviations { + if dev.len() != n { + return Err(format!( + "predict_rating_two: deviations length {} != number of players {n}", + dev.len() + )); + } + if dev.iter().any(|d| d.is_infinite()) { + return Err( + "predict_rating_two: deviations must not be infinite (NaN = missing)".to_string(), + ); + } + } + if ratings.iter().any(|r| r.is_infinite()) { + return Err("predict_rating_two: ratings must not be infinite (NaN = missing)".to_string()); + } + let g = white.len(); + if g == 0 { + return Err("predict_rating_two: at least one game row is required".to_string()); + } + if black.len() != g || gamma.len() != g { + return Err(format!( + "predict_rating_two: white/black/gamma lengths must match, got {}/{}/{}", + g, + black.len(), + gamma.len() + )); + } + if gamma.iter().any(|x| !x.is_finite()) { + return Err("predict_rating_two: gamma must be finite".to_string()); + } + if let Some((t1, t2)) = trat { + if !t1.is_finite() { + return Err("predict_rating_two: trat rating must be finite".to_string()); + } + if deviations.is_some() && !t2.is_finite() { + return Err("predict_rating_two: trat deviation must be finite".to_string()); + } + } + if let Some(t) = thresh { + if !t.is_finite() { + return Err("predict_rating_two: thresh must be finite".to_string()); + } + } + for (w, b) in white.iter().zip(black.iter()) { + for idx in [*w, *b] { + if idx < -1 || idx >= n as i64 { + return Err(format!( + "predict_rating_two: player index {idx} out of range (-1 = unmatched, else 0..{n})" + )); + } + } + if *w >= 0 && w == b { + return Err("predict_rating_two: self-play rows are not allowed".to_string()); + } + } + // R 1088-1093: low-games rows lose their stored values before extraction. + let effective = |p: i64, table: &[f64]| -> f64 { + if p < 0 { + f64::NAN + } else if games[p as usize] < tng { + f64::NAN + } else { + table[p as usize] + } + }; + let qv = std::f64::consts::LN_10 / 400.0; + let qip3 = 3.0 * (qv / std::f64::consts::PI) * (qv / std::f64::consts::PI); + let mut preds = Vec::with_capacity(g); + for k in 0..g { + let mut wrat = effective(white[k], ratings); + let mut brat = effective(black[k], ratings); + if let Some((t1, _)) = trat { + // R 1097-1098: trat replaces ALL extracted NA ratings. + if wrat.is_nan() { + wrat = t1; + } + if brat.is_nan() { + brat = t1; + } + } + let pred = if let Some(dev) = deviations { + let mut wdev = effective(white[k], dev); + let mut bdev = effective(black[k], dev); + if let Some((_, t2)) = trat { + if wdev.is_nan() { + wdev = t2; + } + if bdev.is_nan() { + bdev = t2; + } + } + if wrat.is_nan() || brat.is_nan() || wdev.is_nan() || bdev.is_nan() { + f64::NAN + } else { + let vec = 1.0 / (1.0 + qip3 * (wdev * wdev + bdev * bdev)).sqrt(); + 1.0 / (1.0 + 10f64.powf(vec * (brat - wrat - gamma[k]) / 400.0)) + } + } else if wrat.is_nan() || brat.is_nan() { + f64::NAN + } else { + 1.0 / (1.0 + 10f64.powf((brat - wrat - gamma[k]) / 400.0)) + }; + preds.push(pred); + } + if let Some(t) = thresh { + // R 1127-1128: as.numeric(preds >= thresh); NA stays NA. + for p in preds.iter_mut() { + if !p.is_nan() { + *p = if *p >= t { 1.0 } else { 0.0 }; + } + } + } + Ok(preds) +} + +/// Multi-player (EloM) branch of `predict.rating` from CRAN PlayerRatings +/// 1.1-0 (`R/ratings.R` lines 1103-1105, 1123-1125, 1129-1130; source +/// READ; provenance as for `predict_rating_two`). +/// +/// Per event row: `pred = (rat - rowmean)/40` where `rowmean` is the mean +/// over the row's non-NaN seat ratings (`rowMeans(rats, na.rm=TRUE)`, +/// all-NaN row -> NaN row). Pre-processing (tng/trat) as in +/// `predict_rating_two`. With `placing = true` the per-row predictions +/// are replaced by `rank(-preds, na.last="keep", ties.method="min")`: +/// rank 1 = highest prediction, ties share the minimum rank, NaN seats +/// keep NaN. +/// +/// `players` is row-major `nr x np` with `-1` = empty/unmatched seat. +/// All formulas were verified against an executed oracle transcribed +/// from the R source. +pub fn predict_rating_multi( + ratings: &[f64], + games: &[u64], + players: &[i64], + nr: usize, + np: usize, + tng: u64, + trat: Option, + placing: bool, +) -> Result, String> { + let n = ratings.len(); + if n < 2 || n > 10_000 { + return Err(format!( + "predict_rating_multi: number of players must be in 2..=10000, got {n}" + )); + } + if games.len() != n { + return Err(format!( + "predict_rating_multi: games length {} != number of players {n}", + games.len() + )); + } + if ratings.iter().any(|r| r.is_infinite()) { + return Err( + "predict_rating_multi: ratings must not be infinite (NaN = missing)".to_string(), + ); + } + if nr == 0 { + return Err("predict_rating_multi: at least one event row is required".to_string()); + } + if np < 2 || np > 1000 { + return Err(format!( + "predict_rating_multi: seats per event must be in 2..=1000, got {np}" + )); + } + let total = nr.checked_mul(np).ok_or_else(|| { + format!("predict_rating_multi: nr * np overflows usize (nr = {nr}, np = {np})") + })?; + if players.len() != total { + return Err(format!( + "predict_rating_multi: players length {} != nr*np = {}", + players.len(), + total + )); + } + if let Some(t) = trat { + if !t.is_finite() { + return Err("predict_rating_multi: trat must be finite".to_string()); + } + } + for idx in players { + if *idx < -1 || *idx >= n as i64 { + return Err(format!( + "predict_rating_multi: player index {idx} out of range (-1 = empty, else 0..{n})" + )); + } + } + let mut out = vec![f64::NAN; total]; + let mut row_rat = vec![f64::NAN; np]; + for i in 0..nr { + for s in 0..np { + let p = players[i * np + s]; + let mut v = if p < 0 || games[p as usize] < tng { + f64::NAN + } else { + ratings[p as usize] + }; + if let Some(t) = trat { + if v.is_nan() { + v = t; + } + } + row_rat[s] = v; + } + // rowMeans(rats, na.rm = TRUE) — R 1123. + let mut sum = 0.0; + let mut cnt = 0usize; + for v in row_rat.iter() { + if !v.is_nan() { + sum += v; + cnt += 1; + } + } + let mean = if cnt > 0 { sum / cnt as f64 } else { f64::NAN }; + for s in 0..np { + if !row_rat[s].is_nan() { + out[i * np + s] = (row_rat[s] - mean) / 40.0; + } + } + if placing { + // rank(-preds, na.last="keep", ties.method="min") — R 1129-1130: + // rank of v = 1 + count of non-NaN row values strictly greater. + for s in 0..np { + row_rat[s] = out[i * np + s]; + } + for s in 0..np { + if !row_rat[s].is_nan() { + let greater = row_rat + .iter() + .filter(|u| !u.is_nan() && **u > row_rat[s]) + .count(); + out[i * np + s] = 1.0 + greater as f64; + } + } + } + } + Ok(out) +} + +/// Result of a Bradley-Terry-with-ties (VGAM `bratt`) MM fit. +#[derive(Clone, Debug)] +pub struct BrattResult { + /// Worth parameters `alpha_i > 0` with `alpha[ref_index] == ref_value`. + pub alpha: Vec, + /// Tie parameter `alpha0 > 0` (same joint scale as `alpha`). + pub alpha0: f64, + /// Number of MM updates performed when convergence fired. + pub iterations: usize, + /// Log-likelihood at the returned parameters. + pub log_likelihood: f64, +} + +/// Bradley-Terry model with ties (additive `alpha0` tie parameter, the +/// VGAM `bratt` family) fitted by maximum likelihood via an MM algorithm. +/// +/// Model (VGAM `bratt@linkinv`): for contestants `i != j` with worths +/// `alpha_i > 0` and tie parameter `alpha0 > 0`, writing +/// `D_ij = alpha_i + alpha_j + alpha0`, +/// +/// ```text +/// P(i beats j) = alpha_i / D_ij +/// P(i ties j) = alpha0 / D_ij +/// ``` +/// +/// Log-likelihood (VGAM `bratt@loglikelihood`, +/// `y*log(mu) + 0.5*ties*log(probtie)` with the symmetric ties matrix +/// collapsing the 0.5 double count to one term per unordered pair): +/// +/// ```text +/// LL = sum_{i != j} y_ij ln alpha_i + T ln alpha0 +/// - sum_{i= -ln D_k - (D - D_k)/D_k`), +/// whose separable surrogate is maximized by +/// +/// ```text +/// alpha_i' = W_i / sum_{j != i} n_ij / D_ij_k (W_i = sum_{j != i} y_ij) +/// alpha0' = T / sum_{i 0`, outside the +/// parameter space (and the reference rescale would divide by zero). +/// - Start `alpha_i = 1, alpha0 = 1`; convergence fires when the max +/// absolute change across `alpha` and `alpha0` (after rescale) is +/// `<= tol`; non-convergence within `max_iter` updates is an error. +pub fn bratt_mm( + wins: &[f64], + ties: &[f64], + n: usize, + ref_index: usize, + ref_value: f64, + max_iter: usize, + tol: f64, +) -> Result { + if n < 2 { + return Err("bratt_mm needs at least 2 contestants".into()); + } + if n > 10000 { + return Err(format!("n = {n} exceeds the bratt_mm cap of 10000")); + } + let nn = n * n; + if wins.len() != nn { + return Err(format!( + "wins must be a row-major {n}x{n} matrix ({nn} entries), got {}", + wins.len() + )); + } + if ties.len() != nn { + return Err(format!( + "ties must be a row-major {n}x{n} matrix ({nn} entries), got {}", + ties.len() + )); + } + for (name, m) in [("win", wins), ("tie", ties)] { + for (k, &c) in m.iter().enumerate() { + if !c.is_finite() { + return Err(format!("{name} counts must be finite")); + } + if c < 0.0 { + return Err(format!("{name} counts must be nonnegative")); + } + if k / n == k % n && c != 0.0 { + return Err(format!("diagonal of the {name}s matrix must be zero")); + } + } + } + for i in 0..n { + for j in (i + 1)..n { + if ties[i * n + j] != ties[j * n + i] { + return Err(format!( + "ties matrix must be symmetric (ties[{i}][{j}] != ties[{j}][{i}])" + )); + } + } + } + if ref_index >= n { + return Err(format!("ref_index = {ref_index} out of range for n = {n}")); + } + if !ref_value.is_finite() || ref_value <= 0.0 { + return Err("ref_value must be finite and positive".into()); + } + if !tol.is_finite() || tol <= 0.0 { + return Err("tol must be finite and positive".into()); + } + if max_iter == 0 { + return Err("max_iter must be at least 1".into()); + } + let w_tot: Vec = (0..n) + .map(|i| (0..n).filter(|&j| j != i).map(|j| wins[i * n + j]).sum()) + .collect(); + if let Some(i) = w_tot.iter().position(|&w| w == 0.0) { + return Err(format!( + "contestant {i} has no wins: its ML worth is 0, outside the \ + positive parameter space of bratt_mm" + )); + } + let t_tot: f64 = (0..n) + .flat_map(|i| ((i + 1)..n).map(move |j| (i, j))) + .map(|(i, j)| ties[i * n + j]) + .sum(); + if t_tot == 0.0 { + return Err("ties matrix has no ties: alpha0 has no positive MLE; \ + use bradley_terry_mm for tie-free data" + .into()); + } + // Impl-review guard: entries are finite, but derived aggregates + // (row win totals, tie total, pair totals n_ij) can still overflow + // to +inf, which would let the MM update return exactly-zero + // parameters and a NaN log-likelihood as Ok. + if w_tot.iter().any(|w| !w.is_finite()) || !t_tot.is_finite() { + return Err( + "aggregate win/tie counts overflow f64; counts are too large for bratt_mm".into(), + ); + } + for i in 0..n { + for j in (i + 1)..n { + if !(wins[i * n + j] + wins[j * n + i] + ties[i * n + j]).is_finite() { + return Err( + "aggregate win/tie counts overflow f64; counts are too large for bratt_mm" + .into(), + ); + } + } + } + + let mut alpha = vec![1.0f64; n]; + let mut alpha0 = 1.0f64; + for it in 1..=max_iter { + let mut denom = vec![0.0f64; n]; + let mut denom0 = 0.0f64; + for i in 0..n { + for j in (i + 1)..n { + let n_ij = wins[i * n + j] + wins[j * n + i] + ties[i * n + j]; + if n_ij > 0.0 { + let val = n_ij / (alpha[i] + alpha[j] + alpha0); + denom[i] += val; + denom[j] += val; + denom0 += val; + } + } + } + let mut anew: Vec = (0..n).map(|i| w_tot[i] / denom[i]).collect(); + let mut a0new = t_tot / denom0; + let c = ref_value / anew[ref_index]; + for a in anew.iter_mut() { + *a *= c; + } + a0new *= c; + if anew.iter().any(|a| !a.is_finite() || *a <= 0.0) || !a0new.is_finite() || a0new <= 0.0 { + return Err("MM update produced non-finite or non-positive parameters".into()); + } + let mut delta = (a0new - alpha0).abs(); + for i in 0..n { + delta = delta.max((anew[i] - alpha[i]).abs()); + } + alpha = anew; + alpha0 = a0new; + if delta <= tol { + let mut ll = 0.0f64; + for i in 0..n { + for j in 0..n { + if i != j && wins[i * n + j] > 0.0 { + ll += wins[i * n + j] * alpha[i].ln(); + } + } + } + for i in 0..n { + for j in (i + 1)..n { + let t_ij = ties[i * n + j]; + if t_ij > 0.0 { + ll += t_ij * alpha0.ln(); + } + let n_ij = wins[i * n + j] + wins[j * n + i] + t_ij; + if n_ij > 0.0 { + ll -= n_ij * (alpha[i] + alpha[j] + alpha0).ln(); + } + } + } + if !ll.is_finite() { + return Err("bratt_mm log-likelihood is not finite".into()); + } + return Ok(BrattResult { + alpha, + alpha0, + iterations: it, + log_likelihood: ll, + }); + } + } + Err(format!( + "bratt_mm did not converge within max_iter = {max_iter} updates" + )) +} + #[cfg(test)] #[path = "../../../tests/unit/scaling_tests.rs"] mod tests; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index 7bc8786df..2c7a80105 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -24,6 +24,7 @@ from .rt import fit_response_times as fit_response_times, RtFit as RtFit, fit_speed_accuracy as fit_speed_accuracy, rt_person_fit as rt_person_fit from .scaling import thurstone_case_v as thurstone_case_v, ThurstoneResult as ThurstoneResult from .scaling import bradley_terry_mm as bradley_terry_mm, BradleyTerryResult as BradleyTerryResult +from .scaling import bratt_mm as bratt_mm, BrattResult as BrattResult from .scaling import lsr_pairwise as lsr_pairwise, ilsr_pairwise as ilsr_pairwise, LsrResult as LsrResult from .scaling import rank_centrality as rank_centrality from .scaling import lsr_rankings as lsr_rankings, ilsr_rankings as ilsr_rankings @@ -41,6 +42,8 @@ metrics_rating as metrics_rating, fide_rating as fide_rating, FideResult as FideResult, + predict_rating as predict_rating, + predict_rating_multi as predict_rating_multi, ) from .cdm import fit_cdm as fit_cdm, CdmFit as CdmFit, fit_gdina as fit_gdina, GdinaFit as GdinaFit, validate_q_matrix as validate_q_matrix, QMatrixValidation as QMatrixValidation, gdina_wald_selection as gdina_wald_selection, WaldModelSelection as WaldModelSelection, fit_ho_cdm as fit_ho_cdm, HoCdmFit as HoCdmFit, fit_ho_gdina as fit_ho_gdina, HoGdinaFit as HoGdinaFit, fit_seq_gdina as fit_seq_gdina, SeqGdinaFit as SeqGdinaFit, fit_seq_gdina_qr as fit_seq_gdina_qr, SeqGdinaQrFit as SeqGdinaQrFit from .mixture import fit_mixture as fit_mixture, MixtureFit as MixtureFit @@ -172,7 +175,9 @@ from .lltm import fit_lltm as fit_lltm, LltmFit as LltmFit from .testlet import fit_testlet as fit_testlet, TestletFit as TestletFit from .report import render_diagnostics_report as render_diagnostics_report -from .validation import (ValidationVerdict as ValidationVerdict, +from .validation import (FleissKappaResult as FleissKappaResult, + ValidationVerdict as ValidationVerdict, + fleiss_kappa as fleiss_kappa, validate_judge as validate_judge) from .serving import (bank_information as bank_information, cat_next_item as cat_next_item, @@ -217,6 +222,8 @@ "assemble_test_form", "dimensionality_diagnostics", "ValidationVerdict", + "FleissKappaResult", + "fleiss_kappa", "benjamini_hochberg", "chi2_sf", "dif_analysis", @@ -249,6 +256,8 @@ "ThurstoneResult", "bradley_terry_mm", "BradleyTerryResult", + "bratt_mm", + "BrattResult", "lsr_pairwise", "ilsr_pairwise", "rank_centrality", @@ -274,6 +283,8 @@ "metrics_rating", "fide_rating", "FideResult", + "predict_rating", + "predict_rating_multi", "circle_arc_middle_anchor", "CircleArcResult", "fit_response_times", diff --git a/python/fast_mlsirm/scaling.py b/python/fast_mlsirm/scaling.py index efaded88a..7d41a8e3b 100644 --- a/python/fast_mlsirm/scaling.py +++ b/python/fast_mlsirm/scaling.py @@ -182,6 +182,89 @@ def bradley_terry_mm(wins, alpha=0.0, max_iter=10000, tol=1e-8): ) +@dataclass +class BrattResult: + """Bradley-Terry-with-ties (VGAM ``bratt``) MM fit: ``alpha[i]`` is the + worth of contestant i (``alpha[ref_index] == ref_value``), ``alpha0`` the + additive tie parameter on the same joint scale, ``iterations`` the number + of MM updates performed when the max absolute parameter change (across + ``alpha`` AND ``alpha0``, after the reference rescale) first fell to + ``tol``, and ``log_likelihood`` the model log-likelihood at the returned + parameters.""" + + alpha: np.ndarray + alpha0: float + iterations: int + log_likelihood: float + + +def bratt_mm(wins, ties, ref_index=0, ref_value=1.0, max_iter=10000, tol=1e-10): + """Fit the Bradley-Terry model with ties (additive ``alpha0``) by MM. + + Model (VGAM 1.1-14 ``bratt()`` family, R source READ; Bradley & Terry + 1952 NOT READ, cited as the model origin): ``P(i beats j) = + alpha_i / (alpha_i + alpha_j + alpha0)`` and ``P(i ties j) = + alpha0 / (alpha_i + alpha_j + alpha0)``. This additive-``alpha0`` ties + model is NOT the Rao-Kupper or Davidson ties model (neither read; named + only to disambiguate). The MM ascent is hand-derived using the same + supporting-hyperplane pattern as :func:`bradley_terry_mm`. + + ``wins[i, j]`` is the (possibly fractional, nonnegative) count of wins + of *i* over *j*; ``ties[i, j]`` the tie count of the unordered pair, + stored symmetrically. Diagonals must be zero. Data with no ties at all + are rejected (``alpha0`` has no positive MLE -- use + :func:`bradley_terry_mm`; this is a contract of this API, not VGAM + behavior), as is any contestant with zero wins (its ML worth is 0, + outside the positive parameter space). + """ + from .fitstats import _core_module + + mats = [] + for name, x in (("wins", wins), ("ties", ties)): + if isinstance(x, np.ma.MaskedArray): + raise ValueError(f"bratt_mm: {name} must not be a masked array") + arr = np.asarray(x) + if np.iscomplexobj(arr): + raise ValueError(f"bratt_mm: {name} must be real-valued") + if arr.dtype == object: + if any(isinstance(v, (bool, np.bool_)) for v in arr.flat): + raise ValueError(f"bratt_mm: {name} must be numeric, not boolean") + try: + arr = arr.astype(np.float64) + except (TypeError, ValueError) as exc: + raise ValueError(f"bratt_mm: {name} must be numeric") from exc + if arr.dtype.kind not in "fiu": + raise ValueError(f"bratt_mm: {name} must be numeric (got {arr.dtype})") + arr = np.ascontiguousarray(arr, dtype=np.float64) + if arr.ndim != 2 or arr.shape[0] != arr.shape[1]: + raise ValueError(f"bratt_mm: {name} must be a square 2-D matrix") + mats.append(arr) + w, t = mats + if w.shape != t.shape: + raise ValueError("bratt_mm: wins and ties must have the same shape") + if isinstance(ref_index, bool) or not isinstance(ref_index, (int, np.integer)): + raise ValueError("bratt_mm: ref_index must be an integer") + if ref_index < 0: + raise ValueError("bratt_mm: ref_index must be nonnegative") + n = w.shape[0] + core = _core_module() + res = core.bratt_mm( + w.ravel(), + t.ravel(), + n, + int(ref_index), + float(ref_value), + int(max_iter), + float(tol), + ) + return BrattResult( + alpha=np.asarray(res["alpha"], dtype=np.float64), + alpha0=float(res["alpha0"]), + iterations=int(res["iterations"]), + log_likelihood=float(res["log_likelihood"]), + ) + + @dataclass class LsrResult: """Luce Spectral Ranking fit: ``params[i]`` is the centered log-worth of @@ -1760,3 +1843,312 @@ def fide_rating(games, n_players, init=2200.0, kv=(10.0, 15.0, 30.0), gamma=None elite=np.asarray(res["elite"]), opponent=np.asarray(res["opponent"]), ) + + +def _predict_int_index_array(x, name, fname): + """Validate a player-index array: integers >= -1 (-1 = unmatched).""" + if isinstance(x, np.ma.MaskedArray): + raise ValueError(f"{fname}: masked arrays are not supported for {name}") + raw = x if isinstance(x, np.ndarray) else np.asarray(x, dtype=object) + if np.iscomplexobj(raw) or raw.dtype.kind == "b": + raise ValueError(f"{fname}: {name} must be integer indices, not complex/bool") + if raw.dtype == object: + for v in np.ravel(raw): + if v is None or isinstance( + v, (bool, np.bool_, str, bytes, np.datetime64, np.timedelta64) + ): + raise ValueError(f"{fname}: {name} contains a non-numeric value") + try: + arr = np.asarray(x, dtype=float) + except (TypeError, ValueError) as exc: + raise ValueError(f"{fname}: {name} is not numeric: {exc}") from None + if arr.ndim != 1: + raise ValueError(f"{fname}: {name} must be one-dimensional") + if not np.all(np.isfinite(arr)): + raise ValueError(f"{fname}: {name} must be finite") + if not np.all(arr == np.floor(arr)): + raise ValueError(f"{fname}: {name} must contain integers") + if np.any(arr < -1): + raise ValueError(f"{fname}: {name} indices must be >= -1 (-1 = unmatched)") + return arr.astype(np.int64) + + +def _predict_float_array(x, name, fname, allow_nan): + """Validate a float array; NaN optionally allowed (R NA), Inf rejected.""" + if isinstance(x, np.ma.MaskedArray): + raise ValueError(f"{fname}: masked arrays are not supported for {name}") + raw = x if isinstance(x, np.ndarray) else np.asarray(x, dtype=object) + if np.iscomplexobj(raw) or raw.dtype.kind == "b": + raise ValueError(f"{fname}: {name} must be real numeric, not complex/bool") + if raw.dtype == object: + for v in np.ravel(raw): + if isinstance( + v, (bool, np.bool_, str, bytes, np.datetime64, np.timedelta64) + ) or v is None: + raise ValueError(f"{fname}: {name} contains a non-numeric value") + try: + arr = np.asarray(x, dtype=float) + except (TypeError, ValueError) as exc: + raise ValueError(f"{fname}: {name} is not numeric: {exc}") from None + if arr.ndim != 1: + raise ValueError(f"{fname}: {name} must be one-dimensional") + if np.any(np.isinf(arr)): + raise ValueError(f"{fname}: {name} must not contain infinities") + if not allow_nan and np.any(np.isnan(arr)): + raise ValueError(f"{fname}: {name} must not contain NaN") + return arr + + +def _predict_scalar(x, name, fname): + """Validate a finite real scalar parameter.""" + import math + if isinstance(x, (bool, np.bool_)): + raise ValueError(f"{fname}: {name} must be real numeric, not bool") + raw = np.asarray(x) + if np.iscomplexobj(raw) or raw.dtype.kind == "b": + raise ValueError(f"{fname}: {name} must be real numeric, not complex/bool") + if raw.dtype == object and raw.ndim == 0 and isinstance( + raw.item(), (bool, np.bool_) + ): + raise ValueError(f"{fname}: {name} must be real numeric, not bool") + try: + v = float(x) + except (TypeError, ValueError) as exc: + raise ValueError(f"{fname}: {name} is not numeric: {exc}") from None + if not math.isfinite(v): + raise ValueError(f"{fname}: {name} must be finite") + return v + + +def _predict_games_u64(x, fname): + """Validate games counts losslessly into u64. + + Integer-dtype ndarrays (and Python-int sequences, which numpy keeps in + an integer dtype) pass through without a float round-trip, preserving + counts at or above 2**53. Float inputs are bounded by the source + dtype's exact-integer limit (np.finfo(...).nmant + 1 bits: float64 + 2**53, float32 2**24) so a rounded count can never silently shift the + strict games < tng cutoff.""" + name = "games" + if isinstance(x, np.ma.MaskedArray): + raise ValueError(f"{fname}: masked arrays are not supported for {name}") + if isinstance(x, np.ndarray): + raw = x + else: + probe = np.asarray(x, dtype=object) + if np.iscomplexobj(probe): + raise ValueError(f"{fname}: {name} must be real numeric, not complex/bool") + for v in np.ravel(probe): + if v is None or isinstance( + v, (bool, np.bool_, str, bytes, np.datetime64, np.timedelta64) + ): + raise ValueError(f"{fname}: {name} contains a non-numeric value") + try: + raw = np.asarray(x) + except (TypeError, ValueError, OverflowError) as exc: + raise ValueError(f"{fname}: {name} is not numeric: {exc}") from None + if np.iscomplexobj(raw) or raw.dtype.kind == "b": + raise ValueError(f"{fname}: {name} must be real numeric, not complex/bool") + if raw.ndim != 1: + raise ValueError(f"{fname}: {name} must be one-dimensional") + if raw.dtype.kind in "iu": + if raw.dtype.kind == "i" and np.any(raw < 0): + raise ValueError(f"{fname}: {name} must be nonnegative integers") + return raw.astype(np.uint64) + arr = _predict_float_array(x, name, fname, allow_nan=False) + if np.any(arr < 0) or not np.all(arr == np.floor(arr)): + raise ValueError(f"{fname}: {name} must be nonnegative integers") + if raw.dtype.kind == "f": + fidelity = 2.0 ** (np.finfo(raw.dtype).nmant + 1) + else: + fidelity = 2.0**53 + if np.any(arr >= fidelity): + raise ValueError( + f"{fname}: {name} values at or above {int(fidelity)} are not " + "reliably representable in this float dtype; pass games as an " + "integer array" + ) + return arr.astype(np.uint64) + + +def _predict_tng_u64(tng, fname): + """Validate tng losslessly into u64 (no float round-trip for ints).""" + import math + + if isinstance(tng, (bool, np.bool_)): + raise ValueError(f"{fname}: tng must be real numeric, not bool") + if isinstance(tng, (int, np.integer)): + t = int(tng) + if t < 0: + raise ValueError(f"{fname}: tng must be a nonnegative integer") + if t > 2**64 - 1: + raise ValueError(f"{fname}: tng must fit in an unsigned 64-bit integer") + return t + if ( + isinstance(tng, np.ndarray) + and tng.ndim == 0 + and tng.dtype.kind in "iu" + ): + return _predict_tng_u64(int(tng), fname) + v = _predict_scalar(tng, "tng", fname) + if v < 0 or v != math.floor(v): + raise ValueError(f"{fname}: tng must be a nonnegative integer") + if v >= 2.0**53: + raise ValueError( + f"{fname}: tng values at or above 2**53 are not reliably " + "representable as float; pass tng as an int" + ) + return int(v) + + +def predict_rating( + ratings, + games, + white, + black, + deviations=None, + gamma=30.0, + tng=15, + trat=None, + thresh=None, +): + """Predicted win probabilities for two-player games from fitted ratings. + + Implements the two-player branches of ``predict.rating`` from CRAN + PlayerRatings 1.1-0 (R/ratings.R lines 1056-1133; source read). Without + ``deviations`` this is the Elo branch (used for ``elo()``/``fide()`` + fits); with ``deviations`` it is the deviation-shrunk branch shared by + Glicko, Glicko-2, and Stephenson fits. Players are addressed by index; + ``-1`` marks an unmatched player (R's ``match()`` NA). Players with + ``games < tng`` (strict) are treated as unrated; ``trat`` (a scalar, or + a ``(rating, deviation)`` pair when ``deviations`` is supplied) + replaces ALL missing extracted values, matching R. ``thresh`` maps + predictions to 1 when ``pred >= thresh`` else 0 (NaN kept). ``gamma`` + is the per-game (or scalar, broadcast) first-player advantage. + """ + from .fitstats import _core_module + + fname = "predict_rating" + rat = _predict_float_array(ratings, "ratings", fname, allow_nan=True) + n = rat.shape[0] + if n < 2 or n > 10000: + raise ValueError(f"{fname}: number of players must be in 2..=10000, got {n}") + games_u64 = _predict_games_u64(games, fname) + if games_u64.shape[0] != n: + raise ValueError( + f"{fname}: games must have one entry per player ({n}), " + f"got {games_u64.shape[0]}" + ) + w = _predict_int_index_array(white, "white", fname) + b = _predict_int_index_array(black, "black", fname) + dev = None + if deviations is not None: + dev = _predict_float_array(deviations, "deviations", fname, allow_nan=True) + ng = w.shape[0] + gam_raw = np.asarray(gamma) + if np.iscomplexobj(gam_raw) or gam_raw.dtype.kind == "b": + raise ValueError(f"{fname}: gamma must be real numeric, not complex/bool") + if gam_raw.ndim == 0: + gam = np.full(ng, _predict_scalar(gamma, "gamma", fname)) + else: + gam = _predict_float_array(gamma, "gamma", fname, allow_nan=False) + if gam.shape != (ng,): + raise ValueError( + f"{fname}: gamma must be a scalar or length-{ng} array, got {gam.shape}" + ) + tng_v = _predict_tng_u64(tng, fname) + trat_rating = None + trat_deviation = None + if trat is not None: + if dev is not None: + if not (isinstance(trat, (tuple, list)) and len(trat) == 2): + raise ValueError( + f"{fname}: trat must be a (rating, deviation) pair when " + "deviations are supplied" + ) + trat_rating = _predict_scalar(trat[0], "trat rating", fname) + trat_deviation = _predict_scalar(trat[1], "trat deviation", fname) + else: + if isinstance(trat, (tuple, list)): + if len(trat) != 1: + raise ValueError( + f"{fname}: trat must be a scalar (length 1) without deviations" + ) + trat = trat[0] + trat_rating = _predict_scalar(trat, "trat", fname) + thresh_v = None if thresh is None else _predict_scalar(thresh, "thresh", fname) + out = _core_module().predict_rating_two( + np.ascontiguousarray(rat), + None if dev is None else np.ascontiguousarray(dev), + np.ascontiguousarray(games_u64), + np.ascontiguousarray(w), + np.ascontiguousarray(b), + np.ascontiguousarray(gam), + int(tng_v), + trat_rating, + trat_deviation, + thresh_v, + ) + return np.asarray(out) + + +def predict_rating_multi( + ratings, + games, + players, + tng=15, + trat=None, + placing=False, +): + """Predicted expected scores (or placings) for multi-player EloM events. + + Implements the EloM branch of ``predict.rating`` from CRAN + PlayerRatings 1.1-0 (R/ratings.R lines 1103-1105, 1123-1125, + 1129-1130; source read): per event row, + ``pred = (rating - rowmean) / 40`` with the row mean over non-missing + seats (``rowMeans(rats, na.rm=TRUE)``). ``players`` is an ``(nr, np)`` + index matrix with ``-1`` = empty/unmatched seat. Players with + ``games < tng`` are treated as unrated; scalar ``trat`` replaces all + missing extracted ratings. With ``placing=True`` each row is replaced + by min-tie ranks of the predictions (rank 1 = highest; NaN kept), + matching R's ``rank(-preds, na.last="keep", ties.method="min")``. + """ + from .fitstats import _core_module + + fname = "predict_rating_multi" + rat = _predict_float_array(ratings, "ratings", fname, allow_nan=True) + n = rat.shape[0] + if n < 2 or n > 10000: + raise ValueError(f"{fname}: number of players must be in 2..=10000, got {n}") + games_u64 = _predict_games_u64(games, fname) + if games_u64.shape[0] != n: + raise ValueError( + f"{fname}: games must have one entry per player ({n}), " + f"got {games_u64.shape[0]}" + ) + if isinstance(players, np.ma.MaskedArray): + raise ValueError(f"{fname}: masked arrays are not supported for players") + p_raw = players if isinstance(players, np.ndarray) else np.asarray(players, dtype=object) + if p_raw.ndim != 2: + raise ValueError(f"{fname}: players must be a 2-D (events, seats) matrix") + nr, np_seats = p_raw.shape + flat = _predict_int_index_array(np.ravel(p_raw), "players", fname) + if not (2 <= np_seats <= 1000): + raise ValueError( + f"{fname}: seats per event must be in 2..=1000, got {np_seats}" + ) + tng_v = _predict_tng_u64(tng, fname) + trat_v = None if trat is None else _predict_scalar(trat, "trat", fname) + if not isinstance(placing, (bool, np.bool_)): + raise ValueError(f"{fname}: placing must be a bool") + out = _core_module().predict_rating_multi( + np.ascontiguousarray(rat), + np.ascontiguousarray(games_u64), + np.ascontiguousarray(flat), + int(nr), + int(np_seats), + int(tng_v), + trat_v, + bool(placing), + ) + return np.asarray(out).reshape(nr, np_seats) diff --git a/python/fast_mlsirm/validation.py b/python/fast_mlsirm/validation.py index 61803de66..99d108c7f 100644 --- a/python/fast_mlsirm/validation.py +++ b/python/fast_mlsirm/validation.py @@ -101,3 +101,95 @@ def validate_judge( passed=bool(res["pass"]), failed_gates=[g["name"] for g in gates if not g["pass"]], ) + + +@dataclass +class FleissKappaResult: + """Result of :func:`fleiss_kappa`. In exact (Conger) mode ``z`` and + ``p_value`` are NaN and the category arrays are empty, mirroring irr's + ``kappam.fleiss(exact=TRUE)`` which returns neither.""" + + kappa: float + subjects_used: int + z: float + p_value: float + category_kappa: np.ndarray + category_z: np.ndarray + category_p: np.ndarray + + +def fleiss_kappa( + ratings: np.ndarray, + k: int | None = None, + exact: bool = False, +) -> FleissKappaResult: + """Fleiss' kappa for nominal agreement among multiple raters, with the + exact (Conger) chance-agreement variant. + + Reimplements ``kappam.fleiss()`` from CRAN irr 0.85 (R source READ in + full; algorithm source of truth). Model origins — cited as origins only, + NOT READ: Fleiss, J. L. (1971). Measuring nominal scale agreement among + many raters. *Psychological Bulletin, 76*(5), 378-382; Conger, A. J. + (1980). Integration and generalization of kappas for multiple raters. + *Psychological Bulletin, 88*(2), 322-328. Computation runs in the Rust + core (``mlsirm_core::agreement::fleiss_kappa``). + + ``ratings`` is a 2-D ``(n_subjects, n_raters)`` array of integer category + codes ``0..k-1``. NaN or any negative value marks a missing rating and + drops the whole subject row (listwise, as in irr). ``k=None`` infers + ``max(code)+1``; pass ``k`` explicitly to include trailing empty + categories in the category-wise detail (their kappas are NaN, matching + R's 0/0). + """ + from . import _core # computation lives in the Rust core + + if isinstance(ratings, np.ma.MaskedArray): + raise ValueError("masked arrays are not supported; use NaN for missing") + arr = np.asarray(ratings) + if arr.ndim != 2: + raise ValueError("ratings must be a 2-D (subjects x raters) array") + if np.iscomplexobj(arr): + raise ValueError("ratings must be real-valued") + if arr.dtype == object: + for v in arr.flat: + if v is None or isinstance(v, (bool, np.bool_, str, bytes)): + raise ValueError("ratings must be numeric, not boolean/str/None") + arr = arr.astype(np.float64) + if arr.dtype.kind == "b": + raise ValueError("ratings must be integer codes, not booleans") + if arr.dtype.kind not in "fiu": + raise ValueError(f"ratings dtype {arr.dtype} is not numeric") + ns, nr = arr.shape + if arr.dtype.kind == "f": + finite = np.isfinite(arr) + if np.any(np.isinf(arr)): + raise ValueError("ratings must not contain infinities") + if np.any(arr[finite] != np.floor(arr[finite])): + raise ValueError("ratings must be integer category codes") + if np.any(np.abs(arr[finite]) > 2.0**53): + raise ValueError("ratings exceed exact float64 integer range") + codes = np.where(finite, arr, -1.0).astype(np.int64) + else: + if arr.dtype.kind == "u" and arr.size and int(arr.max()) > np.iinfo(np.int64).max: + raise ValueError("ratings values must fit in int64") + codes = arr.astype(np.int64) + if k is not None: + if isinstance(k, (bool, np.bool_)) or not isinstance(k, (int, np.integer)): + raise ValueError("k must be an integer") + k = int(k) + if k is None: + if codes.size == 0 or int(codes.max()) < 0: + raise ValueError("cannot infer k: no observed category codes") + k = int(codes.max()) + 1 + res = _core.fleiss_kappa( + np.ascontiguousarray(codes.reshape(-1)), int(ns), int(nr), int(k), bool(exact) + ) + return FleissKappaResult( + kappa=float(res["kappa"]), + subjects_used=int(res["subjects_used"]), + z=float(res["z"]), + p_value=float(res["p_value"]), + category_kappa=np.asarray(res["category_kappa"]), + category_z=np.asarray(res["category_z"]), + category_p=np.asarray(res["category_p"]), + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index e45cdc783..dd77ec801 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10067,3 +10067,371 @@ def test_elo_reduction(self): assert f.ratings.tolist() == e.ratings.tolist() assert f.games.tolist() == e.games.tolist() assert f.lag.tolist() == e.lag.tolist() + + +class TestPredict: + """predict_rating / predict_rating_multi (PlayerRatings predict.rating, + R lines 1056-1133). Pins from the executed oracle + (files/predict_oracle.py); every assert reads wrapper/crate outputs.""" + + def test_elo_branch_pins(self): + import numpy as np + from fast_mlsirm import predict_rating + + p1 = predict_rating([2200.0, 2200.0], [20, 20], [0], [1], gamma=0.0) + assert p1.shape == (1,) + assert p1[0] == 0.5 + p2 = predict_rating([2200.0, 2000.0], [20, 20], [0], [1], gamma=30.0) + assert abs(p2[0] - 0.7898441797581306) < 1e-14 + # per-game gamma array and scalar broadcast agree (crate values) + p2b = predict_rating( + [2200.0, 2000.0], [20, 20], [0, 0], [1, 1], gamma=np.array([30.0, 30.0]) + ) + assert p2b[0] == p2[0] and p2b[1] == p2[0] + + def test_deviation_branch(self): + from fast_mlsirm import predict_rating + + p3 = predict_rating( + [2200.0, 2000.0], [20, 20], [0], [1], + deviations=[50.0, 100.0], gamma=30.0, + ) + assert abs(p3[0] - 0.776912664201114) < 1e-14 + p3z = predict_rating( + [2200.0, 2000.0], [20, 20], [0], [1], + deviations=[0.0, 0.0], gamma=30.0, + ) + p2 = predict_rating([2200.0, 2000.0], [20, 20], [0], [1], gamma=30.0) + assert p3z[0] == p2[0] + + def test_tng_trat_unmatched(self): + import math + from fast_mlsirm import predict_rating + + kept = predict_rating([2200.0, 2000.0], [15, 15], [0], [1], gamma=0.0) + assert abs(kept[0] - 0.7597469266479578) < 1e-14 + dropped = predict_rating([2200.0, 2000.0], [14, 15], [0], [1], gamma=0.0) + assert math.isnan(dropped[0]) + replaced = predict_rating( + [2200.0, 2000.0], [14, 15], [0], [1], gamma=0.0, trat=2000.0 + ) + assert replaced[0] == 0.5 + un = predict_rating([2200.0, 2000.0], [20, 20], [-1], [1], gamma=0.0) + assert math.isnan(un[0]) + un_t = predict_rating( + [2200.0, 2000.0], [20, 20], [-1], [1], gamma=0.0, trat=2000.0 + ) + assert un_t[0] == 0.5 + # stored-NaN + pair trat with deviations + p9c = predict_rating( + [2200.0, 2000.0], [20, 20], [0], [1], + deviations=[float("nan"), 50.0], gamma=30.0, trat=(2200.0, 50.0), + ) + assert abs(p9c[0] - 0.7844611342833985) < 1e-14 + + def test_thresh(self): + import math + from fast_mlsirm import predict_rating + + p = predict_rating( + [2200.0, 2200.0], [20, 20], [0, -1], [1, 1], gamma=0.0, thresh=0.5 + ) + assert p[0] == 1.0 # exact equality maps to 1 (>=) + assert math.isnan(p[1]) + + def test_multi_and_placing(self): + import math + import numpy as np + from fast_mlsirm import predict_rating_multi + + p7 = predict_rating_multi( + [2300.0, 2200.0, 2100.0, 2000.0], [20] * 4, + np.array([[0, 1, 2], [0, 3, -1]]), + ) + assert p7.shape == (2, 3) + assert p7[0].tolist() == [2.5, 0.0, -2.5] + assert p7[1][0] == 3.75 and p7[1][1] == -3.75 and math.isnan(p7[1][2]) + p8 = predict_rating_multi( + [2300.0, 2300.0, 2100.0, 2000.0], [20] * 4, + np.array([[0, 1, 2, -1]]), placing=True, + ) + assert p8[0][0] == 1.0 and p8[0][1] == 1.0 and p8[0][2] == 3.0 + assert math.isnan(p8[0][3]) + + def test_validation(self): + import numpy as np + import pytest + from fast_mlsirm import predict_rating, predict_rating_multi + + with pytest.raises(ValueError, match="2..=10000"): + predict_rating([2200.0], [20], [0], [0], gamma=0.0) + with pytest.raises(ValueError, match="complex"): + predict_rating([2200.0 + 1j, 2000.0], [20, 20], [0], [1]) + with pytest.raises(ValueError, match="complex|bool"): + predict_rating(np.array([True, False]), [20, 20], [0], [1]) + with pytest.raises(ValueError, match="non-numeric"): + predict_rating([2200.0, "x"], [20, 20], [0], [1]) + with pytest.raises(ValueError, match="masked"): + predict_rating( + np.ma.masked_array([2200.0, 2000.0], mask=[False, True]), + [20, 20], [0], [1], + ) + with pytest.raises(ValueError, match="infinit"): + predict_rating([np.inf, 2000.0], [20, 20], [0], [1]) + with pytest.raises(ValueError, match="NaN"): + predict_rating([2200.0, 2000.0], [np.nan, 20], [0], [1]) + with pytest.raises(ValueError, match=">= -1"): + predict_rating([2200.0, 2000.0], [20, 20], [-2], [1]) + with pytest.raises(ValueError, match="self-play"): + predict_rating([2200.0, 2000.0], [20, 20], [1], [1]) + with pytest.raises(ValueError, match="pair"): + predict_rating( + [2200.0, 2000.0], [20, 20], [0], [1], + deviations=[50.0, 50.0], trat=2000.0, + ) + with pytest.raises(ValueError, match="bool"): + predict_rating([2200.0, 2000.0], [20, 20], [0], [1], thresh=True) + with pytest.raises(ValueError, match="2-D"): + predict_rating_multi([2200.0, 2000.0], [20, 20], [0, 1]) + with pytest.raises(ValueError, match="bool"): + predict_rating_multi( + [2200.0, 2000.0], [20, 20], np.array([[0, 1]]), trat=True + ) + # NaN ratings are allowed (R NA), not an error + import math + out = predict_rating([float("nan"), 2000.0], [20, 20], [0], [1], gamma=0.0) + assert math.isnan(out[0]) + + def test_u64_fidelity_and_overflow(self): + import math + import numpy as np + import pytest + from fast_mlsirm import predict_rating, predict_rating_multi + + big = 2**53 + # int-dtype games keep exact counts above 2**53: games == 2**53 is + # strictly below tng == 2**53 + 1, so both players are unrated -> NaN. + out = predict_rating( + [2200.0, 2000.0], [big, big], [0], [1], gamma=0.0, tng=big + 1 + ) + assert math.isnan(out[0]) + outm = predict_rating_multi( + [2300.0, 2200.0], [big, big], np.array([[0, 1]]), tng=big + 1 + ) + assert math.isnan(outm[0][0]) and math.isnan(outm[0][1]) + # float games at/above the dtype exact-integer bound are rejected + with pytest.raises(ValueError, match='integer array'): + predict_rating( + [2200.0, 2000.0], np.array([float(big), 20.0]), [0], [1] + ) + with pytest.raises(ValueError, match='integer array'): + predict_rating_multi( + [2200.0, 2000.0], np.array([float(big), 20.0]), np.array([[0, 1]]) + ) + # float tng at/above 2**53 rejected; huge int tng still exact + with pytest.raises(ValueError, match='pass tng as an int'): + predict_rating([2200.0, 2000.0], [20, 20], [0], [1], tng=float(big)) + kept = predict_rating( + [2200.0, 2000.0], [big + 2, big + 2], [0], [1], gamma=0.0, tng=big + 1 + ) + assert abs(kept[0] - 0.7597469266479578) < 1e-14 + # core-level nr*np overflow is a checked ValueError, not a panic + import fast_mlsirm._core as core + with pytest.raises(ValueError, match='overflows'): + core.predict_rating_multi( + np.array([1.0, 2.0]), + np.array([1, 1], dtype=np.uint64), + np.array([], dtype=np.int64), + 1 << 61, + 8, + 15, + None, + False, + ) + + +class TestBratt: + """bratt_mm: Bradley-Terry with additive-alpha0 ties (VGAM bratt). + + Oracle: exact-Fraction MM iterates (files/bratt_oracle.py, session + evidence). Every assert reads values returned by the crate through the + binding. + """ + + @staticmethod + def _fixture(): + import numpy as np + + y = np.array([[0, 3, 1], [1, 0, 2], [2, 1, 0]], dtype=float) + t = np.array([[0, 1, 1], [1, 0, 2], [1, 2, 0]], dtype=float) + return y, t + + def test_b1_iter1_exact(self): + import numpy as np + from fast_mlsirm import bratt_mm + + y, t = self._fixture() + r = bratt_mm(y, t, ref_index=0, ref_value=1.0, max_iter=100, tol=0.5) + assert r.iterations == 1 + np.testing.assert_allclose(r.alpha, [1.0, 27.0 / 40.0, 3.0 / 4.0], rtol=1e-15) + np.testing.assert_allclose(r.alpha0, 9.0 / 14.0, rtol=1e-15) + + def test_b2_converged(self): + import numpy as np + from fast_mlsirm import bratt_mm + + y, t = self._fixture() + r = bratt_mm(y, t, tol=1e-13) + np.testing.assert_allclose( + r.alpha, [1.0, 0.6150318884241122, 0.686344995662376], rtol=1e-10 + ) + np.testing.assert_allclose(r.alpha0, 0.6038293879270596, rtol=1e-10) + np.testing.assert_allclose(r.log_likelihood, -15.12765635227613, rtol=1e-12) + + def test_b4_reference_rescale(self): + import numpy as np + from fast_mlsirm import bratt_mm + + y, t = self._fixture() + r = bratt_mm(y, t, tol=1e-13) + r2 = bratt_mm(y, t, ref_index=1, ref_value=2.0, tol=1e-13) + c = 2.0 / r.alpha[1] + np.testing.assert_allclose(r2.alpha, r.alpha * c, rtol=1e-12) + np.testing.assert_allclose(r2.alpha0, r.alpha0 * c, rtol=1e-12) + np.testing.assert_allclose(r2.log_likelihood, r.log_likelihood, rtol=1e-12) + + def test_validation(self): + import numpy as np + import pytest + from fast_mlsirm import bratt_mm + + y, t = self._fixture() + with pytest.raises(ValueError, match="masked"): + bratt_mm(np.ma.masked_array(y), t) + with pytest.raises(ValueError, match="real-valued"): + bratt_mm(y.astype(complex), t) + with pytest.raises(ValueError, match="numeric"): + bratt_mm(y.astype(bool), t) + with pytest.raises(ValueError, match="numeric"): + bratt_mm(np.array([["a"] * 3] * 3, dtype=object), t) + with pytest.raises(ValueError, match="square"): + bratt_mm(y[:2], t) + with pytest.raises(ValueError, match="same shape"): + bratt_mm(y, np.zeros((4, 4))) + with pytest.raises(ValueError, match="ref_index must be an integer"): + bratt_mm(y, t, ref_index=0.5) + with pytest.raises(ValueError, match="ref_index must be an integer"): + bratt_mm(y, t, ref_index=True) + with pytest.raises(ValueError, match="nonnegative"): + bratt_mm(y, t, ref_index=-1) + # Crate error contract surfaces through the binding: + with pytest.raises(ValueError, match="use bradley_terry_mm"): + bratt_mm(y, np.zeros((3, 3))) + with pytest.raises(ValueError, match="symmetric"): + tt = t.copy() + tt[0, 1] = 5.0 + bratt_mm(y, tt) + with pytest.raises(ValueError, match="no wins"): + y0 = y.copy() + y0[1] = 0.0 + bratt_mm(y0, t) + with pytest.raises(ValueError, match="did not converge"): + bratt_mm(y, t, max_iter=1, tol=1e-15) + # Impl-review regressions: + with pytest.raises(ValueError, match="not boolean"): + yb = np.full((3, 3), True, dtype=object) + np.fill_diagonal(yb, False) + tb = np.full((3, 3), True, dtype=object) + np.fill_diagonal(tb, False) + bratt_mm(yb, tb) + with pytest.raises(ValueError, match="overflow"): + h = 9e307 + yh = np.array([[0, 1, 1], [1, 0, h], [1, h, 0]], dtype=float) + bratt_mm(yh, t, max_iter=3, tol=2.0) + + +class TestFleiss: + """fleiss_kappa binding + wrapper (irr 0.85 kappam.fleiss oracle anchors).""" + + @staticmethod + def _fk(): + import numpy as np + + return np.array( + [[0, 0, 0, 1], [0, 1, 1, 1], [2, 2, 2, 2], [0, 0, 2, 2], [1, 1, 1, 0]] + ) + + def test_classic_anchor(self): + import numpy as np + + from fast_mlsirm import fleiss_kappa + + res = fleiss_kappa(self._fk()) + assert abs(res.kappa - 139.0 / 399.0) < 1e-14 + assert res.subjects_used == 5 + assert abs(res.z - 2.694739854085488) < 1e-10 + assert abs(res.p_value - 0.007044360582468963) < 5e-7 + np.testing.assert_allclose( + res.category_kappa, [1 / 21, 31 / 91, 43 / 63], rtol=1e-14 + ) + + def test_exact_and_missing(self): + import numpy as np + + from fast_mlsirm import fleiss_kappa + + res = fleiss_kappa(self._fk(), exact=True) + assert abs(res.kappa - 37.0 / 102.0) < 1e-14 + assert np.isnan(res.z) and np.isnan(res.p_value) + assert res.category_kappa.size == 0 + # NaN row drops listwise and reproduces the classic anchor. + withnan = np.vstack([[0.0, np.nan, 1.0, 2.0], self._fk().astype(float)]) + res2 = fleiss_kappa(withnan) + assert res2.subjects_used == 5 + assert abs(res2.kappa - 139.0 / 399.0) < 1e-14 + # Explicit k adds an empty category with NaN detail. + res3 = fleiss_kappa(self._fk(), k=4) + assert np.isnan(res3.category_kappa[3]) + assert abs(res3.kappa - 139.0 / 399.0) < 1e-14 + + def test_validation(self): + import numpy as np + import pytest + + from fast_mlsirm import fleiss_kappa + + fk = self._fk() + with pytest.raises(ValueError, match="2-D"): + fleiss_kappa(fk.reshape(-1)) + with pytest.raises(ValueError, match="masked"): + fleiss_kappa(np.ma.masked_array(fk, mask=False)) + with pytest.raises(ValueError, match="complex|real"): + fleiss_kappa(fk.astype(complex)) + with pytest.raises(ValueError, match="boolean"): + fleiss_kappa(np.array([[True, False], [False, True]], dtype=object)) + with pytest.raises(ValueError, match="boolean"): + fleiss_kappa(np.array([[True, False], [False, True]])) + with pytest.raises(ValueError, match="integer"): + fleiss_kappa(fk + 0.5) + with pytest.raises(ValueError, match="infinit"): + fleiss_kappa(np.array([[np.inf, 0.0], [1.0, 0.0]])) + with pytest.raises(ValueError, match="raters"): + fleiss_kappa(fk[:, :1]) + with pytest.raises(ValueError, match="infer"): + fleiss_kappa(np.full((2, 2), np.nan)) + with pytest.raises(ValueError, match="dropped"): + fleiss_kappa(np.full((2, 2), np.nan), k=2) + with pytest.raises(ValueError, match="degenerate"): + fleiss_kappa(np.ones((3, 2), dtype=int), k=2) + # uint64 above i64::MAX must be rejected, not wrapped negative + # (silent listwise drop). Reads the wrapper's guard, killed by + # removing the unsigned-range check before astype(int64). + big = np.array([[2**63, 0], [0, 1], [1, 1]], dtype=np.uint64) + with pytest.raises(ValueError, match="int64"): + fleiss_kappa(big, k=2) + # Explicit k must be a true integer, not lossy-coerced. + for bad_k in (3.9, "3", np.float64(3.0), True): + with pytest.raises(ValueError, match="k must be an integer"): + fleiss_kappa(fk, k=bad_k) + # np.integer k still accepted. + assert fleiss_kappa(fk, k=np.int64(3)).subjects_used == 5 diff --git a/tests/unit/agreement_tests.rs b/tests/unit/agreement_tests.rs index 0b3ae850f..43dd61743 100644 --- a/tests/unit/agreement_tests.rs +++ b/tests/unit/agreement_tests.rs @@ -134,3 +134,238 @@ fn rejects_degenerate_inputs() { .iter() .any(|gate| gate.name == "subgroup_smd")); } + +// --------------------------------------------------------------------------- +// fleiss_kappa (irr 0.85 kappam.fleiss; oracle anchors FK1-FK5 with exact +// Fractions). Every assert below reads crate outputs (FleissKappaResult +// fields returned by fleiss_kappa); pins are hand-derived rationals. +// z pins are pure arithmetic (rel 1e-12); p-value pins go through the +// crate's Numerical-Recipes erfc (|error| < 1.2e-7), hence abs 5e-7. +// --------------------------------------------------------------------------- + +/// FK fixture: ns=5, nr=4, k=3 (asymmetric; 5x4 also breaks stride +/// transposition and row-vs-column chanceP confusion). +fn fk_fixture() -> Vec { + vec![ + 0, 0, 0, 1, // S1 + 0, 1, 1, 1, // S2 + 2, 2, 2, 2, // S3 + 0, 0, 2, 2, // S4 + 1, 1, 1, 0, // S5 + ] +} + +fn rel_eq(a: f64, b: f64, tol: f64) -> bool { + (a - b).abs() <= tol * b.abs().max(1.0) +} + +#[test] +fn fk_anchor_fk1_classic() { + // Kills MU1 (agreeP drops -nr centering: kappa would be 113/133), + // MU2 (row-sum chanceP: 1/5 -> kappa 11/24), and via the z pin MU4 + // (variance second-term sign: Sum p q (q-p) = 441/2000 != 0). + let res = fleiss_kappa(&fk_fixture(), 5, 4, 3, false).unwrap(); + assert!( + rel_eq(res.kappa, 139.0 / 399.0, 1e-15), + "kappa {}", + res.kappa + ); + assert_eq!(res.subjects_used, 5); + // var = 181/10830; z = kappa/sqrt(var) (oracle 2.694739854085488). + assert!( + rel_eq(res.z, (139.0 / 399.0) / (181.0_f64 / 10830.0).sqrt(), 1e-12), + "z {}", + res.z + ); + assert!( + (res.p_value - 0.007044360582468963).abs() < 5e-7, + "p {}", + res.p_value + ); +} + +#[test] +fn fk_exact_fk2() { + // Kills MU3 (exact == classic): Conger correction Sum s2_j = 3/50 gives + // chanceP 8/25 and kappa 37/102 != 139/399. + let res = fleiss_kappa(&fk_fixture(), 5, 4, 3, true).unwrap(); + assert!( + rel_eq(res.kappa, 37.0 / 102.0, 1e-15), + "exact kappa {}", + res.kappa + ); + assert_eq!(res.subjects_used, 5); + // Exact mode returns no test statistic or detail (irr returns neither). + assert!(res.z.is_nan() && res.p_value.is_nan()); + assert!(res.category_kappa.is_empty()); + assert!(res.category_z.is_empty()); + assert!(res.category_p.is_empty()); +} + +#[test] +fn fk_missing_drop_fk3() { + // Kills MU5 (missing code counted as a category instead of listwise row + // drop): the prepended row must vanish, leaving FK1 exactly. + let mut ratings = vec![0, -1, 1, 2]; + ratings.extend(fk_fixture()); + let res = fleiss_kappa(&ratings, 6, 4, 3, false).unwrap(); + assert_eq!(res.subjects_used, 5); + let base = fleiss_kappa(&fk_fixture(), 5, 4, 3, false).unwrap(); + assert_eq!(res.kappa, base.kappa, "drop must reproduce FK1 bitwise"); + assert!(rel_eq(res.kappa, 139.0 / 399.0, 1e-15)); +} + +#[test] +fn fk_category_detail_fk4() { + // Kills MU6 (pjk drops the m*nr*p_j subtraction: category kappas would + // be [51/91, 233/273, 73/63]). + let res = fleiss_kappa(&fk_fixture(), 5, 4, 3, false).unwrap(); + let expect = [1.0 / 21.0, 31.0 / 91.0, 43.0 / 63.0]; + assert_eq!(res.category_kappa.len(), 3); + for j in 0..3 { + assert!( + rel_eq(res.category_kappa[j], expect[j], 1e-15), + "kappa_{j} {}", + res.category_kappa[j] + ); + // var_j = 1/30 for all j -> z_j = kappa_j * sqrt(30). + assert!( + rel_eq(res.category_z[j], expect[j] * 30.0_f64.sqrt(), 1e-12), + "z_{j} {}", + res.category_z[j] + ); + } + // Oracle p-values (through math.erfc; crate erfc abs 5e-7). + let expect_p = [ + 0.7942311156261253, + 0.062059828219847915, + 0.00018517759699332675, + ]; + for j in 0..3 { + assert!( + (res.category_p[j] - expect_p[j]).abs() < 5e-7, + "p_{j} {}", + res.category_p[j] + ); + } +} + +#[test] +fn fk_empty_category_fk5() { + // k=4 with no code 3: R's 0/0 -> NaN preserved; overall kappa unchanged + // (C_3 = 0 contributes nothing to chanceP). + let res = fleiss_kappa(&fk_fixture(), 5, 4, 4, false).unwrap(); + assert!(rel_eq(res.kappa, 139.0 / 399.0, 1e-15)); + assert_eq!(res.category_kappa.len(), 4); + assert!(res.category_kappa[3].is_nan()); + assert!(res.category_z[3].is_nan()); + assert!(res.category_p[3].is_nan()); + // Non-empty categories unaffected by the extra level. + assert!(rel_eq(res.category_kappa[0], 1.0 / 21.0, 1e-15)); +} + +#[test] +fn fk_error_contract() { + let fk = fk_fixture(); + assert!(fleiss_kappa(&fk, 0, 4, 3, false).is_err(), "ns == 0"); + assert!(fleiss_kappa(&fk, 5, 1, 3, false).is_err(), "nr < 2"); + assert!(fleiss_kappa(&fk, 5, 4, 1, false).is_err(), "k < 2"); + assert!( + fleiss_kappa(&[], 2_000_000, 2, 2, false).is_err(), + "ns cap before ns*nr" + ); + assert!(fleiss_kappa(&fk, 5, 4, 4, false).is_ok()); + assert!( + fleiss_kappa(&fk[..19], 5, 4, 3, false).is_err(), + "length mismatch" + ); + assert!( + fleiss_kappa(&[0, 1, 3, 0, 1, 2], 3, 2, 3, false).is_err(), + "code >= k" + ); + assert!( + fleiss_kappa(&[-1, 0, 1, -1], 2, 2, 2, false).is_err(), + "all rows dropped" + ); + assert!( + fleiss_kappa(&[1, 1, 1, 1, 1, 1], 3, 2, 2, false).is_err(), + "degenerate chanceP == 1" + ); +} + +/// 500-rep invariance MC: kappa (classic and exact), z, and category detail +/// are invariant under subject-row permutation; kappa and z are invariant +/// under rater-column permutation (ttab is unchanged by either; rtab is +/// permuted across raters, leaving Sum_j s2_j unchanged). Asserts compare +/// two crate outputs, so any asymmetry introduced into the aggregation +/// breaks them. Disclosure: this cannot detect a wrong-but-symmetric +/// formula; the FK1-FK5 value pins above are the discriminating anchors. +#[test] +#[ignore] +fn fk_mc_500_permutation_invariance() { + struct Lcg(u64); + impl Lcg { + fn next_u64(&mut self) -> u64 { + self.0 = self + .0 + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + self.0 + } + fn below(&mut self, n: usize) -> usize { + (self.next_u64() >> 33) as usize % n + } + } + let mut rng = Lcg(0x5eed_f1e5_5000_0001); + let (ns, nr, k) = (8usize, 5usize, 4usize); + for rep in 0..500 { + let ratings: Vec = (0..ns * nr).map(|_| rng.below(k) as i64).collect(); + let Ok(base) = fleiss_kappa(&ratings, ns, nr, k, false) else { + continue; // degenerate draw + }; + let base_ex = fleiss_kappa(&ratings, ns, nr, k, true).unwrap(); + // Fisher-Yates over subject rows. + let mut order: Vec = (0..ns).collect(); + for i in (1..ns).rev() { + order.swap(i, rng.below(i + 1)); + } + let by_subj: Vec = order + .iter() + .flat_map(|&i| ratings[i * nr..(i + 1) * nr].iter().copied()) + .collect(); + let perm = fleiss_kappa(&by_subj, ns, nr, k, false).unwrap(); + let perm_ex = fleiss_kappa(&by_subj, ns, nr, k, true).unwrap(); + assert!( + rel_eq(perm.kappa, base.kappa, 1e-12), + "rep {rep} subject-perm kappa" + ); + assert!(rel_eq(perm.z, base.z, 1e-12), "rep {rep} subject-perm z"); + assert!( + rel_eq(perm_ex.kappa, base_ex.kappa, 1e-12), + "rep {rep} subject-perm exact" + ); + for j in 0..k { + let (a, b) = (perm.category_kappa[j], base.category_kappa[j]); + assert!( + (a.is_nan() && b.is_nan()) || rel_eq(a, b, 1e-12), + "rep {rep} cat {j}" + ); + } + // Fisher-Yates over rater columns. + let mut rorder: Vec = (0..nr).collect(); + for i in (1..nr).rev() { + rorder.swap(i, rng.below(i + 1)); + } + let by_rater: Vec = (0..ns) + .flat_map(|i| { + let ratings = &ratings; + rorder.iter().map(move |&r| ratings[i * nr + r]) + }) + .collect(); + let rperm = fleiss_kappa(&by_rater, ns, nr, k, true).unwrap(); + assert!( + rel_eq(rperm.kappa, base_ex.kappa, 1e-12), + "rep {rep} rater-perm exact kappa" + ); + } +} diff --git a/tests/unit/scaling_tests.rs b/tests/unit/scaling_tests.rs index 1371cc86d..4cc6ccd84 100644 --- a/tests/unit/scaling_tests.rs +++ b/tests/unit/scaling_tests.rs @@ -5615,3 +5615,968 @@ fn fd_mc_500_elo_reduction() { assert_eq!(f.lag, e.lag, "rep {rep}"); } } + +// --------------------------------------------------------------------------- +// predict_rating_two / predict_rating_multi (PlayerRatings predict.rating, +// R lines 1056-1133 READ). Pins from the executed oracle +// (files/predict_oracle.py / predict_oracle_output.txt). +// --------------------------------------------------------------------------- + +/// P1+P2: Elo branch exact half + gamma-sign pin. +/// Asserts read `predict_rating_two` return values. Killing mutants: +/// gamma sign flip (MU1) changes the P2 pin; equal-rating case pins 1/2. +#[test] +fn pr_anchor_elo_gamma_sign() { + let p1 = predict_rating_two( + &[2200.0, 2200.0], + None, + &[20, 20], + &[0], + &[1], + &[0.0], + 15, + None, + None, + ) + .unwrap(); + assert_eq!(p1, vec![0.5]); + let p2 = predict_rating_two( + &[2200.0, 2000.0], + None, + &[20, 20], + &[0], + &[1], + &[30.0], + 15, + None, + None, + ) + .unwrap(); + assert!((p2[0] - 0.7898441797581306).abs() < 1e-14, "{}", p2[0]); + // MU1 (gamma sign flip) would give ~0.7276; assert distance from it. + let mu1 = 1.0 / (1.0 + 10f64.powf((2000.0 - 2200.0 + 30.0) / 400.0)); + assert!((p2[0] - mu1).abs() > 1e-3); +} + +/// P3: deviation-branch shrink pin + dev=0 reduction to the Elo branch. +/// Asserts read `predict_rating_two` values from BOTH branches. Killing +/// mutants: qip3 factor 3->2 (MU2) changes the nonzero-dev pin; the dev=0 +/// crate-vs-crate identity kills stray vec offsets. +#[test] +fn pr_deviation_shrink() { + let p3 = predict_rating_two( + &[2200.0, 2000.0], + Some(&[50.0, 100.0]), + &[20, 20], + &[0], + &[1], + &[30.0], + 15, + None, + None, + ) + .unwrap(); + assert!((p3[0] - 0.776912664201114).abs() < 1e-14, "{}", p3[0]); + let p3z = predict_rating_two( + &[2200.0, 2000.0], + Some(&[0.0, 0.0]), + &[20, 20], + &[0], + &[1], + &[30.0], + 15, + None, + None, + ) + .unwrap(); + let p2 = predict_rating_two( + &[2200.0, 2000.0], + None, + &[20, 20], + &[0], + &[1], + &[30.0], + 15, + None, + None, + ) + .unwrap(); + assert_eq!(p3z[0], p2[0], "dev=0 must reduce to the Elo branch"); +} + +/// P4: tng boundary — games == tng is KEPT (strict <); games < tng is +/// replaced by trat or NaN. Asserts read `predict_rating_two` values. +/// Killing mutant: `<` -> `<=` (MU3) turns the kept pin into NaN. +#[test] +fn pr_tng_boundary() { + let kept = predict_rating_two( + &[2200.0, 2000.0], + None, + &[15, 15], + &[0], + &[1], + &[0.0], + 15, + None, + None, + ) + .unwrap(); + assert!((kept[0] - 0.7597469266479578).abs() < 1e-14, "{}", kept[0]); + let dropped = predict_rating_two( + &[2200.0, 2000.0], + None, + &[14, 15], + &[0], + &[1], + &[0.0], + 15, + None, + None, + ) + .unwrap(); + assert!(dropped[0].is_nan()); + let replaced = predict_rating_two( + &[2200.0, 2000.0], + None, + &[14, 15], + &[0], + &[1], + &[0.0], + 15, + Some((2000.0, 0.0)), + None, + ) + .unwrap(); + assert_eq!(replaced[0], 0.5); +} + +/// P5+P9: unmatched (-1) and matched-but-stored-NaN players — trat +/// replaces ALL extracted NAs; without trat they propagate. Asserts read +/// `predict_rating_two` values (incl. a crate-vs-crate deviation anchor). +#[test] +fn pr_unmatched_and_stored_na_trat() { + let una = predict_rating_two( + &[2200.0, 2000.0], + None, + &[20, 20], + &[-1], + &[1], + &[0.0], + 15, + None, + None, + ) + .unwrap(); + assert!(una[0].is_nan()); + let unb = predict_rating_two( + &[2200.0, 2000.0], + None, + &[20, 20], + &[-1], + &[1], + &[0.0], + 15, + Some((2000.0, 0.0)), + None, + ) + .unwrap(); + assert_eq!(unb[0], 0.5); + let p9a = predict_rating_two( + &[f64::NAN, 2000.0], + None, + &[20, 20], + &[0], + &[1], + &[0.0], + 15, + Some((2000.0, 0.0)), + None, + ) + .unwrap(); + assert_eq!(p9a[0], 0.5); + let p9b = predict_rating_two( + &[f64::NAN, 2000.0], + None, + &[20, 20], + &[0], + &[1], + &[0.0], + 15, + None, + None, + ) + .unwrap(); + assert!(p9b[0].is_nan()); + // stored-NaN deviation replaced by trat.1 == crate value with real dev. + let p9c = predict_rating_two( + &[2200.0, 2000.0], + Some(&[f64::NAN, 50.0]), + &[20, 20], + &[0], + &[1], + &[30.0], + 15, + Some((2200.0, 50.0)), + None, + ) + .unwrap(); + let p9c_ref = predict_rating_two( + &[2200.0, 2000.0], + Some(&[50.0, 50.0]), + &[20, 20], + &[0], + &[1], + &[30.0], + 15, + None, + None, + ) + .unwrap(); + assert_eq!(p9c[0], p9c_ref[0]); + assert!((p9c[0] - 0.7844611342833985).abs() < 1e-14, "{}", p9c[0]); +} + +/// P6: thresh uses >= (exact equality -> 1) and NaN preds stay NaN. +/// Asserts read `predict_rating_two` values. Killing mutant: `>=` -> `>` +/// (MU6) turns the exact-equality 1.0 into 0.0. +#[test] +fn pr_thresh_ge_and_nan() { + let p = predict_rating_two( + &[2200.0, 2200.0], + None, + &[20, 20], + &[0, -1], + &[1, 1], + &[0.0, 0.0], + 15, + None, + Some(0.5), + ) + .unwrap(); + assert_eq!(p[0], 1.0, "pred exactly == thresh must map to 1 (>=)"); + assert!(p[1].is_nan()); + let below = predict_rating_two( + &[2000.0, 2200.0], + None, + &[20, 20], + &[0], + &[1], + &[0.0], + 15, + None, + Some(0.5), + ) + .unwrap(); + assert_eq!(below[0], 0.0); +} + +/// P7: EloM (rat - rowmean)/40 pins, na.rm rowmean, all-NaN row. +/// Asserts read `predict_rating_multi` values. Killing mutant: divisor +/// 40 -> 400 (MU4) changes the exact 2.5/3.75 pins. +#[test] +fn pr_elom_rowmean() { + let p7 = predict_rating_multi( + &[2300.0, 2200.0, 2100.0, 2000.0], + &[20, 20, 20, 20], + &[0, 1, 2, 0, 3, -1], + 2, + 3, + 15, + None, + false, + ) + .unwrap(); + assert_eq!(&p7[0..3], &[2.5, 0.0, -2.5]); + assert_eq!(p7[3], 3.75); + assert_eq!(p7[4], -3.75); + assert!(p7[5].is_nan()); + let allna = predict_rating_multi( + &[2300.0, 2200.0], + &[20, 20], + &[-1, -1], + 1, + 2, + 15, + None, + false, + ) + .unwrap(); + assert!(allna.iter().all(|v| v.is_nan())); + // P9d: stored-NaN seat replaced by trat. + let p9d = predict_rating_multi( + &[f64::NAN, 2200.0, 2100.0], + &[20, 20, 20], + &[0, 1, 2], + 1, + 3, + 15, + Some(2300.0), + false, + ) + .unwrap(); + // tng boundary in the multi branch: games == tng KEPT (strict <). + let boundary = predict_rating_multi( + &[2300.0, 2200.0, 2100.0], + &[15, 20, 20], + &[0, 1, 2], + 1, + 3, + 15, + None, + false, + ) + .unwrap(); + assert_eq!(&boundary[0..3], &[2.5, 0.0, -2.5]); + let dropped = predict_rating_multi( + &[2300.0, 2200.0, 2100.0], + &[14, 20, 20], + &[0, 1, 2], + 1, + 3, + 15, + None, + false, + ) + .unwrap(); + assert!(dropped[0].is_nan()); + assert_eq!(dropped[1], 1.25); + assert_eq!(dropped[2], -1.25); +} + +/// P8: placing ranks — ties share the MINIMUM rank, NaN kept. +/// Asserts read `predict_rating_multi` values. Killing mutant: min -> +/// average/max tie handling (MU5) changes the (1,1,3) pattern. +#[test] +fn pr_placing_min_ties() { + let p8 = predict_rating_multi( + &[2300.0, 2300.0, 2100.0, 2000.0], + &[20, 20, 20, 20], + &[0, 1, 2, -1], + 1, + 4, + 15, + None, + true, + ) + .unwrap(); + assert_eq!(p8[0], 1.0); + assert_eq!(p8[1], 1.0); + assert_eq!(p8[2], 3.0); + assert!(p8[3].is_nan()); +} + +/// Error contract for both functions. Asserts read Err values. +#[test] +fn pr_error_contract() { + let r = |x: Result, String>| x.unwrap_err(); + // player count bounds + assert!(r(predict_rating_two( + &[1500.0], + None, + &[0], + &[0], + &[0], + &[0.0], + 15, + None, + None + )) + .contains("2..=10000")); + // games length mismatch + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0], + &[0], + &[1], + &[0.0], + 15, + None, + None + )) + .contains("games length")); + // empty game rows + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[], + &[], + &[], + 15, + None, + None + )) + .contains("at least one game")); + // infinite rating rejected (NaN allowed elsewhere) + assert!(r(predict_rating_two( + &[f64::INFINITY, 1500.0], + None, + &[0, 0], + &[0], + &[1], + &[0.0], + 15, + None, + None + )) + .contains("infinite")); + // non-finite gamma + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[0], + &[1], + &[f64::NAN], + 15, + None, + None + )) + .contains("gamma")); + // non-finite trat + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[0], + &[1], + &[0.0], + 15, + Some((f64::NAN, 0.0)), + None + )) + .contains("trat")); + // trat deviation checked only when deviations supplied + assert!(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[0], + &[1], + &[0.0], + 15, + Some((1500.0, f64::NAN)), + None + ) + .is_ok()); + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + Some(&[50.0, 50.0]), + &[0, 0], + &[0], + &[1], + &[0.0], + 15, + Some((1500.0, f64::NAN)), + None + )) + .contains("trat deviation")); + // non-finite thresh + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[0], + &[1], + &[0.0], + 15, + None, + Some(f64::NAN) + )) + .contains("thresh")); + // out-of-range index (only -1 sentinel allowed) + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[-2], + &[1], + &[0.0], + 15, + None, + None + )) + .contains("out of range")); + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[0], + &[2], + &[0.0], + 15, + None, + None + )) + .contains("out of range")); + // self-play + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + None, + &[0, 0], + &[1], + &[1], + &[0.0], + 15, + None, + None + )) + .contains("self-play")); + // deviations length mismatch + assert!(r(predict_rating_two( + &[1500.0, 1500.0], + Some(&[50.0]), + &[0, 0], + &[0], + &[1], + &[0.0], + 15, + None, + None + )) + .contains("deviations length")); + // multi: seat bounds, players length, nr=0, index range, trat + assert!(r(predict_rating_multi( + &[1500.0, 1500.0], + &[0, 0], + &[0], + 1, + 1, + 15, + None, + false + )) + .contains("2..=1000")); + assert!(r(predict_rating_multi( + &[1500.0, 1500.0], + &[0, 0], + &[0, 1, 0], + 1, + 2, + 15, + None, + false + )) + .contains("players length")); + // nr*np usize overflow must be a checked error, not a wrap + panic + // (kills: replacing checked_mul with wrapping `nr * np`, which lets an + // empty players slice pass the length check and then index OOB). + assert!(r(predict_rating_multi( + &[1500.0, 1500.0], + &[0, 0], + &[], + 1usize << 61, + 8, + 15, + None, + false + )) + .contains("overflows")); + assert!(r(predict_rating_multi( + &[1500.0, 1500.0], + &[0, 0], + &[], + 0, + 2, + 15, + None, + false + )) + .contains("at least one event")); + assert!(r(predict_rating_multi( + &[1500.0, 1500.0], + &[0, 0], + &[0, 5], + 1, + 2, + 15, + None, + false + )) + .contains("out of range")); + assert!(r(predict_rating_multi( + &[1500.0, 1500.0], + &[0, 0], + &[0, 1], + 1, + 2, + 15, + Some(f64::INFINITY), + false + )) + .contains("trat")); +} + +/// MC-500: structural invariants over random inputs. Every assert reads +/// crate return values: complement symmetry pred(w,b)+pred(b,w) ~= 1 at +/// gamma=0 (both branches), preds in (0,1), thresh output consistent with +/// the crate's own unthresholded preds, placing ranks are a valid +/// min-tie ranking of the crate's own EloM preds. +#[test] +#[ignore] +fn pr_mc_500_invariants() { + let mut rng = Lcg(0x5eed_cafe_1234_0001); + for rep in 0..500 { + let n = 3 + (rng.next_f64() * 8.0) as usize; + let ratings: Vec = (0..n).map(|_| 1200.0 + 1600.0 * rng.next_f64()).collect(); + let devs: Vec = (0..n).map(|_| 30.0 + 300.0 * rng.next_f64()).collect(); + let games: Vec = (0..n) + .map(|_| 15 + (rng.next_f64() * 40.0) as u64) + .collect(); + let w = (rng.next_f64() * n as f64) as i64; + let mut b = (rng.next_f64() * n as f64) as i64; + if b == w { + b = (b + 1) % n as i64; + } + for dev_opt in [None, Some(devs.as_slice())] { + let fwd = predict_rating_two( + &ratings, + dev_opt, + &games, + &[w], + &[b], + &[0.0], + 15, + None, + None, + ) + .unwrap(); + let bwd = predict_rating_two( + &ratings, + dev_opt, + &games, + &[b], + &[w], + &[0.0], + 15, + None, + None, + ) + .unwrap(); + assert!(fwd[0] > 0.0 && fwd[0] < 1.0, "rep {rep}"); + assert!( + (fwd[0] + bwd[0] - 1.0).abs() < 1e-12, + "rep {rep}: complement symmetry" + ); + let th = predict_rating_two( + &ratings, + dev_opt, + &games, + &[w], + &[b], + &[0.0], + 15, + None, + Some(0.5), + ) + .unwrap(); + let expect = if fwd[0] >= 0.5 { 1.0 } else { 0.0 }; + assert_eq!(th[0], expect, "rep {rep}: thresh vs crate pred"); + } + // EloM: placing must be the min-tie ranking of the crate's preds. + let np = 3.min(n); + let seats: Vec = (0..np as i64).collect(); + let preds = predict_rating_multi(&ratings, &games, &seats, 1, np, 15, None, false).unwrap(); + let ranks = predict_rating_multi(&ratings, &games, &seats, 1, np, 15, None, true).unwrap(); + for s in 0..np { + let expected = 1.0 + + preds + .iter() + .filter(|u| !u.is_nan() && **u > preds[s]) + .count() as f64; + assert_eq!(ranks[s], expected, "rep {rep} seat {s}"); + } + // Row mean of preds is 0 (crate values; na.rm mean identity). + let sum: f64 = preds.iter().sum(); + assert!(sum.abs() < 1e-9, "rep {rep}: pred row sum {sum}"); + } +} + +// --------------------------------------------------------------------------- +// bratt_mm (VGAM bratt: Bradley-Terry with additive-alpha0 ties) tests. +// Oracle: files/bratt_oracle.py (exact Fractions, EXECUTED). Every assert +// reads crate outputs (BrattResult fields or returned Err strings). +// --------------------------------------------------------------------------- + +fn bt2_fixture() -> (Vec, Vec) { + // B1: n = 3. wins y_ij and symmetric ties t_ij. + let y = vec![0.0, 3.0, 1.0, 1.0, 0.0, 2.0, 2.0, 1.0, 0.0]; + let t = vec![0.0, 1.0, 1.0, 1.0, 0.0, 2.0, 1.0, 2.0, 0.0]; + (y, t) +} + +fn bt2_rel(a: f64, b: f64) -> f64 { + (a - b).abs() / b.abs().max(1e-300) +} + +#[test] +fn bt2_anchor_b1_exact_iter1() { + // Oracle B1 iteration 1 (exact Fractions): alpha = [1, 27/40, 3/4], + // alpha0 = 9/14. tol = 0.5 > delta_1 = 0.357... so convergence fires at + // the first update and the crate returns the iteration-1 parameters. + // Kills MU1 (W_i includes ties), MU2 (alpha0 denominator double-counts + // ordered pairs), MU3 (rescale skips alpha0). + let (y, t) = bt2_fixture(); + let r = bratt_mm(&y, &t, 3, 0, 1.0, 100, 0.5).unwrap(); + assert_eq!(r.iterations, 1); + assert!(bt2_rel(r.alpha[0], 1.0) < 1e-15, "alpha0={}", r.alpha[0]); + assert!( + bt2_rel(r.alpha[1], 27.0 / 40.0) < 1e-15, + "alpha1={}", + r.alpha[1] + ); + assert!( + bt2_rel(r.alpha[2], 3.0 / 4.0) < 1e-15, + "alpha2={}", + r.alpha[2] + ); + assert!(bt2_rel(r.alpha0, 9.0 / 14.0) < 1e-15, "a0={}", r.alpha0); +} + +#[test] +fn bt2_anchor_b1_iter2() { + // Oracle B1 iteration 2 (exact Fractions -> float): alpha = + // [1, 0.6276558170061743, 0.7021719314027881], alpha0 = + // 0.6129259568116415. delta_2 = 0.0478 <= tol = 0.1 < delta_1, so + // convergence fires at update 2. Kills MU4 (D_ij missing alpha0): + // iteration 1 is blind to MU4 (constant D shift cancels in the + // rescale from the all-one start) but iteration 2 discriminates it. + let (y, t) = bt2_fixture(); + let r = bratt_mm(&y, &t, 3, 0, 1.0, 100, 0.1).unwrap(); + assert_eq!(r.iterations, 2); + assert!(bt2_rel(r.alpha[1], 0.6276558170061743) < 1e-13); + assert!(bt2_rel(r.alpha[2], 0.7021719314027881) < 1e-13); + assert!(bt2_rel(r.alpha0, 0.6129259568116415) < 1e-13); +} + +#[test] +fn bt2_converged_b2() { + // Oracle B2: tol = 1e-13 converges (oracle: 22 updates) to + // alpha = [1, 0.6150318884241122, 0.686344995662376], + // alpha0 = 0.6038293879270596, LL = -15.12765635227613 with + // stationarity gradient < 1e-13. Kills MU5 (convergence ignores the + // alpha0 delta) together with the iteration-count window: alpha0 + // still moves by more than alpha near the fixed point on this + // fixture, so dropping it from the delta stops too early. + let (y, t) = bt2_fixture(); + let r = bratt_mm(&y, &t, 3, 0, 1.0, 1000, 1e-13).unwrap(); + assert!((r.alpha[0] - 1.0).abs() < 1e-15); + assert!(bt2_rel(r.alpha[1], 0.6150318884241122) < 1e-10); + assert!(bt2_rel(r.alpha[2], 0.686344995662376) < 1e-10); + assert!(bt2_rel(r.alpha0, 0.6038293879270596) < 1e-10); + assert!(bt2_rel(r.log_likelihood, -15.12765635227613) < 1e-12); + assert!( + (18..=26).contains(&r.iterations), + "iterations = {}", + r.iterations + ); + // Spec gradient (independent formula) evaluated at the crate's + // returned parameters must vanish: W_i/a_i = sum_j n_ij/D_ij. + for i in 0..3 { + let w_i: f64 = (0..3).filter(|&j| j != i).map(|j| y[i * 3 + j]).sum(); + let s: f64 = (0..3) + .filter(|&j| j != i) + .map(|j| { + (y[i * 3 + j] + y[j * 3 + i] + t[i * 3 + j]) / (r.alpha[i] + r.alpha[j] + r.alpha0) + }) + .sum(); + assert!((w_i / r.alpha[i] - s).abs() < 1e-9, "grad[{i}]"); + } +} + +#[test] +fn bt2_permutation_b3() { + // Relabeling contestants by the permutation 0<->2 permutes alpha and + // preserves alpha0 and the log-likelihood (both fits read from crate). + let (y, t) = bt2_fixture(); + let perm = [2usize, 1, 0]; + let mut yp = vec![0.0; 9]; + let mut tp = vec![0.0; 9]; + for i in 0..3 { + for j in 0..3 { + yp[perm[i] * 3 + perm[j]] = y[i * 3 + j]; + tp[perm[i] * 3 + perm[j]] = t[i * 3 + j]; + } + } + let r = bratt_mm(&y, &t, 3, 0, 1.0, 1000, 1e-13).unwrap(); + let rp = bratt_mm(&yp, &tp, 3, perm[0], 1.0, 1000, 1e-13).unwrap(); + for i in 0..3 { + assert!(bt2_rel(rp.alpha[perm[i]], r.alpha[i]) < 1e-12, "alpha[{i}]"); + } + assert!(bt2_rel(rp.alpha0, r.alpha0) < 1e-12); + assert!(bt2_rel(rp.log_likelihood, r.log_likelihood) < 1e-12); +} + +#[test] +fn bt2_rescale_b4() { + // Oracle B4: refitting with ref_index = 1, ref_value = 2 returns the + // B2 solution jointly rescaled by c = 2/alpha_B2[1]; the LL is + // invariant under the joint rescale (crate-vs-crate cross pin). + let (y, t) = bt2_fixture(); + let r = bratt_mm(&y, &t, 3, 0, 1.0, 1000, 1e-13).unwrap(); + let r2 = bratt_mm(&y, &t, 3, 1, 2.0, 1000, 1e-13).unwrap(); + assert!(bt2_rel(r2.alpha[0], 3.2518639076171687) < 1e-10); + assert!((r2.alpha[1] - 2.0).abs() < 1e-14); + assert!(bt2_rel(r2.alpha[2], 2.2319005195681427) < 1e-10); + assert!(bt2_rel(r2.alpha0, 1.9635709929585714) < 1e-10); + assert!(bt2_rel(r2.log_likelihood, r.log_likelihood) < 1e-12); + let c = 2.0 / r.alpha[1]; + for i in 0..3 { + assert!(bt2_rel(r2.alpha[i], r.alpha[i] * c) < 1e-12, "alpha[{i}]"); + } + assert!(bt2_rel(r2.alpha0, r.alpha0 * c) < 1e-12); +} + +#[test] +fn bt2_error_contract() { + let (y, t) = bt2_fixture(); + let e = bratt_mm(&y, &t, 1, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("at least 2"), "{e}"); + // n cap is validated before any n*n length arithmetic. + let e = bratt_mm(&[], &[], 10001, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("cap of 10000"), "{e}"); + let e = bratt_mm(&y[..8], &t, 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("wins must be"), "{e}"); + let e = bratt_mm(&y, &t[..8], 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("ties must be"), "{e}"); + let mut bad = y.clone(); + bad[1] = f64::NAN; + let e = bratt_mm(&bad, &t, 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("win counts must be finite"), "{e}"); + let mut bad = t.clone(); + bad[1] = -1.0; + bad[3] = -1.0; + let e = bratt_mm(&y, &bad, 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("tie counts must be nonnegative"), "{e}"); + let mut bad = y.clone(); + bad[4] = 1.0; + let e = bratt_mm(&bad, &t, 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("diagonal of the wins"), "{e}"); + let mut bad = t.clone(); + bad[1] = 5.0; + let e = bratt_mm(&y, &bad, 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("symmetric"), "{e}"); + let e = bratt_mm(&y, &t, 3, 3, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("ref_index"), "{e}"); + let e = bratt_mm(&y, &t, 3, 0, 0.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("ref_value"), "{e}"); + let e = bratt_mm(&y, &t, 3, 0, 1.0, 10, 0.0).unwrap_err(); + assert!(e.contains("tol"), "{e}"); + let e = bratt_mm(&y, &t, 3, 0, 1.0, 0, 1e-6).unwrap_err(); + assert!(e.contains("max_iter"), "{e}"); + // W_i == 0: contestant 1 loses everything. + let y0 = vec![0.0, 3.0, 1.0, 0.0, 0.0, 0.0, 2.0, 1.0, 0.0]; + let e = bratt_mm(&y0, &t, 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("contestant 1 has no wins"), "{e}"); + // T == 0: no ties anywhere. + let t0 = vec![0.0; 9]; + let e = bratt_mm(&y, &t0, 3, 0, 1.0, 10, 1e-6).unwrap_err(); + assert!(e.contains("use bradley_terry_mm"), "{e}"); + // Non-convergence. + let e = bratt_mm(&y, &t, 3, 0, 1.0, 1, 1e-15).unwrap_err(); + assert!(e.contains("did not converge"), "{e}"); +} + +#[test] +#[ignore] +fn bt2_mc_500_dominance() { + // 500 random 4-player tournaments. The crate's fitted (alpha, alpha0) + // must dominate 20 random non-scale perturbations under the spec LL + // formula computed IN THE TEST from crate outputs (a wrong crate + // optimum fails the dominance; a wrong log_likelihood field fails the + // recomputation cross-check). Fractional weighted counts included. + let n = 4usize; + let spec_ll = |y: &[f64], t: &[f64], a: &[f64], a0: f64| -> f64 { + let mut ll = 0.0; + for i in 0..n { + for j in 0..n { + if i != j && y[i * n + j] > 0.0 { + ll += y[i * n + j] * a[i].ln(); + } + } + } + for i in 0..n { + for j in (i + 1)..n { + if t[i * n + j] > 0.0 { + ll += t[i * n + j] * a0.ln(); + } + let n_ij = y[i * n + j] + y[j * n + i] + t[i * n + j]; + if n_ij > 0.0 { + ll -= n_ij * (a[i] + a[j] + a0).ln(); + } + } + } + ll + }; + let mut rng = Lcg(0x5eed_b2a7_7001_u64); + for rep in 0..500 { + let mut y = vec![0.0f64; n * n]; + let mut t = vec![0.0f64; n * n]; + for i in 0..n { + for j in 0..n { + if i != j { + y[i * n + j] = 1.0 + (rng.next_f64() * 9.0).floor() + 0.5 * rng.next_f64(); + } + } + } + for i in 0..n { + for j in (i + 1)..n { + let v = 1.0 + (rng.next_f64() * 5.0).floor(); + t[i * n + j] = v; + t[j * n + i] = v; + } + } + let r = bratt_mm(&y, &t, n, 0, 1.0, 20000, 1e-12).unwrap(); + let base = spec_ll(&y, &t, &r.alpha, r.alpha0); + assert!( + bt2_rel(base, r.log_likelihood) < 1e-10, + "rep {rep}: LL field vs spec recomputation" + ); + for _ in 0..20 { + let mut ap = r.alpha.clone(); + for a in ap.iter_mut().skip(1) { + *a *= (0.1 * rng.normal()).exp(); + } + let a0p = r.alpha0 * (0.1 * rng.normal()).exp(); + assert!( + spec_ll(&y, &t, &ap, a0p) <= base + 1e-9, + "rep {rep}: dominance" + ); + } + } +} + +#[test] +fn bt2_alpha0_delta_convergence() { + // MU5 killer (convergence check dropping the alpha0 delta). On the B1 + // fixture the update-1 deltas are max|alpha' - alpha| = 0.325 and + // |alpha0' - alpha0| = 5/14 = 0.357 (oracle per-iteration trace), so + // with tol = 0.34 the correct check (max over alpha AND alpha0) does + // NOT fire at update 1 and fires at update 2 (delta_2 = 0.0478), + // returning the iteration-2 parameters; a check that ignores alpha0 + // fires at update 1 and returns alpha[1] = 27/40 = 0.675 instead of + // 0.6276558170061743. Both asserts read crate outputs. + let (y, t) = bt2_fixture(); + let r = bratt_mm(&y, &t, 3, 0, 1.0, 100, 0.34).unwrap(); + assert_eq!(r.iterations, 2); + assert!(bt2_rel(r.alpha[1], 0.6276558170061743) < 1e-13); + assert!(bt2_rel(r.alpha0, 0.6129259568116415) < 1e-13); +} + +#[test] +fn bt2_huge_counts_overflow_rejected() { + // Impl-review regression: entries are individually finite, but the + // pair total n_12 = 9e307 + 9e307 + 1 overflows to +inf. Before the + // guard this returned Ok with alpha = [1, 0, 0], alpha0 = 0 and a NaN + // log-likelihood. Assert reads the crate Err. + let h = 9e307; + let y = vec![0.0, 1.0, 1.0, 1.0, 0.0, h, 1.0, h, 0.0]; + let t = vec![0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0]; + let err = bratt_mm(&y, &t, 3, 0, 1.0, 100, 2.0).unwrap_err(); + assert!(err.contains("overflow"), "unexpected error: {err}"); + // Row-total overflow (w_tot) is also rejected, not just pair totals. + let y2 = vec![0.0, h, h, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0]; + let err2 = bratt_mm(&y2, &t, 3, 0, 1.0, 100, 2.0).unwrap_err(); + assert!(err2.contains("overflow"), "unexpected error: {err2}"); +}