helena 0.1.0

//! Metric-vs-human meta-evaluation: does an automated metric track human judgment?
//!
//! [`pairwise`](super::pairwise) turns human comparisons into a system ranking; this
//! module asks the complementary question the small-sample brief makes central, *is a
//! cheap automated metric a trustworthy stand-in for that ranking?* The answer is a
//! correlation, not a benchmark score: on the same conditioning items humans saw, does
//! the metric's per-item score *delta* between two systems move with the human
//! *preference margin*?
//!
//! [`rank_agreement`] reports the brief's secondary statistic, item-blocked Kendall's
//! tau-b and Spearman's rho between those matched margins and deltas. Blocking by
//! conditioning item is the load-bearing choice: correlating raw human scores against
//! raw metric scores across items would reward a metric that merely tracks item
//! difficulty, so the unit of correlation is a within-item `(margin, delta)` pair, which
//! differences item difficulty out. [`soft_pairwise_accuracy`] reports the primary
//! system-level read: do the metric's pairwise system preferences match the human
//! preferences once both are softened by paired permutation tests? The interval helpers
//! use a non-mean clustered bootstrap that resamples conditioning items and recomputes
//! the whole statistic for each draw.
//!
//! [`kendall_tau_b`] and [`spearman_rho`] are the underlying rank-correlation
//! primitives, exposed because they are well-specified on their own and reused by the
//! bootstrap. Total, deterministic, dependency-free host-side folds in the
//! [`pairwise`](super::pairwise) / [`ClusterBootstrap`](super::ClusterBootstrap) mould
//! (no Burn, no FFT, no new dependency).

use std::cmp::Ordering;
use std::collections::{BTreeMap, BTreeSet};
use std::num::NonZeroUsize;

use super::uncertainty::{ClusterBootstrap, ConfidenceInterval, SplitMix64};
use super::{Judgment, Outcome};
use crate::{Error, Result};

/// One automated-metric score: on conditioning `item`, the metric scored the output
/// from `system` at `value`.
///
/// `value` is oriented so that *higher means more preferred*; a lower-is-better metric
/// (a distance) is negated by the caller before it gets here, so a positive
/// [`rank_agreement`] always reads as "the metric tracks humans". Finiteness is a
/// whole-dataset invariant [`MetricScores::new`] checks, so a single record is a
/// transparent carrier with no constructor of its own.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct MetricScore {
    /// The conditioning item the output was generated for.
    pub item: usize,
    /// The system that produced the output.
    pub system: usize,
    /// The metric's scalar score, higher meaning more preferred.
    pub value: f64,
}

/// A metric's scores over `(item, system)` outputs, the automated signal
/// [`rank_agreement`] meta-evaluates against human judgment.
///
/// Constructed once and queried by [`get`](Self::get); the internal `BTreeMap` keys the
/// lookup on `(item, system)` and keeps iteration deterministic. The invariants
/// ([`new`](Self::new) enforces them) are that every `value` is finite and no
/// `(item, system)` is scored twice, so a metric cannot silently contribute a `NaN`
/// delta or a doubled observation.
#[derive(Clone, Debug, PartialEq)]
pub struct MetricScores {
    by_cell: BTreeMap<(usize, usize), f64>,
}

impl MetricScores {
    /// Collect `scores` into a queryable table.
    ///
    /// Returns [`Error::Validation`] if `scores` is empty, if any `value` is non-finite,
    /// or if two records share an `(item, system)` cell.
    pub fn new(scores: Vec<MetricScore>) -> Result<Self> {
        if scores.is_empty() {
            return Err(Error::validation("no metric scores"));
        }
        let mut by_cell = BTreeMap::new();
        for score in scores {
            if !score.value.is_finite() {
                return Err(Error::validation(format!(
                    "metric score for (item {}, system {}) is not finite: {}",
                    score.item, score.system, score.value
                )));
            }
            if by_cell
                .insert((score.item, score.system), score.value)
                .is_some()
            {
                return Err(Error::validation(format!(
                    "metric scores (item {}, system {}) twice",
                    score.item, score.system
                )));
            }
        }
        Ok(Self { by_cell })
    }

    /// The score for `system` on `item`, or `None` if that cell was not scored.
    pub fn get(&self, item: usize, system: usize) -> Option<f64> {
        self.by_cell.get(&(item, system)).copied()
    }
}

/// Deterministic paired sign-flip permutation test used by [`soft_pairwise_accuracy`].
///
/// Each input is an item-level paired difference, oriented so positive means the
/// lower-indexed system is better than the higher-indexed system. The one-tailed
/// p-value estimates `P(null_diff >= observed_diff)`: values below `0.5` mean the
/// first system is preferred, values above `0.5` mean the second is preferred, and
/// values near `0.5` mean the pair is not separated. Tied permutation draws get half
/// credit, so all-zero or exactly tied differences return `0.5` rather than a
/// misleading edge p-value.
#[derive(Clone, Copy, Debug)]
pub struct PairedPermutationTest {
    resamples: NonZeroUsize,
    seed: u64,
}

impl PairedPermutationTest {
    /// Configure a deterministic paired permutation test with `resamples` random
    /// sign-flip draws and a `seed`.
    pub fn new(resamples: NonZeroUsize, seed: u64) -> Self {
        Self { resamples, seed }
    }

    /// One-tailed paired sign-flip p-value for `differences`.
    ///
    /// Returns [`Error::Validation`] if `differences` is empty or contains a
    /// non-finite value.
    pub fn p_value(&self, differences: &[f64]) -> Result<f64> {
        self.p_value_with_salt(differences, 0)
    }

