helena 0.1.0

//! Pairwise human-judgment aggregation: the meta-evaluation core that turns
//! blocked pairwise comparisons into a system ranking.
//!
//! The small-sample regime has no stable absolute scale, so the trustworthy
//! human signal is *relative*: raters compare two systems on a shared
//! conditioning item and say which they prefer, ties allowed. [`BradleyTerry`]
//! fits a latent strength per system by maximum likelihood, the system-level
//! ranking the metric-vs-human meta-evaluation (a later slice) validates an
//! automated metric against. [`Reliability`] reports Krippendorff alpha over the
//! same records, because a confident ranking over judgments the raters disagreed
//! on is false precision.
//!
//! Total, deterministic, dependency-free host-side folds in the
//! [`sample_dispersion`](super::sample_dispersion) mould (no Burn, no FFT, no new
//! dependency). Unit tests use synthetic records; the feature-gated CLI ingests named
//! judgment files and maps them onto this dense-index core.

use std::collections::BTreeMap;

use crate::{Error, Result};

/// The outcome of one pairwise comparison between two systems `a` and `b`.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Outcome {
    /// System `a` was preferred.
    AWins,
    /// System `b` was preferred.
    BWins,
    /// Neither was preferred: a tie, which the two-question protocol allows.
    Tie,
}

/// One pairwise judgment: on conditioning `item`, rater `rater` compared systems
/// `a` and `b` and returned [`outcome`](Self::outcome).
///
/// Systems, items, and raters are caller-assigned dense indices (the CLI maps
/// names to indices when it ingests a judgments file). The one record-level
/// invariant is `a != b` (a system is not compared with itself), enforced by
/// [`new`](Self::new); index bounds against the system count are a whole-dataset
/// property [`BradleyTerry::fit`] checks. `item` and `rater` are carried for the
/// reliability and item-blocked meta-evaluation reads and are not consulted by
/// the strength fit.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Judgment {
    item: usize,
    rater: usize,
    a: usize,
    b: usize,
    outcome: Outcome,
}

impl Judgment {
    /// A judgment comparing the distinct systems `a` and `b`. Returns
    /// [`Error::Validation`] if `a == b`.
    pub fn new(item: usize, rater: usize, a: usize, b: usize, outcome: Outcome) -> Result<Self> {
        if a == b {
            return Err(Error::validation(format!(
                "a judgment compares two distinct systems; got a == b == {a}"
            )));
        }
        Ok(Self {
            item,
            rater,
            a,
            b,
            outcome,
        })
    }

    /// The conditioning item the comparison was made on (the meta-evaluation block).
    pub fn item(&self) -> usize {
        self.item
    }

    /// The rater who made the comparison (the reliability unit).
    pub fn rater(&self) -> usize {
        self.rater
    }

    /// The first compared system.
    pub fn a(&self) -> usize {
        self.a
    }

    /// The second compared system.
    pub fn b(&self) -> usize {
        self.b
    }

    /// Which system was preferred.
    pub fn outcome(&self) -> Outcome {
        self.outcome
    }
}

/// Bradley-Terry strengths fit from pairwise judgments by maximum likelihood.
///
/// The Bradley-Terry model gives system `i` a positive strength `p_i` with
/// `P(i beats j) = p_i / (p_i + p_j)`; [`fit`](Self::fit) recovers the strengths
/// from the observed win counts by the MM (minorization-maximization) iteration
/// (Hunter, 2004), the system-level ranking the meta-evaluation validates a
/// metric against. Strengths are identified only up to a common scale, so they
/// are normalized to sum to one and read as relative shares.
///
/// **Ties** split as half a win to each side (the standard reduction when the
/// tie rate is not itself modeled); a Davidson-style explicit tie term is a
/// documented refinement.
///
/// **Identifiability is required, not patched.** The MLE is finite and unique
/// iff the win graph (an edge `i -> j` whenever `i` took a positive win share
/// from `j`) is strongly connected (Ford, 1957): otherwise some system never
/// loses (its strength diverges) or the comparisons split into groups with no
/// cross-result to scale against. [`fit`](Self::fit) checks this and returns
/// [`Error::Validation`] rather than report a regularized strength whose
/// magnitude is set by a prior, not the data, the false precision the
/// small-sample evaluation brief warns against.
#[derive(Clone, Debug, PartialEq)]
pub struct BradleyTerry {
    strengths: Vec<f64>,
}

