//! Criterion benchmarks for nexus-stats hot-path update methods.
//!
//! Each benchmark primes the type with 1000 samples, then measures the
//! steady-state update cost. Uses `black_box` to prevent dead code elimination.
//!
//! Run on controlled system for meaningful results:
//! ```bash
//! echo 1 | sudo tee /sys/devices/system/cpu/intel_pstate/no_turbo
//! taskset -c 0 cargo bench -p nexus-stats --bench update_bench
//! echo 0 | sudo tee /sys/devices/system/cpu/intel_pstate/no_turbo
//! ```

use criterion::{Criterion, black_box, criterion_group, criterion_main};
use nexus_stats::{
    control::{DeadBandF64, FirstDiffF64, HysteresisF64, LevelCrossingF64, SecondDiffF64},
    detection::{
        AdaptiveThresholdF64, CusumF64, RobustZScoreF64, ShiryaevRobertsF64, TrendAlertF64,
    },
    estimation::{Kalman2dF64, Kalman3dF64},
    learning::{
        AdaGradF64, AdamF64, LmsFilterF64, NlmsFilterF64, OnlineGdF64, OnlineKMeansF64,
        RlsFilterF64,
    },
    monitoring::{
        DrawdownF64, ErrorRateF64, JitterF64, RunningMaxF64, RunningMinF64, SaturationF64,
    },
    regression::{
        EwLinearRegressionF64, LinearRegressionF64, LogisticRegressionF64, PolynomialRegressionF64,
    },
    signal::{AutocorrelationF64, CrossCorrelationF64, EntropyF64, TransferEntropyF64},
    smoothing::{AsymEmaF64, EmaF64, HoltF64, Kalman1dF64, SlewF64, SpringF64},
    statistics::{
        CovarianceF64, EwmaVarF64, HarmonicMeanF64, MomentsF64, PercentileF64, WelfordF64,
    },
};

// Simple LCG for deterministic pseudo-random data in benchmarks.
// Not crypto-quality — just needs to defeat constant propagation.
struct Lcg(u64);

impl Lcg {
    fn new(seed: u64) -> Self {
        Self(seed)
    }
    fn next_f64(&mut self) -> f64 {
        self.0 = self.0.wrapping_mul(6364136223846793005).wrapping_add(1);
        // Map to [0, 100) range — realistic for most stats
        (self.0 >> 33) as f64 * (100.0 / (1u64 << 31) as f64)
    }
    fn next_feature(&mut self) -> f64 {
        self.0 = self.0.wrapping_mul(6364136223846793005).wrapping_add(1);
        // Map to [-1, 1) range — realistic for feature vectors
        ((self.0 >> 33) as f64 / (1u64 << 31) as f64) * 2.0 - 1.0
    }
}

// ============================================================
// Group 1: Core statistics
// ============================================================

