anofox-forecast 0.15.9

//! PELT (Pruned Exact Linear Time) algorithm for changepoint detection.
//!
//! An exact method for detecting multiple changepoints with O(n) average complexity.

use super::cost::{segment_cost, CostFunction};

/// Configuration for PELT algorithm.
#[derive(Debug, Clone)]
pub struct PeltConfig {
    /// Cost function to use
    pub cost_fn: CostFunction,
    /// Penalty for each changepoint (controls number of changepoints)
    pub penalty: f64,
    /// Minimum segment length
    pub min_segment_length: usize,
}

impl Default for PeltConfig {
    fn default() -> Self {
        Self {
            cost_fn: CostFunction::L2,
            penalty: 1.0,
            min_segment_length: 2,
        }
    }
}

impl PeltConfig {
    /// Create a new config with BIC penalty.
    ///
    /// BIC penalty = log(n) where n is the series length.
    pub fn with_bic_penalty(n: usize) -> Self {
        Self {
            penalty: (n as f64).ln(),
            ..Default::default()
        }
    }

    /// Create a new config with AIC penalty.
    ///
    /// AIC penalty = 2.
    pub fn with_aic_penalty() -> Self {
        Self {
            penalty: 2.0,
            ..Default::default()
        }
    }

    /// Set the cost function.
    pub fn cost_function(mut self, cost_fn: CostFunction) -> Self {
        self.cost_fn = cost_fn;
        self
    }

    /// Set the penalty.
    pub fn penalty(mut self, penalty: f64) -> Self {
        self.penalty = penalty;
        self
    }

    /// Set minimum segment length.
    pub fn min_segment_length(mut self, min_len: usize) -> Self {
        self.min_segment_length = min_len.max(1);
        self
    }
}

/// Convenience wrapper for PELT changepoint detection with a builder API.
///
/// # Example
///
/// ```
/// use anofox_forecast::changepoint::{Pelt, CostFunction};
///
/// let mut series = vec![0.0; 50];
/// series.extend(vec![10.0; 50]);
///
/// let result = Pelt::new(CostFunction::L2)
///     .min_size(5)
///     .penalty(5.0)
///     .detect(&series);
///
/// assert!(result.n_changepoints >= 1);
/// ```
#[derive(Debug, Clone)]
pub struct Pelt {
    config: PeltConfig,
}

impl Pelt {
    /// Create a new PELT detector with the given cost function.
    pub fn new(cost_fn: CostFunction) -> Self {
        Self {
            config: PeltConfig {
                cost_fn,
                ..PeltConfig::default()
            },
        }
    }

    /// Set the minimum segment size.
    pub fn min_size(mut self, min_size: usize) -> Self {
        self.config.min_segment_length = min_size.max(1);
        self
    }

    /// Set the penalty for each changepoint.
    pub fn penalty(mut self, penalty: f64) -> Self {
        self.config.penalty = penalty;
        self
    }

    /// Run changepoint detection on the given series.
    pub fn detect(&self, series: &[f64]) -> PeltResult {
        pelt_detect(series, &self.config)
    }

    /// Get a reference to the underlying configuration.
    pub fn config(&self) -> &PeltConfig {
        &self.config
    }

    /// Automatically select the penalty and detect changepoints.
    ///
    /// Uses CROPS (Changepoints for a Range of Penalties) to evaluate
    /// PELT over a geometric range of penalties, then selects the "elbow"
    /// where adding more changepoints yields diminishing cost reduction.
    ///
    /// This is the recommended entry point when you don't know the penalty.
    ///
    /// # Example
    ///
    /// ```
    /// use anofox_forecast::changepoint::{Pelt, CostFunction};
    ///
    /// let mut series = vec![0.0; 50];
    /// series.extend(vec![10.0; 50]);
    /// series.extend(vec![5.0; 50]);
    ///
    /// let result = Pelt::new(CostFunction::L2)
    ///     .min_size(5)
    ///     .auto_detect(&series);
    ///
    /// assert!(result.result.n_changepoints >= 1);
    /// println!("Optimal penalty: {:.2}", result.penalty);
    /// ```
    pub fn auto_detect(&self, series: &[f64]) -> AutoPeltResult {
        let n = series.len();
        if n < 4 {
            return AutoPeltResult {
                result: PeltResult {
                    changepoints: vec![],
                    segments: vec![(0, n)],
                    cost: 0.0,
                    n_changepoints: 0,
                },
                penalty: 0.0,
                crops: vec![],
            };
        }

        // Run CROPS over a geometric range of penalties
        let crops = self.crops(series, None, None);

        // Find the elbow: the penalty where the marginal cost reduction
        // per additional changepoint drops below threshold.
        let best = Self::select_elbow(&crops);

        AutoPeltResult {
            result: best.1.clone(),
            penalty: best.0,
            crops,
        }
    }

