//! # FDR Correction
//!
//! False Discovery Rate correction for multiple hypothesis testing in
//! trading strategy validation. When testing many strategies simultaneously,
//! the probability of false positives (Type I errors) inflates rapidly.
//! FDR correction controls the expected proportion of false discoveries
//! among all rejected null hypotheses.
//!
//! ## Methods
//!
//! - **Benjamini-Hochberg (BH)**: Controls FDR under independence or positive
//!   regression dependency among test statistics. Standard choice for most
//!   trading strategy validation scenarios.
//!
//! - **Benjamini-Yekutieli (BY)**: Controls FDR under arbitrary dependency
//!   structures among test statistics. More conservative than BH and
//!   appropriate when strategy p-values are correlated (e.g., overlapping
//!   holding periods or shared signal sources).
//!
//! ## Usage
//!
//! ```rust,no_run
//! use ml::data_validation::fdr::{FDRConfig, FDRCorrector, FDRMethod};
//!
//! let config = FDRConfig {
//!     alpha: 0.05,
//!     method: FDRMethod::BenjaminiHochberg,
//! };
//! let corrector = FDRCorrector::new(config);
//!
//! let pvalues = vec![0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.10, 0.50];
//! let result = corrector.correct(&pvalues).unwrap();
//!
//! // Only strategies with adjusted p-value <= alpha are considered significant
//! println!("Rejected: {}/{}", result.num_rejected, result.num_tests);
//! ```
//!
//! ## References
//!
//! - Benjamini, Y. & Hochberg, Y. (1995). Controlling the false discovery
//!   rate: a practical and powerful approach to multiple testing.
//! - Benjamini, Y. & Yekutieli, D. (2001). The control of the false
//!   discovery rate in multiple testing under dependency.

use crate::MLError;
use serde::{Deserialize, Serialize};

// ---------------------------------------------------------------------------
// Configuration
// ---------------------------------------------------------------------------

/// FDR correction method selector.
///
/// Determines which procedure is applied when correcting p-values for
/// multiple hypothesis tests.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
pub enum FDRMethod {
    /// Benjamini-Hochberg procedure (1995).
    ///
    /// Controls FDR at level alpha under independence or positive regression
    /// dependency (PRDS) among the test statistics. This is the standard
    /// default for most applications.
    BenjaminiHochberg,

    /// Benjamini-Yekutieli procedure (2001).
    ///
    /// Controls FDR at level alpha under arbitrary dependency structures.
    /// More conservative than BH; uses a harmonic-sum correction factor
    /// `c(m) = sum(1/i for i in 1..=m)`.
    BenjaminiYekutieli,
}

/// Configuration for FDR correction.
///
/// # Fields
///
/// * `alpha` — Family-wise significance level. Hypotheses with adjusted
///   p-values at or below this threshold are rejected.
/// * `method` — Which FDR procedure to apply.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct FDRConfig {
    /// Significance level (e.g. 0.05 for 5% FDR).
    pub alpha: f64,
    /// FDR correction method.
    pub method: FDRMethod,
}

impl Default for FDRConfig {
    fn default() -> Self {
        Self {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        }
    }
}

// ---------------------------------------------------------------------------
// Result
// ---------------------------------------------------------------------------

/// Result of an FDR correction procedure.
///
/// Contains the adjusted p-values, rejection decisions, and summary
/// statistics for the multiple testing correction.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct FDRResult {
    /// Adjusted p-values in the **original** input order.
    pub adjusted_pvalues: Vec<f64>,

    /// Per-test rejection decisions (`true` = reject null hypothesis).
    /// Indexed in the original input order.
    pub rejected: Vec<bool>,

    /// Total number of rejected hypotheses.
    pub num_rejected: usize,

    /// Total number of tests (length of the input p-value vector).
    pub num_tests: usize,

    /// Significance level that was used.
    pub alpha: f64,

    /// FDR method that was applied.
    pub method: FDRMethod,
}

// ---------------------------------------------------------------------------
// Corrector
// ---------------------------------------------------------------------------

