anofox-forecast 0.7.2

//! Conformal prediction for distribution-free prediction intervals.
//!
//! Conformal prediction provides prediction intervals with finite-sample
//! coverage guarantees without distributional assumptions.
//!
//! # Methods
//!
//! - **Split Conformal**: Fast method using a holdout calibration set
//! - **Cross-Validation Conformal**: Uses all data via cross-validation
//! - **Jackknife+**: Leave-one-out with finite sample validity
//!
//! # Example
//!
//! ```ignore
//! use anofox_forecast::postprocess::{ConformalPredictor, ConformalMethod};
//!
//! let mut predictor = ConformalPredictor::new(0.90, ConformalMethod::Split { cal_fraction: 0.2 });
//!
//! // Fit on historical forecast errors
//! let result = predictor.fit(&historical_forecasts, &historical_actuals).unwrap();
//!
//! // Generate prediction intervals
//! let intervals = predictor.predict(&result, &new_forecasts);
//! ```

use crate::error::{ForecastError, Result};
use crate::postprocess::{PointForecasts, PredictionIntervals};

/// Conformal prediction method.
#[derive(Debug, Clone, PartialEq)]
pub enum ConformalMethod {
    /// Split conformal: use a fraction of data for calibration.
    Split {
        /// Fraction of data to use for calibration (0, 1).
        cal_fraction: f64,
    },
    /// Cross-validation conformal: use all data via k-fold CV.
    CrossVal {
        /// Number of folds for cross-validation.
        n_folds: usize,
    },
    /// Jackknife+ (leave-one-out with finite sample validity).
    JackknifePlus,
}

impl Default for ConformalMethod {
    fn default() -> Self {
        Self::Split { cal_fraction: 0.2 }
    }
}

/// Result of fitting a conformal predictor.
#[derive(Debug, Clone)]
pub struct ConformalResult {
    /// Sorted nonconformity scores (absolute residuals).
    scores: Vec<f64>,
    /// The quantile value for the specified coverage.
    quantile_value: f64,
    /// The coverage level used.
    coverage: f64,
    /// The method used.
    method: ConformalMethod,
}

impl ConformalResult {
    /// Get the nonconformity scores.
    pub fn scores(&self) -> &[f64] {
        &self.scores
    }

    /// Get the quantile value (interval half-width).
    pub fn quantile_value(&self) -> f64 {
        self.quantile_value
    }

    /// Get the coverage level.
    pub fn coverage(&self) -> f64 {
        self.coverage
    }

    /// Get the method used.
    pub fn method(&self) -> &ConformalMethod {
        &self.method
    }
}

/// Result of fitting a per-horizon-step conformal predictor.
///
/// Each horizon step gets its own interval half-width, so step h=1 has a tighter
/// interval than step h=12 when forecast error grows with horizon.
#[derive(Debug, Clone)]
pub struct PerStepConformalResult {
    /// Per-step interval half-widths (one per horizon step).
    half_widths: Vec<f64>,
    /// Per-step nonconformity scores (one vec per horizon step).
    scores: Vec<Vec<f64>>,
    /// Coverage level.
    coverage: f64,
    /// Method used.
    method: ConformalMethod,
}

impl PerStepConformalResult {
    /// Get the per-step half-widths.
    pub fn half_widths(&self) -> &[f64] {
        &self.half_widths
    }

    /// Get the per-step nonconformity scores.
    pub fn scores(&self) -> &[Vec<f64>] {
        &self.scores
    }

    /// Get the coverage level.
    pub fn coverage(&self) -> f64 {
        self.coverage
    }

    /// Get the method used.
    pub fn method(&self) -> &ConformalMethod {
        &self.method
    }

    /// Get the number of horizon steps.
    pub fn horizon(&self) -> usize {
        self.half_widths.len()
    }

    /// Apply per-step intervals to a point forecast.
    ///
    /// Returns `(lower, upper)` where each bound has one value per horizon step.
    /// The forecast length must match the horizon used during fitting.
    pub fn predict(&self, point_forecast: &[f64]) -> (Vec<f64>, Vec<f64>) {
        let lower: Vec<f64> = point_forecast
            .iter()
            .zip(self.half_widths.iter())
            .map(|(&p, &hw)| p - hw)
            .collect();
        let upper: Vec<f64> = point_forecast
            .iter()
            .zip(self.half_widths.iter())
            .map(|(&p, &hw)| p + hw)
            .collect();
        (lower, upper)
    }

    /// Apply per-step intervals and return `PredictionIntervals`.
    pub fn predict_intervals(&self, point_forecast: &[f64]) -> PredictionIntervals {
        let (lower, upper) = self.predict(point_forecast);
        PredictionIntervals::from_bounds(lower, upper, self.coverage)
            .expect("Valid prediction intervals")
    }
}

/// Conformal predictor for distribution-free prediction intervals.
///
/// Provides prediction intervals with coverage guarantees based on
/// the conformal prediction framework.
#[derive(Debug, Clone)]
pub struct ConformalPredictor {
    /// Target coverage level (e.g., 0.90 for 90% intervals).
    coverage: f64,
    /// Conformal method to use.
    method: ConformalMethod,
}

impl ConformalPredictor {
    /// Create a new conformal predictor.
    ///
    /// # Arguments
    ///
    /// * `coverage` - Target coverage level in (0, 1)
    /// * `method` - Conformal method to use
    ///
    /// # Panics
    ///
    /// Panics if coverage is not in (0, 1).
    pub fn new(coverage: f64, method: ConformalMethod) -> Self {
        assert!(
            coverage > 0.0 && coverage < 1.0,
            "coverage must be in (0, 1)"
        );
        Self { coverage, method }
    }

    /// Create a split conformal predictor with default calibration fraction.
    pub fn split(coverage: f64) -> Self {
        Self::new(coverage, ConformalMethod::Split { cal_fraction: 0.2 })
    }

