//! V10 Block R — True Hawkes self-exciting process MLE estimator.
//!
//! ## Mathematical model
//!
//! Univariate Hawkes process with single-exponential excitation kernel:
//!
//! ```text
//!   λ(t) = μ + Σ_{t_i < t} α · exp(-β · (t − t_i))
//! ```
//!
//! Three parameters:
//! - `μ ≥ 0` — baseline intensity (events/sec)
//! - `α ≥ 0` — excitation strength per event
//! - `β > 0` — decay rate (1/sec); higher β = shorter memory
//!
//! **Branching ratio** `n = α/β` measures self-excitation criticality:
//! - `n = 0`: pure Poisson, no excitation
//! - `0 < n < 1`: stationary self-exciting
//! - `n → 1`: near-critical, long-memory
//! - `n ≥ 1`: non-stationary (event rate explodes)
//!
//! Per Bacry et al. (2015), Filimonov & Sornette (2012, 2015), markets near
//! `n* ≈ 0.7-0.9` are observationally "near criticality" — high reflexivity,
//! short bursts of correlated activity, characteristic of HFT regimes.
//!
//! ## Estimator
//!
//! Maximum likelihood via **3D grid search** over (μ, n, β) where α = n·β.
//! The log-likelihood
//!
//! ```text
//!   ln L = Σ_i ln(λ(t_i)) − ∫_0^T λ(t) dt
//! ```
//!
//! has a closed-form **recursive** structure for the exponential kernel:
//!
//! ```text
//!   A_i = (A_{i-1} + 1) · exp(-β · Δt_i)
//!   λ(t_i) = μ + α · A_i
//!   ∫_0^T λ(t) dt = μ·T + (α/β) · Σ_j (1 − exp(-β · (T − t_j)))
//! ```
//!
//! Total cost: `|grid| × N` per refit. With 7×7×5 = 245 grid points × 500
//! arrivals = ~120k ops per refit. Refit every 100 arrivals.
//!
//! ## Output features (4 dims)
//!
//! 1. `mu` — baseline intensity μ
//! 2. `alpha` — excitation strength α
//! 3. `beta` — decay rate β
//! 4. `branching_ratio` — n = α/β (the "reflexivity" scalar)

use std::collections::VecDeque;

/// Grid for baseline intensity μ search (events/sec).
pub const HAWKES_GRID_MU: &[f64] = &[0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0];
/// Grid for branching ratio n = α/β search (must be in [0, 1)).
pub const HAWKES_GRID_N: &[f64] = &[0.0, 0.1, 0.3, 0.5, 0.7, 0.85, 0.95];
/// Grid for decay rate β search (1/sec).
pub const HAWKES_GRID_BETA: &[f64] = &[0.1, 0.5, 1.0, 2.0, 5.0];

/// V10 Block R — Hawkes MLE estimator (online + periodic refit).
#[derive(Debug, Clone)]
pub struct HawkesEstimator {
    max_history: usize,
    refit_interval: usize,
    first_arrival_ts_ns: Option<u64>,
    /// Relative event times in seconds (from `first_arrival_ts_ns`).
    times: VecDeque<f64>,
    cached_mu: f64,
    cached_alpha: f64,
    cached_beta: f64,
    obs_since_refit: usize,
}

impl HawkesEstimator {
    pub fn new(max_history: usize, refit_interval: usize) -> Self {
        Self {
            max_history,
            refit_interval: refit_interval.max(1),
            first_arrival_ts_ns: None,
            times: VecDeque::with_capacity(max_history),
            cached_mu: 0.0,
            cached_alpha: 0.0,
            cached_beta: 1.0,
            obs_since_refit: 0,
        }
    }

    /// V10 defaults: 500 events in history, refit every 100 events.
    pub fn with_defaults() -> Self {
        Self::new(500, 100)
    }

    /// Register a new trade arrival. Returns
    /// `[mu, alpha, beta, branching_ratio]` from the most recent fit.
    pub fn update(&mut self, arrival_ts_ns: u64) -> [f64; 4] {
        let first = *self.first_arrival_ts_ns.get_or_insert(arrival_ts_ns);
        let t = (arrival_ts_ns.saturating_sub(first)) as f64 / 1.0e9;
        self.times.push_back(t);
        if self.times.len() > self.max_history {
            self.times.pop_front();
        }
        self.obs_since_refit += 1;
        if self.obs_since_refit >= self.refit_interval && self.times.len() >= 30 {
            self.refit();
            self.obs_since_refit = 0;
        }
        self.features()
    }