impl BradleyTerry {
    /// Fit strengths for `systems` systems from `judgments`.
    ///
    /// Returns [`Error::Validation`] if `systems < 2`, if `judgments` is empty,
    /// if any judgment references a system index `>= systems`, if the win graph
    /// is not strongly connected (see the type docs), or if the iteration fails
    /// to converge.
    pub fn fit(systems: usize, judgments: &[Judgment]) -> Result<Self> {
        if systems < 2 {
            return Err(Error::validation(format!(
                "a Bradley-Terry fit needs at least two systems; got {systems}"
            )));
        }
        if judgments.is_empty() {
            return Err(Error::validation("no judgments to fit"));
        }

        // wins[i][j] is i's (possibly fractional, with ties) win share over j.
        let mut wins = vec![vec![0.0f64; systems]; systems];
        for (k, judgment) in judgments.iter().enumerate() {
            let (a, b) = (judgment.a, judgment.b);
            if a >= systems || b >= systems {
                return Err(Error::validation(format!(
                    "judgment {k} references system {a} / {b} but there are only {systems} systems"
                )));
            }
            match judgment.outcome {
                Outcome::AWins => wins[a][b] += 1.0,
                Outcome::BWins => wins[b][a] += 1.0,
                Outcome::Tie => {
                    wins[a][b] += 0.5;
                    wins[b][a] += 0.5;
                }
            }
        }

        if !is_strongly_connected(&wins, systems) {
            return Err(Error::validation(
                "the comparison design does not identify a ranking: the win graph is not \
                 strongly connected (a system never loses or the comparisons are disconnected)",
            ));
        }
        // Strong connectivity gives every system at least one win and one loss,
        // so no numerator or denominator in the MM step below is zero.

        let wins_total: Vec<f64> = wins.iter().map(|row| row.iter().sum()).collect();
        let counts = pair_counts(&wins, systems);
        let strengths = mm_fit(&wins_total, &counts, systems)?;
        Ok(Self { strengths })
    }

    /// The fitted strengths in system-index order, each positive and summing to one.
    pub fn strengths(&self) -> &[f64] {
        &self.strengths
    }

    /// System indices ordered best to worst by strength, ties broken by index.
    pub fn ranking(&self) -> Vec<usize> {
        let mut order: Vec<usize> = (0..self.strengths.len()).collect();
        order.sort_by(|&i, &j| {
            self.strengths[j]
                .total_cmp(&self.strengths[i])
                .then(i.cmp(&j))
        });
        order
    }
}

/// Symmetric per-pair comparison counts `n[i][j] = wins[i][j] + wins[j][i]`, the
/// total games played between each pair however they were split.
fn pair_counts(wins: &[Vec<f64>], n: usize) -> Vec<Vec<f64>> {
    let mut counts = vec![vec![0.0f64; n]; n];
    for (i, row) in counts.iter_mut().enumerate() {
        for (j, c) in row.iter_mut().enumerate() {
            *c = wins[i][j] + wins[j][i];
        }
    }
    counts
}

/// Strong connectivity of the win digraph (edge `i -> j` iff `wins[i][j] > 0`):
/// node 0 reaches every node and every node reaches node 0, which together imply
/// every pair is mutually reachable. Ford's condition for a unique finite MLE.
fn is_strongly_connected(wins: &[Vec<f64>], n: usize) -> bool {
    reach(n, |i, j| wins[i][j] > 0.0) == n && reach(n, |i, j| wins[j][i] > 0.0) == n
}

