catalogue/statistics

/// Every argument's numbers by the aggregate rule `held`.
fn all(site: &Site<'_>, a: &[CalcValue], sv: &mut EvalServices<'_>, held: Held) -> Step<Vec<f64>> {
    numbers(site, a, 0, a.len(), sv, held)
}

/// Argument `index`'s numbers by the aggregate rule for numbers.
fn list(site: &Site<'_>, a: &[CalcValue], index: usize, sv: &mut EvalServices<'_>) -> Step<Vec<f64>> {
    numbers(site, a, index, index + 1, sv, Held::Numbers)
}

/// Reads the direct arguments through number coercion first, so that their error comes before
/// any held error, as MEDIAN, DEVSQ and HARMEAN order them.
fn direct_first(site: &Site<'_>, a: &[CalcValue]) -> Step<()> {
    for (index, argument) in a.iter().enumerate() {
        if !matches!(argument.core_ref(), CoreRef::Reference(_) | CoreRef::Array(_)) {
            site.number_argument(index, argument)?;
        }
    }
    Ok(())
}

fn sorted(mut values: Vec<f64>) -> Vec<f64> {
    values.sort_by(f64::total_cmp);
    values
}

/// The k-th largest (`largest`) or smallest number of argument 0.
fn nth(site: &Site<'_>, a: &[CalcValue], sv: &mut EvalServices<'_>, largest: bool) -> Step<CalcValue> {
    let values = list(site, a, 0, sv)?;
    let k = num(site, a, 1)?.trunc();
    if k < 1.0 || k > values.len() as f64 {
        return fail(WorksheetErrorCode::Num);
    }
    let mut values = sorted(values);
    if largest {
        values.reverse();
    }
    site.number(values[k as usize - 1])
}

/// RANK's number: a number as it is, an error as itself, anything else #N/A.
fn rank_number(a: &[CalcValue]) -> Step<f64> {
    match a[0].core_ref() {
        CoreRef::Number(n) => Ok(n),
        CoreRef::Error(code) => fail(code),
        _ => fail(WorksheetErrorCode::NA),
    }
}

fn rank(site: &Site<'_>, a: &[CalcValue], sv: &mut EvalServices<'_>, average: bool) -> Step<CalcValue> {
    let number = rank_number(a)?;
    let values = list(site, a, 1, sv)?;
    let order = if given(a, 2) && num(site, a, 2)?.trunc() != 0.0 {
        RankOrder::Ascending
    } else {
        RankOrder::Descending
    };
    let r = if average {
        rank_avg(number, &values, order)
    } else {
        rank_eq(number, &values, order)
    };
    out(site, r)
}

fn percentile(site: &Site<'_>, a: &[CalcValue], sv: &mut EvalServices<'_>, exclusive: bool) -> Step<CalcValue> {
    let mut values = list(site, a, 0, sv)?;
    let p = num(site, a, 1)?;
    out(
        site,
        if exclusive {
            k::percentile_common::percentile_exc_kernel(&mut values, p)
        } else {
            k::percentile_common::percentile_inc_kernel(&mut values, p)
        },
    )
}

fn quartile(site: &Site<'_>, a: &[CalcValue], sv: &mut EvalServices<'_>, exclusive: bool) -> Step<CalcValue> {
    let mut values = list(site, a, 0, sv)?;
    let q = num(site, a, 1)?.trunc();
    let admitted = if exclusive { 1.0..=3.0 } else { 0.0..=4.0 };
    if !admitted.contains(&q) {
        return fail(WorksheetErrorCode::Num);
    }
    out(
        site,
        if exclusive {
            k::percentile_common::percentile_exc_kernel(&mut values, q / 4.0)
        } else {
            k::percentile_common::percentile_inc_kernel(&mut values, q / 4.0)
        },
    )
}

fn percent_rank(site: &Site<'_>, a: &[CalcValue], sv: &mut EvalServices<'_>, mode: PercentRankMode) -> Step<CalcValue> {
    let mut values = list(site, a, 0, sv)?;
    let x = num(site, a, 1)?;
    let significance = opt(site, a, 2, 3.0)?.trunc() as i64;
    out(site, percentrank(&mut values, x, significance, mode))
}

fn variance(site: &Site<'_>, a: &[CalcValue], sv: &mut EvalServices<'_>, held: Held, divisor: VarianceDivisor, root: bool) -> Step<CalcValue> {
    let values = all(site, a, sv, held)?;
    out(
        site,
        if root {
            stdev_from_values(&values, divisor)
        } else {
            variance_from_values(&values, divisor)
        },
    )
}