    /// Force a refit immediately (e.g. before producing v10 fxcache).
    pub fn force_refit(&mut self) {
        if self.times.len() >= 30 {
            self.refit();
            self.obs_since_refit = 0;
        }
    }

    /// Return the current cached features without refit.
    pub fn features(&self) -> [f64; 4] {
        let n = if self.cached_beta > 1e-12 {
            self.cached_alpha / self.cached_beta
        } else {
            0.0
        };
        [self.cached_mu, self.cached_alpha, self.cached_beta, n]
    }

    /// Run 3D grid-search MLE over the cached times.
    fn refit(&mut self) {
        // Time origin: shift to make times[0] = 0 for likelihood numerics
        let t0 = self.times[0];
        let times: Vec<f64> = self.times.iter().map(|t| t - t0).collect();
        let t_end = *times.last().expect("non-empty");
        if t_end <= 0.0 {
            return; // degenerate
        }

        let mut best_ll = f64::NEG_INFINITY;
        let mut best = (self.cached_mu, self.cached_alpha, self.cached_beta);

        for &mu in HAWKES_GRID_MU {
            for &beta in HAWKES_GRID_BETA {
                for &n in HAWKES_GRID_N {
                    let alpha = n * beta;
                    if alpha >= beta && n >= 1.0 {
                        continue; // non-stationary
                    }
                    let ll = hawkes_log_likelihood(&times, t_end, mu, alpha, beta);
                    if ll.is_finite() && ll > best_ll {
                        best_ll = ll;
                        best = (mu, alpha, beta);
                    }
                }
            }
        }

        self.cached_mu = best.0;
        self.cached_alpha = best.1;
        self.cached_beta = best.2;
    }

    pub fn reset(&mut self) {
        self.times.clear();
        self.first_arrival_ts_ns = None;
        self.cached_mu = 0.0;
        self.cached_alpha = 0.0;
        self.cached_beta = 1.0;
        self.obs_since_refit = 0;
    }
}

impl Default for HawkesEstimator {
    fn default() -> Self {
        Self::with_defaults()
    }
}