    fn p_value_with_salt(&self, differences: &[f64], salt: u64) -> Result<f64> {
        if differences.is_empty() {
            return Err(Error::validation(
                "paired permutation test needs at least one difference",
            ));
        }
        if let Some(bad) = differences.iter().find(|v| !v.is_finite()) {
            return Err(Error::validation(format!(
                "paired permutation differences must be finite; found {bad}"
            )));
        }
        let observed: f64 = differences.iter().sum();
        let mut greater = 0usize;
        let mut equal = 1usize; // The observed sign assignment, counted as a mid-p tie.
        let mut rng = SplitMix64::new(self.seed ^ mix_salt(salt));
        for _ in 0..self.resamples.get() {
            let sample = differences
                .iter()
                .map(|&d| if rng.next_u64() & 1 == 0 { d } else { -d })
                .sum::<f64>();
            if sample > observed {
                greater += 1;
            } else if sample == observed {
                equal += 1;
            }
        }
        Ok((greater as f64 + 0.5 * equal as f64) / (self.resamples.get() + 1) as f64)
    }
}

/// Item-blocked rank agreement between a metric and human judgment: how well the
/// metric's per-item deltas track human preference margins.
///
/// Produced by [`rank_agreement`]. [`kendall_tau`](Self::kendall_tau) and
/// [`spearman_rho`](Self::spearman_rho) are both in `[-1, 1]`; `+1` is perfect rank
/// agreement, `0` no monotone association, negative an inverted metric. They are read
/// against [`units`](Self::units): a high correlation over a handful of matched
/// comparisons is not yet evidence, which is why the deferred bootstrap will bracket
/// both with confidence intervals.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct RankAgreement {
    kendall: f64,
    spearman: f64,
    units: usize,
}

impl RankAgreement {
    /// Kendall's tau-b between human margins and metric deltas, in `[-1, 1]`.
    pub fn kendall_tau(&self) -> f64 {
        self.kendall
    }

    /// Spearman's rho between human margins and metric deltas, in `[-1, 1]`.
    pub fn spearman_rho(&self) -> f64 {
        self.spearman
    }

    /// The number of matched `(item, system-pair)` comparison units the correlations
    /// were computed over; the sample size behind the estimate.
    pub fn units(&self) -> usize {
        self.units
    }
}

/// Clustered-bootstrap intervals for [`rank_agreement`].
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct RankAgreementIntervals {
    kendall: ConfidenceInterval,
    spearman: ConfidenceInterval,
    units: usize,
}

impl RankAgreementIntervals {
    /// Interval for Kendall's tau-b.
    pub fn kendall_tau(&self) -> ConfidenceInterval {
        self.kendall
    }

    /// Interval for Spearman's rho.
    pub fn spearman_rho(&self) -> ConfidenceInterval {
        self.spearman
    }

    /// Number of matched comparison units in the observed data.
    pub fn units(&self) -> usize {
        self.units
    }
}

/// System-level Soft Pairwise Accuracy between human judgment and metric scores.
///
/// For every unordered system pair, the human margins and metric deltas are converted
/// into one-tailed paired permutation p-values. The score is the average
/// `1 - abs(p_human - p_metric)`, so `1` means the metric has the same pairwise
/// system preferences and confidence as humans, while `0` means the preferences are
/// maximally opposed. `pairs` is the number of system pairs in the average and `units`
/// is the number of matched `(item, system-pair)` observations that fed those pair
/// p-values.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct SoftPairwiseAccuracy {
    score: f64,
    pairs: usize,
    units: usize,
}

impl SoftPairwiseAccuracy {
    /// SPA score in `[0, 1]`.
    pub fn score(&self) -> f64 {
        self.score
    }

    /// Number of system pairs compared.
    pub fn pairs(&self) -> usize {
        self.pairs
    }

    /// Number of matched item/system-pair units used across all pairs.
    pub fn units(&self) -> usize {
        self.units
    }
}

/// Clustered-bootstrap interval for [`soft_pairwise_accuracy`].
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct SoftPairwiseAccuracyInterval {
    score: ConfidenceInterval,
    pairs: usize,
    units: usize,
}

impl SoftPairwiseAccuracyInterval {
    /// Interval for the SPA score.
    pub fn score(&self) -> ConfidenceInterval {
        self.score
    }

    /// Number of system pairs compared in the observed data.
    pub fn pairs(&self) -> usize {
        self.pairs
    }

    /// Number of matched item/system-pair units in the observed data.
    pub fn units(&self) -> usize {
        self.units
    }
}

/// Meta-evaluate `scores` against `judgments`: the item-blocked rank agreement between
/// the metric's per-item deltas and human preference margins.
///
/// For each conditioning item and unordered system pair, the *human margin* is the net
/// preference (a win for the lower-indexed system `+1`, a loss `-1`, a tie `0`, averaged
/// over the raters who judged that pair, one vote each) and the *metric delta* is the
/// difference of the two systems' scores on that item. A pair contributes a matched unit
/// only when it was both judged and scored for both systems; pairs missing either are
/// dropped. Both quantities are within-item, so item difficulty is differenced out (the
/// blocking the brief requires).
///
/// A margin is one vote per rater, so a rater who judged the same comparison twice is
/// rejected (the same shape [`Reliability::krippendorff_alpha`](super::Reliability)
/// rejects), rather than let one rater silently skew it. This is deliberately *not* a
/// [`BradleyTerry`](super::BradleyTerry) fit, which treats repeated comparisons as extra
/// games; a preference margin is an average of independent opinions.
///
/// Returns [`Error::Validation`] if `judgments` is empty, if a rater judged the same
/// comparison more than once, if fewer than two matched units survive (a correlation
/// needs at least two points), or if either the margins or the deltas have no
/// variability to correlate (see [`kendall_tau_b`] / [`spearman_rho`]).
pub fn rank_agreement(judgments: &[Judgment], scores: &MetricScores) -> Result<RankAgreement> {
    let clusters = matched_unit_clusters(judgments, scores)?;
    rank_agreement_from_clusters(&clusters)
}