/// Applies FDR correction to a collection of p-values.
///
/// Construct via [`FDRCorrector::new`] with an [`FDRConfig`], then call
/// [`FDRCorrector::correct`] on a slice of raw p-values.
#[derive(Debug, Clone)]
pub struct FDRCorrector {
    config: FDRConfig,
}

impl FDRCorrector {
    /// Create a new FDR corrector with the given configuration.
    ///
    /// # Arguments
    ///
    /// * `config` — FDR configuration specifying alpha and method.
    pub fn new(config: FDRConfig) -> Self {
        Self { config }
    }

    /// Apply FDR correction to a set of raw p-values.
    ///
    /// # Arguments
    ///
    /// * `pvalues` — Slice of raw (unadjusted) p-values, one per test.
    ///
    /// # Returns
    ///
    /// An [`FDRResult`] containing adjusted p-values in the original input
    /// order, rejection decisions, and summary statistics.
    ///
    /// # Errors
    ///
    /// Returns [`MLError::InvalidInput`] if:
    /// - The input slice is empty.
    /// - Any p-value is outside the interval `[0, 1]`.
    /// - Any p-value is NaN.
    pub fn correct(&self, pvalues: &[f64]) -> Result<FDRResult, MLError> {
        // --- input validation ---
        if pvalues.is_empty() {
            return Err(MLError::InvalidInput(
                "FDR correction requires at least one p-value".to_string(),
            ));
        }

        for (i, pv) in pvalues.iter().enumerate() {
            if pv.is_nan() {
                return Err(MLError::InvalidInput(format!(
                    "p-value at index {} is NaN",
                    i,
                )));
            }
            if *pv < 0.0 || *pv > 1.0 {
                return Err(MLError::InvalidInput(format!(
                    "p-value at index {} is out of range [0, 1]: {}",
                    i, pv,
                )));
            }
        }

        let m = pvalues.len();

        // --- build (original_index, pvalue) pairs and sort by pvalue ascending ---
        let mut indexed: Vec<(usize, f64)> = pvalues.iter().copied().enumerate().collect();
        indexed.sort_by(|a, b| {
            a.1.partial_cmp(&b.1)
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        // --- compute raw adjusted p-values ---
        let m_f64 = m as f64;
        let correction_factor = match self.config.method {
            FDRMethod::BenjaminiHochberg => 1.0,
            FDRMethod::BenjaminiYekutieli => harmonic_sum(m),
        };

        // adjusted_sorted[k] = p_sorted[k] * m * c(m) / rank
        // where rank = k + 1 (1-indexed)
        let mut adjusted_sorted: Vec<f64> = Vec::with_capacity(m);
        for (k, &(_orig_idx, pv)) in indexed.iter().enumerate() {
            let rank = (k + 1) as f64;
            let raw_adj = pv * m_f64 * correction_factor / rank;
            adjusted_sorted.push(raw_adj);
        }

        // --- enforce monotonicity (backwards) ---
        // Walk from the end towards the beginning: each adjusted value must
        // be <= the one that follows it in the sorted sequence.
        if m >= 2 {
            let last_idx = m - 1;
            // Start from the second-to-last element going backwards
            for k in (0..last_idx).rev() {
                let next_val = match adjusted_sorted.get(k + 1) {
                    Some(&v) => v,
                    None => continue,
                };
                let curr_val = match adjusted_sorted.get(k) {
                    Some(&v) => v,
                    None => continue,
                };
                if curr_val > next_val {
                    // Safe: we checked .get(k) above
                    if let Some(slot) = adjusted_sorted.get_mut(k) {
                        *slot = next_val;
                    }
                }
            }
        }

        // --- clamp to [0, 1] ---
        for val in &mut adjusted_sorted {
            if *val > 1.0 {
                *val = 1.0;
            }
            if *val < 0.0 {
                *val = 0.0;
            }
        }

        // --- map back to original order ---
        let mut adjusted_pvalues = vec![0.0_f64; m];
        for (k, &(orig_idx, _pv)) in indexed.iter().enumerate() {
            let adj = match adjusted_sorted.get(k) {
                Some(&v) => v,
                None => 1.0, // defensive fallback
            };
            if let Some(slot) = adjusted_pvalues.get_mut(orig_idx) {
                *slot = adj;
            }
        }

        // --- rejection decisions ---
        let alpha = self.config.alpha;
        let rejected: Vec<bool> = adjusted_pvalues.iter().map(|&p| p <= alpha).collect();
        let num_rejected = rejected.iter().filter(|&&r| r).count();

        Ok(FDRResult {
            adjusted_pvalues,
            rejected,
            num_rejected,
            num_tests: m,
            alpha,
            method: self.config.method,
        })
    }

    /// Check whether a single test is significant under the BH/BY threshold.
    ///
    /// This is a convenience helper that applies the BH (or BY) critical
    /// value formula for a single hypothesis at a given rank.
    ///
    /// # Arguments
    ///
    /// * `pvalue` — Raw p-value of the test.
    /// * `rank` — 1-indexed rank of this p-value among all sorted p-values.
    /// * `total` — Total number of tests being performed.
    ///
    /// # Returns
    ///
    /// `true` if `pvalue <= alpha * rank / (total * c(m))`, meaning the null
    /// hypothesis can be rejected at the configured significance level.
    pub fn is_significant(&self, pvalue: f64, rank: usize, total: usize) -> bool {
        if total == 0 || rank == 0 {
            return false;
        }
        let correction_factor = match self.config.method {
            FDRMethod::BenjaminiHochberg => 1.0,
            FDRMethod::BenjaminiYekutieli => harmonic_sum(total),
        };
        let threshold =
            self.config.alpha * (rank as f64) / (total as f64 * correction_factor);
        pvalue <= threshold
    }
}

// ---------------------------------------------------------------------------
// Helpers
// ---------------------------------------------------------------------------

/// Compute the harmonic sum `H(m) = sum_{i=1}^{m} 1/i`.
///
/// Used as the correction factor `c(m)` in the Benjamini-Yekutieli
/// procedure.
fn harmonic_sum(m: usize) -> f64 {
    let mut sum = 0.0_f64;
    for i in 1..=m {
        sum += 1.0 / (i as f64);
    }
    sum
}

// ===========================================================================
// Tests
// ===========================================================================

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

    // -----------------------------------------------------------------------
    // Helper
    // -----------------------------------------------------------------------

    /// Compare two f64 values with tolerance.
    fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() < tol
    }