/// Count of nodes reachable from node 0 (inclusive) by following directed edges
/// `edge(from, to)`. With the reversed predicate it counts nodes that reach 0.
fn reach(n: usize, edge: impl Fn(usize, usize) -> bool) -> usize {
    let mut seen = vec![false; n];
    seen[0] = true;
    let mut stack = vec![0usize];
    let mut count = 1;
    while let Some(i) = stack.pop() {
        for (j, seen_j) in seen.iter_mut().enumerate() {
            if !*seen_j && edge(i, j) {
                *seen_j = true;
                count += 1;
                stack.push(j);
            }
        }
    }
    count
}

/// The MM (minorization-maximization) iteration for Bradley-Terry strengths
/// (Hunter, 2004): `p_i <- W_i / sum_j n_ij / (p_i + p_j)`, renormalized to sum
/// to one each step. The update is scale-equivariant, so the per-step
/// renormalization only fixes the free scale; the ratios converge to the MLE.
/// Geometric convergence under a strongly-connected design.
fn mm_fit(wins_total: &[f64], counts: &[Vec<f64>], n: usize) -> Result<Vec<f64>> {
    const TOL: f64 = 1e-10;
    const MAX_ITERS: usize = 10_000;

    let mut p = vec![1.0 / n as f64; n];
    for _ in 0..MAX_ITERS {
        let mut next = vec![0.0f64; n];
        for (i, slot) in next.iter_mut().enumerate() {
            let mut denom = 0.0f64;
            for j in 0..n {
                if i != j && counts[i][j] > 0.0 {
                    denom += counts[i][j] / (p[i] + p[j]);
                }
            }
            *slot = wins_total[i] / denom;
        }
        let sum: f64 = next.iter().sum();
        for x in &mut next {
            *x /= sum;
        }
        let delta = (0..n).map(|i| (next[i] - p[i]).abs()).fold(0.0, f64::max);
        p = next;
        if delta < TOL {
            return Ok(p);
        }
    }
    Err(Error::validation(
        "Bradley-Terry iteration did not converge",
    ))
}

/// Inter-rater reliability of pairwise judgments by Krippendorff's alpha (nominal).
///
/// A Bradley-Terry ranking is only as trustworthy as the judgments behind it, so
/// the small-sample brief pairs any ranking with this reliability read; a
/// confident order over judgments the raters disagreed on is not a real result.
/// The *unit* of agreement is one comparison question, a (conditioning item,
/// unordered system pair); raters who judged the same unit are compared on the
/// canonical outcome (which system of the pair won, or a tie), so the two
/// presentation orders of one comparison count as the same rated value.
///
/// Alpha is `1` for perfect agreement, `0` at chance, and negative for systematic
/// disagreement. Units judged by fewer than two raters carry no agreement
/// information and are excluded. When every usable rating falls in one category
/// there is no variability to disagree over and alpha is `1` by convention (the
/// raters did agree). Nominal alpha treats the three outcomes as unordered; an
/// ordinal metric that counts a tie-vs-win disagreement as smaller than a
/// win-vs-loss one is a documented refinement.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Reliability {
    alpha: f64,
    units: usize,
    ratings: usize,
}