/// Clustered-bootstrap intervals for [`rank_agreement`].
///
/// Resamples conditioning items, preserving the matched system-pair units inside each
/// item, then recomputes Kendall and Spearman on every resample. If a resample makes
/// either correlation undefined, the interval is rejected; the study needs more varied
/// item clusters before it can support a non-mean CI.
pub fn rank_agreement_intervals(
    judgments: &[Judgment],
    scores: &MetricScores,
    bootstrap: &ClusterBootstrap,
) -> Result<RankAgreementIntervals> {
    let clusters = matched_unit_clusters(judgments, scores)?;
    let point = rank_agreement_from_clusters(&clusters)?;
    let kendall = bootstrap.measure_cluster_statistic(&clusters, |sample| {
        Ok(rank_agreement_from_clusters(sample)?.kendall_tau())
    })?;
    let spearman = bootstrap.measure_cluster_statistic(&clusters, |sample| {
        Ok(rank_agreement_from_clusters(sample)?.spearman_rho())
    })?;
    Ok(RankAgreementIntervals {
        kendall,
        spearman,
        units: point.units(),
    })
}

/// System-level Soft Pairwise Accuracy between human judgments and metric scores.
///
/// `systems` is the closed set of system indices in the study; every pair among
/// `0..systems` must have at least one matched item-level human margin and metric
/// delta. Human and metric p-values are computed by `permutation`, then compared by
/// `1 - abs(p_human - p_metric)` and averaged over system pairs.
pub fn soft_pairwise_accuracy(
    systems: usize,
    judgments: &[Judgment],
    scores: &MetricScores,
    permutation: &PairedPermutationTest,
) -> Result<SoftPairwiseAccuracy> {
    validate_spa_inputs(systems, judgments, scores)?;
    let clusters = matched_unit_clusters(judgments, scores)?;
    soft_pairwise_accuracy_from_clusters(systems, &clusters, permutation)
}

/// Clustered-bootstrap interval for [`soft_pairwise_accuracy`].
///
/// The observed study must cover every system pair. Each bootstrap resample must also
/// cover every pair; otherwise the interval is rejected rather than changing the pair
/// denominator from draw to draw.
pub fn soft_pairwise_accuracy_interval(
    systems: usize,
    judgments: &[Judgment],
    scores: &MetricScores,
    permutation: &PairedPermutationTest,
    bootstrap: &ClusterBootstrap,
) -> Result<SoftPairwiseAccuracyInterval> {
    validate_spa_inputs(systems, judgments, scores)?;
    let clusters = matched_unit_clusters(judgments, scores)?;
    let point = soft_pairwise_accuracy_from_clusters(systems, &clusters, permutation)?;
    let score = bootstrap.measure_cluster_statistic(&clusters, |sample| {
        Ok(soft_pairwise_accuracy_from_clusters(systems, sample, permutation)?.score())
    })?;
    Ok(SoftPairwiseAccuracyInterval {
        score,
        pairs: point.pairs(),
        units: point.units(),
    })
}

/// The canonical `(lo, hi, vote)` of a judgment: the system pair in index order and the
/// preference for `lo` over `hi` (`+1` lo preferred, `-1` hi preferred, `0` tie), so the
/// two presentation orders of one comparison contribute the same signed vote.
fn canonical_vote(judgment: &Judgment) -> (usize, usize, f64) {
    let (lo, hi) = if judgment.a() < judgment.b() {
        (judgment.a(), judgment.b())
    } else {
        (judgment.b(), judgment.a())
    };
    let vote = match judgment.outcome() {
        Outcome::Tie => 0.0,
        Outcome::AWins => f64::from(if judgment.a() == lo { 1 } else { -1 }),
        Outcome::BWins => f64::from(if judgment.b() == lo { 1 } else { -1 }),
    };
    (lo, hi, vote)
}

#[derive(Clone, Copy, Debug, PartialEq)]
struct MatchedUnit {
    lo: usize,
    hi: usize,
    human_margin: f64,
    metric_delta: f64,
}

fn matched_unit_clusters(
    judgments: &[Judgment],
    scores: &MetricScores,
) -> Result<Vec<Vec<MatchedUnit>>> {
    if judgments.is_empty() {
        return Err(Error::validation("no judgments to meta-evaluate"));
    }

    let mut votes: BTreeMap<(usize, usize, usize), (f64, usize)> = BTreeMap::new();
    let mut rated: BTreeSet<(usize, usize, usize, usize)> = BTreeSet::new();
    for judgment in judgments {
        let (lo, hi, vote) = canonical_vote(judgment);
        if !rated.insert((judgment.item(), lo, hi, judgment.rater())) {
            return Err(Error::validation(format!(
                "rater {} judged the comparison (item {}, systems {lo}/{hi}) more than once",
                judgment.rater(),
                judgment.item()
            )));
        }
        let entry = votes.entry((judgment.item(), lo, hi)).or_insert((0.0, 0));
        entry.0 += vote;
        entry.1 += 1;
    }

    let mut by_item: BTreeMap<usize, Vec<MatchedUnit>> = BTreeMap::new();
    for ((item, lo, hi), (vote_sum, count)) in votes {
        if let (Some(score_lo), Some(score_hi)) = (scores.get(item, lo), scores.get(item, hi)) {
            by_item.entry(item).or_default().push(MatchedUnit {
                lo,
                hi,
                human_margin: vote_sum / count as f64,
                metric_delta: score_lo - score_hi,
            });
        }
    }
    let clusters: Vec<Vec<MatchedUnit>> = by_item.into_values().filter(|c| !c.is_empty()).collect();
    if clusters.is_empty() {
        return Err(Error::validation(
            "no matched human judgments and metric scores to meta-evaluate",
        ));
    }
    Ok(clusters)
}

fn rank_agreement_from_clusters(clusters: &[Vec<MatchedUnit>]) -> Result<RankAgreement> {
    let units = flatten_units(clusters);
    if units.len() < 2 {
        return Err(Error::validation(
            "rank agreement needs at least two matched (item, system-pair) comparison units",
        ));
    }
    let margins: Vec<f64> = units.iter().map(|u| u.human_margin).collect();
    let deltas: Vec<f64> = units.iter().map(|u| u.metric_delta).collect();
    Ok(RankAgreement {
        kendall: kendall_tau_b(&margins, &deltas)?,
        spearman: spearman_rho(&margins, &deltas)?,
        units: units.len(),
    })
}