    /// Create a cross-validation conformal predictor.
    pub fn cross_val(coverage: f64, n_folds: usize) -> Self {
        Self::new(coverage, ConformalMethod::CrossVal { n_folds })
    }

    /// Create a Jackknife+ conformal predictor.
    pub fn jackknife_plus(coverage: f64) -> Self {
        Self::new(coverage, ConformalMethod::JackknifePlus)
    }

    /// Get the coverage level.
    pub fn coverage(&self) -> f64 {
        self.coverage
    }

    /// Get the method.
    pub fn method(&self) -> &ConformalMethod {
        &self.method
    }

    /// Fit the conformal predictor on historical forecasts and actuals.
    ///
    /// # Arguments
    ///
    /// * `forecasts` - Historical point forecasts
    /// * `actuals` - Corresponding actual values
    ///
    /// # Returns
    ///
    /// A `ConformalResult` containing the calibration information.
    ///
    /// # Errors
    ///
    /// Returns an error if:
    /// - Forecasts and actuals have different lengths
    /// - Not enough data points for the chosen method
    pub fn fit(&self, forecasts: &[f64], actuals: &[f64]) -> Result<ConformalResult> {
        if forecasts.len() != actuals.len() {
            return Err(ForecastError::DimensionMismatch {
                expected: forecasts.len(),
                got: actuals.len(),
            });
        }

        let n = forecasts.len();
        if n == 0 {
            return Err(ForecastError::EmptyData);
        }

        match &self.method {
            ConformalMethod::Split { cal_fraction } => {
                self.fit_split(forecasts, actuals, *cal_fraction)
            }
            ConformalMethod::CrossVal { n_folds } => {
                self.fit_cross_val(forecasts, actuals, *n_folds)
            }
            ConformalMethod::JackknifePlus => self.fit_jackknife_plus(forecasts, actuals),
        }
    }

    /// Fit using split conformal method.
    fn fit_split(
        &self,
        forecasts: &[f64],
        actuals: &[f64],
        cal_fraction: f64,
    ) -> Result<ConformalResult> {
        let n = forecasts.len();
        let cal_size = ((n as f64) * cal_fraction).ceil() as usize;

        if cal_size < 1 {
            return Err(ForecastError::InsufficientData {
                needed: 1,
                got: cal_size,
                hint: None,
            });
        }

        // Use the last cal_size points for calibration
        let cal_start = n - cal_size;

        // Compute nonconformity scores (absolute residuals) on calibration set
        let mut scores: Vec<f64> = forecasts[cal_start..]
            .iter()
            .zip(actuals[cal_start..].iter())
            .map(|(f, a)| (f - a).abs())
            .collect();

        scores.sort_by(|a, b| a.partial_cmp(b).unwrap());

        // Compute the quantile for the coverage level
        // For conformal prediction: quantile at level (1 + 1/n) * coverage
        let adjusted_level = ((cal_size as f64 + 1.0) * self.coverage / cal_size as f64).min(1.0);
        let quantile_idx = ((cal_size as f64) * adjusted_level).ceil() as usize;
        let quantile_idx = quantile_idx.saturating_sub(1).min(scores.len() - 1);
        let quantile_value = scores[quantile_idx];

        Ok(ConformalResult {
            scores,
            quantile_value,
            coverage: self.coverage,
            method: self.method.clone(),
        })
    }

    /// Fit using cross-validation conformal method.
    fn fit_cross_val(
        &self,
        forecasts: &[f64],
        actuals: &[f64],
        n_folds: usize,
    ) -> Result<ConformalResult> {
        let n = forecasts.len();

        if n_folds < 2 {
            return Err(ForecastError::InvalidParameter(
                "n_folds must be at least 2".to_string(),
            ));
        }

        if n < n_folds {
            return Err(ForecastError::InsufficientData {
                needed: n_folds,
                got: n,
                hint: None,
            });
        }

        // For simplicity, we compute leave-one-out scores for each fold
        // In a full implementation, we'd retrain models for each fold
        // Here we use the absolute residuals as nonconformity scores
        let mut scores: Vec<f64> = forecasts
            .iter()
            .zip(actuals.iter())
            .map(|(f, a)| (f - a).abs())
            .collect();

        scores.sort_by(|a, b| a.partial_cmp(b).unwrap());

        // Compute quantile
        let quantile_idx = ((n as f64) * self.coverage).ceil() as usize;
        let quantile_idx = quantile_idx.saturating_sub(1).min(scores.len() - 1);
        let quantile_value = scores[quantile_idx];

        Ok(ConformalResult {
            scores,
            quantile_value,
            coverage: self.coverage,
            method: self.method.clone(),
        })
    }

    /// Fit using Jackknife+ method.
    fn fit_jackknife_plus(&self, forecasts: &[f64], actuals: &[f64]) -> Result<ConformalResult> {
        let n = forecasts.len();

        if n < 2 {
            return Err(ForecastError::InsufficientData {
                needed: 2,
                got: n,
                hint: None,
            });
        }

        // Jackknife+ uses leave-one-out residuals
        // For this implementation, we use the absolute residuals
        // A full implementation would retrain the model for each LOO subset
        let mut scores: Vec<f64> = forecasts
            .iter()
            .zip(actuals.iter())
            .map(|(f, a)| (f - a).abs())
            .collect();

        scores.sort_by(|a, b| a.partial_cmp(b).unwrap());

        // Jackknife+ quantile: ceil((n+1) * coverage) / n
        let adjusted_level = (((n + 1) as f64) * self.coverage / n as f64).min(1.0);
        let quantile_idx = ((n as f64) * adjusted_level).ceil() as usize;
        let quantile_idx = quantile_idx.saturating_sub(1).min(scores.len() - 1);
        let quantile_value = scores[quantile_idx];

        Ok(ConformalResult {
            scores,
            quantile_value,
            coverage: self.coverage,
            method: self.method.clone(),
        })
    }

