scirs2-special 0.2.0

//! QuickCheck-based property testing for special functions
//!
//! This module provides comprehensive randomized property testing
//! to ensure mathematical correctness across wide parameter ranges.
//!
//! Configure test intensity with environment variables:
//! - ADVANCED_FAST_TESTS=1: Optimized mode for rapid development iteration (10 tests)
//! - QUICK_TESTS=1: Run with reduced test cases for faster compilation (50 tests)
//! - COMPREHENSIVE_TESTS=1: Run full test suite (500 tests, default in release mode)

#![allow(dead_code)]

use quickcheck::{Arbitrary, Gen, QuickCheck, TestResult};
use scirs2_core::numeric::Complex64;
use std::f64;

/// Configuration for test intensity
#[derive(Debug, Clone)]
pub struct TestConfig {
    pub test_count: u64,
    pub max_iterations: u64,
    pub enable_expensive_tests: bool,
}

impl Default for TestConfig {
    fn default() -> Self {
        // Check environment variables for test configuration
        let quick_tests = std::env::var("QUICK_TESTS").is_ok();
        let comprehensive_tests = std::env::var("COMPREHENSIVE_TESTS").is_ok();
        let advanced_fast = std::env::var("ADVANCED_FAST_TESTS").is_ok();

        if advanced_fast {
            // Optimized mode for development iteration
            Self {
                test_count: 10,
                max_iterations: 20,
                enable_expensive_tests: false,
            }
        } else if quick_tests {
            Self {
                test_count: 50,
                max_iterations: 100,
                enable_expensive_tests: false,
            }
        } else if comprehensive_tests || cfg!(not(debug_assertions)) {
            Self {
                test_count: 500,
                max_iterations: 1000,
                enable_expensive_tests: true,
            }
        } else {
            Self {
                test_count: 100,
                max_iterations: 200,
                enable_expensive_tests: false,
            }
        }
    }
}

/// Custom type for positive f64 values
#[derive(Clone, Debug)]
struct PositiveF64(f64);

impl Arbitrary for PositiveF64 {
    fn arbitrary(g: &mut Gen) -> Self {
        let val: f64 = Arbitrary::arbitrary(g);
        // Filter out NaN and infinite values, use smaller range for convergence
        let finite_val = if val.is_finite() { val } else { 1.0 };
        PositiveF64((finite_val.abs() % 20.0) + f64::EPSILON)
    }
}

/// Custom type for small positive integers
#[derive(Clone, Debug)]
struct SmallInt(usize);

impl Arbitrary for SmallInt {
    fn arbitrary(g: &mut Gen) -> Self {
        let val: usize = Arbitrary::arbitrary(g);
        SmallInt(val % 20)
    }
}

/// Custom type for reasonable complex numbers
#[derive(Clone, Debug)]
struct ReasonableComplex(Complex64);

impl Arbitrary for ReasonableComplex {
    fn arbitrary(g: &mut Gen) -> Self {
        let re: f64 = Arbitrary::arbitrary(g);
        let im: f64 = Arbitrary::arbitrary(g);
        // Filter out NaN and infinite values
        let finite_re = if re.is_finite() { re } else { 1.0 };
        let finite_im = if im.is_finite() { im } else { 0.0 };
        ReasonableComplex(Complex64::new(
            (finite_re % 5.0).clamp(-5.0, 5.0),
            (finite_im % 5.0).clamp(-5.0, 5.0),
        ))
    }
}

/// Helper function to run QuickCheck tests with custom configuration
#[allow(dead_code)]
pub fn run_quickcheck_test<F, P>(prop: F, config: &TestConfig) -> bool
where
    F: Fn(P) -> bool + Send + Sync + 'static + quickcheck::Testable,
    P: Arbitrary + Clone + Send + std::fmt::Debug + 'static,
{
    QuickCheck::new()
        .tests(config.test_count)
        .max_tests(config.max_iterations)
        .quickcheck(prop);
    true
}

