//! Bayesian Online Changepoint Detection (BOCD)
//!
//! This module implements the Bayesian Online Changepoint Detection algorithm
//! for probabilistic regime change detection in financial time series.
//!
//! # Algorithm Overview
//!
//! BOCD maintains a distribution over the current "run length" (time since last changepoint)
//! and updates it online as new data arrives. The algorithm computes:
//!
//! 1. **Run-length distribution**: P(rₜ|x₁:ₜ) - probability that the current run length is r
//! 2. **Hazard function**: H(r) = 1/λ - probability of changepoint given run length r
//! 3. **Predictive probability**: P(xₜ|x₁:ₜ₋₁, rₜ) - likelihood of observation given run length
//!
//! # Mathematical Formulation
//!
//! The core update equations are:
//!
//! ```text
//! P(rₜ|x₁:ₜ) ∝ P(xₜ|rₜ, x₁:ₜ₋₁) × [
//!     P(rₜ₋₁ = rₜ - 1|x₁:ₜ₋₁) × (1 - H(rₜ-1))  if rₜ > 0
//!     Σᵣ P(rₜ₋₁ = r|x₁:ₜ₋₁) × H(r)              if rₜ = 0
//! ]
//! ```
//!
//! Where:
//! - H(r) = 1/λ is the constant hazard function (λ = expected run length)
//! - P(xₜ|rₜ, x₁:ₜ₋₁) is computed using a conjugate Gaussian model
//!
//! # Usage
//!
//! ```rust
//! use ml::regime::bayesian_changepoint::{BayesianChangepointDetector, ChangepointInfo};
//!
//! // Create detector with hazard rate λ=100 (expect changepoint every 100 bars)
//! let mut detector = BayesianChangepointDetector::new(100.0, 0.3, 200);
//!
//! // Update with new values
//! for price in &[100.0, 101.0, 102.0, 150.0] {  // Jump at 150.0
//!     if let Some(info) = detector.update(*price) {
//!         println!("Changepoint detected! Probability: {}", info.probability);
//!     }
//! }
//!
//! // Query current state
//! let prob = detector.get_changepoint_probability();
//! let run_length = detector.get_expected_run_length();
//! ```
//!
//! # Performance
//!
//! - Target: <150μs per update (Bayesian computation intensive)
//! - Memory: O(max_run_length) for probability distribution
//! - Online: Constant time per update (no recomputation of history)
//!
//! # References
//!
//! - Adams & MacKay (2007): "Bayesian Online Changepoint Detection"
//! - Used in: Regime detection, structural break identification, adaptive strategies

use serde::{Deserialize, Serialize};

/// Information about a detected changepoint
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ChangepointInfo {
    /// Probability of changepoint (0.0 to 1.0)
    pub probability: f64,

    /// Expected run length since last changepoint
    pub expected_run_length: f64,

    /// Maximum a posteriori (MAP) run length
    pub map_run_length: usize,

    /// Current observation value
    pub value: f64,

    /// Time index of detection
    pub time_index: usize,
}

/// Bayesian Online Changepoint Detector
///
/// Maintains a distribution over the current run length (time since last changepoint)
/// and updates it online using Bayesian inference.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct BayesianChangepointDetector {
    /// Hazard rate parameter (λ): Expected run length = 1/hazard_rate
    /// Higher λ → more frequent changepoints
    /// Lower λ → less frequent changepoints
    hazard_rate: f64,

    /// Probability threshold for changepoint detection (0.0 to 1.0)
    /// Typical values: 0.2-0.4
    changepoint_prob_threshold: f64,

    /// Maximum run length to track (truncation for computational efficiency)
    max_run_length: usize,

    /// Current run-length probability distribution P(rₜ|x₁:ₜ)
    /// Index r represents run length r (time since last changepoint)
    run_length_probs: Vec<f64>,

    /// Sufficient statistics for Gaussian model (online updates)
    /// mean`[r]` = mean of observations for run length r
    means: Vec<f64>,

    /// Sufficient statistics for Gaussian model (online updates)
    /// variance`[r]` = variance of observations for run length r
    variances: Vec<f64>,

    /// Number of observations for each run length
    counts: Vec<f64>,

    /// Current time index (number of observations processed)
    time_index: usize,

    /// Prior hyperparameters for Gaussian model
    /// μ₀: Prior mean
    prior_mean: f64,

    /// Prior hyperparameters for Gaussian model
    /// κ₀: Prior precision (pseudo-count)
    prior_precision: f64,

    /// Prior hyperparameters for Gaussian model
    /// α₀: Prior degrees of freedom
    prior_alpha: f64,

    /// Prior hyperparameters for Gaussian model
    /// β₀: Prior variance scale
    prior_beta: f64,
}