impl Reliability {
    /// Krippendorff's nominal alpha over `judgments`.
    ///
    /// Returns [`Error::Validation`] if `judgments` is empty, if any rater judged
    /// the same comparison unit more than once, or if no unit was judged by at
    /// least two raters (agreement is then undefined).
    pub fn krippendorff_alpha(judgments: &[Judgment]) -> Result<Self> {
        if judgments.is_empty() {
            return Err(Error::validation("no judgments to assess"));
        }

        // unit (item, lo, hi) -> rater -> category (0 lo-wins, 1 hi-wins, 2 tie).
        // BTreeMap keeps a deterministic accumulation order, so the floating-point
        // sums below are reproducible run to run.
        let mut units: BTreeMap<(usize, usize, usize), BTreeMap<usize, usize>> = BTreeMap::new();
        for judgment in judgments {
            let (lo, hi, category) = canonical_category(judgment);
            let raters = units.entry((judgment.item, lo, hi)).or_default();
            if raters.insert(judgment.rater, category).is_some() {
                return Err(Error::validation(format!(
                    "rater {} judged the comparison (item {}, systems {lo}/{hi}) more than once",
                    judgment.rater, judgment.item
                )));
            }
        }

        // Nominal alpha needs only the coincidence diagonal and the value
        // marginals, not the full matrix.
        let mut n = 0.0f64;
        let mut diagonal = 0.0f64;
        let mut marginal = [0.0f64; 3];
        let mut usable_units = 0;
        let mut ratings = 0;
        for raters in units.values() {
            let m = raters.len();
            if m < 2 {
                continue;
            }
            usable_units += 1;
            ratings += m;
            let mut counts = [0usize; 3];
            for &category in raters.values() {
                counts[category] += 1;
            }
            let pairable = (m - 1) as f64;
            for (category, &count) in counts.iter().enumerate() {
                let nc = count as f64;
                n += nc;
                marginal[category] += nc;
                diagonal += nc * (nc - 1.0) / pairable;
            }
        }

        if usable_units == 0 {
            return Err(Error::validation(
                "reliability needs at least one comparison judged by two or more raters",
            ));
        }

        let expected = n * n - marginal.iter().map(|&x| x * x).sum::<f64>();
        // No expected disagreement: every usable rating shares one category, so
        // the raters agreed perfectly and alpha is 1 by convention.
        let alpha = if expected == 0.0 {
            1.0
        } else {
            1.0 - (n - 1.0) * (n - diagonal) / expected
        };
        Ok(Self {
            alpha,
            units: usable_units,
            ratings,
        })
    }

    /// Krippendorff's alpha: `1` perfect agreement, `0` chance, negative worse.
    pub fn alpha(&self) -> f64 {
        self.alpha
    }

    /// Number of comparison units two or more raters judged; the rest carry no
    /// agreement information and are excluded.
    pub fn units(&self) -> usize {
        self.units
    }

    /// Total ratings over the usable units.
    pub fn ratings(&self) -> usize {
        self.ratings
    }
}

/// The canonical `(lo, hi, category)` of a judgment: the system pair in index
/// order and which of them won (`0` lo, `1` hi) or a tie (`2`), so the two
/// presentation orders of one comparison map to the same rated value.
fn canonical_category(judgment: &Judgment) -> (usize, usize, usize) {
    let (lo, hi) = if judgment.a < judgment.b {
        (judgment.a, judgment.b)
    } else {
        (judgment.b, judgment.a)
    };
    let category = match judgment.outcome {
        Outcome::Tie => 2,
        Outcome::AWins => usize::from(judgment.a != lo),
        Outcome::BWins => usize::from(judgment.b != lo),
    };
    (lo, hi, category)
}

#[cfg(test)]
mod tests {
    use super::*;

    /// Build judgments from `(a, b, outcome)` triples; item and rater are
    /// irrelevant to the strength fit, so they stay zero here.
    fn judgments(triples: &[(usize, usize, Outcome)]) -> Vec<Judgment> {
        triples
            .iter()
            .map(|&(a, b, o)| Judgment::new(0, 0, a, b, o).unwrap())
            .collect()
    }

    /// A round-robin where, for every pair `i < j`, the lower index wins `win`
    /// of `win + lose` games. With `win > lose > 0` this both ranks `0 > 1 > ...`
    /// and makes the win graph strongly connected (each pair has results both
    /// ways), the canonical identifiable design.
    fn graded_round_robin(systems: usize, win: usize, lose: usize) -> Vec<Judgment> {
        let mut triples = Vec::new();
        for i in 0..systems {
            for j in (i + 1)..systems {
                for _ in 0..win {
                    triples.push((i, j, Outcome::AWins));
                }
                for _ in 0..lose {
                    triples.push((i, j, Outcome::BWins));
                }
            }
        }
        judgments(&triples)
    }

