scirs2-stats 0.4.0

//! Enhanced property-based testing framework for mathematical invariants (v2)
//!
//! This module provides comprehensive property-based testing capabilities for
//! statistical operations, ensuring mathematical correctness, numerical stability,
//! and consistency across different computational approaches.

use crate::error::{StatsError, StatsResult};
use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2};
use scirs2_core::numeric::{Float, NumCast, One, Zero};
use scirs2_core::random::{rngs::StdRng, rng, Rng, SeedableRng};
use scirs2_core::{
    parallel_ops::*,
    simd_ops::SimdUnifiedOps,
    validation::*,
};
use std::collections::HashMap;
use std::fmt::Debug;
use std::marker::PhantomData;

/// Property-based test configuration
#[derive(Debug, Clone)]
pub struct PropertyTestConfig {
    /// Number of test cases to generate
    pub num_test_cases: usize,
    /// Random seed for reproducibility
    pub seed: Option<u64>,
    /// Tolerance for floating-point comparisons
    pub tolerance: f64,
    /// Maximum data size for generated test cases
    pub maxdatasize: usize,
    /// Minimum data size for generated test cases
    pub mindatasize: usize,
    /// Enable parallel test execution
    pub parallel_execution: bool,
    /// Test timeout in milliseconds
    pub timeout_ms: u64,
    /// Enable detailed failure reporting
    pub detailed_failures: bool,
}

impl Default for PropertyTestConfig {
    fn default() -> Self {
        Self {
            num_test_cases: 1000,
            seed: Some(42),
            tolerance: 1e-10,
            maxdatasize: 10000,
            mindatasize: 5,
            parallel_execution: true,
            timeout_ms: 30000,
            detailed_failures: true,
        }
    }
}

/// Test result status
#[derive(Debug, Clone, PartialEq)]
pub enum TestStatus {
    Pass,
    Fail(String),
    Timeout,
    Error(String),
}

/// Individual property test result
#[derive(Debug, Clone)]
pub struct PropertyTestResult {
    /// Test property name
    pub property_name: String,
    /// Test case identifier
    pub test_case_id: usize,
    /// Test status
    pub status: TestStatus,
    /// Input data that caused failure (if any)
    pub failing_input: Option<TestInput>,
    /// Expected vs actual values (if applicable)
    pub comparison: Option<(f64, f64)>,
    /// Execution time in microseconds
    pub execution_time_us: u64,
}

/// Test input data for reproducibility
#[derive(Debug, Clone)]
pub struct TestInput {
    /// Input arrays
    pub arrays: Vec<Array1<f64>>,
    /// Input matrices
    pub matrices: Vec<Array2<f64>>,
    /// Scalar parameters
    pub scalars: Vec<f64>,
    /// Boolean flags
    pub flags: Vec<bool>,
    /// String parameters
    pub strings: Vec<String>,
}

/// Comprehensive test suite result
#[derive(Debug, Clone)]
pub struct TestSuiteResult {
    /// Total number of tests run
    pub total_tests: usize,
    /// Number of passed tests
    pub passed_tests: usize,
    /// Number of failed tests
    pub failed_tests: usize,
    /// Number of timed out tests
    pub timeout_tests: usize,
    /// Number of error tests
    pub error_tests: usize,
    /// Individual test results
    pub test_results: Vec<PropertyTestResult>,
    /// Summary statistics
    pub summary: TestSummary,
}

/// Test summary statistics
#[derive(Debug, Clone)]
pub struct TestSummary {
    /// Success rate (0.0 to 1.0)
    pub success_rate: f64,
    /// Average execution time in microseconds
    pub avg_execution_time_us: f64,
    /// Maximum execution time in microseconds
    pub max_execution_time_us: u64,
    /// Minimum execution time in microseconds
    pub min_execution_time_us: u64,
    /// Properties tested
    pub properties_tested: Vec<String>,
    /// Most common failure reasons
    pub failure_reasons: HashMap<String, usize>,
}

/// Enhanced property-based test framework
pub struct PropertyBasedTestFramework<F> {
    config: PropertyTestConfig,
    rng: StdRng, phantom: PhantomData<F>,
}

