use cord_sparse::*;
use std::time::Instant;

fn product_fn(fp: &FixedPoint) -> impl Fn(&[i64]) -> i64 + '_ {
    move |pt: &[i64]| {
        let mut result = fp.one();
        for &x in pt {
            result = fp.mul(result, x);
        }
        result
    }
}

fn measure_runtime<F: FnMut() -> f64>(mut f: F, k: usize) -> f64 {
    let n = 2 * k + 1;
    let mut results: Vec<f64> = (0..n).map(|_| f()).collect();
    results.sort_by(|a, b| a.partial_cmp(b).unwrap());
    results[k]
}

fn main() {
    let fp = FixedPoint::new(32);
    let basis = MonomialBasis;
    let pts = GoldenPoints;

    let dimensions = [2usize, 4, 8, 16, 32, 64];
    let max_num_points: usize = 100_000_000;
    let min_quotient = 2.0;
    let max_runtime = 10.0;
    let epsilon = 0.01;
    let k = 2;

    // warm-up
    eprintln!("Warm-up");
    {
        let iter = BoundedSumIter::new(5, 60);
        let f = product_fn(&fp);
        let rhs = evaluate_function(&iter, &f, &pts, &fp);
        let mut op = create_interpolation_operator(fp, &iter, &basis, &pts);
        let _ = op.solve(rhs);
    }
    eprintln!("Warm-up finished");

    let mut csv = String::new();

    for &d in &dimensions {
        eprintln!("Dimension: {d}");
        let mut last_num_points = 0usize;
        let mut last_runtime = 0.0f64;

        for bound in 1.. {
            let num_points = index::binom(bound + d, d);
            if num_points > max_num_points {
                break;
            }
            if num_points as f64 / (last_num_points as f64 + epsilon) * last_runtime > max_runtime {
                break;
            }
            if num_points < (last_num_points as f64 * min_quotient) as usize {
                continue;
            }

            last_num_points = num_points;
            eprintln!("  Bound: {bound}, points: {num_points}");

            let runtime = measure_runtime(
                || {
                    let iter = BoundedSumIter::new(d, bound);
                    let f = product_fn(&fp);
                    let rhs = evaluate_function(&iter, &f, &pts, &fp);
                    let mut op = create_interpolation_operator(fp, &iter, &basis, &pts);

                    let start = Instant::now();
                    let _ = op.solve(rhs);
                    let elapsed = start.elapsed().as_secs_f64();
                    eprintln!("    Time for solve(): {elapsed:.6} s");
                    elapsed
                },
                k,
            );

            csv.push_str(&format!("{d}, {bound}, {num_points}, {runtime}\n"));
            last_runtime = runtime;
        }
    }

    println!("{csv}");
}