impl BayesianChangepointDetector {
    /// Create a new Bayesian changepoint detector
    ///
    /// # Arguments
    ///
    /// * `hazard_rate` - Expected run length = 1/hazard_rate (e.g., 100.0 → expect changepoint every 100 bars)
    /// * `threshold` - Probability threshold for changepoint detection (0.0 to 1.0, typical: 0.2-0.4)
    /// * `max_run_length` - Maximum run length to track (truncation, typical: 200-500)
    ///
    /// # Returns
    ///
    /// New detector instance with default Gaussian priors
    ///
    /// # Examples
    ///
    /// ```rust
    /// use ml::regime::bayesian_changepoint::BayesianChangepointDetector;
    ///
    /// // Expect changepoint every 100 bars, detect at 30% probability, track up to 200 bars
    /// let detector = BayesianChangepointDetector::new(100.0, 0.3, 200);
    /// ```
    pub fn new(hazard_rate: f64, threshold: f64, max_run_length: usize) -> Self {
        // Initialize run-length distribution: P(r₀ = 0) = 1.0 (start with run length 0)
        let mut run_length_probs = vec![0.0; max_run_length + 1];
        run_length_probs[0] = 1.0;

        // Non-informative priors for Gaussian model
        // prior_beta = 1e4 gives prior variance = β/α = 1e4, broad enough for any
        // data scale. This prevents the changepoint hypothesis from being starved
        // by a narrow prior when the first observations set the mean far from 0.
        let prior_mean = 0.0; // μ₀: No prior knowledge of mean
        let prior_precision = 0.01; // κ₀: Low confidence in prior mean
        let prior_alpha = 1.0; // α₀: Minimal degrees of freedom
        let prior_beta = 1e4; // β₀: Broad variance scale (supports any data scale)

        Self {
            hazard_rate,
            changepoint_prob_threshold: threshold,
            max_run_length,
            run_length_probs,
            means: vec![prior_mean; max_run_length + 1],
            variances: vec![1.0; max_run_length + 1],
            counts: vec![0.0; max_run_length + 1],
            time_index: 0,
            prior_mean,
            prior_precision,
            prior_alpha,
            prior_beta,
        }
    }

