atelier_quant 0.0.12

#[cfg(test)]
mod tests {
    use atelier_quant::poisson::estimation;
    use atelier_quant::poisson::PoissonProcess;

    // ── validate_events ──────────────────────────────────────────────

    #[test]
    fn test_validate_events_empty() {
        assert!(estimation::validate_events(&[]).is_err());
    }

    #[test]
    fn test_validate_events_single() {
        assert!(estimation::validate_events(&[1.0]).is_err());
    }

    #[test]
    fn test_validate_events_nan() {
        assert!(estimation::validate_events(&[1.0, f64::NAN, 3.0]).is_err());
    }

    #[test]
    fn test_validate_events_non_increasing() {
        assert!(estimation::validate_events(&[1.0, 3.0, 2.0]).is_err());
    }

    #[test]
    fn test_validate_events_ok() {
        assert!(estimation::validate_events(&[0.1, 0.5, 1.0, 2.0]).is_ok());
    }

    // ── log_likelihood ───────────────────────────────────────────────

    /// Hand-computed: events = [0, 1, 2, 3], λ = 1.0
    /// T = 3, n_gaps = 3
    /// ℓ = 3·ln(1) − 1·3 = 0 − 3 = −3.0
    #[test]
    fn test_log_likelihood_known_value() {
        let events = vec![0.0, 1.0, 2.0, 3.0];
        let ll = estimation::log_likelihood(1.0, &events);
        assert!(
            (ll - (-3.0)).abs() < 1e-12,
            "Expected -3.0, got {}",
            ll
        );
    }

    #[test]
    fn test_log_likelihood_negative_lambda() {
        let events = vec![0.0, 1.0, 2.0];
        assert_eq!(estimation::log_likelihood(-1.0, &events), f64::NEG_INFINITY);
    }

    #[test]
    fn test_log_likelihood_too_few_events() {
        assert_eq!(estimation::log_likelihood(1.0, &[0.5]), f64::NEG_INFINITY);
    }

    /// Verify that the MLE λ̂ maximises the log-likelihood by checking
    /// that ℓ(λ̂) > ℓ(λ̂ ± δ) for small perturbations.
    #[test]
    fn test_log_likelihood_maximised_at_mle() {
        let events = vec![0.0, 0.5, 1.2, 1.8, 3.0, 4.5, 5.0];
        let config = estimation::PoissonEstimationConfig;
        let result = estimation::estimate_poisson_mle(&events, &config).unwrap();

        let ll_mle = result.log_likelihood;
        let delta = 0.01;

        let ll_above = estimation::log_likelihood(result.lambda + delta, &events);
        let ll_below = estimation::log_likelihood(result.lambda - delta, &events);

        assert!(
            ll_mle >= ll_above,
            "ℓ(λ̂) = {} should be >= ℓ(λ̂+δ) = {}",
            ll_mle,
            ll_above
        );
        assert!(
            ll_mle >= ll_below,
            "ℓ(λ̂) = {} should be >= ℓ(λ̂-δ) = {}",
            ll_mle,
            ll_below
        );
    }

    // ── estimate_poisson_mle ─────────────────────────────────────────

    /// 10 equispaced events on [0, 10]: λ̂ = 9/10 = 0.9 exactly.
    #[test]
    fn test_mle_closed_form_exact() {
        let events: Vec<f64> = (0..10).map(|i| i as f64).collect(); // [0,1,...,9]
        let config = estimation::PoissonEstimationConfig;
        let result = estimation::estimate_poisson_mle(&events, &config).unwrap();

        // T = 9, n_gaps = 9, λ̂ = 9/9 = 1.0
        assert!(
            (result.lambda - 1.0).abs() < 1e-12,
            "Expected λ̂ = 1.0, got {}",
            result.lambda
        );
        assert_eq!(result.iterations, 1);
        assert!(result.converged);
    }

    /// Generate 500 events from Poisson(5.0), fit, recover λ̂ ≈ 5.0.
    #[test]
    fn test_mle_roundtrip_recovery() {
        let true_lambda = 5.0;
        let pp = PoissonProcess::new(true_lambda).unwrap();
        let events = pp.generate_values(0.0, 500);

        let config = estimation::PoissonEstimationConfig;
        let result = estimation::estimate_poisson_mle(&events, &config).unwrap();

        let rel_err = ((result.lambda - true_lambda) / true_lambda).abs();
        assert!(
            rel_err < 0.30,
            "λ̂ = {:.4} should be near {:.1} (rel_err = {:.2})",
            result.lambda,
            true_lambda,
            rel_err
        );
    }

