ported/paired_stats_common

fn means(pairs: &[(f64, f64)]) -> Result<(f64, f64), WorksheetErrorCode> {
    if pairs.is_empty() {
        return Err(WorksheetErrorCode::Div0);
    }
    let (sum_x, sum_y) = pairs
        .iter()
        .fold((0.0, 0.0), |(ax, ay), (x, y)| (ax + x, ay + y));
    let n = pairs.len() as f64;
    Ok((sum_x / n, sum_y / n))
}

pub fn covariance_from_pairs(
    pairs: &[(f64, f64)],
    divisor: CovarianceDivisor,
) -> Result<f64, WorksheetErrorCode> {
    if pairs.is_empty() {
        return Err(WorksheetErrorCode::Div0);
    }
    if matches!(divisor, CovarianceDivisor::Sample) && pairs.len() < 2 {
        return Err(WorksheetErrorCode::Div0);
    }

    let (mean_x, mean_y) = means(pairs)?;
    let sum = pairs
        .iter()
        .map(|(x, y)| (x - mean_x) * (y - mean_y))
        .sum::<f64>();
    let denom = match divisor {
        CovarianceDivisor::Population => pairs.len() as f64,
        CovarianceDivisor::Sample => (pairs.len() - 1) as f64,
    };
    Ok(sum / denom)
}

pub fn correlation_from_pairs(pairs: &[(f64, f64)]) -> Result<f64, WorksheetErrorCode> {
    if pairs.len() < 2 {
        return Err(WorksheetErrorCode::Div0);
    }
    let (mean_x, mean_y) = means(pairs)?;
    let mut sum_xy = 0.0;
    let mut sum_x2 = 0.0;
    let mut sum_y2 = 0.0;
    for (x, y) in pairs {
        let dx = *x - mean_x;
        let dy = *y - mean_y;
        sum_xy += dx * dy;
        sum_x2 += dx * dx;
        sum_y2 += dy * dy;
    }
    if sum_x2 == 0.0 || sum_y2 == 0.0 {
        return Err(WorksheetErrorCode::Div0);
    }
    let denom = (pairs.len() - 1) as f64;
    let covariance = sum_xy / denom;
    let variance_x = sum_x2 / denom;
    let variance_y = sum_y2 / denom;
    Ok(covariance / (variance_x.sqrt() * variance_y.sqrt()))
}

pub fn slope_from_pairs(pairs: &[(f64, f64)]) -> Result<f64, WorksheetErrorCode> {
    if pairs.len() < 2 {
        return Err(WorksheetErrorCode::Div0);
    }
    let (mean_x, mean_y) = means(pairs)?;
    let mut sum_xy = 0.0;
    let mut sum_x2 = 0.0;
    for (x, y) in pairs {
        let dx = *x - mean_x;
        sum_xy += dx * (*y - mean_y);
        sum_x2 += dx * dx;
    }
    if sum_x2 == 0.0 {
        return Err(WorksheetErrorCode::Div0);
    }
    Ok(sum_xy / sum_x2)
}