    /// CROPS: Changepoints for a Range of Penalties (Haynes et al., 2017).
    ///
    /// Runs PELT over a geometric range of penalties from `pen_min` to `pen_max`,
    /// returning all distinct segmentations found. Efficient because PELT's
    /// O(n) average complexity means each run is cheap.
    ///
    /// Default range: `[0.01 * log(n), 100 * log(n)]` with ~30 steps.
    pub fn crops(
        &self,
        series: &[f64],
        pen_min: Option<f64>,
        pen_max: Option<f64>,
    ) -> Vec<(f64, PeltResult)> {
        let n = series.len();
        let log_n = (n as f64).ln().max(1.0);

        // Default range: [0.5*log(n), 100*log(n)].
        // Starting at 0.5*log(n) avoids over-segmentation on noisy data
        // while still finding changepoints with strong signal.
        let p_min = pen_min.unwrap_or(0.5 * log_n);
        let p_max = pen_max.unwrap_or(100.0 * log_n);

        let n_steps = 30;
        let ratio = (p_max / p_min.max(1e-6)).powf(1.0 / n_steps as f64);

        let mut results: Vec<(f64, PeltResult)> = Vec::new();
        let mut prev_n_cp = usize::MAX;

        let mut pen = p_min;
        for _ in 0..=n_steps {
            let config = PeltConfig {
                penalty: pen,
                cost_fn: self.config.cost_fn,
                min_segment_length: self.config.min_segment_length,
            };
            let result = pelt_detect(series, &config);

            // Only keep distinct segmentations (different n_changepoints)
            if result.n_changepoints != prev_n_cp {
                prev_n_cp = result.n_changepoints;
                results.push((pen, result));
            }

            pen *= ratio;
            if pen > p_max * 1.01 {
                break;
            }
        }

        // Ensure we also try the exact 0-changepoint case (high penalty)
        if results.last().is_none_or(|r| r.1.n_changepoints > 0) {
            let config = PeltConfig {
                penalty: p_max * 10.0,
                cost_fn: self.config.cost_fn,
                min_segment_length: self.config.min_segment_length,
            };
            let result = pelt_detect(series, &config);
            if result.n_changepoints == 0 {
                results.push((p_max * 10.0, result));
            }
        }

        // Sort by penalty (ascending) = by n_changepoints (descending)
        results.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
        results
    }