/// Helper function to run QuickCheck tests that return TestResult
#[allow(dead_code)]
pub fn run_quickcheck_test_result<F, P>(prop: F, config: &TestConfig) -> bool
where
    F: Fn(P) -> TestResult + Send + Sync + 'static + quickcheck::Testable,
    P: Arbitrary + Clone + Send + std::fmt::Debug + 'static,
{
    QuickCheck::new()
        .tests(config.test_count)
        .max_tests(config.max_iterations)
        .quickcheck(prop);
    true
}

#[cfg(test)]
mod gamma_properties {
    use super::*;
    use quickcheck_macros::quickcheck;

    /// Optimized gamma recurrence relation test with early termination
    #[quickcheck]
    fn gamma_recurrence_relation(x: PositiveF64) -> TestResult {
        let x = x.0;

        // More restrictive bounds for faster testing
        if x >= 50.0 || x <= f64::EPSILON {
            return TestResult::discard();
        }

        let gamma_x = crate::gamma::gamma(x);
        let gamma_x_plus_1 = crate::gamma::gamma(x + 1.0);
        let expected = x * gamma_x;

        if !gamma_x.is_finite() || !gamma_x_plus_1.is_finite() {
            return TestResult::discard();
        }

        let relative_error = (gamma_x_plus_1 - expected).abs() / expected.abs();
        TestResult::from_bool(relative_error < 1e-8)
    }

    /// Optimized log gamma test with bounds checking
    #[quickcheck]
    fn log_gamma_additive_property(x: PositiveF64, n: SmallInt) -> TestResult {
        let x = x.0;
        let n_usize = n.0;
        let n = n_usize as f64;

        // Tighter bounds to avoid numerical issues with large values
        // Discard zero n as it's trivial
        if n_usize == 0 || !(1.0..=14.0).contains(&x) || x + n > 20.0 || n > 10.0 {
            return TestResult::discard();
        }

        let log_gamma_x = crate::gamma::loggamma(x);
        let log_gamma_x_n = crate::gamma::loggamma(x + n);

        // Calculate sum of logarithms
        let mut log_sum = log_gamma_x;
        for i in 0..n_usize {
            log_sum += (x + i as f64).ln();
        }

        if !log_gamma_x_n.is_finite() || !log_sum.is_finite() {
            return TestResult::discard();
        }

        // Use relative error tolerance for better numerical stability
        let abs_error = (log_gamma_x_n - log_sum).abs();
        let rel_error = abs_error / log_gamma_x_n.abs().max(1.0);

        TestResult::from_bool(abs_error < 1e-8 || rel_error < 1e-10)
    }

    /// Regression test for specific log_gamma_additive_property case
    #[test]
    fn log_gamma_additive_property_regression() {
        // Test the exact case that was failing: x=1.5, n=1
        let x = 1.5000000000000002f64;
        let n = 1usize;

        let log_gamma_x = crate::gamma::loggamma(x);
        let log_gamma_x_n = crate::gamma::loggamma(x + n as f64);

        println!("log_gamma({}) = {}", x, log_gamma_x);
        println!("log_gamma({}) = {}", x + n as f64, log_gamma_x_n);

        // Expected values based on Γ(1.5) = sqrt(π)/2
        let gamma_1_5 = 0.88622692545275801364908374167057f64;
        let gamma_2_5 = 1.329340388179137f64;
        println!("Expected log_gamma(1.5) = {}", gamma_1_5.ln());
        println!("Expected log_gamma(2.5) = {}", gamma_2_5.ln());

        let mut log_sum = log_gamma_x;
        for i in 0..n {
            log_sum += (x + i as f64).ln();
        }

        println!("log_sum = {} + {} = {}", log_gamma_x, x.ln(), log_sum);
        println!("Difference: {}", (log_gamma_x_n - log_sum).abs());

        let abs_error = (log_gamma_x_n - log_sum).abs();
        let rel_error = abs_error / log_gamma_x_n.abs().max(1.0);

        assert!(
            abs_error < 1e-8 || rel_error < 1e-10,
            "log_gamma property failed: abs_error={}, rel_error={}\n\
             log_gamma({}) = {} (expected {})\n\
             log_gamma({}) = {} (expected {})",
            abs_error,
            rel_error,
            x,
            log_gamma_x,
            gamma_1_5.ln(),
            x + n as f64,
            log_gamma_x_n,
            gamma_2_5.ln()
        );
    }