    /// Update detector with a new observation
    ///
    /// # Arguments
    ///
    /// * `value` - New observation value
    ///
    /// # Returns
    ///
    /// Some(ChangepointInfo) if changepoint detected (P(r=0) > threshold), None otherwise
    ///
    /// # Algorithm
    ///
    /// 1. Compute predictive probabilities P(xₜ|rₜ, x₁:ₜ₋₁) for all run lengths
    /// 2. Update run-length distribution using Bayes' rule and hazard function
    /// 3. Update sufficient statistics for Gaussian model
    /// 4. Check if P(r=0) exceeds threshold
    ///
    /// # Examples
    ///
    /// ```rust
    /// use ml::regime::bayesian_changepoint::BayesianChangepointDetector;
    ///
    /// let mut detector = BayesianChangepointDetector::new(100.0, 0.3, 200);
    ///
    /// // Stable regime
    /// for i in 0..50 {
    ///     detector.update(100.0 + (i as f64 * 0.1));
    /// }
    ///
    /// // Regime change (sudden jump)
    /// if let Some(info) = detector.update(150.0) {
    ///     println!("Changepoint detected at probability {}", info.probability);
    /// }
    /// ```
    pub fn update(&mut self, value: f64) -> Option<ChangepointInfo> {
        self.time_index += 1;

        // Step 1: Compute predictive probabilities P(xₜ|rₜ, x₁:ₜ₋₁) for all run lengths
        let mut predictive_probs = vec![0.0; self.max_run_length + 1];
        for r in 0..=self.max_run_length {
            if self.run_length_probs[r] > 1e-10 {
                // Skip negligible probabilities
                predictive_probs[r] = self.compute_predictive_probability(value, r);
            }
        }

        // Step 2: Update run-length distribution using Bayes' rule and hazard function
        let mut new_probs = vec![0.0; self.max_run_length + 1];

        // Growth probabilities: P(rₜ = r+1|x₁:ₜ) from P(rₜ₋₁ = r|x₁:ₜ₋₁)
        for r in 0..self.max_run_length {
            if self.run_length_probs[r] > 1e-10 {
                let survival_prob = 1.0 - self.hazard_function(r);
                new_probs[r + 1] += self.run_length_probs[r] * predictive_probs[r] * survival_prob;
            }
        }

        // Changepoint probability: P(rₜ = 0|x₁:ₜ) = Σᵣ P(rₜ₋₁ = r|x₁:ₜ₋₁) × H(r) × P(xₜ|rₜ=0)
        let mut changepoint_prob = 0.0;
        for r in 0..=self.max_run_length {
            if self.run_length_probs[r] > 1e-10 {
                changepoint_prob += self.run_length_probs[r] * self.hazard_function(r);
            }
        }
        new_probs[0] = changepoint_prob * predictive_probs[0];

        // Step 3: Normalize probabilities
        let total: f64 = new_probs.iter().sum();
        if total > 1e-10 {
            for p in &mut new_probs {
                *p /= total;
            }
        } else {
            // Numerical underflow: Reset to initial state
            new_probs = vec![0.0; self.max_run_length + 1];
            new_probs[0] = 1.0;
        }

        // Step 3b: If a changepoint is detected, reset sufficient statistics
        // so the new regime starts with fresh priors. Without this, the shared
        // Welford statistics carry contamination from the old regime, preventing
        // the algorithm from properly tracking subsequent regime changes.
        if new_probs[0] > self.changepoint_prob_threshold {
            self.counts = vec![0.0; self.max_run_length + 1];
            self.means = vec![self.prior_mean; self.max_run_length + 1];
            self.variances = vec![1.0; self.max_run_length + 1];
        }

        // Step 4: Update sufficient statistics for Gaussian model
        self.update_sufficient_statistics(value, &new_probs);

        // Update run-length distribution
        self.run_length_probs = new_probs;

        // Step 5: Check for changepoint detection
        let cp_prob = self.run_length_probs[0];
        (cp_prob > self.changepoint_prob_threshold).then(|| ChangepointInfo {
            probability: cp_prob,
            expected_run_length: self.get_expected_run_length(),
            map_run_length: self.get_map_run_length(),
            value,
            time_index: self.time_index,
        })
    }

    /// Get current changepoint probability P(r=0|x₁:ₜ)
    ///
    /// # Returns
    ///
    /// Probability that a changepoint just occurred (0.0 to 1.0)
    pub fn get_changepoint_probability(&self) -> f64 {
        self.run_length_probs[0]
    }

    /// Get expected run length E[r|x₁:ₜ] = Σᵣ r × P(r|x₁:ₜ)
    ///
    /// # Returns
    ///
    /// Expected number of bars since last changepoint
    pub fn get_expected_run_length(&self) -> f64 {
        self.run_length_probs
            .iter()
            .enumerate()
            .map(|(r, &prob)| r as f64 * prob)
            .sum()
    }

    /// Get maximum a posteriori (MAP) run length
    ///
    /// # Returns
    ///
    /// Most likely run length (arg max P(r|x₁:ₜ))
    pub fn get_map_run_length(&self) -> usize {
        self.run_length_probs
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
            .map(|(r, _)| r)
            .unwrap_or(0)
    }

    /// Compute hazard function H(r) = 1/λ (constant hazard)
    ///
    /// # Arguments
    ///
    /// * `run_length` - Current run length r
    ///
    /// # Returns
    ///
    /// Probability of changepoint given run length r
    fn hazard_function(&self, _run_length: usize) -> f64 {
        // Constant hazard: H(r) = 1/λ for all r
        1.0 / self.hazard_rate
    }