impl<F> PropertyBasedTestFramework<F>
where
    F: Float + NumCast + SimdUnifiedOps + Zero + One + PartialOrd + Copy + Send + Sync + Debug
        + std::fmt::Display,
{
    /// Create new property-based test framework
    pub fn new(config: PropertyTestConfig) -> Self {
        let rng = match config.seed {
            Some(seed) => StdRng::seed_from_u64(seed),
            None => StdRng::from_rng(scirs2_core::random::thread_rng()),
        };

        Self {
            config,
            rng_phantom: PhantomData,
        }
    }

    /// Test mathematical invariants for basic statistics
    pub fn test_descriptive_statistics_invariants(&mut self) -> StatsResult<TestSuiteResult> {
        let mut results = Vec::new();
        let start_time = std::time::Instant::now();

        // Test mean properties
        results.extend(self.test_mean_properties()?);
        
        // Test variance properties
        results.extend(self.test_variance_properties()?);
        
        // Test standard deviation properties
        results.extend(self.test_std_properties()?);
        
        // Test skewness properties
        results.extend(self.test_skewness_properties()?);
        
        // Test kurtosis properties
        results.extend(self.test_kurtosis_properties()?);
        
        // Test quantile properties
        results.extend(self.test_quantile_properties()?);

        Ok(self.compile_test_results(results, start_time.elapsed()))
    }

    /// Test correlation analysis invariants
    pub fn test_correlation_invariants(&mut self) -> StatsResult<TestSuiteResult> {
        let mut results = Vec::new();
        let start_time = std::time::Instant::now();

        // Test Pearson correlation properties
        results.extend(self.test_pearson_correlation_properties()?);
        
        // Test Spearman correlation properties
        results.extend(self.test_spearman_correlation_properties()?);
        
        // Test Kendall tau properties
        results.extend(self.test_kendall_tau_properties()?);
        
        // Test correlation matrix properties
        results.extend(self.test_correlation_matrix_properties()?);

        Ok(self.compile_test_results(results, start_time.elapsed()))
    }

    /// Test regression analysis invariants
    pub fn test_regression_invariants(&mut self) -> StatsResult<TestSuiteResult> {
        let mut results = Vec::new();
        let start_time = std::time::Instant::now();

        // Test linear regression properties
        results.extend(self.test_linear_regression_properties()?);
        
        // Test polynomial regression properties
        results.extend(self.test_polynomial_regression_properties()?);
        
        // Test robust regression properties
        results.extend(self.test_robust_regression_properties()?);

        Ok(self.compile_test_results(results, start_time.elapsed()))
    }

    /// Test statistical test invariants
    pub fn test_statistical_test_invariants(&mut self) -> StatsResult<TestSuiteResult> {
        let mut results = Vec::new();
        let start_time = std::time::Instant::now();

        // Test t-test properties
        results.extend(self.test_ttest_properties()?);
        
        // Test ANOVA properties
        results.extend(self.test_anova_properties()?);
        
        // Test non-parametric test properties
        results.extend(self.test_nonparametric_properties()?);
        
        // Test normality test properties
        results.extend(self.test_normality_test_properties()?);

        Ok(self.compile_test_results(results, start_time.elapsed()))
    }

    /// Test SIMD vs scalar consistency
    pub fn test_simd_scalar_consistency(&mut self) -> StatsResult<TestSuiteResult> {
        let mut results = Vec::new();
        let start_time = std::time::Instant::now();

        // Test SIMD mean vs scalar mean
        results.extend(self.test_simd_vs_scalar_mean()?);
        
        // Test SIMD variance vs scalar variance
        results.extend(self.test_simd_vs_scalar_variance()?);
        
        // Test SIMD correlation vs scalar correlation
        results.extend(self.test_simd_vs_scalar_correlation()?);

        Ok(self.compile_test_results(results, start_time.elapsed()))
    }

    /// Test parallel vs sequential consistency
    pub fn test_parallel_sequential_consistency(&mut self) -> StatsResult<TestSuiteResult> {
        let mut results = Vec::new();
        let start_time = std::time::Instant::now();

        // Test parallel mean vs sequential mean
        results.extend(self.test_parallel_vs_sequential_mean()?);
        
        // Test parallel correlation vs sequential correlation
        results.extend(self.test_parallel_vs_sequential_correlation()?);
        
        // Test parallel bootstrap vs sequential bootstrap
        results.extend(self.test_parallel_vs_sequentialbootstrap()?);

        Ok(self.compile_test_results(results, start_time.elapsed()))
    }

    /// Test numerical stability properties
    pub fn test_numerical_stability(&mut self) -> StatsResult<TestSuiteResult> {
        let mut results = Vec::new();
        let start_time = std::time::Instant::now();

        // Test with extreme values
        results.extend(self.test_extreme_values_stability()?);
        
        // Test with near-zero values
        results.extend(self.test_near_zero_stability()?);
        
        // Test with large values
        results.extend(self.test_large_values_stability()?);
        
        // Test with ill-conditioned data
        results.extend(self.test_ill_conditioned_stability()?);

        Ok(self.compile_test_results(results, start_time.elapsed()))
    }

    // Implementation of individual property tests

    fn test_mean_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        let mut results = Vec::new();

        for test_case_id in 0..self.config.num_test_cases {
            let start_time = std::time::Instant::now();
            
            // Generate test data
            let data = self.generate_random_array()?;
            let input = TestInput {
                arrays: vec![data.clone()],
                matrices: vec![],
                scalars: vec![],
                flags: vec![],
                strings: vec![],
            };

            // Test mean invariant: mean of constants should equal the constant
            let constant_value = 5.0;
            let constantdata = Array1::from_elem(data.len(), constant_value);
            
            let result = match crate::descriptive::mean(&constantdata.view()) {
                Ok(computed_mean) => {
                    let diff = (computed_mean - constant_value).abs();
                    if diff < self.config.tolerance {
                        PropertyTestResult {
                            property_name: "mean_of_constants".to_string(),
                            test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((computed_mean, constant_value)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "mean_of_constants".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "Mean of constants failed: expected {}, got {}, diff: {}",
                                constant_value, computed_mean, diff
                            )),
                            failing_input: Some(input),
                            comparison: Some((computed_mean, constant_value)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    }
                }
                Err(e) => PropertyTestResult {
                    property_name: "mean_of_constants".to_string(),
                    test_case_id,
                    status: TestStatus::Error(format!("Error computing mean: {}", e)),
                    failing_input: Some(input),
                    comparison: None,
                    execution_time_us: start_time.elapsed().as_micros() as u64,
                },
            };
            
            results.push(result);

            // Test linearity: mean(a*X + b) = a*mean(X) + b
            if let Ok(original_mean) = crate::descriptive::mean(&data.view()) {
                let a = self.rng.random_range(0.1..10.0);
                let b = self.rng.random_range(-5.0..5.0);
                
                let transformeddata = data.mapv(|x| a * x + b);
                
                if let Ok(transformed_mean) = crate::descriptive::mean(&transformeddata.view()) {
                    let expected_mean = a * original_mean + b;
                    let diff = (transformed_mean - expected_mean).abs();
                    
                    let result = if diff < self.config.tolerance {
                        PropertyTestResult {
                            property_name: "mean_linearity".to_string()..test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((transformed_mean, expected_mean)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "mean_linearity".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "Mean linearity failed: expected {}, got {}, diff: {}",
                                expected_mean, transformed_mean, diff
                            )),
                            failing_input: Some(TestInput {
                                arrays: vec![data.clone(), transformeddata],
                                matrices: vec![],
                                scalars: vec![a, b, original_mean],
                                flags: vec![],
                                strings: vec![],
                            }),
                            comparison: Some((transformed_mean, expected_mean)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    };
                    
                    results.push(result);
                }
            }
        }

        Ok(results)
    }

    fn test_variance_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        let mut results = Vec::new();

        for test_case_id in 0..self.config.num_test_cases {
            let start_time = std::time::Instant::now();
            
            let data = self.generate_random_array()?;
            
            // Test variance of constants should be zero
            let constant_value = 3.0;
            let constantdata = Array1::from_elem(data.len(), constant_value);
            
            let result = match crate::descriptive::var(&constantdata.view(), 1, None) {
                Ok(computed_variance) => {
                    if computed_variance.abs() < self.config.tolerance {
                        PropertyTestResult {
                            property_name: "variance_of_constants".to_string(),
                            test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((computed_variance, 0.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "variance_of_constants".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "Variance of constants should be zero, got: {}",
                                computed_variance
                            )),
                            failing_input: Some(TestInput {
                                arrays: vec![constantdata],
                                matrices: vec![],
                                scalars: vec![constant_value],
                                flags: vec![],
                                strings: vec![],
                            }),
                            comparison: Some((computed_variance, 0.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    }
                }
                Err(e) => PropertyTestResult {
                    property_name: "variance_of_constants".to_string(),
                    test_case_id,
                    status: TestStatus::Error(format!("Error computing variance: {}", e)),
                    failing_input: None,
                    comparison: None,
                    execution_time_us: start_time.elapsed().as_micros() as u64,
                },
            };
            
            results.push(result);

            // Test variance scaling: var(a*X) = a²*var(X)
            if let Ok(original_var) = crate::descriptive::var(&data.view(), 1, None) {
                let a = self.rng.random_range(0.1..5.0);
                let scaleddata = data.mapv(|x| a * x);
                
                if let Ok(scaled_var) = crate::descriptive::var(&scaleddata.view()..1, None) {
                    let expected_var = a * a * original_var;
                    let diff = (scaled_var - expected_var).abs();
                    
                    let result = if diff < self.config.tolerance * expected_var.abs().max(1.0) {
                        PropertyTestResult {
                            property_name: "variance_scaling".to_string(),
                            test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((scaled_var, expected_var)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "variance_scaling".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "Variance scaling failed: expected {}, got {}, diff: {}",
                                expected_var, scaled_var, diff
                            )),
                            failing_input: Some(TestInput {
                                arrays: vec![data.clone(), scaleddata],
                                matrices: vec![],
                                scalars: vec![a, original_var],
                                flags: vec![],
                                strings: vec![],
                            }),
                            comparison: Some((scaled_var, expected_var)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    };
                    
                    results.push(result);
                }
            }
        }

        Ok(results)
    }

    fn test_std_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        let mut results = Vec::new();

        for test_case_id in 0..self.config.num_test_cases {
            let start_time = std::time::Instant::now();
            
            let data = self.generate_random_array()?;
            
            // Test that std = sqrt(variance)
            let variance_result = crate::descriptive::var(&data.view(), 1, None);
            let std_result = crate::descriptive::std(&data.view(), 1, None);
            
            let result = match (variance_result, std_result) {
                (Ok(variance), Ok(std_dev)) => {
                    let expected_std = variance.sqrt();
                    let diff = (std_dev - expected_std).abs();
                    
                    if diff < self.config.tolerance * expected_std.max(1.0) {
                        PropertyTestResult {
                            property_name: "std_sqrt_variance".to_string(),
                            test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((std_dev, expected_std)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "std_sqrt_variance".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "std ≠ sqrt(variance): std={}, sqrt(var)={}, diff={}",
                                std_dev, expected_std, diff
                            )),
                            failing_input: Some(TestInput {
                                arrays: vec![data.clone()],
                                matrices: vec![],
                                scalars: vec![variance, std_dev],
                                flags: vec![],
                                strings: vec![],
                            }),
                            comparison: Some((std_dev, expected_std)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    }
                }
                _ => PropertyTestResult {
                    property_name: "std_sqrt_variance".to_string(),
                    test_case_id,
                    status: TestStatus::Error("Failed to compute variance or std".to_string()),
                    failing_input: None,
                    comparison: None,
                    execution_time_us: start_time.elapsed().as_micros() as u64,
                },
            };
            
            results.push(result);
        }

        Ok(results)
    }

    fn test_skewness_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        let mut results = Vec::new();

        for test_case_id in 0..self.config.num_test_cases {
            let start_time = std::time::Instant::now();
            
            // Generate symmetric data around zero
            let n = self.rng.random_range(self.config.mindatasize..self.config.maxdatasize + 1);
            let mut data = Vec::new();
            
            for _ in 0..n/2 {
                let value = self.rng.random_range(-5.0..5.0);
                data.push(value);
                data.push(-value); // Add symmetric value
            }
            
            if n % 2 == 1 {
                data.push(0.0); // Add center point for odd sizes
            }
            
            let data_array = Array1::from_vec(data);
            
            // Test that symmetric data should have near-zero skewness
            let result = match crate::descriptive::skew(&data_array.view()..false, None) {
                Ok(skewness) => {
                    if skewness.abs() < self.config.tolerance * 10.0 { // Allow some tolerance for finite samples
                        PropertyTestResult {
                            property_name: "symmetricdata_skewness".to_string(),
                            test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((skewness, 0.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "symmetricdata_skewness".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "Symmetric data should have near-zero skewness, got: {}",
                                skewness
                            )),
                            failing_input: Some(TestInput {
                                arrays: vec![data_array],
                                matrices: vec![],
                                scalars: vec![skewness],
                                flags: vec![],
                                strings: vec![],
                            }),
                            comparison: Some((skewness, 0.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    }
                }
                Err(e) => PropertyTestResult {
                    property_name: "symmetricdata_skewness".to_string(),
                    test_case_id,
                    status: TestStatus::Error(format!("Error computing skewness: {}", e)),
                    failing_input: None,
                    comparison: None,
                    execution_time_us: start_time.elapsed().as_micros() as u64,
                },
            };
            
            results.push(result);
        }

        Ok(results)
    }

    fn test_kurtosis_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        let mut results = Vec::new();

        for test_case_id in 0..self.config.num_test_cases {
            let start_time = std::time::Instant::now();
            
            // Generate normal-like data (should have kurtosis ≈ 0 for Fisher definition)
            let n = self.rng.random_range(self.config.mindatasize..self.config.maxdatasize + 1);
            let data: Vec<f64> = (0..n)
                .map(|_| {
                    // Box-Muller transform for normal distribution
                    let u1: f64 = self.rng.random();
                    let u2: f64 = self.rng.random();
                    (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
                })
                .collect();
            
            let data_array = Array1::from_vec(data);
            
            // Test that normal data has kurtosis ≈ 0 (Fisher definition)
            let result = match crate::descriptive::kurtosis(&data_array.view(), true, false, None) {
                Ok(kurtosis_val) => {
                    // Allow larger tolerance for finite samples of normal distribution
                    if kurtosis_val.abs() < 2.0 { 
                        PropertyTestResult {
                            property_name: "normaldata_kurtosis".to_string(),
                            test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((kurtosis_val, 0.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "normaldata_kurtosis".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "Normal data should have kurtosis ≈ 0, got: {}",
                                kurtosis_val
                            )),
                            failing_input: Some(TestInput {
                                arrays: vec![data_array],
                                matrices: vec![],
                                scalars: vec![kurtosis_val],
                                flags: vec![],
                                strings: vec![],
                            }),
                            comparison: Some((kurtosis_val, 0.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    }
                }
                Err(e) => PropertyTestResult {
                    property_name: "normaldata_kurtosis".to_string(),
                    test_case_id,
                    status: TestStatus::Error(format!("Error computing kurtosis: {}", e)),
                    failing_input: None,
                    comparison: None,
                    execution_time_us: start_time.elapsed().as_micros() as u64,
                },
            };
            
            results.push(result);
        }

        Ok(results)
    }

    fn test_quantile_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        let mut results = Vec::new();

        for test_case_id in 0..self.config.num_test_cases {
            let start_time = std::time::Instant::now();
            
            let data = self.generate_random_array()?;
            
            // Test that quantiles are monotonic
            let quantiles = [0.25, 0.5, 0.75];
            let mut computed_quantiles = Vec::new();
            
            for &q in &quantiles {
                if let Ok(quantile_val) = crate::quantile::quantile(&data.view(), q, crate::quantile::QuantileInterpolation::Linear) {
                    computed_quantiles.push(quantile_val);
                }
            }
            
            let result = if computed_quantiles.len() == 3 {
                let monotonic = computed_quantiles[0] <= computed_quantiles[1] && 
                               computed_quantiles[1] <= computed_quantiles[2];
                
                if monotonic {
                    PropertyTestResult {
                        property_name: "quantiles_monotonic".to_string(),
                        test_case_id,
                        status: TestStatus::Pass,
                        failing_input: None,
                        comparison: None,
                        execution_time_us: start_time.elapsed().as_micros() as u64,
                    }
                } else {
                    PropertyTestResult {
                        property_name: "quantiles_monotonic".to_string(),
                        test_case_id,
                        status: TestStatus::Fail(format!(
                            "Quantiles not monotonic: Q25={}, Q50={}, Q75={}",
                            computed_quantiles[0], computed_quantiles[1], computed_quantiles[2]
                        )),
                        failing_input: Some(TestInput {
                            arrays: vec![data.clone()],
                            matrices: vec![],
                            scalars: computed_quantiles,
                            flags: vec![],
                            strings: vec![],
                        }),
                        comparison: None,
                        execution_time_us: start_time.elapsed().as_micros() as u64,
                    }
                }
            } else {
                PropertyTestResult {
                    property_name: "quantiles_monotonic".to_string(),
                    test_case_id,
                    status: TestStatus::Error("Failed to compute all quantiles".to_string()),
                    failing_input: None,
                    comparison: None,
                    execution_time_us: start_time.elapsed().as_micros() as u64,
                }
            };
            
            results.push(result);
        }

        Ok(results)
    }

    // Placeholder implementations for other test methods
    fn test_pearson_correlation_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        let mut results = Vec::new();

        for test_case_id in 0..self.config.num_test_cases {
            let start_time = std::time::Instant::now();
            
            let data = self.generate_random_array()?;
            
            // Test correlation of data with itself should be 1.0
            let result = match crate::correlation::pearson_r(&data.view(), &data.view()) {
                Ok(correlation) => {
                    let diff = (correlation - 1.0).abs();
                    if diff < self.config.tolerance {
                        PropertyTestResult {
                            property_name: "pearson_self_correlation".to_string(),
                            test_case_id,
                            status: TestStatus::Pass,
                            failing_input: None,
                            comparison: Some((correlation, 1.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    } else {
                        PropertyTestResult {
                            property_name: "pearson_self_correlation".to_string(),
                            test_case_id,
                            status: TestStatus::Fail(format!(
                                "Self-correlation should be 1.0, got: {}",
                                correlation
                            )),
                            failing_input: Some(TestInput {
                                arrays: vec![data.clone()],
                                matrices: vec![],
                                scalars: vec![correlation],
                                flags: vec![],
                                strings: vec![],
                            }),
                            comparison: Some((correlation, 1.0)),
                            execution_time_us: start_time.elapsed().as_micros() as u64,
                        }
                    }
                }
                Err(e) => PropertyTestResult {
                    property_name: "pearson_self_correlation".to_string(),
                    test_case_id,
                    status: TestStatus::Error(format!("Error computing correlation: {}", e)),
                    failing_input: None,
                    comparison: None,
                    execution_time_us: start_time.elapsed().as_micros() as u64,
                },
            };
            
            results.push(result);
        }

        Ok(results)
    }

    // Simplified implementations for other methods
    fn test_spearman_correlation_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_kendall_tau_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_correlation_matrix_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_linear_regression_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_polynomial_regression_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_robust_regression_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_ttest_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_anova_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_nonparametric_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_normality_test_properties(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_simd_vs_scalar_mean(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_simd_vs_scalar_variance(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_simd_vs_scalar_correlation(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_parallel_vs_sequential_mean(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_parallel_vs_sequential_correlation(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_parallel_vs_sequentialbootstrap(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_extreme_values_stability(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_near_zero_stability(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_large_values_stability(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    fn test_ill_conditioned_stability(&mut self) -> StatsResult<Vec<PropertyTestResult>> {
        Ok(vec![]) // Placeholder
    }

    // Helper methods

    fn generate_random_array(&mut self) -> StatsResult<Array1<f64>> {
        let size = self.rng.random_range(self.config.mindatasize..self.config.maxdatasize + 1);
        let data: Vec<f64> = (0..size)
            .map(|_| self.rng.random_range(-100.0..100.0))
            .collect();
        Ok(Array1::from_vec(data))
    }

    fn compile_test_results(
        &self..results: Vec<PropertyTestResult>,
        total_duration: std::time::Duration,) -> TestSuiteResult {
        let total_tests = results.len();
        let passed_tests = results.iter().filter(|r| r.status == TestStatus::Pass).count();
        let failed_tests = results.iter().filter(|r| matches!(r.status, TestStatus::Fail(_))).count();
        let timeout_tests = results.iter().filter(|r| r.status == TestStatus::Timeout).count();
        let error_tests = results.iter().filter(|r| matches!(r.status, TestStatus::Error(_))).count();

        let success_rate = if total_tests > 0 {
            passed_tests as f64 / total_tests as f64
        } else {
            0.0
        };

        let execution_times: Vec<u64> = results.iter().map(|r| r.execution_time_us).collect();
        let avg_execution_time_us = if !execution_times.is_empty() {
            execution_times.iter().sum::<u64>() as f64 / execution_times.len() as f64
        } else {
            0.0
        };
        let max_execution_time_us = execution_times.iter().copied().max().unwrap_or(0);
        let min_execution_time_us = execution_times.iter().copied().min().unwrap_or(0);

        let mut properties_tested = Vec::new();
        let mut failure_reasons = HashMap::new();

        for result in &results {
            if !properties_tested.contains(&result.property_name) {
                properties_tested.push(result.property_name.clone());
            }

            if let TestStatus::Fail(reason) = &result.status {
                *failure_reasons.entry(reason.clone()).or_insert(0) += 1;
            }
        }

        let summary = TestSummary {
            success_rate,
            avg_execution_time_us,
            max_execution_time_us,
            min_execution_time_us,
            properties_tested,
            failure_reasons,
        };

        TestSuiteResult {
            total_tests,
            passed_tests,
            failed_tests,
            timeout_tests,
            error_tests,
            test_results: results,
            summary,
        }
    }
}

/// Convenience functions for property-based testing
#[allow(dead_code)]
pub fn test_basic_statistics_properties() -> StatsResult<TestSuiteResult> {
    let config = PropertyTestConfig::default();
    let mut framework = PropertyBasedTestFramework::<f64>::new(config);
    framework.test_descriptive_statistics_invariants()
}

#[allow(dead_code)]
pub fn test_correlation_properties() -> StatsResult<TestSuiteResult> {
    let config = PropertyTestConfig::default();
    let mut framework = PropertyBasedTestFramework::<f64>::new(config);
    framework.test_correlation_invariants()
}

#[allow(dead_code)]
pub fn test_all_mathematical_invariants() -> StatsResult<Vec<TestSuiteResult>> {
    let config = PropertyTestConfig::default();
    let mut framework = PropertyBasedTestFramework::<f64>::new(config);

    let mut all_results = Vec::new();
    
    all_results.push(framework.test_descriptive_statistics_invariants()?);
    all_results.push(framework.test_correlation_invariants()?);
    all_results.push(framework.test_regression_invariants()?);
    all_results.push(framework.test_statistical_test_invariants()?);
    all_results.push(framework.test_simd_scalar_consistency()?);
    all_results.push(framework.test_parallel_sequential_consistency()?);
    all_results.push(framework.test_numerical_stability()?);

    Ok(all_results)
}

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

    #[test]
    fn test_property_test_framework() {
        let config = PropertyTestConfig {
            num_test_cases: 10,
            ..Default::default()
        };
        
        let mut framework = PropertyBasedTestFramework::<f64>::new(config);
        let result = framework.test_descriptive_statistics_invariants();
        
        assert!(result.is_ok());
        let suite_result = result.expect("Operation failed");
        assert!(suite_result.total_tests > 0);
        assert!(suite_result.summary.success_rate >= 0.0 && suite_result.summary.success_rate <= 1.0);
    }

    #[test]
    fn test_mean_properties() {
        let result = test_basic_statistics_properties();
        assert!(result.is_ok());
        
        let suite_result = result.expect("Operation failed");
        assert!(suite_result.total_tests > 0);
    }

    #[test]
    fn test_correlation_properties_basic() {
        let result = test_correlation_properties();
        assert!(result.is_ok());
        
        let suite_result = result.expect("Operation failed");
        assert!(suite_result.total_tests >= 0); // May be 0 if simplified implementations
    }

    #[test]
    fn test_test_input_creation() {
        let input = TestInput {
            arrays: vec![Array1::from_vec(vec![1.0, 2.0, 3.0])],
            matrices: vec![Array2::eye(3)],
            scalars: vec![42.0],
            flags: vec![true, false],
            strings: vec!["test".to_string()],
        };
        
        assert_eq!(input.arrays.len(), 1);
        assert_eq!(input.matrices.len(), 1);
        assert_eq!(input.scalars.len(), 1);
        assert_eq!(input.flags.len(), 2);
        assert_eq!(input.strings.len(), 1);
    }

    #[test]
    fn test_config_validation() {
        let config = PropertyTestConfig {
            num_test_cases: 0,
            mindatasize: 10,
            maxdatasize: 5, // Invalid: max < min
            ..Default::default()
        };
        
        // The framework should handle invalid configurations gracefully
        let framework = PropertyBasedTestFramework::<f64>::new(config);
        assert!(framework.config.num_test_cases == 0);
    }
}