    /// Beta function symmetry test
    #[quickcheck]
    fn beta_symmetry(x: PositiveF64, y: PositiveF64) -> TestResult {
        let x = x.0.min(20.0); // Reduced range
        let y = y.0.min(20.0);

        let beta_xy = crate::gamma::beta(x, y);
        let beta_yx = crate::gamma::beta(y, x);

        if !beta_xy.is_finite() || !beta_yx.is_finite() {
            return TestResult::discard();
        }

        TestResult::from_bool((beta_xy - beta_yx).abs() < 1e-12 * beta_xy.abs())
    }

    /// Comprehensive test runner for gamma properties
    #[test]
    fn test_gamma_properties_comprehensive() {
        let config = TestConfig::default();
        println!("Running gamma property tests with config: {config:?}");

        // Run tests only if not in quick mode
        if !config.enable_expensive_tests {
            println!(
                "Skipping expensive gamma property tests (set COMPREHENSIVE_TESTS=1 to enable)"
            );
        }

        // Note: Individual quickcheck tests will run with default settings
        // This is mainly for documentation and future custom test runners
    }
}

#[cfg(test)]
mod bessel_properties {
    use super::*;
    use quickcheck_macros::quickcheck;

    #[quickcheck]
    fn bessel_j_derivative_relation(x: PositiveF64) -> bool {
        let x = x.0.min(50.0);

        if x < 0.1 {
            return true; // Skip near zero
        }

        let j0_prime = crate::bessel::j0_prime(x);
        let j1 = crate::bessel::j1(x);

        (j0_prime + j1).abs() < 1e-8
    }

    #[quickcheck]
    fn bessel_recurrence_relation(n: SmallInt, x: PositiveF64) -> bool {
        let n = n.0.max(1);
        let x = x.0.min(50.0);

        if x < 0.1 || n == 0 {
            return true;
        }

        let jnminus_1 = crate::bessel::jn((n - 1) as i32, x);
        let jn = crate::bessel::jn(n as i32, x);
        let jn_plus_1 = crate::bessel::jn((n + 1) as i32, x);

        let expected = (2.0 * n as f64 / x) * jn - jnminus_1;

        if !jn_plus_1.is_finite() || !expected.is_finite() {
            return true;
        }

        (jn_plus_1 - expected).abs() < 1e-8 * expected.abs().max(1.0)
    }

    #[quickcheck]
    fn bessel_wronskian(x: PositiveF64) -> bool {
        let x = x.0.min(50.0);

        if x < 0.1 {
            return true;
        }

        let j0 = crate::bessel::j0(x);
        let y0 = crate::bessel::y0(x);
        let j0_prime = crate::bessel::j0_prime(x);
        let y0_prime = crate::bessel::y0_prime(x);

        let wronskian = j0 * y0_prime - j0_prime * y0;
        let expected = 2.0 / (f64::consts::PI * x);

        if !wronskian.is_finite() || !expected.is_finite() {
            return true;
        }

        (wronskian - expected).abs() < 1e-8 * expected.abs()
    }
}

#[cfg(test)]
mod error_function_properties {
    use super::*;
    use quickcheck_macros::quickcheck;

    #[quickcheck]
    fn erf_odd_function(x: f64) -> bool {
        // Filter out NaN and extreme values
        if !x.is_finite() || x.abs() > 10.0 {
            return true; // Skip invalid/extreme values
        }

        let erf_x = crate::erf::erf(x);
        let erf_neg_x = crate::erf::erf(-x);

        (erf_x + erf_neg_x).abs() < 1e-12
    }