    /// Compute predictive probability P(xₜ|rₜ, x₁:ₜ₋₁) using Student's t-distribution
    ///
    /// # Arguments
    ///
    /// * `value` - Observation value xₜ
    /// * `run_length` - Run length r
    ///
    /// # Returns
    ///
    /// Predictive probability (likelihood of observation)
    ///
    /// # Algorithm
    ///
    /// Uses conjugate Gaussian model with Normal-Inverse-Gamma priors.
    /// Predictive distribution is Student's t with parameters:
    /// - Location: μᵣ (posterior mean)
    /// - Scale: √(βᵣ(κᵣ+1)/(αᵣκᵣ)) (posterior variance)
    /// - Degrees of freedom: 2αᵣ
    fn compute_predictive_probability(&self, value: f64, run_length: usize) -> f64 {
        // r=0 (changepoint hypothesis): always use the prior predictive distribution.
        // In Adams & MacKay (2007), a changepoint means a NEW regime starts with
        // no evidence — the prior parameters describe the new regime's distribution.
        // Without this, the shared sufficient statistics at r=0 get contaminated
        // by observations from the current regime, preventing detection.
        if run_length == 0 {
            let variance = self.prior_beta / self.prior_alpha;
            return self.gaussian_pdf(value, self.prior_mean, variance);
        }

        let n = self.counts[run_length];

        if n < 1e-10 {
            let variance = self.prior_beta / self.prior_alpha;
            return self.gaussian_pdf(value, self.prior_mean, variance);
        }

        // Compute posterior parameters
        let mean = self.means[run_length];
        let variance = self.variances[run_length];

        // Posterior precision: κᵣ = κ₀ + n (accumulates confidence with observations)
        let location = mean;
        let kappa_r = self.prior_precision + n;
        let scale_squared = variance * (kappa_r + 1.0) / kappa_r;

        // Cap at 1.0: gaussian_pdf for narrow distributions can exceed 1.0 at the
        // mode (it's a density, not a probability). In BOCD, growth factor =
        // pred × (1−H). When pred > 1, growth > 1−H per step, causing exponential
        // amplification of hazard-redistributed probability mass across ALL run
        // lengths, drowning the changepoint signal. Capping at 1.0 ensures growth
        // ≤ (1−H) < 1, so hazard-redistributed mass decays geometrically.
        let pred = if n > 10.0 {
            self.gaussian_pdf(value, location, scale_squared)
        } else {
            let df = 2.0 * self.prior_alpha + n;
            self.student_t_pdf(value, location, scale_squared, df)
        };
        pred.min(1.0)
    }

    /// Gaussian PDF: N(x|μ, σ²)
    ///
    /// Uses 1e-300 floor (not 1e-10) to preserve relative probability ordering
    /// between run lengths. A flat 1e-10 floor masks the discrimination between
    /// established regimes (very low prob for outliers) and the prior (moderate
    /// prob for any value).
    fn gaussian_pdf(&self, x: f64, mean: f64, variance: f64) -> f64 {
        if variance < 1e-10 {
            return if (x - mean).abs() < 1e-6 { 1.0 } else { 1e-300 };
        }

        let diff = x - mean;
        let exponent = -0.5 * diff * diff / variance;
        let normalization = 1.0 / (2.0 * std::f64::consts::PI * variance).sqrt();

        (normalization * exponent.exp()).max(1e-300)
    }

    /// Student's t-distribution PDF (approximation)
    fn student_t_pdf(&self, x: f64, location: f64, scale: f64, df: f64) -> f64 {
        if scale < 1e-10 {
            return if (x - location).abs() < 1e-6 {
                1.0
            } else {
                1e-300
            };
        }

        // Simplified Student's t approximation (sufficient for BOCD)
        let diff = (x - location) / scale.sqrt();
        let factor = 1.0 + diff * diff / df;
        let exponent = -(df + 1.0) / 2.0;

        (factor.powf(exponent) / scale.sqrt()).max(1e-300)
    }

    /// Update sufficient statistics for Gaussian model (online updates)
    ///
    /// # Arguments
    ///
    /// * `value` - New observation
    /// * `probs` - New run-length distribution
    fn update_sufficient_statistics(&mut self, value: f64, probs: &[f64]) {
        // Update statistics for each run length using online formulas
        for r in 0..=self.max_run_length {
            if probs[r] > 1e-10 {
                let prev_count = self.counts[r];
                let prev_mean = self.means[r];

                // Update count (weighted by probability)
                self.counts[r] = prev_count + probs[r];

                // Update mean (online Welford's algorithm)
                let delta = value - prev_mean;
                self.means[r] = prev_mean + delta * probs[r] / self.counts[r];

                // Update variance (online Welford's algorithm)
                let delta2 = value - self.means[r];
                self.variances[r] =
                    (prev_count * self.variances[r] + probs[r] * delta * delta2) / self.counts[r];
            }
        }
    }

    /// Reset detector to initial state
    pub fn reset(&mut self) {
        self.run_length_probs = vec![0.0; self.max_run_length + 1];
        self.run_length_probs[0] = 1.0;
        self.means = vec![self.prior_mean; self.max_run_length + 1];
        self.variances = vec![1.0; self.max_run_length + 1];
        self.counts = vec![0.0; self.max_run_length + 1];
        self.time_index = 0;
    }
}

// Tests are in ml/tests/bayesian_changepoint_test.rs