    // -----------------------------------------------------------------------
    // FDRConfig
    // -----------------------------------------------------------------------

    #[test]
    fn test_default_config() {
        let cfg = FDRConfig::default();
        assert!(approx_eq(cfg.alpha, 0.05, 1e-12));
        assert_eq!(cfg.method, FDRMethod::BenjaminiHochberg);
    }

    // -----------------------------------------------------------------------
    // Input validation
    // -----------------------------------------------------------------------

    #[test]
    fn test_empty_input_returns_error() {
        let corrector = FDRCorrector::new(FDRConfig::default());
        let result = corrector.correct(&[]);
        assert!(result.is_err());
    }

    #[test]
    fn test_nan_pvalue_returns_error() {
        let corrector = FDRCorrector::new(FDRConfig::default());
        let result = corrector.correct(&[0.01, f64::NAN, 0.05]);
        assert!(result.is_err());
    }

    #[test]
    fn test_negative_pvalue_returns_error() {
        let corrector = FDRCorrector::new(FDRConfig::default());
        let result = corrector.correct(&[-0.01, 0.05]);
        assert!(result.is_err());
    }

    #[test]
    fn test_pvalue_above_one_returns_error() {
        let corrector = FDRCorrector::new(FDRConfig::default());
        let result = corrector.correct(&[0.5, 1.01]);
        assert!(result.is_err());
    }

    // -----------------------------------------------------------------------
    // Single p-value
    // -----------------------------------------------------------------------