    #[test]
    fn recovers_known_ranking() {
        let bt = BradleyTerry::fit(3, &graded_round_robin(3, 2, 1)).unwrap();
        assert_eq!(bt.ranking(), vec![0, 1, 2]);
        let s = bt.strengths();
        assert!(s[0] > s[1] && s[1] > s[2], "{s:?}");
    }

    #[test]
    fn strengths_are_a_positive_simplex() {
        let bt = BradleyTerry::fit(4, &graded_round_robin(4, 3, 1)).unwrap();
        let s = bt.strengths();
        assert_eq!(s.len(), 4);
        assert!(s.iter().all(|&p| p > 0.0), "{s:?}");
        assert!((s.iter().sum::<f64>() - 1.0).abs() < 1e-9, "{s:?}");
    }

    #[test]
    fn balanced_results_give_equal_strengths() {
        // Every pair splits one win each way: perfectly symmetric, so all
        // strengths equal 1/n.
        let bt = BradleyTerry::fit(3, &graded_round_robin(3, 1, 1)).unwrap();
        for &p in bt.strengths() {
            assert!((p - 1.0 / 3.0).abs() < 1e-9, "{:?}", bt.strengths());
        }
    }

    #[test]
    fn ties_are_symmetric_and_keep_the_graph_connected() {
        // Every pair only ever ties: each side takes half a win, so the design
        // is identifiable (edges both ways) and all strengths are equal.
        let all_ties = judgments(&[
            (0, 1, Outcome::Tie),
            (0, 2, Outcome::Tie),
            (1, 2, Outcome::Tie),
        ]);
        let bt = BradleyTerry::fit(3, &all_ties).unwrap();
        for &p in bt.strengths() {
            assert!((p - 1.0 / 3.0).abs() < 1e-9, "{:?}", bt.strengths());
        }
    }

    #[test]
    fn one_tie_breaks_a_near_dominance() {
        // System 0 beats 1 outright but they also tie once, so 0 does not
        // diverge; 0 still ranks first.
        let js = judgments(&[
            (0, 1, Outcome::AWins),
            (0, 1, Outcome::AWins),
            (0, 1, Outcome::Tie),
        ]);
        let bt = BradleyTerry::fit(2, &js).unwrap();
        assert_eq!(bt.ranking(), vec![0, 1]);
        assert!(bt.strengths()[0] > bt.strengths()[1]);
    }

    #[test]
    fn is_deterministic() {
        let js = graded_round_robin(5, 3, 2);
        let a = BradleyTerry::fit(5, &js).unwrap();
        let b = BradleyTerry::fit(5, &js).unwrap();
        assert_eq!(a, b);
    }

    #[test]
    fn ranking_breaks_strength_ties_by_index() {
        let bt = BradleyTerry::fit(3, &graded_round_robin(3, 1, 1)).unwrap();
        // Equal strengths, so the order is the index order.
        assert_eq!(bt.ranking(), vec![0, 1, 2]);
    }

    #[test]
    fn rejects_self_comparison_at_construction() {
        assert!(Judgment::new(0, 0, 2, 2, Outcome::AWins).is_err());
    }

    #[test]
    fn rejects_too_few_systems_and_empty() {
        assert!(BradleyTerry::fit(1, &graded_round_robin(2, 2, 1)).is_err());
        assert!(BradleyTerry::fit(2, &[]).is_err());
    }

    #[test]
    fn rejects_out_of_range_system() {
        let js = judgments(&[(0, 3, Outcome::AWins)]);
        assert!(BradleyTerry::fit(2, &js).is_err());
    }

    #[test]
    fn rejects_a_system_that_never_loses() {
        // 0 always beats 1: 0's strength diverges, so the design is rejected
        // rather than reporting an arbitrary finite number.
        let js = judgments(&[(0, 1, Outcome::AWins), (0, 1, Outcome::AWins)]);
        assert!(BradleyTerry::fit(2, &js).is_err());
    }