fn soft_pairwise_accuracy_from_clusters(
    systems: usize,
    clusters: &[Vec<MatchedUnit>],
    permutation: &PairedPermutationTest,
) -> Result<SoftPairwiseAccuracy> {
    let units = flatten_units(clusters);
    let mut by_pair: BTreeMap<(usize, usize), (Vec<f64>, Vec<f64>)> = BTreeMap::new();
    for unit in &units {
        let entry = by_pair.entry((unit.lo, unit.hi)).or_default();
        entry.0.push(unit.human_margin);
        entry.1.push(unit.metric_delta);
    }

    let mut sum = 0.0f64;
    let mut pairs = 0usize;
    let mut used_units = 0usize;
    for lo in 0..systems {
        for hi in (lo + 1)..systems {
            let Some((human, metric)) = by_pair.get(&(lo, hi)) else {
                return Err(Error::validation(format!(
                    "soft pairwise accuracy needs at least one matched item for system pair {lo}/{hi}"
                )));
            };
            let salt = pair_salt(lo, hi);
            let human_p = permutation.p_value_with_salt(human, salt)?;
            let metric_p = permutation.p_value_with_salt(metric, salt)?;
            sum += 1.0 - (human_p - metric_p).abs();
            pairs += 1;
            used_units += human.len();
        }
    }
    Ok(SoftPairwiseAccuracy {
        score: sum / pairs as f64,
        pairs,
        units: used_units,
    })
}

fn flatten_units(clusters: &[Vec<MatchedUnit>]) -> Vec<MatchedUnit> {
    clusters
        .iter()
        .flat_map(|cluster| cluster.iter().copied())
        .collect()
}

fn validate_spa_inputs(
    systems: usize,
    judgments: &[Judgment],
    scores: &MetricScores,
) -> Result<()> {
    if systems < 2 {
        return Err(Error::validation(format!(
            "soft pairwise accuracy needs at least two systems; got {systems}"
        )));
    }
    for (idx, judgment) in judgments.iter().enumerate() {
        if judgment.a() >= systems || judgment.b() >= systems {
            return Err(Error::validation(format!(
                "judgment {idx} references system {} / {} but there are only {systems} systems",
                judgment.a(),
                judgment.b()
            )));
        }
    }
    for &(item, system) in scores.by_cell.keys() {
        if system >= systems {
            return Err(Error::validation(format!(
                "metric score (item {item}, system {system}) is outside the {systems}-system study"
            )));
        }
    }
    Ok(())
}

fn pair_salt(lo: usize, hi: usize) -> u64 {
    ((lo as u64) << 33) ^ ((hi as u64) << 1)
}

