//! Descriptive statistics and bootstrap confidence intervals for timing samples.
//! Floats are used here only for statistics (never for field arithmetic).

use rand::{RngCore, SeedableRng};
use rand_chacha::ChaCha20Rng;
use serde::Serialize;

/// Linear-interpolated quantile of a sorted slice (same convention as numpy `linear`).
pub fn quantile_sorted(sorted: &[f64], q: f64) -> f64 {
    if sorted.is_empty() {
        return f64::NAN;
    }
    let pos = q * (sorted.len() as f64 - 1.0);
    let lo = pos.floor() as usize;
    let hi = pos.ceil() as usize;
    if lo == hi {
        sorted[lo]
    } else {
        sorted[lo] + (sorted[hi] - sorted[lo]) * (pos - lo as f64)
    }
}

pub fn median(xs: &[f64]) -> f64 {
    let mut v = xs.to_vec();
    v.sort_by(|a, b| a.partial_cmp(b).unwrap());
    quantile_sorted(&v, 0.5)
}

#[derive(Serialize, Clone, Debug)]
pub struct Timing {
    pub reps: usize,
    pub ops_per_rep: usize,
    pub ns_per_op_samples: Vec<f64>,
    pub median_ns: f64,
    pub q1_ns: f64,
    pub q3_ns: f64,
    pub iqr_ns: f64,
    pub min_ns: f64,
    pub mean_ns: f64,
    pub sd_ns: f64,
    pub ops_per_s: f64,
}

impl Timing {
    pub fn from_samples(ops_per_rep: usize, samples: Vec<f64>) -> Timing {
        let mut v = samples.clone();
        v.sort_by(|a, b| a.partial_cmp(b).unwrap());
        let n = v.len() as f64;
        let mean = v.iter().sum::<f64>() / n;
        let sd = if v.len() > 1 {
            (v.iter().map(|x| (x - mean) * (x - mean)).sum::<f64>() / (n - 1.0)).sqrt()
        } else {
            0.0
        };
        let med = quantile_sorted(&v, 0.5);
        let q1 = quantile_sorted(&v, 0.25);
        let q3 = quantile_sorted(&v, 0.75);
        Timing {
            reps: samples.len(),
            ops_per_rep,
            ns_per_op_samples: samples,
            median_ns: med,
            q1_ns: q1,
            q3_ns: q3,
            iqr_ns: q3 - q1,
            min_ns: v[0],
            mean_ns: mean,
            sd_ns: sd,
            ops_per_s: 1e9 / med,
        }
    }
}

#[derive(Serialize, Clone, Debug)]
pub struct RatioCI {
    pub point: f64,
    pub ci95_low: f64,
    pub ci95_high: f64,
    pub bootstrap_resamples: usize,
}

/// Bootstrap 95 % CI for the ratio median(num) / median(den), resampling each sample set with
/// replacement (percentile method). Deterministic given `seed`.
pub fn bootstrap_ratio_of_medians(num: &[f64], den: &[f64], resamples: usize, seed: u64) -> RatioCI {
    let mut rng = ChaCha20Rng::seed_from_u64(seed);
    let mut ratios = Vec::with_capacity(resamples);
    let mut bn = vec![0.0; num.len()];
    let mut bd = vec![0.0; den.len()];
    for _ in 0..resamples {
        for x in bn.iter_mut() {
            *x = num[(rng.next_u64() % num.len() as u64) as usize];
        }
        for x in bd.iter_mut() {
            *x = den[(rng.next_u64() % den.len() as u64) as usize];
        }
        ratios.push(median(&bn) / median(&bd));
    }
    ratios.sort_by(|a, b| a.partial_cmp(b).unwrap());
    RatioCI {
        point: median(num) / median(den),
        ci95_low: quantile_sorted(&ratios, 0.025),
        ci95_high: quantile_sorted(&ratios, 0.975),
        bootstrap_resamples: resamples,
    }
}