    #[test]
    fn rejects_disconnected_comparisons() {
        // {0, 1} and {2, 3} each compared internally but never across: the two
        // groups cannot be placed on one scale.
        let js = judgments(&[
            (0, 1, Outcome::AWins),
            (1, 0, Outcome::AWins),
            (2, 3, Outcome::AWins),
            (3, 2, Outcome::AWins),
        ]);
        assert!(BradleyTerry::fit(4, &js).is_err());
    }

    #[test]
    fn alpha_is_one_for_perfect_agreement() {
        // Two units, each judged by two agreeing raters (on different outcomes).
        let js = vec![
            Judgment::new(0, 0, 1, 2, Outcome::AWins).unwrap(),
            Judgment::new(0, 1, 1, 2, Outcome::AWins).unwrap(),
            Judgment::new(1, 0, 0, 3, Outcome::BWins).unwrap(),
            Judgment::new(1, 1, 0, 3, Outcome::BWins).unwrap(),
        ];
        let r = Reliability::krippendorff_alpha(&js).unwrap();
        assert!((r.alpha() - 1.0).abs() < 1e-12, "{}", r.alpha());
        assert_eq!(r.units(), 2);
        assert_eq!(r.ratings(), 4);
    }

    #[test]
    fn alpha_canonicalizes_presentation_order() {
        // The same preference (system 2 over system 5) presented both ways agrees.
        let js = vec![
            Judgment::new(0, 0, 2, 5, Outcome::AWins).unwrap(),
            Judgment::new(0, 1, 5, 2, Outcome::BWins).unwrap(),
        ];
        let r = Reliability::krippendorff_alpha(&js).unwrap();
        assert!((r.alpha() - 1.0).abs() < 1e-12, "{}", r.alpha());
        assert_eq!(r.units(), 1);
    }