    /// Fit per-horizon-step conformal intervals.
    ///
    /// Each horizon step gets its own quantile, so early steps (h=1) get tighter
    /// intervals than later steps (h=12) when forecast error grows with horizon.
    ///
    /// # Arguments
    ///
    /// * `fold_forecasts` - Per-fold forecast arrays, each of length `horizon`
    /// * `fold_actuals` - Per-fold actual arrays, each of length `horizon`
    ///
    /// # Errors
    ///
    /// Returns an error if:
    /// - Fewer than 2 folds are provided
    /// - Any fold has a different length than the first
    /// - Forecast/actual lengths don't match within a fold
    pub fn fit_per_step(
        &self,
        fold_forecasts: &[Vec<f64>],
        fold_actuals: &[Vec<f64>],
    ) -> Result<PerStepConformalResult> {
        let n_folds = fold_forecasts.len();
        if n_folds != fold_actuals.len() {
            return Err(ForecastError::DimensionMismatch {
                expected: n_folds,
                got: fold_actuals.len(),
            });
        }
        if n_folds < 2 {
            return Err(ForecastError::InsufficientData {
                needed: 2,
                got: n_folds,
                hint: Some("per-step conformal needs at least 2 folds".into()),
            });
        }

        let horizon = fold_forecasts[0].len();
        if horizon == 0 {
            return Err(ForecastError::EmptyData);
        }

        // Validate all folds have the same horizon length
        for (i, (fc, ac)) in fold_forecasts.iter().zip(fold_actuals.iter()).enumerate() {
            if fc.len() != horizon || ac.len() != horizon {
                return Err(ForecastError::InvalidParameter(format!(
                    "fold {} has forecast len {} and actual len {}, expected {}",
                    i,
                    fc.len(),
                    ac.len(),
                    horizon
                )));
            }
        }

        // Compute pooled quantile as fallback for steps with too few residuals
        let all_scores: Vec<f64> = fold_forecasts
            .iter()
            .zip(fold_actuals.iter())
            .flat_map(|(fc, ac)| fc.iter().zip(ac.iter()).map(|(f, a)| (f - a).abs()))
            .collect();
        let pooled_quantile = Self::compute_quantile(&all_scores, self.coverage, &self.method);

        // Compute per-step quantiles
        let mut half_widths = Vec::with_capacity(horizon);
        let mut per_step_scores = Vec::with_capacity(horizon);

        for t in 0..horizon {
            let step_scores: Vec<f64> = fold_forecasts
                .iter()
                .zip(fold_actuals.iter())
                .map(|(fc, ac)| (fc[t] - ac[t]).abs())
                .collect();

            let hw = if step_scores.len() < 2 {
                // Fall back to pooled quantile when too few residuals
                pooled_quantile
            } else {
                Self::compute_quantile(&step_scores, self.coverage, &self.method)
            };

            half_widths.push(hw);
            per_step_scores.push(step_scores);
        }

        Ok(PerStepConformalResult {
            half_widths,
            scores: per_step_scores,
            coverage: self.coverage,
            method: self.method.clone(),
        })
    }

    /// Compute the conformal quantile from a set of nonconformity scores.
    fn compute_quantile(scores: &[f64], coverage: f64, method: &ConformalMethod) -> f64 {
        let mut sorted = scores.to_vec();
        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
        let n = sorted.len();

        let quantile_idx = match method {
            ConformalMethod::Split { .. } | ConformalMethod::JackknifePlus => {
                let adjusted_level = (((n + 1) as f64) * coverage / n as f64).min(1.0);
                ((n as f64) * adjusted_level).ceil() as usize
            }
            ConformalMethod::CrossVal { .. } => ((n as f64) * coverage).ceil() as usize,
        };

        let idx = quantile_idx.saturating_sub(1).min(n - 1);
        sorted[idx]
    }

    /// Generate quantile forecasts at multiple quantile levels.
    ///
    /// For each quantile level, computes the conformal prediction bound
    /// from the fitted nonconformity scores. Returns a `QuantileForecasts`
    /// with one column per quantile level.
    ///
    /// Quantile levels below 0.5 produce lower bounds (point - hw),
    /// levels above 0.5 produce upper bounds (point + hw).
    pub fn predict_quantiles(
        &self,
        result: &ConformalResult,
        point_forecast: &[f64],
        quantile_levels: &[f64],
    ) -> crate::postprocess::QuantileForecasts {
        let horizon = point_forecast.len();
        let scores = result.scores();

        let mut forecast_values = Vec::with_capacity(horizon);

        for h in 0..horizon {
            let mut row = Vec::with_capacity(quantile_levels.len());
            for &q in quantile_levels {
                // Map quantile to a coverage level and direction
                let hw = if q < 0.5 {
                    // Lower tail: use coverage = 1 - 2*q → half-width, subtract
                    let cov = 1.0 - 2.0 * q;
                    let hw = Self::compute_quantile(scores, cov, &result.method);
                    -hw
                } else if q > 0.5 {
                    // Upper tail: use coverage = 2*q - 1 → half-width, add
                    let cov = 2.0 * q - 1.0;
                    Self::compute_quantile(scores, cov, &result.method)
                } else {
                    // Median: no adjustment
                    0.0
                };
                row.push(point_forecast[h] + hw);
            }
            forecast_values.push(row);
        }

        crate::postprocess::QuantileForecasts::from_values(
            quantile_levels.to_vec(),
            forecast_values,
        )
        .expect("Valid quantile forecasts")
    }

    /// Generate prediction intervals for new point forecasts.
    ///
    /// # Arguments
    ///
    /// * `result` - The fitted conformal result
    /// * `point_forecasts` - New point forecasts to generate intervals for
    ///
    /// # Returns
    ///
    /// Prediction intervals with the target coverage.
    pub fn predict(
        &self,
        result: &ConformalResult,
        point_forecasts: &PointForecasts,
    ) -> PredictionIntervals {
        let values = point_forecasts.values();
        let q = result.quantile_value;

        let lower: Vec<f64> = values.iter().map(|&v| v - q).collect();
        let upper: Vec<f64> = values.iter().map(|&v| v + q).collect();

        // Use unwrap since we control the inputs and they should be valid
        PredictionIntervals::new(
            point_forecasts.timestamps().to_vec(),
            lower,
            upper,
            self.coverage,
        )
        .expect("Valid prediction intervals")
    }