    // ── AIC / BIC ────────────────────────────────────────────────────

    /// Verify AIC = 2k − 2ℓ and BIC = k·ln(n) − 2ℓ with k = 1.
    #[test]
    fn test_aic_bic_formula() {
        let events = vec![0.0, 1.0, 2.0, 3.0, 4.0];
        let config = estimation::PoissonEstimationConfig;
        let result = estimation::estimate_poisson_mle(&events, &config).unwrap();

        let k = 1.0_f64;
        let n_f = events.len() as f64;
        let expected_aic = 2.0 * k - 2.0 * result.log_likelihood;
        let expected_bic = k * n_f.ln() - 2.0 * result.log_likelihood;

        assert!(
            (result.aic - expected_aic).abs() < 1e-12,
            "AIC mismatch: {} vs {}",
            result.aic,
            expected_aic
        );
        assert!(
            (result.bic - expected_bic).abs() < 1e-12,
            "BIC mismatch: {} vs {}",
            result.bic,
            expected_bic
        );
    }

    // ── compensator ──────────────────────────────────────────────────

    /// Λ(t) = λ·(t − t₁) is linear.
    #[test]
    fn test_compensator_linear() {
        let events = vec![0.0, 1.0, 2.0, 3.0];
        let lambda = 2.5;

        // Λ(0) = 0
        assert!((estimation::compensator(lambda, &events, 0.0)).abs() < 1e-12);

        // Λ(2) = 2.5 * 2 = 5.0
        assert!((estimation::compensator(lambda, &events, 2.0) - 5.0).abs() < 1e-12);

        // Λ(3) = 2.5 * 3 = 7.5
        assert!((estimation::compensator(lambda, &events, 3.0) - 7.5).abs() < 1e-12);
    }

    #[test]
    fn test_compensator_monotonic() {
        let events = vec![0.0, 1.0, 2.0, 3.0];
        let lambda = 1.5;

        let mut prev = 0.0;
        for t in [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0] {
            let c = estimation::compensator(lambda, &events, t);
            assert!(c >= prev, "Compensator must be monotonically increasing");
            prev = c;
        }
    }

    // ── time_rescaling_residuals ─────────────────────────────────────

    #[test]
    fn test_time_rescaling_residuals_positive() {
        let pp = PoissonProcess::new(3.0).unwrap();
        let events = pp.generate_values(0.0, 100);

        let residuals = estimation::time_rescaling_residuals(3.0, &events);
        assert_eq!(residuals.len(), 99);

        for &r in &residuals {
            assert!(r > 0.0, "Residuals must be positive, got {}", r);
        }
    }

    /// Under the true λ, the mean of the time-rescaling residuals should
    /// be approximately 1.0 (they are i.i.d. Exp(1)).
    #[test]
    fn test_time_rescaling_residuals_mean_near_one() {
        let true_lambda = 4.0;
        let pp = PoissonProcess::new(true_lambda).unwrap();
        let events = pp.generate_values(0.0, 500);

        let residuals = estimation::time_rescaling_residuals(true_lambda, &events);
        let mean = residuals.iter().sum::<f64>() / residuals.len() as f64;

        assert!(
            (mean - 1.0).abs() < 0.50,
            "Mean residual {:.4} should be near 1.0",
            mean
        );
    }

    /// Exact check on uniform events: events = [0, 1, 2, 3], λ = 2.0
    /// Residuals: [2·1, 2·1, 2·1] = [2, 2, 2]
    #[test]
    fn test_time_rescaling_residuals_exact() {
        let events = vec![0.0, 1.0, 2.0, 3.0];
        let residuals = estimation::time_rescaling_residuals(2.0, &events);

        assert_eq!(residuals.len(), 3);
        for &r in &residuals {
            assert!(
                (r - 2.0).abs() < 1e-12,
                "Expected residual 2.0, got {}",
                r
            );
        }
    }
}