    #[quickcheck]
    fn erf_erfc_complement(x: f64) -> bool {
        // Filter out NaN and extreme values
        if !x.is_finite() || x.abs() > 10.0 {
            return true;
        }

        let erf_x = crate::erf::erf(x);
        let erfc_x = crate::erf::erfc(x);

        (erf_x + erfc_x - 1.0).abs() < 1e-12
    }

    #[quickcheck]
    fn erf_bounds(x: f64) -> bool {
        let erf_x = crate::erf::erf(x);
        // Handle NaN case: if input is NaN, output can be NaN (acceptable)
        if x.is_nan() {
            return erf_x.is_nan();
        }
        (-1.0..=1.0).contains(&erf_x)
    }

    #[quickcheck]
    fn erfinv_inverse_property(x: f64) -> bool {
        // Handle NaN input
        if x.is_nan() {
            return true;
        }

        // Handle edge cases ±1.0 first (before clamping)
        if x.abs() >= 1.0 {
            return true; // Skip boundary cases where erfinv is infinite
        }

        let x = x.clamp(-0.999, 0.999); // Keep within valid range

        let erfinv_x = crate::erf::erfinv(x);
        if !erfinv_x.is_finite() {
            return true;
        }

        let erf_erfinv = crate::erf::erf(erfinv_x);
        (erf_erfinv - x).abs() < 1e-10
    }
}

#[cfg(test)]
mod orthogonal_polynomial_properties {
    use super::*;
    use quickcheck_macros::quickcheck;

    #[quickcheck]
    fn legendre_symmetry(n: SmallInt, x: f64) -> bool {
        let n = n.0;

        // Handle NaN input
        if x.is_nan() {
            return true;
        }

        let x = x.clamp(-1.0, 1.0);

        let p_n_x = crate::orthogonal::legendre(n, x);
        let p_n_neg_x = crate::orthogonal::legendre(n, -x);

        let expected = if n.is_multiple_of(2) { p_n_x } else { -p_n_x };

        (p_n_neg_x - expected).abs() < 1e-10
    }

    #[quickcheck]
    fn chebyshev_t_bounds(n: SmallInt, x: f64) -> bool {
        let n = n.0;

        // Handle NaN input
        if x.is_nan() {
            return true; // Skip NaN cases
        }

        let x = x.clamp(-1.0, 1.0);

        let t_n = crate::orthogonal::chebyshev(n, x, true);

        // Handle NaN output (shouldn't happen but safety check)
        if t_n.is_nan() {
            return false;
        }

        // Chebyshev polynomials are bounded by 1 on [-1, 1]
        t_n.abs() <= 1.0 + 1e-10
    }

    #[quickcheck]
    fn hermite_recurrence(n: SmallInt, x: f64) -> bool {
        // Fixed: operator precedence corrected in orthogonal.rs
        // H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)
        let n = n.0.clamp(1, 10);
        let x = x.clamp(-5.0, 5.0);

        if !x.is_finite() {
            return true;
        }

        let h_n = crate::orthogonal::hermite(n, x);
        let h_n_minus_1 = crate::orthogonal::hermite(n - 1, x);
        let h_n_plus_1 = crate::orthogonal::hermite(n + 1, x);

        // Check recurrence: H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)
        let expected = 2.0 * x * h_n - 2.0 * (n as f64) * h_n_minus_1;

        if !expected.is_finite() || !h_n_plus_1.is_finite() {
            return true;
        }

        (h_n_plus_1 - expected).abs() <= 1e-8 * (1.0 + expected.abs().max(h_n_plus_1.abs()))
    }
}

#[cfg(test)]
mod complex_function_properties {
    use super::*;
    use quickcheck_macros::quickcheck;

    #[quickcheck]
    fn complex_erf_conjugate_symmetry(z: ReasonableComplex) -> bool {
        let z = z.0;

        let erf_z = crate::erf::complex::erf_complex(z);
        let erf_conj_z = crate::erf::complex::erf_complex(z.conj());
        let expected = erf_z.conj();

        (erf_conj_z - expected).norm() < 1e-10
    }