fn mix_salt(mut x: u64) -> u64 {
    x = (x ^ (x >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    x = (x ^ (x >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    x ^ (x >> 31)
}

/// Kendall's tau-b rank correlation between paired observations `xs` and `ys`.
///
/// Tau-b is the tie-corrected count `(concordant - discordant)` over all `n(n-1)/2`
/// pairs, normalized by the geometric mean of the untied-pair counts on each side, so a
/// pair tied in either coordinate weakens neither concordance nor discordance but still
/// shrinks the normalizer. The result is in `[-1, 1]`. The tie correction matters here:
/// human margins are heavily tied (many comparisons land at `+1` or `-1`).
///
/// Returns [`Error::Validation`] if the inputs differ in length, hold fewer than two
/// observations, contain a non-finite value, or leave no untied pairs on a side (one
/// coordinate is constant, so no ranking can be recovered). `O(n^2)`, sized for the
/// hundreds of matched comparisons of a small study.
pub fn kendall_tau_b(xs: &[f64], ys: &[f64]) -> Result<f64> {
    let n = check_pair(xs, ys)?;
    let (mut concordant, mut discordant): (i64, i64) = (0, 0);
    let (mut tied_x, mut tied_y): (i64, i64) = (0, 0);
    for i in 0..n {
        for j in (i + 1)..n {
            // Finite (checked), so partial_cmp is total here.
            let dx = xs[i].partial_cmp(&xs[j]).expect("finite");
            let dy = ys[i].partial_cmp(&ys[j]).expect("finite");
            let x_tie = dx == Ordering::Equal;
            let y_tie = dy == Ordering::Equal;
            if x_tie {
                tied_x += 1;
            }
            if y_tie {
                tied_y += 1;
            }
            if !x_tie && !y_tie {
                if dx == dy {
                    concordant += 1;
                } else {
                    discordant += 1;
                }
            }
        }
    }
    let pairs = n as i64 * (n as i64 - 1) / 2;
    let denom_sq = (pairs - tied_x) as f64 * (pairs - tied_y) as f64;
    if denom_sq <= 0.0 {
        return Err(Error::validation(
            "Kendall tau-b is undefined: one coordinate is constant (no untied pairs)",
        ));
    }
    Ok((concordant - discordant) as f64 / denom_sq.sqrt())
}

/// Spearman's rho rank correlation between paired observations `xs` and `ys`: the
/// Pearson correlation of their ranks, with tied values sharing their average rank.
///
/// In `[-1, 1]`. Returns [`Error::Validation`] under the same conditions as
/// [`kendall_tau_b`], including when one coordinate is constant (its ranks have zero
/// variance, so the correlation is undefined).
pub fn spearman_rho(xs: &[f64], ys: &[f64]) -> Result<f64> {
    check_pair(xs, ys)?;
    pearson(&average_ranks(xs), &average_ranks(ys))
}

/// The shared precondition of the rank correlations: equal length, at least two
/// observations, every value finite. Returns the common length.
fn check_pair(xs: &[f64], ys: &[f64]) -> Result<usize> {
    if xs.len() != ys.len() {
        return Err(Error::validation(format!(
            "paired statistics need equal-length inputs; got {} and {}",
            xs.len(),
            ys.len()
        )));
    }
    if xs.len() < 2 {
        return Err(Error::validation(
            "a rank correlation needs at least two observations",
        ));
    }
    if let Some(bad) = xs.iter().chain(ys).find(|v| !v.is_finite()) {
        return Err(Error::validation(format!(
            "rank-correlation inputs must be finite; found {bad}"
        )));
    }
    Ok(xs.len())
}

/// Fractional ranks of `values` (ascending), tied values sharing the average of the
/// ranks they span.
///
/// The sort uses `total_cmp` for a total, deterministic order, but tie *grouping* uses
/// numeric equality (`partial_cmp`), so `-0.0` and `+0.0` share a rank. That keeps
/// Spearman's notion of a tie identical to [`kendall_tau_b`]'s `partial_cmp`; the two
/// must not disagree on what counts as tied. Finiteness is guaranteed by [`check_pair`],
/// so `partial_cmp` is never `None`.
fn average_ranks(values: &[f64]) -> Vec<f64> {
    let n = values.len();
    let mut order: Vec<usize> = (0..n).collect();
    order.sort_by(|&a, &b| values[a].total_cmp(&values[b]));
    let mut ranks = vec![0.0f64; n];
    let mut i = 0;
    while i < n {
        let mut j = i + 1;
        while j < n && values[order[j]].partial_cmp(&values[order[i]]) == Some(Ordering::Equal) {
            j += 1;
        }
        // 1-based ranks (i+1)..=j over the tie group, whose average is their midpoint.
        let average = (i + 1 + j) as f64 / 2.0;
        for &idx in &order[i..j] {
            ranks[idx] = average;
        }
        i = j;
    }
    ranks
}

/// Pearson correlation of two equal-length finite vectors. Returns
/// [`Error::Validation`] if either has zero variance (an undefined correlation).
fn pearson(xs: &[f64], ys: &[f64]) -> Result<f64> {
    let n = xs.len() as f64;
    let mean_x = xs.iter().sum::<f64>() / n;
    let mean_y = ys.iter().sum::<f64>() / n;
    let (mut cov, mut var_x, mut var_y) = (0.0f64, 0.0f64, 0.0f64);
    for (&x, &y) in xs.iter().zip(ys) {
        let (dx, dy) = (x - mean_x, y - mean_y);
        cov += dx * dy;
        var_x += dx * dx;
        var_y += dy * dy;
    }
    let denom = (var_x * var_y).sqrt();
    if denom <= 0.0 {
        return Err(Error::validation(
            "Spearman rho is undefined: one coordinate has zero rank variance",
        ));
    }
    Ok(cov / denom)
}

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

    fn score(item: usize, system: usize, value: f64) -> MetricScore {
        MetricScore {
            item,
            system,
            value,
        }
    }

    fn permutation() -> PairedPermutationTest {
        PairedPermutationTest::new(NonZeroUsize::new(512).unwrap(), 17)
    }

    fn bootstrap() -> ClusterBootstrap {
        ClusterBootstrap::new(NonZeroUsize::new(64).unwrap(), 0.90, 23).unwrap()
    }

    #[test]
    fn kendall_perfect_concordance_and_discordance() {
        let xs = [1.0, 2.0, 3.0, 4.0];
        let asc = [10.0, 20.0, 30.0, 40.0];
        let desc = [40.0, 30.0, 20.0, 10.0];
        assert!((kendall_tau_b(&xs, &asc).unwrap() - 1.0).abs() < 1e-12);
        assert!((kendall_tau_b(&xs, &desc).unwrap() + 1.0).abs() < 1e-12);
    }

    #[test]
    fn kendall_matches_hand_count() {
        // 4 points, one discordant pair: (3,4) vs their y order is swapped.
        // pairs = 6, concordant = 5, discordant = 1 => tau = 4/6.
        let xs = [1.0, 2.0, 3.0, 4.0];
        let ys = [1.0, 2.0, 4.0, 3.0];
        assert!((kendall_tau_b(&xs, &ys).unwrap() - (4.0 / 6.0)).abs() < 1e-12);
    }

    #[test]
    fn spearman_perfect_and_reversed() {
        let xs = [1.0, 2.0, 3.0, 4.0, 5.0];
        let asc = [2.0, 4.0, 6.0, 8.0, 10.0];
        let desc = [10.0, 8.0, 6.0, 4.0, 2.0];
        assert!((spearman_rho(&xs, &asc).unwrap() - 1.0).abs() < 1e-12);
        assert!((spearman_rho(&xs, &desc).unwrap() + 1.0).abs() < 1e-12);
    }

    #[test]
    fn matching_ties_stay_perfect() {
        // The tie block lines up on both coordinates, so the tie correction recovers a
        // perfect +1: the tied pair leaves the numerator and the normalizer together.
        let xs = [1.0, 1.0, 3.0];
        let ys = [5.0, 5.0, 9.0];
        assert!((kendall_tau_b(&xs, &ys).unwrap() - 1.0).abs() < 1e-12);
        assert!((spearman_rho(&xs, &ys).unwrap() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn one_sided_ties_reduce_the_magnitude() {
        // A tie on y that x resolves is a genuine disagreement in rank, so neither
        // statistic reaches 1: tau-b = 2/sqrt(6), rho = 1.5/sqrt(3) by hand.
        let xs = [1.0, 2.0, 3.0];
        let ys = [5.0, 5.0, 9.0];
        assert!((kendall_tau_b(&xs, &ys).unwrap() - 2.0 / 6.0_f64.sqrt()).abs() < 1e-12);
        assert!((spearman_rho(&xs, &ys).unwrap() - 1.5 / 3.0_f64.sqrt()).abs() < 1e-12);
    }

    #[test]
    fn signed_zero_ties_consistently_across_both() {
        // -0.0 and +0.0 are numerically equal, so both statistics must treat them as one
        // tied value. Spearman sorts with total_cmp, which orders -0.0 before +0.0; were
        // tie grouping to follow that order it would split the pair and disagree with
        // Kendall's partial_cmp. Both must read +1 here.
        let xs = [-0.0, 0.0, 1.0];
        let ys = [5.0, 5.0, 9.0];
        assert!((kendall_tau_b(&xs, &ys).unwrap() - 1.0).abs() < 1e-12);
        assert!((spearman_rho(&xs, &ys).unwrap() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn correlations_reject_degenerate_inputs() {
        assert!(kendall_tau_b(&[1.0], &[2.0]).is_err()); // < 2 points
        assert!(kendall_tau_b(&[1.0, 2.0], &[3.0]).is_err()); // length mismatch
        assert!(kendall_tau_b(&[1.0, 2.0], &[f64::NAN, 3.0]).is_err()); // non-finite
        assert!(kendall_tau_b(&[7.0, 7.0, 7.0], &[1.0, 2.0, 3.0]).is_err()); // constant x
        assert!(spearman_rho(&[1.0, 2.0, 3.0], &[4.0, 4.0, 4.0]).is_err()); // constant y
    }

    /// A study comparing systems 0 and 1 over four items whose human preference margins
    /// vary (+1, +0.5, 0, -0.5) and whose metric deltas track them (+4, +2, 0, -1). The
    /// varying margins are what make the rank agreement well-defined: a metric cannot be
    /// validated against a human signal that never moves.
    fn agreeing_study() -> (Vec<Judgment>, MetricScores) {
        let judgments = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(0, 1, 0, 1, Outcome::AWins).unwrap(), // item 0 margin +1.0
            Judgment::new(1, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(1, 1, 0, 1, Outcome::Tie).unwrap(), // item 1 margin +0.5
            Judgment::new(2, 0, 0, 1, Outcome::Tie).unwrap(),
            Judgment::new(2, 1, 0, 1, Outcome::Tie).unwrap(), // item 2 margin 0.0
            Judgment::new(3, 0, 0, 1, Outcome::BWins).unwrap(),
            Judgment::new(3, 1, 0, 1, Outcome::Tie).unwrap(), // item 3 margin -0.5
        ];
        let scores = MetricScores::new(vec![
            score(0, 0, 5.0),
            score(0, 1, 1.0), // delta +4
            score(1, 0, 3.0),
            score(1, 1, 1.0), // delta +2
            score(2, 0, 1.0),
            score(2, 1, 1.0), // delta 0
            score(3, 0, 1.0),
            score(3, 1, 2.0), // delta -1
        ])
        .unwrap();
        (judgments, scores)
    }

    #[test]
    fn rank_agreement_is_perfect_when_metric_matches_humans() {
        let (judgments, scores) = agreeing_study();
        let agreement = rank_agreement(&judgments, &scores).unwrap();
        assert!((agreement.kendall_tau() - 1.0).abs() < 1e-12);
        assert!((agreement.spearman_rho() - 1.0).abs() < 1e-12);
        // One matched unit per item.
        assert_eq!(agreement.units(), 4);
    }

    #[test]
    fn rank_agreement_inverts_for_a_reversed_metric() {
        let (judgments, agreeing) = agreeing_study();
        // Negate every score: the same human margins, a metric ordered backwards.
        let reversed = MetricScores::new(
            [
                (0, 0, -5.0),
                (0, 1, -1.0),
                (1, 0, -3.0),
                (1, 1, -1.0),
                (2, 0, -1.0),
                (2, 1, -1.0),
                (3, 0, -1.0),
                (3, 1, -2.0),
            ]
            .map(|(i, s, v)| score(i, s, v))
            .to_vec(),
        )
        .unwrap();
        assert_eq!(agreeing.get(0, 0), Some(5.0));
        let agreement = rank_agreement(&judgments, &reversed).unwrap();
        assert!((agreement.kendall_tau() + 1.0).abs() < 1e-12);
        assert!((agreement.spearman_rho() + 1.0).abs() < 1e-12);
    }

    #[test]
    fn rank_agreement_drops_unmatched_units() {
        // Each item judges 0-vs-1 and 1-vs-2, but the metric never scored system 2, so
        // only the 0-vs-1 unit survives per item. Their margins differ (+1, +0.5), so
        // the two surviving units still define a correlation.
        let judgments = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(0, 0, 1, 2, Outcome::AWins).unwrap(),
            Judgment::new(1, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(1, 1, 0, 1, Outcome::Tie).unwrap(),
            Judgment::new(1, 0, 1, 2, Outcome::AWins).unwrap(),
        ];
        let scores = MetricScores::new(vec![
            score(0, 0, 2.0),
            score(0, 1, 1.0),
            score(1, 0, 4.0),
            score(1, 1, 1.0),
        ])
        .unwrap();
        let agreement = rank_agreement(&judgments, &scores).unwrap();
        assert_eq!(agreement.units(), 2);
    }

    #[test]
    fn rank_agreement_canonicalizes_presentation_order() {
        // The same item/pair judged both ways with consistent preference is one unit;
        // the flipped presentation must not change the margin's sign.
        let judgments = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(0, 1, 1, 0, Outcome::BWins).unwrap(),
            Judgment::new(1, 0, 0, 1, Outcome::BWins).unwrap(),
        ];
        let scores = MetricScores::new(vec![
            score(0, 0, 5.0),
            score(0, 1, 1.0),
            score(1, 0, 1.0),
            score(1, 1, 5.0),
        ])
        .unwrap();
        // Item 0: both judgments prefer system 0 (margin +1, delta +4). Item 1: prefer
        // system 1 (margin -1, delta -4). Perfect agreement.
        let agreement = rank_agreement(&judgments, &scores).unwrap();
        assert_eq!(agreement.units(), 2);
        assert!((agreement.kendall_tau() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn rank_agreement_rejects_empty_and_too_few_units() {
        let scores = MetricScores::new(vec![score(0, 0, 1.0), score(0, 1, 2.0)]).unwrap();
        assert!(rank_agreement(&[], &scores).is_err());
        // One judged-and-scored pair: a single matched unit cannot be correlated.
        let one = vec![Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap()];
        assert!(rank_agreement(&one, &scores).is_err());
    }

    #[test]
    fn rank_agreement_rejects_duplicate_rater_on_a_unit() {
        // Same rater, same comparison, twice (here via the flipped presentation order):
        // a margin is one vote per rater, so this is rejected, not silently doubled. The
        // sibling reliability read rejects the same shape.
        let scores = MetricScores::new(vec![score(0, 0, 2.0), score(0, 1, 1.0)]).unwrap();
        let js = vec![
            Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap(),
            Judgment::new(0, 0, 1, 0, Outcome::BWins).unwrap(),
        ];
        assert!(rank_agreement(&js, &scores).is_err());
    }

    #[test]
    fn metric_scores_reject_duplicate_empty_and_nonfinite() {
        assert!(MetricScores::new(vec![]).is_err());
        assert!(MetricScores::new(vec![score(0, 0, f64::NAN)]).is_err());
        assert!(MetricScores::new(vec![score(0, 0, 1.0), score(0, 0, 2.0)]).is_err());
    }

    #[test]
    fn paired_permutation_orients_preferences() {
        let p = permutation();
        assert!(p.p_value(&[1.0; 8]).unwrap() < 0.05);
        assert!(p.p_value(&[-1.0; 8]).unwrap() > 0.95);
        assert!((p.p_value(&[0.0; 8]).unwrap() - 0.5).abs() < 1e-12);
        assert!(p.p_value(&[]).is_err());
        assert!(p.p_value(&[1.0, f64::NAN]).is_err());
    }

    fn complete_three_system_study(reversed_metric: bool) -> (Vec<Judgment>, MetricScores) {
        let mut judgments = Vec::new();
        let mut scores = Vec::new();
        for item in 0..8 {
            judgments.push(Judgment::new(item, 0, 0, 1, Outcome::AWins).unwrap());
            judgments.push(Judgment::new(item, 0, 0, 2, Outcome::AWins).unwrap());
            judgments.push(Judgment::new(item, 0, 1, 2, Outcome::AWins).unwrap());
            let values = if reversed_metric {
                [1.0, 2.0, 3.0]
            } else {
                [3.0, 2.0, 1.0]
            };
            for (system, value) in values.into_iter().enumerate() {
                scores.push(score(item, system, value));
            }
        }
        (judgments, MetricScores::new(scores).unwrap())
    }

    #[test]
    fn soft_pairwise_accuracy_rewards_matching_system_preferences() {
        let (judgments, scores) = complete_three_system_study(false);
        let spa = soft_pairwise_accuracy(3, &judgments, &scores, &permutation()).unwrap();
        assert_eq!(spa.pairs(), 3);
        assert_eq!(spa.units(), 24);
        assert!(spa.score() > 0.98, "{}", spa.score());
    }

    #[test]
    fn soft_pairwise_accuracy_penalizes_reversed_system_preferences() {
        let (judgments, scores) = complete_three_system_study(true);
        let spa = soft_pairwise_accuracy(3, &judgments, &scores, &permutation()).unwrap();
        assert!(spa.score() < 0.05, "{}", spa.score());
    }

    #[test]
    fn soft_pairwise_accuracy_is_exact_for_identical_margins_and_deltas() {
        let margins = [
            (Outcome::AWins, Outcome::AWins, 1.0),
            (Outcome::AWins, Outcome::AWins, 1.0),
            (Outcome::AWins, Outcome::AWins, 1.0),
            (Outcome::BWins, Outcome::BWins, -1.0),
            (Outcome::BWins, Outcome::BWins, -1.0),
            (Outcome::AWins, Outcome::Tie, 0.5),
            (Outcome::BWins, Outcome::Tie, -0.5),
            (Outcome::Tie, Outcome::Tie, 0.0),
        ];
        let mut judgments = Vec::new();
        let mut scores = Vec::new();
        for (item, (r0, r1, delta)) in margins.into_iter().enumerate() {
            judgments.push(Judgment::new(item, 0, 0, 1, r0).unwrap());
            judgments.push(Judgment::new(item, 1, 0, 1, r1).unwrap());
            scores.push(score(item, 0, delta));
            scores.push(score(item, 1, 0.0));
        }
        let scores = MetricScores::new(scores).unwrap();
        let permutation = PairedPermutationTest::new(NonZeroUsize::new(64).unwrap(), 17);
        let spa = soft_pairwise_accuracy(2, &judgments, &scores, &permutation).unwrap();
        assert_eq!(spa.score(), 1.0);
    }

    #[test]
    fn soft_pairwise_accuracy_requires_declared_pair_coverage() {
        let judgments = vec![Judgment::new(0, 0, 0, 1, Outcome::AWins).unwrap()];
        let scores = MetricScores::new(vec![score(0, 0, 2.0), score(0, 1, 1.0)]).unwrap();
        assert!(soft_pairwise_accuracy(1, &judgments, &scores, &permutation()).is_err());
        assert!(soft_pairwise_accuracy(3, &judgments, &scores, &permutation()).is_err());
        let bad_scores = MetricScores::new(vec![score(0, 0, 2.0), score(0, 3, 1.0)]).unwrap();
        assert!(soft_pairwise_accuracy(3, &judgments, &bad_scores, &permutation()).is_err());
    }

    fn varied_two_system_study() -> (Vec<Judgment>, MetricScores) {
        let patterns = [
            (Outcome::AWins, Outcome::AWins, 4.0, 1.0),
            (Outcome::AWins, Outcome::Tie, 3.0, 1.0),
            (Outcome::Tie, Outcome::Tie, 2.0, 2.0),
            (Outcome::BWins, Outcome::Tie, 1.0, 2.0),
            (Outcome::BWins, Outcome::BWins, 1.0, 4.0),
            (Outcome::AWins, Outcome::AWins, 5.0, 0.0),
            (Outcome::Tie, Outcome::AWins, 3.0, 2.0),
            (Outcome::BWins, Outcome::BWins, 0.0, 5.0),
        ];
        let mut judgments = Vec::new();
        let mut scores = Vec::new();
        for (item, (r0, r1, s0, s1)) in patterns.into_iter().enumerate() {
            judgments.push(Judgment::new(item, 0, 0, 1, r0).unwrap());
            judgments.push(Judgment::new(item, 1, 0, 1, r1).unwrap());
            scores.push(score(item, 0, s0));
            scores.push(score(item, 1, s1));
        }
        (judgments, MetricScores::new(scores).unwrap())
    }

    #[test]
    fn rank_agreement_intervals_resample_items() {
        let (judgments, scores) = varied_two_system_study();
        let point = rank_agreement(&judgments, &scores).unwrap();
        let intervals = rank_agreement_intervals(&judgments, &scores, &bootstrap()).unwrap();
        assert_eq!(intervals.units(), point.units());
        assert!((intervals.kendall_tau().point() - point.kendall_tau()).abs() < 1e-12);
        assert!((intervals.spearman_rho().point() - point.spearman_rho()).abs() < 1e-12);
        assert!(intervals.kendall_tau().width() >= 0.0);
        assert!(intervals.spearman_rho().width() >= 0.0);
    }

    #[test]
    fn soft_pairwise_accuracy_interval_resamples_items() {
        let (judgments, scores) = complete_three_system_study(false);
        let point = soft_pairwise_accuracy(3, &judgments, &scores, &permutation()).unwrap();
        let interval =
            soft_pairwise_accuracy_interval(3, &judgments, &scores, &permutation(), &bootstrap())
                .unwrap();
        assert_eq!(interval.pairs(), point.pairs());
        assert_eq!(interval.units(), point.units());
        assert!((interval.score().point() - point.score()).abs() < 1e-12);
        assert!(interval.score().lo() >= 0.0);
        assert!(interval.score().hi() <= 1.0);
    }

    use proptest::prelude::*;

    /// Paired integer-valued vectors (as `f64`), 2..12 long. Integers make ties common,
    /// exercising the tie corrections; the value range keeps things small.
    fn arb_pairs() -> impl Strategy<Value = (Vec<f64>, Vec<f64>)> {
        (2usize..12).prop_flat_map(|n| {
            let col = || proptest::collection::vec((-5i8..5).prop_map(f64::from), n);
            (col(), col())
        })
    }

    proptest! {
        #[test]
        fn correlations_stay_in_range((xs, ys) in arb_pairs()) {
            // Only defined when both coordinates vary; when defined, both are in [-1, 1].
            if let Ok(tau) = kendall_tau_b(&xs, &ys) {
                prop_assert!((-1.0 - 1e-9..=1.0 + 1e-9).contains(&tau), "tau {tau}");
            }
            if let Ok(rho) = spearman_rho(&xs, &ys) {
                prop_assert!((-1.0 - 1e-9..=1.0 + 1e-9).contains(&rho), "rho {rho}");
            }
        }

        #[test]
        fn correlations_are_symmetric((xs, ys) in arb_pairs()) {
            // corr(x, y) == corr(y, x); both branches agree on definedness.
            match (kendall_tau_b(&xs, &ys), kendall_tau_b(&ys, &xs)) {
                (Ok(a), Ok(b)) => prop_assert!((a - b).abs() < 1e-9),
                (Err(_), Err(_)) => {}
                _ => prop_assert!(false, "kendall definedness differs under swap"),
            }
            match (spearman_rho(&xs, &ys), spearman_rho(&ys, &xs)) {
                (Ok(a), Ok(b)) => prop_assert!((a - b).abs() < 1e-9),
                (Err(_), Err(_)) => {}
                _ => prop_assert!(false, "spearman definedness differs under swap"),
            }
        }

        #[test]
        fn correlations_are_permutation_invariant((xs, ys) in arb_pairs(), seed in any::<u64>()) {
            // Reordering the paired observations together cannot change a correlation.
            let mut idx: Vec<usize> = (0..xs.len()).collect();
            // A deterministic shuffle from the seed (Fisher-Yates, SplitMix-style step).
            let mut state = seed;
            for i in (1..idx.len()).rev() {
                state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
                idx.swap(i, (state >> 33) as usize % (i + 1));
            }
            let px: Vec<f64> = idx.iter().map(|&i| xs[i]).collect();
            let py: Vec<f64> = idx.iter().map(|&i| ys[i]).collect();
            if let (Ok(a), Ok(b)) = (kendall_tau_b(&xs, &ys), kendall_tau_b(&px, &py)) {
                prop_assert!((a - b).abs() < 1e-9);
            }
            if let (Ok(a), Ok(b)) = (spearman_rho(&xs, &ys), spearman_rho(&px, &py)) {
                prop_assert!((a - b).abs() < 1e-9);
            }
        }

        #[test]
        fn strictly_monotone_is_perfect(xs in proptest::collection::vec(-50.0f64..50.0, 2..10)) {
            // Distinct, strictly increasing x and y => both correlations are exactly +1.
            let mut sorted: Vec<f64> = xs;
            sorted.sort_by(f64::total_cmp);
            sorted.dedup();
            prop_assume!(sorted.len() >= 2);
            let ys: Vec<f64> = sorted.iter().map(|v| v * 2.0 + 1.0).collect();
            prop_assert!((kendall_tau_b(&sorted, &ys).unwrap() - 1.0).abs() < 1e-9);
            prop_assert!((spearman_rho(&sorted, &ys).unwrap() - 1.0).abs() < 1e-9);
        }

        #[test]
        fn correlations_are_reproducible((xs, ys) in arb_pairs()) {
            prop_assert_eq!(kendall_tau_b(&xs, &ys).ok(), kendall_tau_b(&xs, &ys).ok());
            prop_assert_eq!(spearman_rho(&xs, &ys).ok(), spearman_rho(&xs, &ys).ok());
        }
    }
}