fn bench_welford(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut w = WelfordF64::new();
    for _ in 0..1000 {
        let _ = w.update(rng.next_f64());
    }
    c.bench_function("WelfordF64::update", |b| {
        b.iter(|| w.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_moments(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut m = MomentsF64::new();
    for _ in 0..1000 {
        let _ = m.update(rng.next_f64());
    }
    c.bench_function("MomentsF64::update", |b| {
        b.iter(|| m.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_ema(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut e = EmaF64::builder().span(20).min_samples(1).build().unwrap();
    for _ in 0..1000 {
        let _ = e.update(rng.next_f64());
    }
    c.bench_function("EmaF64::update (span=20)", |b| {
        b.iter(|| e.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_asym_ema(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut e = AsymEmaF64::builder()
        .alpha_up(0.3)
        .alpha_down(0.1)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = e.update(rng.next_f64());
    }
    c.bench_function("AsymEmaF64::update", |b| {
        b.iter(|| e.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_holt(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut h = HoltF64::builder()
        .alpha(0.3)
        .beta(0.1)
        .min_samples(2)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = h.update(rng.next_f64());
    }
    c.bench_function("HoltF64::update", |b| {
        b.iter(|| h.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_ewma_var(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut e = EwmaVarF64::builder().alpha(0.1).build().unwrap();
    for _ in 0..1000 {
        let _ = e.update(rng.next_f64());
    }
    c.bench_function("EwmaVarF64::update", |b| {
        b.iter(|| e.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_covariance(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut cov = CovarianceF64::new();
    for _ in 0..1000 {
        let _ = cov.update(rng.next_f64(), rng.next_f64());
    }
    c.bench_function("CovarianceF64::update", |b| {
        b.iter(|| {
            cov.update(black_box(rng.next_f64()), black_box(rng.next_f64()))
                .unwrap()
        })
    });
}

fn bench_harmonic_mean(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut h = HarmonicMeanF64::new();
    for _ in 0..1000 {
        let _ = h.update(rng.next_f64() + 1.0); // must be positive
    }
    c.bench_function("HarmonicMeanF64::update", |b| {
        b.iter(|| h.update(black_box(rng.next_f64() + 1.0)).unwrap())
    });
}

fn bench_percentile_p50(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut p = PercentileF64::builder().percentile(0.5).build().unwrap();
    for _ in 0..1000 {
        let _ = p.update(rng.next_f64());
    }
    c.bench_function("PercentileF64::update (p50)", |b| {
        b.iter(|| p.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_percentile_p99(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut p = PercentileF64::builder().percentile(0.99).build().unwrap();
    for _ in 0..1000 {
        let _ = p.update(rng.next_f64());
    }
    c.bench_function("PercentileF64::update (p99)", |b| {
        b.iter(|| p.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_jitter(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut j = JitterF64::builder().alpha(0.1).build().unwrap();
    for _ in 0..1000 {
        let _ = j.update(rng.next_f64());
    }
    c.bench_function("JitterF64::update", |b| {
        b.iter(|| j.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_drawdown(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut d = DrawdownF64::new();
    for _ in 0..1000 {
        let _ = d.update(rng.next_f64());
    }
    c.bench_function("DrawdownF64::update", |b| {
        b.iter(|| d.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_running_min(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut r = RunningMinF64::new();
    for _ in 0..1000 {
        r.update(rng.next_f64());
    }
    c.bench_function("RunningMinF64::update", |b| {
        b.iter(|| r.update(black_box(rng.next_f64())))
    });
}

fn bench_running_max(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut r = RunningMaxF64::new();
    for _ in 0..1000 {
        r.update(rng.next_f64());
    }
    c.bench_function("RunningMaxF64::update", |b| {
        b.iter(|| r.update(black_box(rng.next_f64())))
    });
}

fn bench_slew(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut s = SlewF64::new(5.0).unwrap();
    for _ in 0..1000 {
        s.update(rng.next_f64());
    }
    c.bench_function("SlewF64::update", |b| {
        b.iter(|| s.update(black_box(rng.next_f64())))
    });
}

fn bench_spring(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut s = SpringF64::new(0.1).unwrap();
    for _ in 0..1000 {
        let _ = s.update(rng.next_f64(), 1.0);
    }
    c.bench_function("SpringF64::update", |b| {
        b.iter(|| s.update(black_box(rng.next_f64()), black_box(1.0)))
    });
}

fn bench_dead_band(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut d = DeadBandF64::new(1.0);
    for _ in 0..1000 {
        d.update(rng.next_f64());
    }
    c.bench_function("DeadBandF64::update", |b| {
        b.iter(|| d.update(black_box(rng.next_f64())))
    });
}

fn bench_hysteresis(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut h = HysteresisF64::new(40.0, 60.0).unwrap();
    for _ in 0..1000 {
        h.update(rng.next_f64());
    }
    c.bench_function("HysteresisF64::update", |b| {
        b.iter(|| h.update(black_box(rng.next_f64())))
    });
}

fn bench_level_crossing(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut l = LevelCrossingF64::new(50.0);
    for _ in 0..1000 {
        l.update(rng.next_f64());
    }
    c.bench_function("LevelCrossingF64::update", |b| {
        b.iter(|| l.update(black_box(rng.next_f64())))
    });
}

fn bench_first_diff(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut d = FirstDiffF64::new();
    for _ in 0..1000 {
        d.update(rng.next_f64());
    }
    c.bench_function("FirstDiffF64::update", |b| {
        b.iter(|| d.update(black_box(rng.next_f64())))
    });
}

fn bench_second_diff(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut d = SecondDiffF64::new();
    for _ in 0..1000 {
        d.update(rng.next_f64());
    }
    c.bench_function("SecondDiffF64::update", |b| {
        b.iter(|| d.update(black_box(rng.next_f64())))
    });
}

fn bench_cusum(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut cs = CusumF64::builder(50.0)
        .threshold(10.0)
        .slack(0.5)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = cs.update(rng.next_f64());
    }
    c.bench_function("CusumF64::update", |b| {
        b.iter(|| cs.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_shiryaev_roberts(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut sr = ShiryaevRobertsF64::builder()
        .pre_change_mean(50.0)
        .post_change_mean(51.0)
        .variance(1.0)
        .threshold(100.0)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = sr.update(rng.next_f64());
    }
    c.bench_function("ShiryaevRobertsF64::update", |b| {
        b.iter(|| sr.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_trend_alert(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut ta = TrendAlertF64::builder()
        .alpha(0.3)
        .beta(0.1)
        .trend_threshold(5.0)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = ta.update(rng.next_f64());
    }
    c.bench_function("TrendAlertF64::update", |b| {
        b.iter(|| ta.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_adaptive_threshold(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut at = AdaptiveThresholdF64::builder()
        .alpha(0.1)
        .z_threshold(3.0)
        .min_samples(20)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = at.update(rng.next_f64());
    }
    c.bench_function("AdaptiveThresholdF64::update", |b| {
        b.iter(|| at.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_robust_z(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut rz = RobustZScoreF64::builder()
        .alpha(0.1)
        .reject_threshold(3.0)
        .min_samples(5)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = rz.update(rng.next_f64());
    }
    c.bench_function("RobustZScoreF64::update", |b| {
        b.iter(|| rz.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_error_rate(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut er = ErrorRateF64::builder()
        .alpha(0.1)
        .threshold(0.5)
        .build()
        .unwrap();
    for _ in 0..1000 {
        er.update(rng.next_f64() > 50.0);
    }
    c.bench_function("ErrorRateF64::update", |b| {
        b.iter(|| er.update(black_box(rng.next_f64() > 50.0)))
    });
}

fn bench_saturation(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut s = SaturationF64::builder()
        .alpha(0.1)
        .threshold(80.0)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = s.update(rng.next_f64());
    }
    c.bench_function("SaturationF64::update", |b| {
        b.iter(|| s.update(black_box(rng.next_f64())).unwrap())
    });
}

// ============================================================
// Group 2: Const-generic types
// ============================================================

fn bench_autocorrelation(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    let mut ac1 = AutocorrelationF64::builder().lag(1).build().unwrap();
    for _ in 0..1000 {
        let _ = ac1.update(rng.next_f64());
    }
    c.bench_function("AutocorrelationF64::update (lag=1)", |b| {
        b.iter(|| ac1.update(black_box(rng.next_f64())).unwrap())
    });

    let mut ac10 = AutocorrelationF64::builder().lag(10).build().unwrap();
    for _ in 0..1000 {
        let _ = ac10.update(rng.next_f64());
    }
    c.bench_function("AutocorrelationF64::update (lag=10)", |b| {
        b.iter(|| ac10.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_cross_correlation(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    let mut cc1 = CrossCorrelationF64::builder().lag(1).build().unwrap();
    for _ in 0..1000 {
        let _ = cc1.update(rng.next_f64(), rng.next_f64());
    }
    c.bench_function("CrossCorrelationF64::update (lag=1)", |b| {
        b.iter(|| {
            cc1.update(black_box(rng.next_f64()), black_box(rng.next_f64()))
                .unwrap()
        })
    });

    let mut cc10 = CrossCorrelationF64::builder().lag(10).build().unwrap();
    for _ in 0..1000 {
        let _ = cc10.update(rng.next_f64(), rng.next_f64());
    }
    c.bench_function("CrossCorrelationF64::update (lag=10)", |b| {
        b.iter(|| {
            cc10.update(black_box(rng.next_f64()), black_box(rng.next_f64()))
                .unwrap()
        })
    });
}

fn bench_entropy(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    let mut e4 = EntropyF64::builder().bins(4).build().unwrap();
    for _ in 0..1000 {
        e4.update((rng.next_f64() as usize) % 4);
    }
    c.bench_function("EntropyF64::update (bins=4)", |b| {
        b.iter(|| e4.update(black_box((rng.next_f64() as usize) % 4)))
    });

    let mut e16 = EntropyF64::builder().bins(16).build().unwrap();
    for _ in 0..1000 {
        e16.update((rng.next_f64() as usize) % 16);
    }
    c.bench_function("EntropyF64::update (bins=16)", |b| {
        b.iter(|| e16.update(black_box((rng.next_f64() as usize) % 16)))
    });
}

// ============================================================
// Group 3: Runtime-configured types
// ============================================================

fn bench_linear_regression(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    let mut lr = LinearRegressionF64::builder()
        .intercept(true)
        .build()
        .unwrap();
    for i in 0..1000 {
        let _ = lr.update(i as f64, rng.next_f64());
    }
    c.bench_function("LinearRegressionF64::update (intercept)", |b| {
        let mut x = 1000.0;
        b.iter(|| {
            x += 1.0;
            lr.update(black_box(x), black_box(rng.next_f64())).unwrap()
        })
    });

    let mut lr_no = LinearRegressionF64::builder()
        .intercept(false)
        .build()
        .unwrap();
    for i in 0..1000 {
        let _ = lr_no.update(i as f64, rng.next_f64());
    }
    c.bench_function("LinearRegressionF64::update (no intercept)", |b| {
        let mut x = 1000.0;
        b.iter(|| {
            x += 1.0;
            lr_no
                .update(black_box(x), black_box(rng.next_f64()))
                .unwrap()
        })
    });
}

fn bench_ew_linear_regression(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut lr = EwLinearRegressionF64::builder()
        .alpha(0.05)
        .intercept(true)
        .build()
        .unwrap();
    for i in 0..1000 {
        let _ = lr.update(i as f64, rng.next_f64());
    }
    c.bench_function("EwLinearRegressionF64::update (alpha=0.05)", |b| {
        let mut x = 1000.0;
        b.iter(|| {
            x += 1.0;
            lr.update(black_box(x), black_box(rng.next_f64())).unwrap()
        })
    });
}

fn bench_polynomial_regression(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for degree in [2, 3, 5] {
        let mut pr = PolynomialRegressionF64::builder()
            .degree(degree)
            .intercept(true)
            .build()
            .unwrap();
        for i in 0..1000 {
            let _ = pr.update(i as f64 * 0.01, rng.next_f64());
        }
        c.bench_function(
            &format!("PolynomialRegressionF64::update (deg={degree})"),
            |b| {
                let mut x = 10.0;
                b.iter(|| {
                    x += 0.01;
                    pr.update(black_box(x), black_box(rng.next_f64())).unwrap()
                })
            },
        );
    }
}

fn bench_transfer_entropy(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for bins in [4usize, 8] {
        let mut te = TransferEntropyF64::builder()
            .bins(bins)
            .lag(1)
            .build()
            .unwrap();
        for _ in 0..1000 {
            te.update(
                (rng.next_f64() as usize) % bins,
                (rng.next_f64() as usize) % bins,
            );
        }
        c.bench_function(&format!("TransferEntropyF64::update (bins={bins})"), |b| {
            b.iter(|| {
                te.update(
                    black_box((rng.next_f64() as usize) % bins),
                    black_box((rng.next_f64() as usize) % bins),
                )
            })
        });
    }
}

// ============================================================
// Group 4: Adaptive filters
// ============================================================

fn bench_lms_filter(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for dims in [3, 10] {
        let mut f = LmsFilterF64::builder()
            .dimensions(dims)
            .learning_rate(0.01)
            .build()
            .unwrap();
        let mut features: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
        for _ in 0..1000 {
            for v in features.iter_mut() {
                *v = rng.next_feature();
            }
            let _ = f.update(&features, rng.next_f64());
        }
        c.bench_function(&format!("LmsFilterF64::update (d={dims})"), |b| {
            b.iter(|| {
                for v in features.iter_mut() {
                    *v = rng.next_feature();
                }
                f.update(black_box(&features), black_box(rng.next_f64()))
                    .unwrap()
            })
        });
    }
}

fn bench_nlms_filter(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for dims in [3, 10] {
        let mut f = NlmsFilterF64::builder()
            .dimensions(dims)
            .learning_rate(0.5)
            .build()
            .unwrap();
        let mut features: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
        for _ in 0..1000 {
            for v in features.iter_mut() {
                *v = rng.next_feature();
            }
            let _ = f.update(&features, rng.next_f64());
        }
        c.bench_function(&format!("NlmsFilterF64::update (d={dims})"), |b| {
            b.iter(|| {
                for v in features.iter_mut() {
                    *v = rng.next_feature();
                }
                f.update(black_box(&features), black_box(rng.next_f64()))
                    .unwrap()
            })
        });
    }
}

fn bench_rls_filter(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for dims in [3, 5, 10] {
        let mut f = RlsFilterF64::builder()
            .dimensions(dims)
            .forgetting_factor(0.99)
            .build()
            .unwrap();
        let mut features: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
        for _ in 0..1000 {
            for v in features.iter_mut() {
                *v = rng.next_feature();
            }
            let _ = f.update(&features, rng.next_f64());
        }
        c.bench_function(&format!("RlsFilterF64::update (d={dims})"), |b| {
            b.iter(|| {
                for v in features.iter_mut() {
                    *v = rng.next_feature();
                }
                f.update(black_box(&features), black_box(rng.next_f64()))
                    .unwrap()
            })
        });
    }
}

fn bench_logistic_regression(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let dims = 5;
    let mut lr = LogisticRegressionF64::builder()
        .dimensions(dims)
        .learning_rate(0.1)
        .build()
        .unwrap();
    let mut features: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
    for _ in 0..1000 {
        for v in features.iter_mut() {
            *v = rng.next_feature();
        }
        lr.update(&features, rng.next_f64() > 50.0);
    }
    c.bench_function("LogisticRegressionF64::update (d=5)", |b| {
        b.iter(|| {
            for v in features.iter_mut() {
                *v = rng.next_feature();
            }
            lr.update(black_box(&features), black_box(rng.next_f64() > 50.0))
        })
    });
}

fn bench_online_kmeans(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let dims = 5;
    let mut km = OnlineKMeansF64::builder()
        .dimensions(dims)
        .clusters(3)
        .learning_rate(0.1)
        .build()
        .unwrap();
    let mut features: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
    for _ in 0..1000 {
        for v in features.iter_mut() {
            *v = rng.next_feature();
        }
        km.update(&features);
    }
    c.bench_function("OnlineKMeansF64::update (k=3, d=5)", |b| {
        b.iter(|| {
            for v in features.iter_mut() {
                *v = rng.next_feature();
            }
            km.update(black_box(&features))
        })
    });
}

// ============================================================
// Group 5: State estimation
// ============================================================

fn bench_kalman1d(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut k = Kalman1dF64::builder()
        .process_noise(0.01)
        .measurement_noise(1.0)
        .build()
        .unwrap();
    for _ in 0..1000 {
        let _ = k.update(rng.next_f64());
    }
    c.bench_function("Kalman1dF64::update", |b| {
        b.iter(|| k.update(black_box(rng.next_f64())).unwrap())
    });
}

fn bench_kalman2d(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut k = Kalman2dF64::builder()
        .process_noise([[0.01, 0.0], [0.0, 0.01]])
        .measurement_noise(1.0)
        .build()
        .unwrap();
    for _ in 0..1000 {
        k.predict();
        let _ = k.update(rng.next_f64(), [1.0, 0.0]);
    }

    c.bench_function("Kalman2dF64::predict+update", |b| {
        b.iter(|| {
            k.predict();
            k.update(black_box(rng.next_f64()), black_box([1.0, 0.0]))
                .unwrap()
        })
    });

    c.bench_function("Kalman2dF64::predict_with_dynamics", |b| {
        b.iter(|| {
            k.predict_with_dynamics(black_box([[1.0, 1.0], [0.0, 1.0]]));
            k.update(black_box(rng.next_f64()), black_box([1.0, 0.0]))
                .unwrap()
        })
    });
}

fn bench_kalman3d(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut k = Kalman3dF64::builder()
        .process_noise([[0.01, 0.0, 0.0], [0.0, 0.01, 0.0], [0.0, 0.0, 0.01]])
        .measurement_noise(1.0)
        .build()
        .unwrap();
    for _ in 0..1000 {
        k.predict();
        let _ = k.update(rng.next_f64(), [1.0, 0.0, 0.0]);
    }
    c.bench_function("Kalman3dF64::predict+update", |b| {
        b.iter(|| {
            k.predict();
            k.update(black_box(rng.next_f64()), black_box([1.0, 0.0, 0.0]))
                .unwrap()
        })
    });
}

// ============================================================
// Group 6: Optimizers
// ============================================================

fn bench_online_gd(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for dims in [5, 20] {
        let mut opt = OnlineGdF64::builder()
            .dimensions(dims)
            .learning_rate(0.01)
            .build()
            .unwrap();
        let mut grad: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
        for _ in 0..1000 {
            for v in grad.iter_mut() {
                *v = rng.next_feature();
            }
            let _ = opt.step(&grad);
        }
        c.bench_function(&format!("OnlineGdF64::step (d={dims})"), |b| {
            b.iter(|| {
                for v in grad.iter_mut() {
                    *v = rng.next_feature();
                }
                opt.step(black_box(&grad)).unwrap()
            })
        });
    }
}

fn bench_adagrad(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for dims in [5, 20] {
        let mut opt = AdaGradF64::builder()
            .dimensions(dims)
            .learning_rate(0.1)
            .build()
            .unwrap();
        let mut grad: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
        for _ in 0..1000 {
            for v in grad.iter_mut() {
                *v = rng.next_feature();
            }
            let _ = opt.step(&grad);
        }
        c.bench_function(&format!("AdaGradF64::step (d={dims})"), |b| {
            b.iter(|| {
                for v in grad.iter_mut() {
                    *v = rng.next_feature();
                }
                opt.step(black_box(&grad)).unwrap()
            })
        });
    }
}

fn bench_adam(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for dims in [5, 20] {
        let mut opt = AdamF64::builder()
            .dimensions(dims)
            .learning_rate(0.001)
            .build()
            .unwrap();
        let mut grad: Vec<f64> = (0..dims).map(|_| rng.next_feature()).collect();
        for _ in 0..1000 {
            for v in grad.iter_mut() {
                *v = rng.next_feature();
            }
            let _ = opt.step(&grad);
        }
        c.bench_function(&format!("AdamF64::step (d={dims})"), |b| {
            b.iter(|| {
                for v in grad.iter_mut() {
                    *v = rng.next_feature();
                }
                opt.step(black_box(&grad)).unwrap()
            })
        });
    }
}

// ============================================================
// Group 7: Query methods (cold path, reference)
// ============================================================

fn bench_polynomial_coefficients(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for degree in [2, 3, 5] {
        let mut pr = PolynomialRegressionF64::builder()
            .degree(degree)
            .intercept(true)
            .build()
            .unwrap();
        for i in 0..1000 {
            let _ = pr.update(i as f64 * 0.01, rng.next_f64());
        }
        c.bench_function(
            &format!("PolynomialRegressionF64::coefficients (deg={degree})"),
            |b| b.iter(|| pr.coefficients()),
        );
    }
}

fn bench_transfer_entropy_query(c: &mut Criterion) {
    let mut rng = Lcg::new(42);

    for bins in [4usize, 8] {
        let mut te = TransferEntropyF64::builder()
            .bins(bins)
            .lag(1)
            .build()
            .unwrap();
        for _ in 0..1000 {
            te.update(
                (rng.next_f64() as usize) % bins,
                (rng.next_f64() as usize) % bins,
            );
        }
        c.bench_function(
            &format!("TransferEntropyF64::te_x_to_y (bins={bins})"),
            |b| b.iter(|| te.te_x_to_y()),
        );
    }
}

fn bench_entropy_query(c: &mut Criterion) {
    let mut rng = Lcg::new(42);
    let mut e = EntropyF64::builder().bins(8).build().unwrap();
    for _ in 0..1000 {
        e.update((rng.next_f64() as usize) % 8);
    }
    c.bench_function("EntropyF64::entropy (bins=8)", |b| b.iter(|| e.entropy()));
}

// ============================================================
// Criterion groups and main
// ============================================================

criterion_group!(
    core_stats,
    bench_welford,
    bench_moments,
    bench_ema,
    bench_asym_ema,
    bench_holt,
    bench_ewma_var,
    bench_covariance,
    bench_harmonic_mean,
    bench_percentile_p50,
    bench_percentile_p99,
    bench_jitter,
    bench_drawdown,
    bench_running_min,
    bench_running_max,
    bench_slew,
    bench_spring,
    bench_dead_band,
    bench_hysteresis,
    bench_level_crossing,
    bench_first_diff,
    bench_second_diff,
    bench_cusum,
    bench_shiryaev_roberts,
    bench_trend_alert,
    bench_adaptive_threshold,
    bench_robust_z,
    bench_error_rate,
    bench_saturation,
);

criterion_group!(
    const_generic,
    bench_autocorrelation,
    bench_cross_correlation,
    bench_entropy,
);

criterion_group!(
    runtime_configured,
    bench_linear_regression,
    bench_ew_linear_regression,
    bench_polynomial_regression,
    bench_transfer_entropy,
);

criterion_group!(
    adaptive_filters,
    bench_lms_filter,
    bench_nlms_filter,
    bench_rls_filter,
    bench_logistic_regression,
    bench_online_kmeans,
);

criterion_group!(
    state_estimation,
    bench_kalman1d,
    bench_kalman2d,
    bench_kalman3d,
);

criterion_group!(optimizers, bench_online_gd, bench_adagrad, bench_adam,);

criterion_group!(
    queries,
    bench_polynomial_coefficients,
    bench_transfer_entropy_query,
    bench_entropy_query,
);

criterion_main!(
    core_stats,
    const_generic,
    runtime_configured,
    adaptive_filters,
    state_estimation,
    optimizers,
    queries,
);