    #[quickcheck]
    fn complex_gamma_conjugate_symmetry(z: ReasonableComplex) -> bool {
        let z = z.0;

        // Skip near poles
        if z.re <= 0.0 && (z.re.fract().abs() < 0.1 || z.im.abs() < 0.1) {
            return true;
        }

        let gamma_z = crate::gamma::complex::gamma_complex(z);
        let gamma_conj_z = crate::gamma::complex::gamma_complex(z.conj());
        let expected = gamma_z.conj();

        if !gamma_z.is_finite() || !gamma_conj_z.is_finite() {
            return true;
        }

        (gamma_conj_z - expected).norm() < 1e-8 * gamma_z.norm()
    }
}

#[cfg(test)]
mod statistical_function_properties {
    use super::*;
    use quickcheck_macros::quickcheck;

    #[quickcheck]
    fn logistic_bounds(x: f64) -> bool {
        // Handle NaN input
        if x.is_nan() || x.abs() > 100.0 {
            return true;
        }

        let sigma = crate::statistical::logistic(x);
        (0.0..=1.0).contains(&sigma)
    }

    #[quickcheck]
    fn logistic_symmetry(x: f64) -> bool {
        // Handle NaN input
        if x.is_nan() || x.abs() > 100.0 {
            return true;
        }

        let sigma_x = crate::statistical::logistic(x);
        let sigma_neg_x = crate::statistical::logistic(-x);

        (sigma_x + sigma_neg_x - 1.0).abs() < 1e-12
    }

    #[quickcheck]
    fn softmax_sum_to_one(xs: Vec<f64>) -> bool {
        if xs.is_empty() || xs.len() > 100 {
            return true;
        }

        // Clamp values to reasonable range
        let _xs: Vec<f64> = xs.iter().map(|&x| x.clamp(-50.0, 50.0)).collect();

        let xs_array = scirs2_core::ndarray::Array1::from(_xs.clone());
        let softmax_result = crate::statistical::softmax(xs_array.view());
        let sum: f64 = match softmax_result {
            Ok(arr) => arr.iter().sum(),
            Err(_) => return true, // Skip failed computations
        };

        (sum - 1.0).abs() < 1e-10
    }

    #[quickcheck]
    fn logsumexp_accuracy(xs: Vec<f64>) -> bool {
        if xs.is_empty() || xs.len() > 100 {
            return true;
        }

        // Clamp to reasonable range to avoid overflow in direct calculation
        let clamped_xs: Vec<f64> = xs.iter().map(|&x| x.clamp(-100.0, 100.0)).collect();

        let xs_array = scirs2_core::ndarray::Array1::from(clamped_xs.clone());
        let lse_result = crate::statistical::logsumexp(xs_array.view());
        let lse = lse_result.unwrap_or(f64::NAN);

        // Direct calculation using clamped values (same as logsumexp input)
        let direct: f64 = clamped_xs.iter().map(|&x| x.exp()).sum::<f64>().ln();

        if !lse.is_finite() || !direct.is_finite() {
            // If direct overflows but logsumexp doesn't, that's good
            return lse.is_finite() || !direct.is_finite();
        }

        (lse - direct).abs() < 1e-8 * direct.abs().max(1.0)
    }
}

/// Run all QuickCheck property tests
#[allow(dead_code)]
pub fn run_all_quickcheck_tests() {
    println!("Running QuickCheck property tests...");

    // The tests are automatically run by cargo test
    // This function is for documentation purposes
}

#[cfg(test)]
mod integration {

    #[test]
    fn test_quickcheck_infrastructure() {
        // Basic test to ensure QuickCheck is working
        fn prop_reversing_twice_is_identity(xs: Vec<i32>) -> bool {
            let mut rev = xs.clone();
            rev.reverse();
            rev.reverse();
            xs == rev
        }

        quickcheck::quickcheck(prop_reversing_twice_is_identity as fn(Vec<i32>) -> bool);
    }
}