/// Log-likelihood of exponential Hawkes via recursive `A_i` accumulation.
///
/// Returns `NEG_INFINITY` for invalid parameters (μ ≤ 0, β ≤ 0, or any
/// negative intensity).
pub fn hawkes_log_likelihood(times: &[f64], t_end: f64, mu: f64, alpha: f64, beta: f64) -> f64 {
    let n_events = times.len();
    if n_events == 0 || mu <= 0.0 || beta <= 0.0 || alpha < 0.0 {
        return f64::NEG_INFINITY;
    }

    let mut log_l = 0.0_f64;
    let mut excitation = 0.0_f64; // A_i recursion: starts at 0 before first event

    // First event: λ(t_0) = μ (no prior events to excite)
    log_l += mu.ln();

    let mut prev_t = times[0];
    for i in 1..n_events {
        let dt = times[i] - prev_t;
        if dt < 0.0 {
            return f64::NEG_INFINITY;
        }
        // Recursive update: A_i = (A_{i-1} + 1) · exp(-β · Δt)
        excitation = (excitation + 1.0) * (-beta * dt).exp();
        let lambda_i = mu + alpha * excitation;
        if lambda_i <= 0.0 {
            return f64::NEG_INFINITY;
        }
        log_l += lambda_i.ln();
        prev_t = times[i];
    }

    // Compensator integral: ∫_0^T λ(t) dt = μ·T + (α/β) · Σ_j (1 − exp(-β·(T−t_j)))
    let mu_term = mu * t_end;
    let alpha_over_beta = alpha / beta;
    let excitation_sum: f64 = times
        .iter()
        .map(|&tj| 1.0 - (-beta * (t_end - tj)).exp())
        .sum();
    let comp = mu_term + alpha_over_beta * excitation_sum;

    log_l - comp
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_cold_start_features_zero() {
        let h = HawkesEstimator::with_defaults();
        let features = h.features();
        assert_eq!(features[0], 0.0, "cold start: μ = 0");
        assert_eq!(features[1], 0.0, "cold start: α = 0");
        assert_eq!(features[3], 0.0, "cold start: n = 0");
    }

    #[test]
    fn test_log_likelihood_zero_events_is_neg_infinity() {
        let ll = hawkes_log_likelihood(&[], 1.0, 1.0, 0.5, 1.0);
        assert_eq!(ll, f64::NEG_INFINITY);
    }

    #[test]
    fn test_log_likelihood_invalid_params_neg_infinity() {
        assert_eq!(
            hawkes_log_likelihood(&[0.0, 0.5, 1.0], 1.0, -0.5, 0.5, 1.0),
            f64::NEG_INFINITY,
            "negative μ → -inf"
        );
        assert_eq!(
            hawkes_log_likelihood(&[0.0, 0.5, 1.0], 1.0, 1.0, 0.5, -0.5),
            f64::NEG_INFINITY,
            "negative β → -inf"
        );
    }

    #[test]
    fn test_log_likelihood_poisson_prefers_zero_alpha() {
        // Generate regular Poisson-like times. MLE should prefer α ≈ 0
        // (= no self-excitation; just Poisson).
        // With 100 events spaced uniformly Δt = 1, the rate is 1.0 events/sec.
        let times: Vec<f64> = (0..100).map(|i| i as f64).collect();
        let t_end = 99.0;

        let ll_no_excite = hawkes_log_likelihood(&times, t_end, 1.0, 0.0, 1.0);
        let ll_high_excite = hawkes_log_likelihood(&times, t_end, 1.0, 0.5, 1.0);

        // Regular spacing has no self-excitation signal → low-α should win
        assert!(
            ll_no_excite > ll_high_excite,
            "regular spacing should prefer α=0 over α=0.5; got ll_no={ll_no_excite} ll_high={ll_high_excite}"
        );
    }

    #[test]
    fn test_log_likelihood_clustered_prefers_high_alpha() {
        // Heavily clustered: 50 events in [0, 1] then 50 events in [10, 11].
        // High α → high λ during clusters → high likelihood per event.
        let mut times = Vec::with_capacity(100);
        for i in 0..50 {
            times.push((i as f64) * 0.02); // [0, 1)
        }
        for i in 0..50 {
            times.push(10.0 + (i as f64) * 0.02); // [10, 11)
        }
        let t_end = 11.0;

        let ll_no_excite = hawkes_log_likelihood(&times, t_end, 5.0, 0.0, 1.0);
        let ll_high_excite = hawkes_log_likelihood(&times, t_end, 0.5, 0.8, 1.0);

        // Clustered data should prefer high α (μ low, α high, branching ratio near 1)
        assert!(
            ll_high_excite > ll_no_excite,
            "clustered events should prefer α>0; got ll_no={ll_no_excite} ll_high={ll_high_excite}"
        );
    }

    #[test]
    fn test_estimator_refit_on_regular_poisson_picks_low_branching() {
        // Feed 200 regular events at 1Hz. MLE should pick small branching ratio.
        let mut h = HawkesEstimator::new(500, 50);
        let start_ns = 1_700_000_000_000_000_000_u64;
        for i in 0..200 {
            h.update(start_ns + (i as u64) * 1_000_000_000);
        }
        let [_, _, _, n] = h.features();
        assert!(n < 0.5, "regular Poisson → low branching ratio, got {n}");
    }

    #[test]
    fn test_estimator_refit_on_clustered_picks_high_branching() {
        // Heavily clustered events: 100 events in 1s, then nothing for 10s, repeat
        let mut h = HawkesEstimator::new(500, 50);
        let mut t = 1_700_000_000_000_000_000_u64;
        for _cluster in 0..3 {
            for _ in 0..50 {
                h.update(t);
                t += 20_000_000; // 20ms apart
            }
            t += 5_000_000_000; // 5s gap
        }
        let [_, _, _, n] = h.features();
        assert!(n >= 0.3, "clustered events → branching ratio ≥ 0.3, got {n}");
    }

    #[test]
    fn test_reset_clears_state() {
        let mut h = HawkesEstimator::new(500, 50);
        let start_ns = 1_700_000_000_000_000_000_u64;
        for i in 0..100 {
            h.update(start_ns + (i as u64) * 1_000_000_000);
        }
        h.reset();
        let features = h.features();
        assert_eq!(features[0], 0.0);
        assert_eq!(features[1], 0.0);
        assert_eq!(features[3], 0.0);
    }

    #[test]
    fn test_force_refit_does_nothing_without_enough_data() {
        let mut h = HawkesEstimator::new(500, 50);
        // Fewer than 30 arrivals
        for i in 0..10 {
            h.update(1_700_000_000_000_000_000_u64 + (i as u64) * 1_000_000_000);
        }
        let before = h.features();
        h.force_refit();
        let after = h.features();
        assert_eq!(before, after, "<30 events → force_refit is no-op");
    }
}