    /// Generate prediction intervals from raw values.
    pub fn predict_values(&self, result: &ConformalResult, values: &[f64]) -> PredictionIntervals {
        let q = result.quantile_value;

        let lower: Vec<f64> = values.iter().map(|&v| v - q).collect();
        let upper: Vec<f64> = values.iter().map(|&v| v + q).collect();

        PredictionIntervals::from_bounds(lower, upper, self.coverage)
            .expect("Valid prediction intervals")
    }
}

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

    // =========================================================================
    // ConformalMethod tests
    // =========================================================================

    mod conformal_method {
        use super::*;

        #[test]
        fn default_is_split_with_20_percent() {
            let method = ConformalMethod::default();
            match method {
                ConformalMethod::Split { cal_fraction } => {
                    assert!((cal_fraction - 0.2).abs() < 1e-10);
                }
                _ => panic!("Expected Split method"),
            }
        }

        #[test]
        fn split_stores_cal_fraction() {
            let method = ConformalMethod::Split { cal_fraction: 0.3 };
            if let ConformalMethod::Split { cal_fraction } = method {
                assert!((cal_fraction - 0.3).abs() < 1e-10);
            } else {
                panic!("Expected Split method");
            }
        }

        #[test]
        fn cross_val_stores_n_folds() {
            let method = ConformalMethod::CrossVal { n_folds: 5 };
            if let ConformalMethod::CrossVal { n_folds } = method {
                assert_eq!(n_folds, 5);
            } else {
                panic!("Expected CrossVal method");
            }
        }

        #[test]
        fn jackknife_plus_variant_exists() {
            let method = ConformalMethod::JackknifePlus;
            assert_eq!(method, ConformalMethod::JackknifePlus);
        }

        #[test]
        fn methods_are_clonable() {
            let method = ConformalMethod::Split { cal_fraction: 0.25 };
            let cloned = method.clone();
            assert_eq!(method, cloned);
        }
    }

    // =========================================================================
    // ConformalPredictor construction tests
    // =========================================================================

    mod construction {
        use super::*;

        #[test]
        fn new_creates_predictor() {
            let predictor =
                ConformalPredictor::new(0.90, ConformalMethod::Split { cal_fraction: 0.2 });
            assert!((predictor.coverage() - 0.90).abs() < 1e-10);
        }

        #[test]
        fn split_creates_split_predictor() {
            let predictor = ConformalPredictor::split(0.95);
            assert!((predictor.coverage() - 0.95).abs() < 1e-10);
            match predictor.method() {
                ConformalMethod::Split { cal_fraction } => {
                    assert!((cal_fraction - 0.2).abs() < 1e-10);
                }
                _ => panic!("Expected Split method"),
            }
        }

        #[test]
        fn cross_val_creates_cv_predictor() {
            let predictor = ConformalPredictor::cross_val(0.90, 5);
            assert!((predictor.coverage() - 0.90).abs() < 1e-10);
            match predictor.method() {
                ConformalMethod::CrossVal { n_folds } => {
                    assert_eq!(*n_folds, 5);
                }
                _ => panic!("Expected CrossVal method"),
            }
        }

        #[test]
        fn jackknife_plus_creates_jackknife_predictor() {
            let predictor = ConformalPredictor::jackknife_plus(0.90);
            assert!((predictor.coverage() - 0.90).abs() < 1e-10);
            assert_eq!(predictor.method(), &ConformalMethod::JackknifePlus);
        }

        #[test]
        #[should_panic(expected = "coverage must be in (0, 1)")]
        fn new_panics_on_zero_coverage() {
            ConformalPredictor::new(0.0, ConformalMethod::default());
        }

        #[test]
        #[should_panic(expected = "coverage must be in (0, 1)")]
        fn new_panics_on_one_coverage() {
            ConformalPredictor::new(1.0, ConformalMethod::default());
        }

        #[test]
        #[should_panic(expected = "coverage must be in (0, 1)")]
        fn new_panics_on_negative_coverage() {
            ConformalPredictor::new(-0.1, ConformalMethod::default());
        }

        #[test]
        fn predictor_is_clonable() {
            let predictor = ConformalPredictor::split(0.90);
            let cloned = predictor.clone();
            assert!((cloned.coverage() - 0.90).abs() < 1e-10);
        }
    }

    // =========================================================================
    // Split conformal fit tests
    // =========================================================================

    mod fit_split {
        use super::*;

        #[test]
        fn fit_returns_result() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();

            assert!((result.coverage() - 0.90).abs() < 1e-10);
            assert!(!result.scores().is_empty());
        }

        #[test]
        fn fit_fails_on_length_mismatch() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0];
            let actuals = vec![10.5, 10.5];

            let result = predictor.fit(&forecasts, &actuals);
            assert!(result.is_err());
        }

        #[test]
        fn fit_fails_on_empty_data() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts: Vec<f64> = vec![];
            let actuals: Vec<f64> = vec![];

            let result = predictor.fit(&forecasts, &actuals);
            assert!(result.is_err());
        }

        #[test]
        fn scores_are_sorted() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0];
            let actuals = vec![10.5, 10.0, 12.5, 13.5, 14.0, 15.5, 16.0, 17.5, 18.0, 19.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            let scores = result.scores();

            for i in 1..scores.len() {
                assert!(scores[i] >= scores[i - 1], "Scores should be sorted");
            }
        }

        #[test]
        fn quantile_value_is_positive() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            assert!(result.quantile_value() >= 0.0);
        }

        #[test]
        fn higher_coverage_gives_larger_quantile() {
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0];
            let actuals = vec![9.0, 12.0, 11.0, 14.0, 13.0, 16.0, 15.0, 18.0, 17.0, 20.0];

            let predictor_90 = ConformalPredictor::split(0.50);
            let predictor_95 = ConformalPredictor::split(0.90);

            let result_90 = predictor_90.fit(&forecasts, &actuals).unwrap();
            let result_95 = predictor_95.fit(&forecasts, &actuals).unwrap();

            assert!(result_95.quantile_value() >= result_90.quantile_value());
        }
    }

    // =========================================================================
    // Cross-validation conformal fit tests
    // =========================================================================

    mod fit_cross_val {
        use super::*;

        #[test]
        fn fit_returns_result() {
            let predictor = ConformalPredictor::cross_val(0.90, 5);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            assert!((result.coverage() - 0.90).abs() < 1e-10);
        }

        #[test]
        fn fit_fails_on_insufficient_folds() {
            let predictor = ConformalPredictor::cross_val(0.90, 1);
            let forecasts = vec![10.0, 11.0, 12.0];
            let actuals = vec![10.5, 10.5, 12.5];

            let result = predictor.fit(&forecasts, &actuals);
            assert!(result.is_err());
        }

        #[test]
        fn fit_fails_when_n_less_than_folds() {
            let predictor = ConformalPredictor::cross_val(0.90, 10);
            let forecasts = vec![10.0, 11.0, 12.0];
            let actuals = vec![10.5, 10.5, 12.5];

            let result = predictor.fit(&forecasts, &actuals);
            assert!(result.is_err());
        }

        #[test]
        fn uses_all_data_for_scores() {
            let predictor = ConformalPredictor::cross_val(0.90, 5);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5, 15.5, 15.5, 17.5, 18.5, 19.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            // CV uses all data points for scores
            assert_eq!(result.scores().len(), 10);
        }
    }

    // =========================================================================
    // Jackknife+ fit tests
    // =========================================================================

    mod fit_jackknife_plus {
        use super::*;

        #[test]
        fn fit_returns_result() {
            let predictor = ConformalPredictor::jackknife_plus(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            assert!((result.coverage() - 0.90).abs() < 1e-10);
        }

        #[test]
        fn fit_fails_on_single_point() {
            let predictor = ConformalPredictor::jackknife_plus(0.90);
            let forecasts = vec![10.0];
            let actuals = vec![10.5];

            let result = predictor.fit(&forecasts, &actuals);
            assert!(result.is_err());
        }

        #[test]
        fn works_with_two_points() {
            let predictor = ConformalPredictor::jackknife_plus(0.90);
            let forecasts = vec![10.0, 11.0];
            let actuals = vec![10.5, 10.5];

            let result = predictor.fit(&forecasts, &actuals);
            assert!(result.is_ok());
        }

        #[test]
        fn uses_all_data_for_scores() {
            let predictor = ConformalPredictor::jackknife_plus(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            assert_eq!(result.scores().len(), 5);
        }
    }

    // =========================================================================
    // Predict tests
    // =========================================================================

    mod predict {
        use super::*;
        use chrono::{TimeZone, Utc};

        fn make_timestamps(n: usize) -> Vec<chrono::DateTime<Utc>> {
            (0..n)
                .map(|i| {
                    Utc.with_ymd_and_hms(2024, 1, 1, 0, 0, 0).unwrap()
                        + chrono::Duration::days(i as i64)
                })
                .collect()
        }

        #[test]
        fn predict_returns_intervals() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5, 15.5, 15.5, 17.5, 18.5, 19.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();

            let new_forecasts = PointForecasts::from_values(vec![20.0, 21.0, 22.0]);
            let intervals = predictor.predict(&result, &new_forecasts);

            assert_eq!(intervals.len(), 3);
            assert!((intervals.coverage() - 0.90).abs() < 1e-10);
        }

        #[test]
        fn predict_with_timestamps() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5, 15.5, 15.5, 17.5, 18.5, 19.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();

            let timestamps = make_timestamps(3);
            let new_forecasts =
                PointForecasts::new(timestamps.clone(), vec![20.0, 21.0, 22.0]).unwrap();
            let intervals = predictor.predict(&result, &new_forecasts);

            assert!(intervals.has_timestamps());
            assert_eq!(intervals.timestamps(), &timestamps);
        }

        #[test]
        fn intervals_are_symmetric() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5, 15.5, 15.5, 17.5, 18.5, 19.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();

            let new_forecasts = PointForecasts::from_values(vec![20.0]);
            let intervals = predictor.predict(&result, &new_forecasts);

            let point = 20.0;
            let lower = intervals.lower()[0];
            let upper = intervals.upper()[0];

            let lower_diff = point - lower;
            let upper_diff = upper - point;

            assert!(
                (lower_diff - upper_diff).abs() < 1e-10,
                "Intervals should be symmetric"
            );
        }

        #[test]
        fn predict_values_works() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5, 15.5, 15.5, 17.5, 18.5, 19.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            let intervals = predictor.predict_values(&result, &[20.0, 21.0]);

            assert_eq!(intervals.len(), 2);
            assert!(!intervals.has_timestamps());
        }

        #[test]
        fn larger_errors_give_wider_intervals() {
            // Small errors
            let forecasts_small = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals_small = vec![10.1, 11.1, 12.1, 13.1, 14.1];

            // Large errors
            let forecasts_large = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals_large = vec![8.0, 13.0, 10.0, 15.0, 12.0];

            let predictor = ConformalPredictor::split(0.90);

            let result_small = predictor.fit(&forecasts_small, &actuals_small).unwrap();
            let result_large = predictor.fit(&forecasts_large, &actuals_large).unwrap();

            assert!(result_large.quantile_value() > result_small.quantile_value());
        }
    }

    // =========================================================================
    // Coverage validation tests
    // =========================================================================

    mod coverage_validation {
        use super::*;

        #[test]
        fn empirical_coverage_approximately_matches_target() {
            // Generate synthetic data with known error distribution
            let n = 100;
            let forecasts: Vec<f64> = (0..n).map(|i| i as f64).collect();
            // Errors uniformly distributed in [-1, 1]
            let errors: Vec<f64> = (0..n)
                .map(|i| ((i * 7 + 3) % 21) as f64 / 10.0 - 1.0)
                .collect();
            let actuals: Vec<f64> = forecasts
                .iter()
                .zip(errors.iter())
                .map(|(f, e)| f + e)
                .collect();

            let predictor = ConformalPredictor::split(0.90);
            let result = predictor.fit(&forecasts, &actuals).unwrap();

            // Test on new data with similar error distribution
            let new_forecasts: Vec<f64> = (100..150).map(|i| i as f64).collect();
            let new_errors: Vec<f64> = (100..150)
                .map(|i| ((i * 7 + 3) % 21) as f64 / 10.0 - 1.0)
                .collect();
            let new_actuals: Vec<f64> = new_forecasts
                .iter()
                .zip(new_errors.iter())
                .map(|(f, e)| f + e)
                .collect();

            let intervals = predictor.predict_values(&result, &new_forecasts);
            let empirical = intervals.empirical_coverage(&new_actuals).unwrap();

            // Allow some slack due to finite sample
            assert!(
                empirical >= 0.70,
                "Empirical coverage {} should be reasonably high",
                empirical
            );
        }
    }

    // =========================================================================
    // ConformalResult tests
    // =========================================================================

    mod conformal_result {
        use super::*;

        #[test]
        fn accessors_return_correct_values() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();

            assert!(!result.scores().is_empty());
            assert!(result.quantile_value() >= 0.0);
            assert!((result.coverage() - 0.90).abs() < 1e-10);
            match result.method() {
                ConformalMethod::Split { .. } => {}
                _ => panic!("Expected Split method"),
            }
        }

        #[test]
        fn result_is_clonable() {
            let predictor = ConformalPredictor::split(0.90);
            let forecasts = vec![10.0, 11.0, 12.0, 13.0, 14.0];
            let actuals = vec![10.5, 10.5, 12.5, 12.5, 14.5];

            let result = predictor.fit(&forecasts, &actuals).unwrap();
            let cloned = result.clone();

            assert_eq!(result.scores(), cloned.scores());
            assert!((result.quantile_value() - cloned.quantile_value()).abs() < 1e-10);
        }
    }

    // =========================================================================
    // Per-step conformal tests
    // =========================================================================

    mod per_step {
        use super::*;

        fn make_folds(n_folds: usize, horizon: usize) -> (Vec<Vec<f64>>, Vec<Vec<f64>>) {
            // Errors grow with horizon step: step t has error ~ t * 0.5
            let mut forecasts = Vec::new();
            let mut actuals = Vec::new();
            for fold in 0..n_folds {
                let fc: Vec<f64> = (0..horizon).map(|t| 100.0 + t as f64).collect();
                let ac: Vec<f64> = (0..horizon)
                    .map(|t| {
                        let error = (t as f64 + 1.0) * 0.5 * if fold % 2 == 0 { 1.0 } else { -1.0 };
                        fc[t] + error
                    })
                    .collect();
                forecasts.push(fc);
                actuals.push(ac);
            }
            (forecasts, actuals)
        }

        #[test]
        fn fit_per_step_returns_result() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(10, 5);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            assert_eq!(result.horizon(), 5);
            assert_eq!(result.half_widths().len(), 5);
            assert_eq!(result.scores().len(), 5);
            assert!((result.coverage() - 0.90).abs() < 1e-10);
        }

        #[test]
        fn later_steps_have_wider_intervals() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(20, 6);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            let hw = result.half_widths();
            // Errors grow with horizon, so half-widths should generally increase
            assert!(
                hw[5] > hw[0],
                "Last step hw ({}) should be > first step hw ({})",
                hw[5],
                hw[0]
            );
        }

        #[test]
        fn predict_returns_correct_bounds() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(10, 3);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            let point = vec![50.0, 51.0, 52.0];
            let (lower, upper) = result.predict(&point);

            assert_eq!(lower.len(), 3);
            assert_eq!(upper.len(), 3);

            for t in 0..3 {
                assert!(lower[t] < point[t]);
                assert!(upper[t] > point[t]);
                // Symmetric
                let diff_low = point[t] - lower[t];
                let diff_high = upper[t] - point[t];
                assert!(
                    (diff_low - diff_high).abs() < 1e-10,
                    "Step {} should be symmetric",
                    t
                );
            }
        }

        #[test]
        fn predict_intervals_returns_prediction_intervals() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(10, 4);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            let intervals = result.predict_intervals(&[10.0, 20.0, 30.0, 40.0]);
            assert_eq!(intervals.len(), 4);
            assert!((intervals.coverage() - 0.90).abs() < 1e-10);
        }

        #[test]
        fn fails_with_fewer_than_2_folds() {
            let predictor = ConformalPredictor::split(0.90);
            let fc = vec![vec![1.0, 2.0]];
            let ac = vec![vec![1.5, 2.5]];
            assert!(predictor.fit_per_step(&fc, &ac).is_err());
        }

        #[test]
        fn fails_with_mismatched_fold_lengths() {
            let predictor = ConformalPredictor::split(0.90);
            let fc = vec![vec![1.0, 2.0], vec![1.0, 2.0, 3.0]];
            let ac = vec![vec![1.5, 2.5], vec![1.5, 2.5, 3.5]];
            assert!(predictor.fit_per_step(&fc, &ac).is_err());
        }

        #[test]
        fn fails_with_mismatched_forecast_actual_count() {
            let predictor = ConformalPredictor::split(0.90);
            let fc = vec![vec![1.0, 2.0], vec![1.0, 2.0]];
            let ac = vec![vec![1.5, 2.5]];
            assert!(predictor.fit_per_step(&fc, &ac).is_err());
        }

        #[test]
        fn works_with_jackknife_plus() {
            let predictor = ConformalPredictor::jackknife_plus(0.90);
            let (fc, ac) = make_folds(10, 3);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();
            assert_eq!(result.horizon(), 3);
            assert!(result.half_widths().iter().all(|&hw| hw > 0.0));
        }

        #[test]
        fn works_with_cross_val() {
            let predictor = ConformalPredictor::cross_val(0.90, 5);
            let (fc, ac) = make_folds(10, 3);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();
            assert_eq!(result.horizon(), 3);
        }

        #[test]
        fn each_step_has_n_folds_scores() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(8, 4);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            for step_scores in result.scores() {
                assert_eq!(step_scores.len(), 8);
            }
        }

        #[test]
        fn result_is_clonable() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(5, 3);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();
            let cloned = result.clone();
            assert_eq!(result.half_widths(), cloned.half_widths());
            assert_eq!(result.horizon(), cloned.horizon());
        }

        // =================================================================
        // All three methods work: split, cross_val, jackknife_plus
        // =================================================================

        #[test]
        fn all_three_methods_produce_valid_per_step_results() {
            let (fc, ac) = make_folds(10, 4);

            let split = ConformalPredictor::split(0.90);
            let cv = ConformalPredictor::cross_val(0.90, 5);
            let jk = ConformalPredictor::jackknife_plus(0.90);

            let r_split = split.fit_per_step(&fc, &ac).unwrap();
            let r_cv = cv.fit_per_step(&fc, &ac).unwrap();
            let r_jk = jk.fit_per_step(&fc, &ac).unwrap();

            for r in [&r_split, &r_cv, &r_jk] {
                assert_eq!(r.horizon(), 4);
                assert_eq!(r.half_widths().len(), 4);
                assert!(r.half_widths().iter().all(|&hw| hw >= 0.0));
                assert!((r.coverage() - 0.90).abs() < 1e-10);
            }
        }

        // =================================================================
        // Fallback to pooled quantile when step has < 2 residuals
        // (use exactly 1 fold — triggers the InsufficientData guard so
        //  we need exactly 2 folds but a step that has few residuals;
        //  since each fold contributes 1 residual per step, 2 folds
        //  gives exactly 2 per step which is the minimum. To test the
        //  fallback we call compute_quantile directly isn't possible, so
        //  instead we verify it doesn't panic with exactly 2 folds.)
        // =================================================================

        #[test]
        fn two_folds_works_uses_per_step_or_pooled() {
            let predictor = ConformalPredictor::split(0.90);
            let fc = vec![vec![10.0, 20.0, 30.0], vec![11.0, 21.0, 31.0]];
            let ac = vec![vec![10.5, 20.5, 30.5], vec![11.5, 21.5, 31.5]];

            let result = predictor.fit_per_step(&fc, &ac).unwrap();
            assert_eq!(result.horizon(), 3);
            // With only 2 scores per step, the quantile must still be computed
            for &hw in result.half_widths() {
                assert!(hw >= 0.0, "half-width should be non-negative");
                assert!(hw.is_finite(), "half-width should be finite");
            }
        }

        // =================================================================
        // Half-widths are non-negative
        // =================================================================

        #[test]
        fn half_widths_are_non_negative() {
            let predictor = ConformalPredictor::split(0.95);
            let (fc, ac) = make_folds(15, 8);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            for (t, &hw) in result.half_widths().iter().enumerate() {
                assert!(
                    hw >= 0.0,
                    "Step {}: half-width {} should be non-negative",
                    t,
                    hw
                );
            }
        }

        // =================================================================
        // predict() returns correct length
        // =================================================================

        #[test]
        fn predict_returns_correct_length() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(10, 7);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            let point: Vec<f64> = (0..7).map(|i| 50.0 + i as f64).collect();
            let (lower, upper) = result.predict(&point);

            assert_eq!(lower.len(), 7);
            assert_eq!(upper.len(), 7);
        }

        // =================================================================
        // predict_intervals() returns valid PredictionIntervals
        // =================================================================

        #[test]
        fn predict_intervals_returns_valid_prediction_intervals() {
            let predictor = ConformalPredictor::split(0.90);
            let (fc, ac) = make_folds(10, 5);
            let result = predictor.fit_per_step(&fc, &ac).unwrap();

            let point = vec![100.0, 200.0, 300.0, 400.0, 500.0];
            let intervals = result.predict_intervals(&point);

            assert_eq!(intervals.len(), 5);
            assert!((intervals.coverage() - 0.90).abs() < 1e-10);
            for i in 0..5 {
                assert!(
                    intervals.lower()[i] <= intervals.upper()[i],
                    "lower > upper at step {}",
                    i
                );
                assert!(intervals.lower()[i].is_finite());
                assert!(intervals.upper()[i].is_finite());
            }
        }

        // =================================================================
        // Fold with constant errors → all half-widths should be equal
        // =================================================================

        #[test]
        fn constant_errors_give_equal_half_widths() {
            let predictor = ConformalPredictor::split(0.90);

            // Every fold, every step has the same absolute error = 2.0
            let n_folds = 10;
            let horizon = 5;
            let fc: Vec<Vec<f64>> = (0..n_folds).map(|_| vec![10.0; horizon]).collect();
            let ac: Vec<Vec<f64>> = (0..n_folds)
                .map(|fold| vec![if fold % 2 == 0 { 12.0 } else { 8.0 }; horizon])
                .collect();

            let result = predictor.fit_per_step(&fc, &ac).unwrap();
            let hw = result.half_widths();

            let first = hw[0];
            for (t, &h) in hw.iter().enumerate() {
                assert!(
                    (h - first).abs() < 1e-10,
                    "Step {} half-width {} should equal step 0 half-width {}",
                    t,
                    h,
                    first
                );
            }
        }

        // =================================================================
        // Fold with growing errors → later half-widths should be larger
        // =================================================================

        #[test]
        fn growing_errors_give_increasing_half_widths() {
            let predictor = ConformalPredictor::split(0.90);

            // Error at step t is proportional to (t+1)
            let n_folds = 20;
            let horizon = 6;
            let fc: Vec<Vec<f64>> = (0..n_folds)
                .map(|_| (0..horizon).map(|t| 100.0 + t as f64).collect())
                .collect();
            let ac: Vec<Vec<f64>> = (0..n_folds)
                .map(|fold| {
                    let sign = if fold % 2 == 0 { 1.0 } else { -1.0 };
                    (0..horizon)
                        .map(|t| 100.0 + t as f64 + sign * (t as f64 + 1.0) * 2.0)
                        .collect()
                })
                .collect();

            let result = predictor.fit_per_step(&fc, &ac).unwrap();
            let hw = result.half_widths();

            assert!(
                hw[horizon - 1] > hw[0],
                "Last step hw ({}) should exceed first step hw ({})",
                hw[horizon - 1],
                hw[0]
            );
        }

        // =================================================================
        // Empty fold_forecasts → error
        // =================================================================

        #[test]
        fn empty_fold_forecasts_errors() {
            let predictor = ConformalPredictor::split(0.90);
            let fc: Vec<Vec<f64>> = vec![];
            let ac: Vec<Vec<f64>> = vec![];
            let result = predictor.fit_per_step(&fc, &ac);
            assert!(result.is_err(), "Empty folds should produce an error");
        }

        // =================================================================
        // Mismatched fold lengths → error
        // =================================================================

        #[test]
        fn mismatched_fold_lengths_errors() {
            let predictor = ConformalPredictor::split(0.90);
            let fc = vec![vec![1.0, 2.0, 3.0], vec![1.0, 2.0]]; // horizon 3 vs 2
            let ac = vec![vec![1.5, 2.5, 3.5], vec![1.5, 2.5]];
            let result = predictor.fit_per_step(&fc, &ac);
            assert!(
                result.is_err(),
                "Mismatched fold lengths should produce an error"
            );
        }

        // =================================================================
        // Mismatched fold_forecasts vs fold_actuals count → error
        // =================================================================

        #[test]
        fn mismatched_fold_forecast_actual_count_errors() {
            let predictor = ConformalPredictor::split(0.90);
            let fc = vec![vec![1.0, 2.0], vec![1.0, 2.0], vec![1.0, 2.0]];
            let ac = vec![vec![1.5, 2.5], vec![1.5, 2.5]]; // 3 vs 2
            let result = predictor.fit_per_step(&fc, &ac);
            assert!(
                result.is_err(),
                "Mismatched forecast/actual fold count should produce an error"
            );
        }
    }

    // =========================================================================
    // predict_quantiles tests
    // =========================================================================

    mod quantile_forecasts {
        use super::*;

        fn make_calibration_data() -> (Vec<f64>, Vec<f64>) {
            let forecasts: Vec<f64> = (0..50).map(|i| i as f64).collect();
            let actuals: Vec<f64> = (0..50)
                .map(|i| i as f64 + ((i * 7 + 3) % 11) as f64 * 0.3 - 1.5)
                .collect();
            (forecasts, actuals)
        }

        #[test]
        fn returns_correct_shape() {
            let (fc, ac) = make_calibration_data();
            let cp = ConformalPredictor::split(0.90);
            let result = cp.fit(&fc, &ac).unwrap();

            let levels = vec![0.10, 0.25, 0.50, 0.75, 0.90];
            let qf = cp.predict_quantiles(&result, &[50.0, 51.0, 52.0], &levels);

            assert_eq!(qf.n_times(), 3);
            assert_eq!(qf.quantiles().len(), 5);
        }

        #[test]
        fn quantiles_are_monotonically_ordered() {
            let (fc, ac) = make_calibration_data();
            let cp = ConformalPredictor::split(0.90);
            let result = cp.fit(&fc, &ac).unwrap();

            let levels = vec![0.05, 0.25, 0.50, 0.75, 0.95];
            let qf = cp.predict_quantiles(&result, &[50.0], &levels);

            let row = qf.at_time(0).unwrap();
            for i in 1..row.len() {
                assert!(
                    row[i] >= row[i - 1],
                    "q[{}]={} < q[{}]={}",
                    i,
                    row[i],
                    i - 1,
                    row[i - 1]
                );
            }
        }

        #[test]
        fn median_equals_point_forecast() {
            let (fc, ac) = make_calibration_data();
            let cp = ConformalPredictor::split(0.90);
            let result = cp.fit(&fc, &ac).unwrap();

            let qf = cp.predict_quantiles(&result, &[100.0], &[0.50]);
            let median = qf.at_time(0).unwrap()[0];
            assert!(
                (median - 100.0).abs() < 1e-10,
                "Median {} should equal point forecast 100.0",
                median
            );
        }

        #[test]
        fn symmetric_quantiles_are_symmetric() {
            let (fc, ac) = make_calibration_data();
            let cp = ConformalPredictor::split(0.90);
            let result = cp.fit(&fc, &ac).unwrap();

            let point = 100.0;
            let qf = cp.predict_quantiles(&result, &[point], &[0.10, 0.50, 0.90]);
            let row = qf.at_time(0).unwrap();

            let lower_dist = point - row[0]; // distance from point to 10th percentile
            let upper_dist = row[2] - point; // distance from point to 90th percentile
            assert!(
                (lower_dist - upper_dist).abs() < 1e-10,
                "Symmetric quantiles: lower_dist={}, upper_dist={}",
                lower_dist,
                upper_dist
            );
        }

        #[test]
        fn works_with_jackknife_plus() {
            let (fc, ac) = make_calibration_data();
            let cp = ConformalPredictor::jackknife_plus(0.90);
            let result = cp.fit(&fc, &ac).unwrap();

            let qf = cp.predict_quantiles(&result, &[50.0, 51.0], &[0.25, 0.75]);
            assert_eq!(qf.n_times(), 2);
            assert_eq!(qf.quantiles().len(), 2);
        }
    }
}