    /// Select the optimal penalty from CROPS results using the largest-gap method.
    ///
    /// Computes the marginal cost per changepoint at each penalty level.
    /// The "elbow" is at the largest jump in marginal cost — where adding
    /// another changepoint suddenly becomes much more expensive.
    /// Falls back to BIC penalty when no clear elbow is found.
    fn select_elbow(crops: &[(f64, PeltResult)]) -> (f64, &PeltResult) {
        if crops.is_empty() {
            unreachable!("crops should never be empty");
        }
        if crops.len() == 1 {
            return (crops[0].0, &crops[0].1);
        }

        // Compute marginal cost reduction for each adjacent pair.
        // crops is sorted by penalty ascending = n_changepoints descending.
        let mut marginals: Vec<(usize, f64, f64)> = Vec::new(); // (index, marginal, penalty)
        for i in 0..crops.len() - 1 {
            let cp_diff = crops[i].1.n_changepoints as f64 - crops[i + 1].1.n_changepoints as f64;
            if cp_diff > 0.0 {
                let cost_diff = crops[i + 1].1.cost - crops[i].1.cost;
                let marginal = cost_diff / cp_diff;
                marginals.push((i, marginal, crops[i + 1].0));
            }
        }

        if marginals.is_empty() {
            let mid = crops.len() / 2;
            return (crops[mid].0, &crops[mid].1);
        }

        // Find the largest gap in marginal costs.
        // The transition from "cheap to add CPs" → "expensive to add CPs"
        // indicates the natural number of changepoints.
        if marginals.len() >= 2 {
            let mut best_gap = 0.0f64;
            let mut best_idx = 0;

            for i in 0..marginals.len() - 1 {
                let gap = marginals[i + 1].1 - marginals[i].1;
                if gap > best_gap {
                    best_gap = gap;
                    best_idx = i + 1;
                }
            }

            if best_gap > 0.0 {
                // Select the penalty at the gap transition
                let idx = marginals[best_idx].0;
                return (crops[idx].0, &crops[idx].1);
            }
        }

        // Fallback: select the result closest to BIC penalty = 2*log(n)
        let n_total = crops[0]
            .1
            .segments
            .iter()
            .map(|(_, e)| *e)
            .max()
            .unwrap_or(100);
        let bic_pen = 2.0 * (n_total as f64).ln();

        let bic_idx = crops
            .iter()
            .enumerate()
            .min_by_key(|(_, (pen, _))| ((pen - bic_pen).abs() * 1000.0) as i64)
            .map(|(i, _)| i)
            .unwrap_or(0);

        (crops[bic_idx].0, &crops[bic_idx].1)
    }
}

/// Result of automatic changepoint detection.
#[derive(Debug, Clone)]
pub struct AutoPeltResult {
    /// The best changepoint detection result.
    pub result: PeltResult,
    /// The automatically selected penalty.
    pub penalty: f64,
    /// All CROPS results: (penalty, result) pairs sorted by penalty.
    pub crops: Vec<(f64, PeltResult)>,
}

/// Result of PELT changepoint detection.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct PeltResult {
    /// Detected changepoint indices
    pub changepoints: Vec<usize>,
    /// Segment boundaries (start, end) pairs
    pub segments: Vec<(usize, usize)>,
    /// Total cost (excluding penalty)
    pub cost: f64,
    /// Number of changepoints
    pub n_changepoints: usize,
}

impl PeltResult {
    /// Get the segment containing a specific index.
    pub fn segment_for_index(&self, index: usize) -> Option<(usize, usize)> {
        self.segments
            .iter()
            .find(|&&(start, end)| index >= start && index < end)
            .copied()
    }

    /// Get segment means.
    pub fn segment_means(&self, series: &[f64]) -> Vec<f64> {
        self.segments
            .iter()
            .map(|&(start, end)| {
                let segment = &series[start..end];
                if segment.is_empty() {
                    f64::NAN
                } else {
                    segment.iter().sum::<f64>() / segment.len() as f64
                }
            })
            .collect()
    }
}