/// The most frequent numbers, each with its count, in order of first appearance.
fn modes(values: &[f64]) -> Vec<f64> {
    let mut counts: BTreeMap<u64, (usize, usize)> = BTreeMap::new();
    for (position, value) in values.iter().enumerate() {
        let key = (if *value == 0.0 { 0.0 } else { *value }).to_bits();
        counts.entry(key).or_insert((position, 0)).1 += 1;
    }
    let best = counts.values().map(|(_, c)| *c).max().unwrap_or(0);
    if best < 2 {
        return Vec::new();
    }
    let mut found: Vec<(usize, f64)> = counts
        .iter()
        .filter(|(_, (_, c))| *c == best)
        .map(|(key, (first, _))| (*first, f64::from_bits(*key)))
        .collect();
    found.sort_by_key(|(first, _)| *first);
    found.into_iter().map(|(_, v)| v).collect()
}

/// FORECAST's line through the pairs (x, y), evaluated at `x`.
fn forecast(site: &Site<'_>, x: f64, pairs: &[(f64, f64)]) -> Step<CalcValue> {
    let n = pairs.len();
    if n == 0 {
        return fail(WorksheetErrorCode::Div0);
    }
    let nf = n as f64;
    let mean_x = pairs.iter().map(|(x, _)| *x).sum::<f64>() / nf;
    let mean_y = pairs.iter().map(|(_, y)| *y).sum::<f64>() / nf;
    let (mut num, mut den) = (0.0, 0.0);
    for (xv, yv) in pairs {
        let dx = xv - mean_x;
        num += dx * (yv - mean_y);
        den += dx * dx;
    }
    if den == 0.0 {
        return fail(WorksheetErrorCode::Div0);
    }
    let b = num / den;
    let a = mean_y - b * mean_x;
    site.number(a + b * x)
}

/// LINEST, LOGEST, TREND and GROWTH: OxFunc's evaluation over the arguments with each reference
/// read into an array (a blank cell is Empty, which those functions refuse with #VALUE!), its
/// number results published by the number rule.
fn regression(
    site: &Site<'_>,
    a: &[CalcValue],
    sv: &mut EvalServices<'_>,
    evaluate: fn(&[CalcValue]) -> Result<CalcValue, WorksheetErrorCode>,
    publication: RegressionNumbers,
) -> Step<CalcValue> {
    let mut read = Vec::with_capacity(a.len());
    for (index, argument) in a.iter().enumerate() {
        read.push(match argument.core_ref() {
            CoreRef::Reference(_) => {
                let g = grid(site, a, index, sv)?;
                let rows: Vec<Vec<CalcValue>> =
                    g.cells.chunks(g.cols).map(<[CalcValue]>::to_vec).collect();
                CalcValue::array(
                    CalcArray::from_rows(rows).ok_or_else(|| site.not_modelled(Some(index)))?,
                )
            }
            _ => argument.clone(),
        });
    }
    let result = evaluate(&read).map_err(Exit::Error)?;
    // The retained coefficient-publication observations distinguish LOGEST's
    // subnormal coefficient row from GROWTH's zero forecasts. Statistics rows
    // keep the existing number rule; normal and nonfinite outputs do too.
    let number = |n: f64, first_row: bool| match publication {
        RegressionNumbers::Coefficients if first_row && n.is_subnormal() => out_with_subnormal(site, Ok(n)),
        RegressionNumbers::Forecast if n.is_subnormal() => site.number(0.0),
        _ => site.number(n),
    };
    match result.core_ref() {
        CoreRef::Number(n) => number(n, true),
        CoreRef::Array(array) => {
            let mut rows = Vec::with_capacity(array.shape().rows);
            for r in 0..array.shape().rows {
                let mut row = Vec::with_capacity(array.shape().cols);
                for cell in array.row_slice(r).unwrap_or(&[]) {
                    row.push(match cell.core_ref() {
                        CoreRef::Number(n) => match number(n, r == 0) {
                            Ok(v) => v,
                            Err(Exit::Error(code)) => CalcValue::error(code),
                            Err(fault) => return Err(fault),
                        },
                        _ => cell.clone(),
                    });
                }
                rows.push(row);
            }
            CalcArray::from_rows(rows)
                .map(CalcValue::array)
                .ok_or_else(|| site.not_modelled(None))
        }
        _ => Ok(result),
    }
}