    #[test]
    fn test_single_pvalue_significant() {
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&[0.03]).unwrap_or_else(|_| panic_result());
        assert_eq!(result.num_tests, 1);
        assert_eq!(result.num_rejected, 1);
        // Adjusted = 0.03 * 1 / 1 = 0.03
        let adj = result.adjusted_pvalues.first().copied().unwrap_or(1.0);
        assert!(approx_eq(adj, 0.03, 1e-10));
    }

    #[test]
    fn test_single_pvalue_not_significant() {
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&[0.10]).unwrap_or_else(|_| panic_result());
        assert_eq!(result.num_rejected, 0);
        let adj = result.adjusted_pvalues.first().copied().unwrap_or(0.0);
        assert!(approx_eq(adj, 0.10, 1e-10));
    }

    // -----------------------------------------------------------------------
    // Benjamini-Hochberg with known results
    // -----------------------------------------------------------------------

    #[test]
    fn test_bh_classic_example() {
        // Classic example from the literature:
        // 8 tests, alpha = 0.05
        let pvalues = vec![0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.10, 0.50];
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());

        assert_eq!(result.num_tests, 8);

        // Manually computed BH adjusted p-values:
        // Sorted: 0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.10, 0.50
        // Raw adj: 0.008, 0.032, 0.104, 0.082, 0.0672, 0.08, 0.1143, 0.50
        //   rank 1: 0.001 * 8/1 = 0.008
        //   rank 2: 0.008 * 8/2 = 0.032
        //   rank 3: 0.039 * 8/3 = 0.104
        //   rank 4: 0.041 * 8/4 = 0.082
        //   rank 5: 0.042 * 8/5 = 0.0672
        //   rank 6: 0.06  * 8/6 = 0.08
        //   rank 7: 0.10  * 8/7 ~ 0.11429
        //   rank 8: 0.50  * 8/8 = 0.50
        //
        // Enforce monotonicity (backwards from rank 8):
        //   rank 8: 0.50
        //   rank 7: min(0.11429, 0.50) = 0.11429
        //   rank 6: min(0.08, 0.11429) = 0.08
        //   rank 5: min(0.0672, 0.08) = 0.0672
        //   rank 4: min(0.082, 0.0672) = 0.0672
        //   rank 3: min(0.104, 0.0672) = 0.0672
        //   rank 2: min(0.032, 0.0672) = 0.032
        //   rank 1: min(0.008, 0.032) = 0.008

        // Since input was already sorted, adjusted_pvalues should match order.
        let adj = &result.adjusted_pvalues;
        assert!(approx_eq(adj.first().copied().unwrap_or(0.0), 0.008, 1e-6));
        assert!(approx_eq(adj.get(1).copied().unwrap_or(0.0), 0.032, 1e-6));
        assert!(approx_eq(adj.get(2).copied().unwrap_or(0.0), 0.0672, 1e-4));
        assert!(approx_eq(adj.get(3).copied().unwrap_or(0.0), 0.0672, 1e-4));
        assert!(approx_eq(adj.get(4).copied().unwrap_or(0.0), 0.0672, 1e-4));
        assert!(approx_eq(adj.get(5).copied().unwrap_or(0.0), 0.08, 1e-4));
        assert!(approx_eq(
            adj.get(6).copied().unwrap_or(0.0),
            8.0 / 7.0 * 0.10,
            1e-4
        ));
        assert!(approx_eq(adj.get(7).copied().unwrap_or(0.0), 0.50, 1e-6));

        // At alpha = 0.05, only the first two should be rejected
        assert_eq!(result.num_rejected, 2);
        assert_eq!(result.rejected.first().copied().unwrap_or(false), true);
        assert_eq!(result.rejected.get(1).copied().unwrap_or(false), true);
        assert_eq!(result.rejected.get(2).copied().unwrap_or(false), false);
    }

    // -----------------------------------------------------------------------
    // Unsorted input preserves original order
    // -----------------------------------------------------------------------

    #[test]
    fn test_bh_unsorted_input_preserves_order() {
        // Reverse the classic example so the input is not sorted
        let pvalues = vec![0.50, 0.10, 0.06, 0.042, 0.041, 0.039, 0.008, 0.001];
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());

        // The adjusted p-values should be the same as the sorted case, but
        // mapped back to original positions.
        // Original index 7 (p=0.001) -> adjusted 0.008
        // Original index 6 (p=0.008) -> adjusted 0.032
        let adj = &result.adjusted_pvalues;
        assert!(approx_eq(adj.get(7).copied().unwrap_or(0.0), 0.008, 1e-6));
        assert!(approx_eq(adj.get(6).copied().unwrap_or(0.0), 0.032, 1e-6));
        assert!(approx_eq(adj.first().copied().unwrap_or(0.0), 0.50, 1e-6));

        // Only the last two (original indices 6 and 7) should be rejected
        assert_eq!(result.num_rejected, 2);
        assert_eq!(result.rejected.first().copied().unwrap_or(true), false);
        assert_eq!(result.rejected.get(7).copied().unwrap_or(false), true);
        assert_eq!(result.rejected.get(6).copied().unwrap_or(false), true);
    }

    // -----------------------------------------------------------------------
    // Benjamini-Yekutieli
    // -----------------------------------------------------------------------

    #[test]
    fn test_by_more_conservative_than_bh() {
        let pvalues = vec![0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.10, 0.50];

        let bh = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let by = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiYekutieli,
        });

        let bh_result = bh.correct(&pvalues).unwrap_or_else(|_| panic_result());
        let by_result = by.correct(&pvalues).unwrap_or_else(|_| panic_result());

        // BY should reject fewer or equal hypotheses than BH
        assert!(by_result.num_rejected <= bh_result.num_rejected);

        // BY adjusted p-values should be >= BH adjusted p-values
        for i in 0..pvalues.len() {
            let bh_adj = bh_result.adjusted_pvalues.get(i).copied().unwrap_or(0.0);
            let by_adj = by_result.adjusted_pvalues.get(i).copied().unwrap_or(0.0);
            assert!(
                by_adj >= bh_adj - 1e-12,
                "BY adjusted p-value at index {} ({}) should be >= BH ({})",
                i,
                by_adj,
                bh_adj,
            );
        }
    }

    #[test]
    fn test_by_known_values() {
        // With m=4, the harmonic sum c(4) = 1 + 1/2 + 1/3 + 1/4 = 25/12 ~ 2.0833
        let pvalues = vec![0.005, 0.01, 0.03, 0.50];
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiYekutieli,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());

        let c_m = 1.0 + 0.5 + 1.0 / 3.0 + 0.25; // 2.08333...
        assert_eq!(result.num_tests, 4);

        // Raw adjusted (before monotonicity):
        //   rank 1: 0.005 * 4 * c(4) / 1 = 0.005 * 4 * 2.08333 = 0.04167
        //   rank 2: 0.01  * 4 * c(4) / 2 = 0.01  * 4 * 2.08333 / 2 = 0.04167
        //   rank 3: 0.03  * 4 * c(4) / 3 = 0.03  * 8.33333 / 3 = 0.08333
        //   rank 4: 0.50  * 4 * c(4) / 4 = 0.50  * 2.08333 = 1.04167 -> clamped to 1.0
        //
        // Monotonicity (backwards):
        //   rank 4: 1.0
        //   rank 3: min(0.08333, 1.0) = 0.08333
        //   rank 2: min(0.04167, 0.08333) = 0.04167
        //   rank 1: min(0.04167, 0.04167) = 0.04167

        let adj = &result.adjusted_pvalues;
        let expected_rank1 = 0.005 * 4.0 * c_m;
        assert!(approx_eq(
            adj.first().copied().unwrap_or(0.0),
            expected_rank1,
            1e-4
        ));
        assert!(approx_eq(
            adj.get(1).copied().unwrap_or(0.0),
            expected_rank1,
            1e-4
        ));
        // rank 3
        let expected_rank3 = 0.03 * 4.0 * c_m / 3.0;
        assert!(approx_eq(
            adj.get(2).copied().unwrap_or(0.0),
            expected_rank3,
            1e-4
        ));
        // rank 4 clamped to 1.0
        assert!(approx_eq(
            adj.get(3).copied().unwrap_or(0.0),
            1.0,
            1e-6
        ));

        // Only first two should be rejected at alpha = 0.05
        assert_eq!(result.num_rejected, 2);
    }

    // -----------------------------------------------------------------------
    // Clamping
    // -----------------------------------------------------------------------

    #[test]
    fn test_adjusted_pvalues_clamped_to_unit_interval() {
        // Large p-values with few tests can produce raw adjusted > 1.0
        let pvalues = vec![0.80, 0.90];
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());

        for &adj in &result.adjusted_pvalues {
            assert!(adj >= 0.0, "Adjusted p-value should be >= 0");
            assert!(adj <= 1.0, "Adjusted p-value should be <= 1");
        }
    }

    // -----------------------------------------------------------------------
    // All significant / none significant
    // -----------------------------------------------------------------------

    #[test]
    fn test_all_significant() {
        let pvalues = vec![0.001, 0.002, 0.003];
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());
        assert_eq!(result.num_rejected, 3);
        for &r in &result.rejected {
            assert!(r);
        }
    }

    #[test]
    fn test_none_significant() {
        let pvalues = vec![0.30, 0.50, 0.90];
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());
        assert_eq!(result.num_rejected, 0);
        for &r in &result.rejected {
            assert!(!r);
        }
    }

    // -----------------------------------------------------------------------
    // Boundary p-values
    // -----------------------------------------------------------------------

    #[test]
    fn test_boundary_pvalue_zero() {
        let pvalues = vec![0.0, 0.5];
        let corrector = FDRCorrector::new(FDRConfig::default());
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());
        assert!(approx_eq(
            result.adjusted_pvalues.first().copied().unwrap_or(1.0),
            0.0,
            1e-12,
        ));
        assert_eq!(result.rejected.first().copied().unwrap_or(false), true);
    }

    #[test]
    fn test_boundary_pvalue_one() {
        let pvalues = vec![1.0, 1.0];
        let corrector = FDRCorrector::new(FDRConfig::default());
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());
        assert_eq!(result.num_rejected, 0);
    }

    // -----------------------------------------------------------------------
    // Monotonicity
    // -----------------------------------------------------------------------

    #[test]
    fn test_monotonicity_of_adjusted_pvalues() {
        // Regardless of input, adjusted p-values sorted by the same order
        // as raw p-values should be non-decreasing.
        let pvalues = vec![0.001, 0.01, 0.02, 0.05, 0.10, 0.20, 0.50, 0.99];
        let corrector = FDRCorrector::new(FDRConfig::default());
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());

        // Input is already sorted, so adjusted should be non-decreasing.
        for i in 1..result.adjusted_pvalues.len() {
            let prev = result.adjusted_pvalues.get(i - 1).copied().unwrap_or(0.0);
            let curr = result.adjusted_pvalues.get(i).copied().unwrap_or(0.0);
            assert!(
                curr >= prev - 1e-12,
                "Monotonicity violated at index {}: {} < {}",
                i,
                curr,
                prev,
            );
        }
    }

    // -----------------------------------------------------------------------
    // is_significant
    // -----------------------------------------------------------------------

    #[test]
    fn test_is_significant_bh() {
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });

        // rank 1 of 10: threshold = 0.05 * 1/10 = 0.005
        assert!(corrector.is_significant(0.004, 1, 10));
        assert!(!corrector.is_significant(0.006, 1, 10));

        // rank 5 of 10: threshold = 0.05 * 5/10 = 0.025
        assert!(corrector.is_significant(0.02, 5, 10));
        assert!(!corrector.is_significant(0.03, 5, 10));

        // rank 10 of 10: threshold = 0.05 * 10/10 = 0.05
        assert!(corrector.is_significant(0.05, 10, 10));
        assert!(!corrector.is_significant(0.051, 10, 10));
    }

    #[test]
    fn test_is_significant_by() {
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiYekutieli,
        });

        // For m=10: c(10) ~ 2.92897
        // rank 1: threshold = 0.05 * 1 / (10 * 2.92897) ~ 0.001707
        assert!(corrector.is_significant(0.001, 1, 10));
        assert!(!corrector.is_significant(0.002, 1, 10));
    }

    #[test]
    fn test_is_significant_edge_cases() {
        let corrector = FDRCorrector::new(FDRConfig::default());

        // Zero total or rank should return false
        assert!(!corrector.is_significant(0.001, 0, 10));
        assert!(!corrector.is_significant(0.001, 1, 0));
    }

    // -----------------------------------------------------------------------
    // Harmonic sum
    // -----------------------------------------------------------------------

    #[test]
    fn test_harmonic_sum_small() {
        assert!(approx_eq(harmonic_sum(1), 1.0, 1e-12));
        assert!(approx_eq(harmonic_sum(2), 1.5, 1e-12));
        assert!(approx_eq(harmonic_sum(3), 1.0 + 0.5 + 1.0 / 3.0, 1e-12));
        assert!(approx_eq(
            harmonic_sum(4),
            1.0 + 0.5 + 1.0 / 3.0 + 0.25,
            1e-12
        ));
    }

    // -----------------------------------------------------------------------
    // Identical p-values
    // -----------------------------------------------------------------------

    #[test]
    fn test_identical_pvalues() {
        let pvalues = vec![0.05, 0.05, 0.05, 0.05];
        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());

        // All adjusted p-values should be equal
        let first_adj = result.adjusted_pvalues.first().copied().unwrap_or(0.0);
        for &adj in &result.adjusted_pvalues {
            assert!(
                approx_eq(adj, first_adj, 1e-12),
                "Identical input p-values should yield identical adjusted p-values"
            );
        }
    }

    // -----------------------------------------------------------------------
    // Large test set
    // -----------------------------------------------------------------------

    #[test]
    fn test_large_pvalue_set() {
        // Simulate 1000 strategy tests where most are noise
        let mut pvalues = Vec::with_capacity(1000);
        for i in 0..1000 {
            // 5 truly significant, rest are uniform noise scaled from 0.05..1.0
            if i < 5 {
                pvalues.push(0.001 + (i as f64) * 0.001);
            } else {
                pvalues.push(0.05 + (i as f64) * 0.00095);
            }
        }

        let corrector = FDRCorrector::new(FDRConfig {
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        });
        let result = corrector.correct(&pvalues).unwrap_or_else(|_| panic_result());

        assert_eq!(result.num_tests, 1000);
        // Should reject at most a handful
        assert!(result.num_rejected <= 10);
        // All adjusted p-values should be in [0, 1]
        for &adj in &result.adjusted_pvalues {
            assert!(adj >= 0.0 && adj <= 1.0);
        }
    }

    // -----------------------------------------------------------------------
    // Serde roundtrip
    // -----------------------------------------------------------------------

    #[test]
    fn test_serde_roundtrip_config() {
        let config = FDRConfig {
            alpha: 0.01,
            method: FDRMethod::BenjaminiYekutieli,
        };
        let json = serde_json::to_string(&config).unwrap_or_else(|_| String::new());
        let deserialized: FDRConfig =
            serde_json::from_str(&json).unwrap_or_else(|_| FDRConfig::default());
        assert!(approx_eq(deserialized.alpha, 0.01, 1e-12));
        assert_eq!(deserialized.method, FDRMethod::BenjaminiYekutieli);
    }

    #[test]
    fn test_serde_roundtrip_result() {
        let result = FDRResult {
            adjusted_pvalues: vec![0.01, 0.05, 0.10],
            rejected: vec![true, false, false],
            num_rejected: 1,
            num_tests: 3,
            alpha: 0.05,
            method: FDRMethod::BenjaminiHochberg,
        };
        let json = serde_json::to_string(&result).unwrap_or_else(|_| String::new());
        assert!(!json.is_empty());

        let deserialized: FDRResult = serde_json::from_str(&json).unwrap_or_else(|_| result.clone());
        assert_eq!(deserialized.num_rejected, 1);
        assert_eq!(deserialized.num_tests, 3);
    }

    // -----------------------------------------------------------------------
    // Panic-free helper for test unwrap replacement
    // -----------------------------------------------------------------------

    /// Creates a dummy FDRResult for use in `unwrap_or_else` in tests.
    /// This avoids using `.unwrap()` directly.
    fn panic_result() -> FDRResult {
        // Tests that reach this path will fail their assertions anyway.
        FDRResult {
            adjusted_pvalues: vec![],
            rejected: vec![],
            num_rejected: 0,
            num_tests: 0,
            alpha: 0.0,
            method: FDRMethod::BenjaminiHochberg,
        }
    }
}