/// Detect changepoints using the PELT algorithm.
///
/// # Arguments
/// * `series` - Input time series
/// * `config` - PELT configuration
///
/// # Returns
/// PELT result containing detected changepoints
pub fn pelt_detect(series: &[f64], config: &PeltConfig) -> PeltResult {
    let n = series.len();

    if n < 2 * config.min_segment_length {
        return PeltResult {
            changepoints: Vec::new(),
            segments: vec![(0, n)],
            cost: if n > 0 {
                segment_cost(series, config.cost_fn)
            } else {
                0.0
            },
            n_changepoints: 0,
        };
    }

    // F[t] = minimum cost of segmenting series[0..t]
    let mut f = vec![f64::INFINITY; n + 1];
    f[0] = -config.penalty; // So first segment doesn't get penalized twice

    // cp[t] = optimal last changepoint for series[0..t]
    let mut cp: Vec<usize> = vec![0; n + 1];

    // R = set of candidate changepoints (pruned)
    let mut candidates: Vec<usize> = vec![0];

    // Precompute cumulative sums for efficient cost calculation
    let cum_sum: Vec<f64> = std::iter::once(0.0)
        .chain(series.iter().scan(0.0, |acc, &x| {
            *acc += x;
            Some(*acc)
        }))
        .collect();

    let cum_sum_sq: Vec<f64> = std::iter::once(0.0)
        .chain(series.iter().scan(0.0, |acc, &x| {
            *acc += x * x;
            Some(*acc)
        }))
        .collect();

    // Precompute cum_ixy for LinearTrend: cumulative sum of i * x[i]
    let cum_ixy: Vec<f64> = std::iter::once(0.0)
        .chain(series.iter().enumerate().scan(0.0, |acc, (i, &x)| {
            *acc += i as f64 * x;
            Some(*acc)
        }))
        .collect();

    for t in config.min_segment_length..=n {
        let (best_cost, best_cp) = find_best_candidate(
            &candidates,
            t,
            config,
            &f,
            &cum_sum,
            &cum_sum_sq,
            &cum_ixy,
            series,
        );

        f[t] = best_cost;
        cp[t] = best_cp;

        // Pruning: remove candidates that can never be optimal
        candidates.retain(|&s| {
            if t - s < config.min_segment_length {
                return true;
            }
            let seg_cost = compute_segment_cost_fast(
                s,
                t,
                &cum_sum,
                &cum_sum_sq,
                &cum_ixy,
                config.cost_fn,
                series,
            );
            f[s] + seg_cost <= f[t]
        });

        candidates.push(t);
    }

    let changepoints = backtrack_changepoints(&cp, n);
    let segments = build_segments(&changepoints, n);

    // Compute total cost
    let total_cost: f64 = segments
        .iter()
        .map(|&(s, e)| segment_cost(&series[s..e], config.cost_fn))
        .sum();

    PeltResult {
        n_changepoints: changepoints.len(),
        changepoints,
        segments,
        cost: total_cost,
    }
}

/// Find the best candidate changepoint for position t.
#[inline]
fn find_best_candidate(
    candidates: &[usize],
    t: usize,
    config: &PeltConfig,
    f: &[f64],
    cum_sum: &[f64],
    cum_sum_sq: &[f64],
    cum_ixy: &[f64],
    series: &[f64],
) -> (f64, usize) {
    let mut best_cost = f64::INFINITY;
    let mut best_cp = 0;

    for &s in candidates {
        if t - s >= config.min_segment_length {
            let seg_cost = compute_segment_cost_fast(
                s,
                t,
                cum_sum,
                cum_sum_sq,
                cum_ixy,
                config.cost_fn,
                series,
            );
            let total = f[s] + seg_cost + config.penalty;
            if total < best_cost {
                best_cost = total;
                best_cp = s;
            }
        }
    }

    (best_cost, best_cp)
}

/// Backtrack through the changepoint array to recover all changepoints.
#[inline]
fn backtrack_changepoints(cp: &[usize], n: usize) -> Vec<usize> {
    let mut changepoints = Vec::new();
    let mut t = n;
    while t > 0 {
        let prev = cp[t];
        if prev > 0 {
            changepoints.push(prev);
        }
        t = prev;
    }
    changepoints.reverse();
    changepoints
}

/// Build segment boundaries from changepoints.
#[inline]
fn build_segments(changepoints: &[usize], n: usize) -> Vec<(usize, usize)> {
    let mut segments = Vec::with_capacity(changepoints.len() + 1);
    let mut start = 0;
    for &cp_idx in changepoints {
        segments.push((start, cp_idx));
        start = cp_idx;
    }
    segments.push((start, n));
    segments
}

/// Fast segment cost computation using precomputed cumulative sums.
///
/// For L2, Normal, MeanVariance, and LinearTrend cost functions, uses O(1) computation.
/// For other cost functions, falls back to direct computation.
fn compute_segment_cost_fast(
    start: usize,
    end: usize,
    cum_sum: &[f64],
    cum_sum_sq: &[f64],
    cum_ixy: &[f64],
    cost_fn: CostFunction,
    series: &[f64],
) -> f64 {
    let n = end - start;
    if n == 0 {
        return 0.0;
    }

    match cost_fn {
        CostFunction::L2 | CostFunction::Normal | CostFunction::MeanVariance => {
            compute_l2_family_cost(start, end, n, cum_sum, cum_sum_sq, cost_fn)
        }
        CostFunction::LinearTrend => {
            compute_linear_trend_cost(start, end, n, cum_sum, cum_sum_sq, cum_ixy)
        }
        _ => segment_cost(&series[start..end], cost_fn),
    }
}