    #[test]
    fn alpha_is_zero_at_chance() {
        // One unit split one-each: chance-level agreement.
        let js = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(0, 1, 0, 1, Outcome::BWins).unwrap(),
        ];
        let r = Reliability::krippendorff_alpha(&js).unwrap();
        assert!(r.alpha().abs() < 1e-12, "{}", r.alpha());
    }

    #[test]
    fn alpha_is_negative_for_systematic_disagreement() {
        // Two units, each split the same way: worse than chance.
        let js = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(0, 1, 0, 1, Outcome::BWins).unwrap(),
            Judgment::new(1, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(1, 1, 0, 1, Outcome::BWins).unwrap(),
        ];
        let r = Reliability::krippendorff_alpha(&js).unwrap();
        assert!(r.alpha() < 0.0, "{}", r.alpha());
    }

    #[test]
    fn alpha_rejects_duplicate_rater_on_a_unit() {
        // Same rater, same unit, twice (here via the flipped presentation order).
        let js = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(0, 0, 1, 0, Outcome::AWins).unwrap(),
        ];
        assert!(Reliability::krippendorff_alpha(&js).is_err());
    }

    #[test]
    fn alpha_rejects_no_multi_rater_unit() {
        // Two distinct units, one rater each: no agreement to measure.
        let js = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(1, 0, 0, 1, Outcome::AWins).unwrap(),
        ];
        assert!(Reliability::krippendorff_alpha(&js).is_err());
    }

    #[test]
    fn alpha_rejects_empty() {
        assert!(Reliability::krippendorff_alpha(&[]).is_err());
    }

    use proptest::prelude::*;

    /// A strongly-connected design: a balanced base round-robin (one win each
    /// way for every pair, guaranteeing identifiability) plus random extra
    /// decisive games and ties on top.
    fn arb_connected_design() -> impl Strategy<Value = (usize, Vec<Judgment>)> {
        (2usize..6)
            .prop_flat_map(|systems| {
                let pairs = systems * (systems - 1) / 2;
                (Just(systems), proptest::collection::vec(0u8..3, pairs * 3))
            })
            .prop_map(|(systems, extra)| {
                let mut triples = Vec::new();
                let mut e = extra.into_iter();
                for i in 0..systems {
                    for j in (i + 1)..systems {
                        triples.push((i, j, Outcome::AWins));
                        triples.push((j, i, Outcome::AWins));
                        for _ in 0..3 {
                            match e.next() {
                                Some(0) => triples.push((i, j, Outcome::AWins)),
                                Some(1) => triples.push((i, j, Outcome::BWins)),
                                _ => triples.push((i, j, Outcome::Tie)),
                            }
                        }
                    }
                }
                (systems, judgments(&triples))
            })
    }

    proptest! {
        #[test]
        fn fit_is_a_positive_simplex_permutation((systems, js) in arb_connected_design()) {
            let bt = BradleyTerry::fit(systems, &js).unwrap();
            let s = bt.strengths();
            prop_assert_eq!(s.len(), systems);
            prop_assert!(s.iter().all(|&p| p > 0.0 && p.is_finite()));
            prop_assert!((s.iter().sum::<f64>() - 1.0).abs() < 1e-9);

            let mut sorted = bt.ranking();
            sorted.sort_unstable();
            prop_assert_eq!(sorted, (0..systems).collect::<Vec<_>>());
        }

        #[test]
        fn fit_is_reproducible((systems, js) in arb_connected_design()) {
            let a = BradleyTerry::fit(systems, &js).unwrap();
            let b = BradleyTerry::fit(systems, &js).unwrap();
            prop_assert_eq!(a, b);
        }

        #[test]
        fn strength_order_matches_win_share((systems, js) in arb_connected_design()) {
            // The MLE ranks systems by fitted strength; the top-ranked system's
            // strength is at least the bottom-ranked one's (a sanity floor that
            // the fit is monotone in the latent scale).
            let bt = BradleyTerry::fit(systems, &js).unwrap();
            let order = bt.ranking();
            let s = bt.strengths();
            prop_assert!(s[order[0]] >= s[order[systems - 1]]);
        }
    }

    /// A grid of `units x raters` outcome categories (`0` a-wins, `1` b-wins,
    /// `2` tie), each unit a distinct item with at least two raters.
    fn arb_unit_grid() -> impl Strategy<Value = Vec<Vec<u8>>> {
        proptest::collection::vec(proptest::collection::vec(0u8..3, 2..5), 1..6)
    }

    /// Build judgments over the fixed pair `(0, 1)` from a category grid. With
    /// `agree`, every rater of a unit copies the first rater's category, forcing
    /// perfect within-unit agreement. Rater indices are unique per unit, so no
    /// unit is double-rated.
    fn from_grid(grid: &[Vec<u8>], agree: bool) -> Vec<Judgment> {
        let mut js = Vec::new();
        for (item, raters) in grid.iter().enumerate() {
            for (rater, &cat) in raters.iter().enumerate() {
                let category = if agree { raters[0] } else { cat };
                let outcome = match category {
                    0 => Outcome::AWins,
                    1 => Outcome::BWins,
                    _ => Outcome::Tie,
                };
                js.push(Judgment::new(item, rater, 0, 1, outcome).unwrap());
            }
        }
        js
    }

    proptest! {
        #[test]
        fn alpha_at_most_one(grid in arb_unit_grid()) {
            let r = Reliability::krippendorff_alpha(&from_grid(&grid, false)).unwrap();
            prop_assert!(r.alpha() <= 1.0 + 1e-9);
        }

        #[test]
        fn alpha_is_one_under_within_unit_agreement(grid in arb_unit_grid()) {
            let r = Reliability::krippendorff_alpha(&from_grid(&grid, true)).unwrap();
            prop_assert!((r.alpha() - 1.0).abs() < 1e-9, "{}", r.alpha());
        }

        #[test]
        fn alpha_is_reproducible(grid in arb_unit_grid()) {
            let js = from_grid(&grid, false);
            let a = Reliability::krippendorff_alpha(&js).unwrap();
            let b = Reliability::krippendorff_alpha(&js).unwrap();
            prop_assert_eq!(a, b);
        }
    }
}