/// Compute L2, Normal, or MeanVariance cost using cumulative sums.
#[inline]
fn compute_l2_family_cost(
    start: usize,
    end: usize,
    n: usize,
    cum_sum: &[f64],
    cum_sum_sq: &[f64],
    cost_fn: CostFunction,
) -> f64 {
    let n_f64 = n as f64;
    let sum_y = cum_sum[end] - cum_sum[start];
    let sum_y2 = cum_sum_sq[end] - cum_sum_sq[start];
    let mean = sum_y / n_f64;
    let l2 = sum_y2 - n_f64 * mean * mean;

    if n < 2 {
        return l2.max(0.0);
    }

    let var = l2 / n_f64;
    match cost_fn {
        CostFunction::Normal => {
            if var > 1e-10 {
                n_f64 * var.ln()
            } else {
                0.0
            }
        }
        CostFunction::MeanVariance => {
            if var > 1e-10 {
                n_f64 * (1.0 + var.ln())
            } else {
                n_f64
            }
        }
        _ => l2.max(0.0),
    }
}

/// Compute LinearTrend cost (RSS of linear regression) using cumulative sums.
#[inline]
fn compute_linear_trend_cost(
    start: usize,
    end: usize,
    n: usize,
    cum_sum: &[f64],
    cum_sum_sq: &[f64],
    cum_ixy: &[f64],
) -> f64 {
    if n < 2 {
        return 0.0;
    }

    let n_f64 = n as f64;
    let sum_x = n_f64 * (n_f64 - 1.0) / 2.0;
    let sum_x2 = n_f64 * (n_f64 - 1.0) * (2.0 * n_f64 - 1.0) / 6.0;

    let sum_y = cum_sum[end] - cum_sum[start];
    let sum_y2 = cum_sum_sq[end] - cum_sum_sq[start];
    let sum_xy = (cum_ixy[end] - cum_ixy[start]) - (start as f64) * sum_y;

    let ss_xx = sum_x2 - sum_x * sum_x / n_f64;
    let ss_yy = sum_y2 - sum_y * sum_y / n_f64;
    let ss_xy = sum_xy - sum_x * sum_y / n_f64;

    if ss_xx.abs() < 1e-10 {
        ss_yy.max(0.0)
    } else {
        (ss_yy - ss_xy * ss_xy / ss_xx).max(0.0)
    }
}

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

    #[test]
    fn pelt_no_changepoint() {
        // Constant series - no changepoints
        let series = vec![5.0; 20];
        let config = PeltConfig::default().penalty(10.0);
        let result = pelt_detect(&series, &config);

        assert_eq!(result.n_changepoints, 0);
        assert_eq!(result.segments.len(), 1);
        assert_eq!(result.segments[0], (0, 20));
    }

    #[test]
    fn pelt_one_clear_changepoint() {
        // Clear level shift at position 10
        let mut series = vec![0.0; 10];
        series.extend(vec![10.0; 10]);

        let config = PeltConfig::default().penalty(2.0);
        let result = pelt_detect(&series, &config);

        assert_eq!(result.n_changepoints, 1);
        assert_eq!(result.changepoints[0], 10);
        assert_eq!(result.segments, vec![(0, 10), (10, 20)]);
    }

    #[test]
    fn pelt_two_changepoints() {
        // Three distinct levels
        let mut series = vec![0.0; 10];
        series.extend(vec![10.0; 10]);
        series.extend(vec![0.0; 10]);

        let config = PeltConfig::default().penalty(2.0);
        let result = pelt_detect(&series, &config);

        assert_eq!(result.n_changepoints, 2);
        assert!(result.changepoints.contains(&10));
        assert!(result.changepoints.contains(&20));
    }

    #[test]
    fn pelt_short_series() {
        let series = vec![1.0, 2.0, 3.0];
        let config = PeltConfig::default();
        let result = pelt_detect(&series, &config);

        assert_eq!(result.n_changepoints, 0);
    }

    #[test]
    fn pelt_empty_series() {
        let series: Vec<f64> = vec![];
        let config = PeltConfig::default();
        let result = pelt_detect(&series, &config);

        assert_eq!(result.n_changepoints, 0);
        assert!(result.changepoints.is_empty());
    }

    #[test]
    fn pelt_high_penalty_no_changepoints() {
        // Even with clear changepoint, very high penalty prevents detection
        // For series [0]*10 + [100]*10, L2 cost without CP ≈ 50000
        // With CP at 10, cost = 0, so we need penalty > 50000
        let mut series = vec![0.0; 10];
        series.extend(vec![100.0; 10]);

        let config = PeltConfig::default().penalty(100000.0);
        let result = pelt_detect(&series, &config);

        assert_eq!(result.n_changepoints, 0);
    }

    #[test]
    fn pelt_low_penalty_many_changepoints() {
        // Very low penalty may detect spurious changepoints
        let series: Vec<f64> = (0..50).map(|i| i as f64 + ((i * 7) % 3) as f64).collect();
        let config = PeltConfig::default().penalty(0.01);
        let result = pelt_detect(&series, &config);

        // Result should be valid (n_changepoints is usize, always >= 0)
        let _ = result.n_changepoints;
    }

    #[test]
    fn pelt_config_bic() {
        let config = PeltConfig::with_bic_penalty(100);
        assert_relative_eq!(config.penalty, 100.0_f64.ln(), epsilon = 1e-10);
    }

    #[test]
    fn pelt_config_aic() {
        let config = PeltConfig::with_aic_penalty();
        assert_relative_eq!(config.penalty, 2.0, epsilon = 1e-10);
    }

    #[test]
    fn pelt_config_builder() {
        let config = PeltConfig::default()
            .cost_function(CostFunction::L1)
            .penalty(5.0)
            .min_segment_length(5);

        assert_eq!(config.cost_fn, CostFunction::L1);
        assert_relative_eq!(config.penalty, 5.0, epsilon = 1e-10);
        assert_eq!(config.min_segment_length, 5);
    }

    #[test]
    fn pelt_segment_means() {
        let mut series = vec![1.0; 5];
        series.extend(vec![10.0; 5]);

        let config = PeltConfig::default().penalty(1.0);
        let result = pelt_detect(&series, &config);

        let means = result.segment_means(&series);

        // Should have 2 segments
        assert_eq!(means.len(), 2);
        assert_relative_eq!(means[0], 1.0, epsilon = 1e-10);
        assert_relative_eq!(means[1], 10.0, epsilon = 1e-10);
    }

    #[test]
    fn pelt_segment_for_index() {
        let mut series = vec![0.0; 10];
        series.extend(vec![10.0; 10]);

        let config = PeltConfig::default().penalty(1.0);
        let result = pelt_detect(&series, &config);

        // Index 5 should be in first segment
        assert_eq!(result.segment_for_index(5), Some((0, 10)));
        // Index 15 should be in second segment
        assert_eq!(result.segment_for_index(15), Some((10, 20)));
    }

    #[test]
    fn pelt_min_segment_length() {
        // Changepoint at position 2, but min_segment_length = 5
        let mut series = vec![0.0; 2];
        series.extend(vec![100.0; 18]);

        let config = PeltConfig::default().penalty(1.0).min_segment_length(5);
        let result = pelt_detect(&series, &config);

        // Changepoint at 2 should not be detected due to min segment length
        for cp in &result.changepoints {
            assert!(*cp >= 5);
        }
    }

    // ==================== Integration tests for new cost functions ====================

    #[test]
    fn pelt_linear_trend_detects_slope_change() {
        // First segment: slope +1 (y = x)
        // Second segment: slope -1 (y = 100 - x)
        let mut series: Vec<f64> = (0..50).map(|i| i as f64).collect();
        series.extend((0..50).map(|i| 100.0 - i as f64));

        let config = PeltConfig::default()
            .cost_function(CostFunction::LinearTrend)
            .penalty(100.0);
        let result = pelt_detect(&series, &config);

        // Should detect the slope change around index 50
        assert!(result.n_changepoints >= 1);
        let cp = result.changepoints[0];
        assert!(
            (45..=55).contains(&cp),
            "Expected changepoint near 50, got {}",
            cp
        );
    }

    #[test]
    fn pelt_linear_trend_no_change_for_constant_slope() {
        // Constant slope across entire series
        let series: Vec<f64> = (0..100).map(|i| 2.0 * i as f64 + 5.0).collect();

        let config = PeltConfig::default()
            .cost_function(CostFunction::LinearTrend)
            .penalty(50.0);
        let result = pelt_detect(&series, &config);

        // No changepoints in a perfectly linear series
        assert_eq!(result.n_changepoints, 0);
    }

    #[test]
    fn pelt_mean_variance_detects_variance_shift() {
        use rand::{rngs::StdRng, Rng, SeedableRng};
        let mut rng = StdRng::seed_from_u64(42);

        // First segment: low variance (std = 1)
        // Second segment: high variance (std = 10)
        let mut series: Vec<f64> = (0..50).map(|_| 10.0 + rng.gen_range(-1.0..1.0)).collect();
        series.extend((0..50).map(|_| 10.0 + rng.gen_range(-10.0..10.0)));

        let config = PeltConfig::default()
            .cost_function(CostFunction::MeanVariance)
            .penalty(50.0);
        let result = pelt_detect(&series, &config);

        // Should detect variance change
        assert!(result.n_changepoints >= 1);
    }

    #[test]
    fn pelt_mean_variance_detects_joint_change() {
        // First segment: mean=0, variance=1
        // Second segment: mean=10, variance=4
        let mut series: Vec<f64> = vec![
            -0.5, 0.3, -0.2, 0.8, -0.1, 0.4, -0.6, 0.2, -0.3, 0.5, // std ~0.5
            -0.4, 0.1, -0.7, 0.6, -0.2, 0.3, -0.5, 0.4, -0.1, 0.2,
        ];
        series.extend(vec![
            8.0, 12.0, 7.0, 13.0, 9.0, 11.0, 6.0, 14.0, 8.0, 12.0, // mean 10, larger spread
            7.0, 13.0, 8.0, 12.0, 9.0, 11.0, 6.0, 14.0, 7.0, 13.0,
        ]);

        let config = PeltConfig::default()
            .cost_function(CostFunction::MeanVariance)
            .penalty(5.0);
        let result = pelt_detect(&series, &config);

        // Should detect change
        assert!(result.n_changepoints >= 1);
        let cp = result.changepoints[0];
        assert!(
            (15..=25).contains(&cp),
            "Expected changepoint near 20, got {}",
            cp
        );
    }

    #[test]
    fn pelt_cusum_detects_sustained_shift() {
        // First segment: centered around 0
        // Second segment: sustained positive shift to 5
        let mut series: Vec<f64> = vec![
            0.1, -0.2, 0.3, -0.1, 0.2, -0.3, 0.1, -0.2, 0.3, -0.1, 0.0, 0.1, -0.1, 0.2, -0.2, 0.1,
            -0.3, 0.2, -0.1, 0.0,
        ];
        series.extend(vec![
            5.1, 4.9, 5.2, 4.8, 5.0, 5.1, 4.9, 5.2, 4.8, 5.0, 5.1, 4.9, 5.2, 4.8, 5.0, 5.1, 4.9,
            5.2, 4.8, 5.0,
        ]);

        let config = PeltConfig::default()
            .cost_function(CostFunction::Cusum)
            .penalty(2.0);
        let result = pelt_detect(&series, &config);

        // Should detect the sustained shift
        assert!(result.n_changepoints >= 1);
        let cp = result.changepoints[0];
        assert!(
            (15..=25).contains(&cp),
            "Expected changepoint near 20, got {}",
            cp
        );
    }

    #[test]
    fn pelt_cusum_no_change_for_balanced() {
        // Series that oscillates evenly around mean - no sustained shift
        let series: Vec<f64> = (0..40)
            .map(|i| if i % 2 == 0 { 1.0 } else { -1.0 })
            .collect();

        let config = PeltConfig::default()
            .cost_function(CostFunction::Cusum)
            .penalty(5.0);
        let result = pelt_detect(&series, &config);

        // Balanced oscillations shouldn't trigger CUSUM changepoints
        // (might still detect some due to segment boundaries)
        assert!(result.n_changepoints <= 2);
    }
}