scirs2-stats 0.4.1

//! Advanced-parallel statistical computing framework for scirs2-stats v1.0.0
//!
//! This module provides an advanced parallel processing system that goes beyond
//! basic parallelization to offer intelligent work distribution, adaptive load
//! balancing, NUMA-aware processing, and hybrid CPU-GPU computation strategies.

use crate::error::{StatsError, StatsResult};
use scirs2_core::ndarray::{s, Array2, ArrayView1, ArrayView2};
use scirs2_core::numeric::{Float, NumCast, Zero};
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::sync::{Arc, Mutex, RwLock};
use std::thread;
use std::time::{Duration, Instant};

/// Configuration for advanced-parallel statistical operations
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct AdvancedParallelConfig {
    /// Enable adaptive work distribution
    pub adaptive_work_distribution: bool,
    /// Enable NUMA-aware processing
    pub numa_aware: bool,
    /// Enable hybrid CPU-GPU processing
    pub hybrid_processing: bool,
    /// Work stealing strategy
    pub work_stealing: WorkStealingStrategy,
    /// Load balancing algorithm
    pub load_balancing: LoadBalancingAlgorithm,
    /// Thread pool configuration
    pub thread_pool_config: ThreadPoolConfig,
    /// Memory management strategy
    pub memory_strategy: ParallelMemoryStrategy,
    /// Performance monitoring
    pub performance_monitoring: bool,
    /// Minimum work size per thread
    pub min_worksize: usize,
    /// Maximum parallelization depth
    pub max_parallel_depth: usize,
    /// Enable vectorized operations within parallel tasks
    pub enable_simd_in_parallel: bool,
    /// Cache optimization level
    pub cache_optimization: CacheOptimizationLevel,
}

impl Default for AdvancedParallelConfig {
    fn default() -> Self {
        Self {
            adaptive_work_distribution: true,
            numa_aware: true,
            hybrid_processing: false,
            work_stealing: WorkStealingStrategy::Adaptive,
            load_balancing: LoadBalancingAlgorithm::DynamicRoundRobin,
            thread_pool_config: ThreadPoolConfig::default(),
            memory_strategy: ParallelMemoryStrategy::CacheAware,
            performance_monitoring: true,
            min_worksize: 1000,
            max_parallel_depth: 3,
            enable_simd_in_parallel: true,
            cache_optimization: CacheOptimizationLevel::Aggressive,
        }
    }
}

/// Work stealing strategies
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum WorkStealingStrategy {
    /// No work stealing
    None,
    /// Simple random work stealing
    Random,
    /// Locality-aware work stealing
    LocalityAware,
    /// Adaptive work stealing based on performance
    Adaptive,
    /// NUMA-topology aware work stealing
    NumaAware,
}

/// Load balancing algorithms
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum LoadBalancingAlgorithm {
    /// Static round-robin distribution
    StaticRoundRobin,
    /// Dynamic round-robin with feedback
    DynamicRoundRobin,
    /// Work-based load balancing
    WorkBased,
    /// Performance-based load balancing
    PerformanceBased,
    /// Hierarchical load balancing
    Hierarchical,
    /// Machine learning-based load balancing
    MLBased,
}

/// Thread pool configuration
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ThreadPoolConfig {
    /// Number of worker threads
    pub num_workers: Option<usize>,
    /// Thread priority
    pub thread_priority: ThreadPriority,
    /// Thread affinity strategy
    pub affinity_strategy: ThreadAffinityStrategy,
    /// Stack size per thread
    pub stacksize: Option<usize>,
    /// Idle thread timeout
    pub idle_timeout: Duration,
    /// Thread pool scaling strategy
    pub scaling_strategy: ScalingStrategy,
}

impl Default for ThreadPoolConfig {
    fn default() -> Self {
        Self {
            num_workers: None, // Auto-detect
            thread_priority: ThreadPriority::Normal,
            affinity_strategy: ThreadAffinityStrategy::NUMA,
            stacksize: Some(2 * 1024 * 1024), // 2MB
            idle_timeout: Duration::from_secs(60),
            scaling_strategy: ScalingStrategy::Adaptive,
        }
    }
}

/// Thread priority levels
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum ThreadPriority {
    Low,
    Normal,
    High,
    RealTime,
}

/// Thread affinity strategies
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum ThreadAffinityStrategy {
    /// No affinity setting
    None,
    /// NUMA-aware affinity
    NUMA,
    /// Core-based affinity
    CoreBased,
    /// Custom affinity pattern
    Custom(Vec<usize>),
}

/// Thread pool scaling strategies
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum ScalingStrategy {
    /// Fixed number of threads
    Fixed,
    /// Adaptive scaling based on workload
    Adaptive,
    /// Performance-based scaling
    PerformanceBased,
    /// Resource-aware scaling
    ResourceAware,
}

/// Parallel memory management strategies
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum ParallelMemoryStrategy {
    /// Simple memory allocation
    Simple,
    /// Cache-aware memory allocation
    CacheAware,
    /// NUMA-aware memory allocation
    NumaAware,
    /// Memory pool-based allocation
    PoolBased,
    /// Lock-free memory management
    LockFree,
}

/// Cache optimization levels
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
pub enum CacheOptimizationLevel {
    /// No cache optimization
    None,
    /// Basic cache-friendly patterns
    Basic,
    /// Aggressive cache optimization
    Aggressive,
    /// Hardware-specific optimization
    HardwareSpecific,
}

/// Work unit for parallel processing
#[derive(Debug, Clone)]
pub struct WorkUnit<T> {
    /// Work identifier
    pub id: usize,
    /// Work data
    pub data: T,
    /// Estimated computational cost
    pub cost: f64,
    /// Dependencies on other work units
    pub dependencies: Vec<usize>,
    /// Priority level
    pub priority: WorkPriority,
    /// NUMA node preference
    pub numa_node: Option<usize>,
}

/// Work priority levels
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub enum WorkPriority {
    Low = 0,
    Normal = 1,
    High = 2,
    Critical = 3,
}

/// Parallel execution context
#[derive(Debug, Clone)]
pub struct ParallelExecutionContext {
    /// Thread ID
    pub thread_id: usize,
    /// NUMA node ID
    pub numa_node: usize,
    /// Local memory pool
    pub memory_pool: Option<Arc<Mutex<MemoryPool>>>,
    /// Performance counters
    pub counters: PerformanceCounters,
    /// Thread-local cache
    pub cache: ThreadLocalCache,
}

/// Performance counters for monitoring
#[derive(Debug, Clone, Default)]
pub struct PerformanceCounters {
    /// Tasks completed
    pub tasks_completed: usize,
    /// Total execution time
    pub total_time: Duration,
    /// Cache hits
    pub cache_hits: usize,
    /// Cache misses
    pub cache_misses: usize,
    /// Work stolen from other threads
    pub work_stolen: usize,
    /// Work given to other threads
    pub work_given: usize,
    /// Memory allocations
    pub memory_allocations: usize,
    /// Memory deallocations
    pub memory_deallocations: usize,
}

/// Thread-local cache for frequently used data
#[derive(Debug, Clone)]
pub struct ThreadLocalCache {
    /// Cached computation results
    pub results: HashMap<String, CachedResult>,
    /// Cache hit statistics
    pub stats: CacheStatistics,
}

/// Cached computation result
#[derive(Debug, Clone)]
pub struct CachedResult {
    /// Cached value
    pub value: Vec<f64>,
    /// Timestamp when cached
    pub timestamp: Instant,
    /// Access count
    pub access_count: usize,
    /// Cost to recompute
    pub recompute_cost: f64,
}

/// Cache statistics
#[derive(Debug, Clone, Default)]
pub struct CacheStatistics {
    /// Total cache accesses
    pub total_accesses: usize,
    /// Cache hits
    pub hits: usize,
    /// Cache misses
    pub misses: usize,
    /// Evictions
    pub evictions: usize,
}

/// Memory pool for efficient allocation
#[derive(Debug)]
pub struct MemoryPool {
    /// Pool of pre-allocated memory blocks
    blocks: Vec<Vec<u8>>,
    /// Available block indices
    available: Vec<usize>,
    /// Block size
    blocksize: usize,
    /// Total allocations
    total_allocations: usize,
}

impl MemoryPool {
    /// Create new memory pool
    pub fn new(_num_blocks: usize, blocksize: usize) -> Self {
        let mut _blocks = Vec::with_capacity(_num_blocks);
        let mut available = Vec::with_capacity(_num_blocks);

        for i in 0.._num_blocks {
            blocks.push(vec![0u8; blocksize]);
            available.push(i);
        }

        Self {
            blocks,
            available,
            blocksize,
            total_allocations: 0,
        }
    }

    /// Allocate a block
    pub fn allocate(&mut self) -> Option<*mut u8> {
        if let Some(index) = self.available.pop() {
            self.total_allocations += 1;
            Some(self.blocks[index].as_mut_ptr())
        } else {
            None
        }
    }

    /// Deallocate a block
    pub fn deallocate(&mut self, ptr: *mut u8) {
        // Find which block this pointer belongs to
        for (i, block) in self.blocks.iter().enumerate() {
            if ptr == block.as_ptr() as *mut u8 {
                self.available.push(i);
                break;
            }
        }
    }
}

/// Parallel computation result with detailed metrics
#[derive(Debug, Clone)]
pub struct AdvancedParallelResult<T> {
    /// Computation result
    pub result: T,
    /// Execution metrics
    pub metrics: ParallelExecutionMetrics,
    /// Performance analysis
    pub analysis: ParallelPerformanceAnalysis,
    /// Resource utilization
    pub resource_utilization: ResourceUtilization,
}

/// Detailed execution metrics
#[derive(Debug, Clone)]
pub struct ParallelExecutionMetrics {
    /// Total execution time
    pub total_time: Duration,
    /// Time spent in parallel execution
    pub parallel_time: Duration,
    /// Time spent in sequential execution
    pub sequential_time: Duration,
    /// Time spent in synchronization
    pub sync_time: Duration,
    /// Number of threads used
    pub threads_used: usize,
    /// Load balancing efficiency
    pub load_balance_efficiency: f64,
    /// Parallel efficiency
    pub parallel_efficiency: f64,
    /// Speedup achieved
    pub speedup: f64,
    /// Work distribution quality
    pub work_distribution_quality: f64,
}

/// Performance analysis results
#[derive(Debug, Clone)]
pub struct ParallelPerformanceAnalysis {
    /// Bottleneck analysis
    pub bottlenecks: Vec<PerformanceBottleneck>,
    /// Scaling characteristics
    pub scaling_analysis: ScalingAnalysis,
    /// Optimization opportunities
    pub optimization_opportunities: Vec<OptimizationOpportunity>,
    /// Performance rating
    pub performance_rating: PerformanceRating,
}

/// Performance bottleneck identification
#[derive(Debug, Clone)]
pub struct PerformanceBottleneck {
    /// Bottleneck type
    pub bottleneck_type: BottleneckType,
    /// Severity (0.0-1.0)
    pub severity: f64,
    /// Description
    pub description: String,
    /// Suggested mitigation
    pub mitigation: String,
}

/// Types of performance bottlenecks
#[derive(Debug, Clone)]
pub enum BottleneckType {
    /// Memory bandwidth limitation
    MemoryBandwidth,
    /// Cache thrashing
    CacheContention,
    /// Load imbalance
    LoadImbalance,
    /// Synchronization overhead
    SynchronizationOverhead,
    /// NUMA effects
    NumaEffects,
    /// False sharing
    FalseSharing,
    /// Context switching overhead
    ContextSwitching,
    /// Insufficient parallelism
    InsufficientParallelism,
}

/// Scaling analysis results
#[derive(Debug, Clone)]
pub struct ScalingAnalysis {
    /// Theoretical maximum speedup
    pub theoretical_max_speedup: f64,
    /// Achieved speedup
    pub achieved_speedup: f64,
    /// Parallel fraction (Amdahl's law)
    pub parallel_fraction: f64,
    /// Serial bottleneck impact
    pub serial_bottleneck_impact: f64,
    /// Scaling efficiency by thread count
    pub scaling_efficiency: HashMap<usize, f64>,
    /// Optimal thread count
    pub optimal_thread_count: usize,
}

/// Optimization opportunity
#[derive(Debug, Clone)]
pub struct OptimizationOpportunity {
    /// Opportunity type
    pub opportunity_type: OptimizationType,
    /// Potential improvement
    pub potential_improvement: f64,
    /// Implementation complexity
    pub complexity: OptimizationComplexity,
    /// Description
    pub description: String,
    /// Implementation steps
    pub implementation_steps: Vec<String>,
}

/// Types of optimization opportunities
#[derive(Debug, Clone)]
pub enum OptimizationType {
    /// Better work distribution
    WorkDistribution,
    /// Memory layout optimization
    MemoryLayout,
    /// Cache optimization
    CacheOptimization,
    /// SIMD integration
    SimdIntegration,
    /// Algorithm selection
    AlgorithmSelection,
    /// Resource allocation
    ResourceAllocation,
}

/// Optimization complexity levels
#[derive(Debug, Clone)]
pub enum OptimizationComplexity {
    Low,
    Medium,
    High,
    VeryHigh,
}

/// Performance rating
#[derive(Debug, Clone)]
pub enum PerformanceRating {
    Excellent,
    Good,
    Acceptable,
    Poor,
    Unacceptable,
}

/// Performance statistics summary
#[derive(Debug, Clone)]
pub struct PerformanceStatistics {
    /// Total number of operations tracked
    pub total_operations: usize,
    /// Average speedup achieved
    pub average_speedup: f64,
    /// Best performing strategies
    pub best_strategies: Vec<String>,
    /// Hardware utilization metrics
    pub hardware_utilization: HardwareUtilization,
}

/// Hardware utilization metrics
#[derive(Debug, Clone)]
pub struct HardwareUtilization {
    /// SIMD utilization percentage
    pub simd_utilization: f64,
    /// Memory bandwidth utilization
    pub memory_bandwidth_utilization: f64,
    /// Cache efficiency
    pub cache_efficiency: f64,
    /// Energy efficiency (if available)
    pub energy_efficiency: Option<f64>,
}

/// Resource utilization metrics
#[derive(Debug, Clone)]
pub struct ResourceUtilization {
    /// CPU utilization per core
    pub cpu_utilization: Vec<f64>,
    /// Memory utilization
    pub memory_utilization: f64,
    /// Cache utilization
    pub cache_utilization: CacheUtilization,
    /// NUMA node utilization
    pub numa_utilization: Vec<f64>,
    /// Energy consumption (if available)
    pub energy_consumption: Option<f64>,
}

/// Cache utilization metrics
#[derive(Debug, Clone)]
pub struct CacheUtilization {
    /// L1 cache hit rate
    pub l1_hit_rate: f64,
    /// L2 cache hit rate
    pub l2_hit_rate: f64,
    /// L3 cache hit rate
    pub l3_hit_rate: f64,
    /// Cache line utilization
    pub cache_line_utilization: f64,
}

/// Main advanced-parallel statistical processor
pub struct AdvancedParallelStatsProcessor {
    config: AdvancedParallelConfig,
    execution_contexts: Vec<Arc<RwLock<ParallelExecutionContext>>>,
    work_queue: Arc<Mutex<Vec<WorkUnit<Vec<f64>>>>>,
    performance_history: Arc<Mutex<Vec<ParallelExecutionMetrics>>>,
    optimization_cache: Arc<RwLock<HashMap<String, OptimizationStrategy>>>,
}

/// Optimization strategy cache
#[derive(Debug, Clone)]
pub struct OptimizationStrategy {
    /// Strategy name
    pub name: String,
    /// Thread count
    pub thread_count: usize,
    /// Work distribution method
    pub work_distribution: WorkDistributionMethod,
    /// Memory layout
    pub memory_layout: MemoryLayoutStrategy,
    /// Expected performance
    pub expected_performance: f64,
}

/// Work distribution methods
#[derive(Debug, Clone)]
pub enum WorkDistributionMethod {
    /// Equal-sized chunks
    EqualChunks,
    /// Size-based distribution
    SizeBased,
    /// Cost-based distribution
    CostBased,
    /// Adaptive distribution
    Adaptive,
    /// Locality-aware distribution
    LocalityAware,
}

/// Memory layout strategies
#[derive(Debug, Clone)]
pub enum MemoryLayoutStrategy {
    /// Contiguous layout
    Contiguous,
    /// Interleaved layout
    Interleaved,
    /// NUMA-aware layout
    NumaAware,
    /// Cache-optimized layout
    CacheOptimized,
}

impl AdvancedParallelStatsProcessor {
    /// Create new advanced-parallel processor
    pub fn new(config: AdvancedParallelConfig) -> StatsResult<Self> {
        let num_threads = _config
            .thread_pool_config
            .num_workers
            .unwrap_or_else(|| num_threads().max(1));

        let mut execution_contexts = Vec::with_capacity(num_threads);

        for i in 0..num_threads {
            let context = ParallelExecutionContext {
                thread_id: i,
                numa_node: i % 2, // Simplified NUMA assignment
                memory_pool: Some(Arc::new(Mutex::new(MemoryPool::new(100, 4096)))),
                counters: PerformanceCounters::default(),
                cache: ThreadLocalCache {
                    results: HashMap::new(),
                    stats: CacheStatistics::default(),
                },
            };
            execution_contexts.push(Arc::new(RwLock::new(context)));
        }

        Ok(Self {
            config,
            execution_contexts,
            work_queue: Arc::new(Mutex::new(Vec::new())),
            performance_history: Arc::new(Mutex::new(Vec::new())),
            optimization_cache: Arc::new(RwLock::new(HashMap::new())),
        })
    }

    /// Create with default configuration
    pub fn default() -> StatsResult<Self> {
        Self::new(AdvancedParallelConfig::default())
    }

    /// Compute mean using advanced-parallel processing
    pub fn mean_advanced_parallel<F>(&self, data: ArrayView1<F>) -> StatsResult<AdvancedParallelResult<F>>
    where
        F: Float + NumCast + Send + Sync + Zero + std::iter::Sum + std::fmt::Display,
    {
        let start_time = Instant::now();

        // Analyze workload and select optimal strategy
        let strategy = self.select_optimization_strategy("mean", data.len())?;

        // Distribute work among threads
        let work_units = self.create_work_units(&data, &strategy)?;

        // Execute parallel computation
        let partial_results = self.execute_parallel_work(&work_units)?;

        // Combine results
        let result = self.combine_mean_results(&partial_results, data.len())?;

        // Calculate metrics and analysis
        let total_time = start_time.elapsed();
        let metrics = self.calculate_execution_metrics(total_time, &work_units)?;
        let analysis = self.analyze_performance(&metrics)?;
        let resource_utilization = self.measure_resource_utilization()?;

        // Update performance history
        self.update_performance_history(&metrics);

        Ok(AdvancedParallelResult {
            result,
            metrics,
            analysis,
            resource_utilization,
        })
    }

    /// Compute variance using advanced-parallel processing
    pub fn variance_advanced_parallel<F>(
        &self,
        data: ArrayView1<F>,
        ddof: usize,
    ) -> StatsResult<AdvancedParallelResult<F>>
    where
        F: Float + NumCast + Send + Sync + Zero + std::iter::Sum + std::fmt::Display,
    {
        let start_time = Instant::now();

        // First compute mean in parallel
        let mean_result = self.mean_advanced_parallel(data)?;
        let mean_val = mean_result.result;

        // Select strategy for variance computation
        let strategy = self.select_optimization_strategy("variance", data.len())?;

        // Create work units for variance computation
        let work_units = self.create_variance_work_units(&data, mean_val, ddof, &strategy)?;

        // Execute parallel computation
        let partial_results = self.execute_parallel_work(&work_units)?;

        // Combine variance results
        let result = self.combine_variance_results(&partial_results, data.len(), ddof)?;

        let total_time = start_time.elapsed();
        let metrics = self.calculate_execution_metrics(total_time, &work_units)?;
        let analysis = self.analyze_performance(&metrics)?;
        let resource_utilization = self.measure_resource_utilization()?;

        self.update_performance_history(&metrics);

        Ok(AdvancedParallelResult {
            result,
            metrics,
            analysis,
            resource_utilization,
        })
    }

    /// Compute correlation matrix using advanced-parallel processing
    pub fn correlation_matrix_advanced_parallel<F>(
        &self,
        data: ArrayView2<F>,
    ) -> StatsResult<AdvancedParallelResult<Array2<F>>>
    where
        F: Float + NumCast + Send + Sync + Zero + std::iter::Sum + Clone + std::fmt::Display,
    {
        let start_time = Instant::now();
        let (n_rows, n_cols) = data.dim();

        // Create work units for each correlation pair
        let mut correlation_work = Vec::new();
        let mut work_id = 0;

        for i in 0..n_cols {
            for j in i..n_cols {
                let col_i = data.column(i).to_owned();
                let col_j = data.column(j).to_owned();

                correlation_work.push(WorkUnit {
                    id: work_id,
                    data: (
                        col_i.into_raw_vec_and_offset().0,
                        col_j.into_raw_vec_and_offset().0,
                        i,
                        j,
                    ),
                    cost: (n_rows as f64).sqrt(), // Correlation cost is roughly O(sqrt(n))
                    dependencies: Vec::new(),
                    priority: WorkPriority::Normal,
                    numa_node: Some(work_id % 2),
                });
                work_id += 1;
            }
        }

        // Execute parallel correlation computations
        let correlation_results = self.execute_correlation_work(correlation_work.as_slice())?;

        // Assemble correlation matrix
        let mut result_matrix = Array2::zeros((n_cols, n_cols));
        for ((i, j), correlation) in correlation_results {
            result_matrix[[i, j]] = correlation;
            if i != j {
                result_matrix[[j, i]] = correlation; // Symmetric matrix
            }
        }

        let total_time = start_time.elapsed();
        let metrics = self.calculate_matrix_execution_metrics(total_time, &correlation_work)?;
        let analysis = self.analyze_performance(&metrics)?;
        let resource_utilization = self.measure_resource_utilization()?;

        self.update_performance_history(&metrics);

        Ok(AdvancedParallelResult {
            result: result_matrix,
            metrics,
            analysis,
            resource_utilization,
        })
    }

    /// Select optimal strategy for the given operation
    fn select_optimization_strategy(
        &self,
        operation: &str,
        datasize: usize,
    ) -> StatsResult<OptimizationStrategy> {
        let cache_key = format!("{}_{}", operation, datasize / 1000); // Granular caching

        // Check cache first
        if let Ok(cache) = self.optimization_cache.read() {
            if let Some(strategy) = cache.get(&cache_key) {
                return Ok(strategy.clone());
            }
        }

        // Generate new strategy based on data characteristics
        let optimal_threads = self.calculate_optimal_thread_count(datasize);
        let work_distribution = if datasize > 1_000_000 {
            WorkDistributionMethod::CostBased
        } else if datasize > 100_000 {
            WorkDistributionMethod::SizeBased
        } else {
            WorkDistributionMethod::EqualChunks
        };

        let memory_layout = if self.config.numa_aware {
            MemoryLayoutStrategy::NumaAware
        } else if self.config.cache_optimization != CacheOptimizationLevel::None {
            MemoryLayoutStrategy::CacheOptimized
        } else {
            MemoryLayoutStrategy::Contiguous
        };

        let strategy = OptimizationStrategy {
            name: format!("{}_optimized", operation),
            thread_count: optimal_threads,
            work_distribution,
            memory_layout,
            expected_performance: self.estimate_performance(optimal_threads, datasize),
        };

        // Cache the strategy
        if let Ok(mut cache) = self.optimization_cache.write() {
            cache.insert(cache_key, strategy.clone());
        }

        Ok(strategy)
    }

    /// Calculate optimal thread count for given data size
    fn calculate_optimal_thread_count(&self, datasize: usize) -> usize {
        let available_threads = self.execution_contexts.len();
        let min_work_per_thread = self.config.min_worksize;

        // Don't use more threads than can be effectively utilized
        let max_useful_threads = (datasize / min_work_per_thread).max(1);

        // Consider NUMA topology
        let numa_optimal = if self.config.numa_aware {
            // Simplified: assume 2 NUMA nodes
            (available_threads / 2) * 2
        } else {
            available_threads
        };

        max_useful_threads.min(numa_optimal).min(available_threads)
    }

    /// Estimate performance for given configuration
    fn estimate_performance(&self, thread_count: usize, datasize: usize) -> f64 {
        // Simplified performance model
        let sequential_time = datasize as f64;
        let parallel_efficiency = 0.8; // Account for overhead
        let parallel_time = sequential_time / (thread_count as f64 * parallel_efficiency);

        sequential_time / parallel_time
    }

    /// Create work units for mean computation
    fn create_work_units<F>(
        &self,
        data: &ArrayView1<F>,
        strategy: &OptimizationStrategy,
    ) -> StatsResult<Vec<WorkUnit<Vec<f64>>>>
    where
        F: Float + NumCast + std::fmt::Display,
    {
        let mut work_units = Vec::new();
        let datasize = data.len();
        let chunksize = datasize / strategy.thread_count;

        for i in 0..strategy.thread_count {
            let start = i * chunksize;
            let end = if i == strategy.thread_count - 1 {
                datasize
            } else {
                (i + 1) * chunksize
            };

            let chunkdata: Vec<f64> = data
                .slice(s![start..end])
                .iter()
                .map(|&x| x.to_f64().unwrap_or(0.0))
                .collect();

            work_units.push(WorkUnit {
                id: i,
                data: chunkdata,
                cost: (end - start) as f64,
                dependencies: Vec::new(),
                priority: WorkPriority::Normal,
                numa_node: if self.config.numa_aware {
                    Some(i % 2) // Simplified NUMA assignment
                } else {
                    None
                },
            });
        }

        Ok(work_units)
    }

    /// Create work units for variance computation
    fn create_variance_work_units<F>(
        &self,
        data: &ArrayView1<F>,
        mean_val: F, _ddof: usize,
        strategy: &OptimizationStrategy,
    ) -> StatsResult<Vec<WorkUnit<Vec<f64>>>>
    where
        F: Float + NumCast + std::fmt::Display,
    {
        let mut work_units = Vec::new();
        let datasize = data.len();
        let chunksize = datasize / strategy.thread_count;
        let mean_f64 = mean_val.to_f64().unwrap_or(0.0);

        for i in 0..strategy.thread_count {
            let start = i * chunksize;
            let end = if i == strategy.thread_count - 1 {
                datasize
            } else {
                (i + 1) * chunksize
            };

            let chunkdata: Vec<f64> = data
                .slice(s![start..end])
                .iter()
                .map(|&x| {
                    let _val = x.to_f64().unwrap_or(0.0);
                    let diff = _val - mean_f64;
                    diff * diff
                })
                .collect();

            work_units.push(WorkUnit {
                id: i,
                data: chunkdata,
                cost: (end - start) as f64,
                dependencies: Vec::new(),
                priority: WorkPriority::Normal,
                numa_node: if self.config.numa_aware {
                    Some(i % 2)
                } else {
                    None
                },
            });
        }

        Ok(work_units)
    }

    /// Execute parallel work units
    fn execute_parallel_work(&self, workunits: &[WorkUnit<Vec<f64>>]) -> StatsResult<Vec<f64>> {
        let num_threads = work_units.len();
        let results = Arc::new(Mutex::new(vec![0.0; num_threads]));
        let work_units = Arc::new(work_units.to_vec());

        thread::scope(|s| {
            let handles: Vec<_> = (0..num_threads)
                .map(|thread_id| {
                    let results = Arc::clone(&results);
                    let work_units = Arc::clone(&work_units);

                    s.spawn(move || {
                        let work_unit = &work_units[thread_id];
                        let sum: f64 = work_unit.data.iter().sum();

                        if let Ok(mut results) = results.lock() {
                            results[thread_id] = sum;
                        }
                    })
                })
                .collect();

            // Wait for all threads to complete
            for handle in handles {
                let _ = handle.join();
            }
        });

        let results = results.lock().expect("Operation failed").clone();
        Ok(results)
    }

    /// Execute correlation work units
    fn execute_correlation_work<F>(
        &self,
        work_units: &[WorkUnit<(Vec<F>, Vec<F>, usize, usize)>],
    ) -> StatsResult<Vec<((usize, usize), F)>>
    where
        F: Float + NumCast + Send + Sync + Clone + std::iter::Sum + std::fmt::Display,
    {
        let num_work_units = work_units.len();
        let results = Arc::new(Mutex::new(Vec::with_capacity(num_work_units)));
        let work_units = Arc::new(work_units.to_vec());

        thread::scope(|s| {
            let handles: Vec<_> = (0..num_work_units)
                .map(|work_id| {
                    let results = Arc::clone(&results);
                    let work_units = Arc::clone(&work_units);

                    s.spawn(move || {
                        let work_unit = &work_units[work_id];
                        let (ref x, ref y, i, j) = work_unit.data;

                        // Compute Pearson correlation
                        let correlation = self.compute_correlation(x, y).unwrap_or(F::zero());

                        if let Ok(mut results) = results.lock() {
                            results.push(((i, j), correlation));
                        }
                    })
                })
                .collect();

            for handle in handles {
                let _ = handle.join();
            }
        });

        let results = results.lock().expect("Operation failed").clone();
        Ok(results)
    }

    /// Compute Pearson correlation between two vectors
    fn compute_correlation<F>(&self, x: &[F], y: &[F]) -> StatsResult<F>
    where
        F: Float + NumCast + Send + Sync + Clone + std::iter::Sum + std::fmt::Display,
    {
        if x.len() != y.len() || x.is_empty() {
            return Ok(F::zero());
        }

        let n = F::from(x.len()).expect("Operation failed");
        let sum_x: F = x.iter().cloned().sum();
        let sum_y: F = y.iter().cloned().sum();
        let sum_xx: F = x.iter().map(|&xi| xi * xi).sum();
        let sum_yy: F = y.iter().map(|&yi| yi * yi).sum();
        let sum_xy: F = x.iter().zip(y).map(|(&xi, &yi)| xi * yi).sum();

        let numerator = n * sum_xy - sum_x * sum_y;
        let denominator = ((n * sum_xx - sum_x * sum_x) * (n * sum_yy - sum_y * sum_y)).sqrt();

        if denominator == F::zero() {
            Ok(F::zero())
        } else {
            Ok(numerator / denominator)
        }
    }

    /// Combine mean results from parallel computation
    fn combine_mean_results<F>(&self, partial_results: &[f64], totalcount: usize) -> StatsResult<F>
    where
        F: Float + NumCast + std::fmt::Display,
    {
        let total_sum: f64 = partial_results.iter().sum();
        let mean = total_sum / total_count as f64;
        F::from(mean).ok_or_else(|| {
            StatsError::ComputationError("Failed to convert mean result".to_string())
        })
    }

    /// Combine variance results from parallel computation
    fn combine_variance_results<F>(
        &self,
        partial_results: &[f64],
        total_count: usize,
        ddof: usize,
    ) -> StatsResult<F>
    where
        F: Float + NumCast + std::fmt::Display,
    {
        let total_sum_sq_dev: f64 = partial_results.iter().sum();
        let variance = total_sum_sq_dev / (total_count - ddof) as f64;
        F::from(variance).ok_or_else(|| {
            StatsError::ComputationError("Failed to convert variance result".to_string())
        })
    }

    /// Calculate execution metrics
    fn calculate_execution_metrics(
        &self,
        total_time: Duration,
        work_units: &[WorkUnit<Vec<f64>>],
    ) -> StatsResult<ParallelExecutionMetrics> {
        let threads_used = work_units.len();
        let total_work: f64 = work_units.iter().map(|wu| wu.cost).sum();
        let avg_work_per_thread = total_work / threads_used as f64;

        // Estimate load balance efficiency
        let work_variance = work_units
            .iter()
            .map(|wu| (wu.cost - avg_work_per_thread).powi(2))
            .sum::<f64>()
            / threads_used as f64;
        let load_balance_efficiency = 1.0 - (work_variance.sqrt() / avg_work_per_thread).min(1.0);

        // Estimate parallel efficiency (simplified)
        let sequential_time_estimate = total_time.mul_f64(threads_used as f64);
        let parallel_efficiency = total_time.as_secs_f64() / sequential_time_estimate.as_secs_f64();

        let speedup = threads_used as f64 * parallel_efficiency;

        Ok(ParallelExecutionMetrics {
            total_time,
            parallel_time: total_time.mul_f64(0.9), // Estimate
            sequential_time: total_time.mul_f64(0.1), // Estimate
            sync_time: total_time.mul_f64(0.05),    // Estimate
            threads_used,
            load_balance_efficiency,
            parallel_efficiency,
            speedup,
            work_distribution_quality: load_balance_efficiency,
        })
    }

    /// Calculate matrix operation execution metrics
    fn calculate_matrix_execution_metrics<F>(
        &self,
        total_time: Duration,
        work_units: &[WorkUnit<(Vec<F>, Vec<F>, usize, usize)>],
    ) -> StatsResult<ParallelExecutionMetrics>
    where
        F: Float + NumCast + Send + Sync + Clone + std::iter::Sum + std::fmt::Display,
    {
        let threads_used = work_units.len();
        let _total_work: f64 = work_units.iter().map(|wu| wu.cost).sum();
        let load_balance_efficiency = 0.85; // Estimate for matrix operations

        Ok(ParallelExecutionMetrics {
            total_time,
            parallel_time: total_time.mul_f64(0.85),
            sequential_time: total_time.mul_f64(0.15),
            sync_time: total_time.mul_f64(0.08),
            threads_used,
            load_balance_efficiency,
            parallel_efficiency: 0.8, // Estimate
            speedup: threads_used as f64 * 0.8,
            work_distribution_quality: load_balance_efficiency,
        })
    }

    /// Analyze performance characteristics
    fn analyze_performance(
        &self,
        metrics: &ParallelExecutionMetrics,
    ) -> StatsResult<ParallelPerformanceAnalysis> {
        let mut bottlenecks = Vec::new();

        // Detect load imbalance
        if metrics.load_balance_efficiency < 0.8 {
            bottlenecks.push(PerformanceBottleneck {
                bottleneck_type: BottleneckType::LoadImbalance,
                severity: 1.0 - metrics.load_balance_efficiency,
                description: "Load imbalance detected among threads".to_string(),
                mitigation: "Consider dynamic work distribution".to_string(),
            });
        }

        // Detect synchronization overhead
        if metrics.sync_time.as_secs_f64() / metrics.total_time.as_secs_f64() > 0.1 {
            bottlenecks.push(PerformanceBottleneck {
                bottleneck_type: BottleneckType::SynchronizationOverhead,
                severity: metrics.sync_time.as_secs_f64() / metrics.total_time.as_secs_f64(),
                description: "High synchronization overhead".to_string(),
                mitigation: "Reduce synchronization points or use lock-free algorithms".to_string(),
            });
        }

        // Performance rating
        let performance_rating = if metrics.parallel_efficiency > 0.9 {
            PerformanceRating::Excellent
        } else if metrics.parallel_efficiency > 0.7 {
            PerformanceRating::Good
        } else if metrics.parallel_efficiency > 0.5 {
            PerformanceRating::Acceptable
        } else if metrics.parallel_efficiency > 0.3 {
            PerformanceRating::Poor
        } else {
            PerformanceRating::Unacceptable
        };

        Ok(ParallelPerformanceAnalysis {
            bottlenecks,
            scaling_analysis: ScalingAnalysis {
                theoretical_max_speedup: metrics.threads_used as f64,
                achieved_speedup: metrics.speedup,
                parallel_fraction: 0.9,             // Estimate
                serial_bottleneck_impact: 0.1,      // Estimate
                scaling_efficiency: HashMap::new(), // Would implement proper analysis
                optimal_thread_count: metrics.threads_used,
            },
            optimization_opportunities: Vec::new(), // Would implement opportunity detection
            performance_rating,
        })
    }

    /// Measure resource utilization
    fn measure_resource_utilization(&self) -> StatsResult<ResourceUtilization> {
        // Simplified resource measurement
        Ok(ResourceUtilization {
            cpu_utilization: vec![0.8; self.execution_contexts.len()], // 80% utilization estimate
            memory_utilization: 0.6, // 60% memory utilization estimate
            cache_utilization: CacheUtilization {
                l1_hit_rate: 0.95,
                l2_hit_rate: 0.85,
                l3_hit_rate: 0.75,
                cache_line_utilization: 0.8,
            },
            numa_utilization: vec![0.8, 0.8], // Balanced NUMA utilization
            energy_consumption: None,         // Would require hardware monitoring
        })
    }

    /// Update performance history for learning
    fn update_performance_history(&self, metrics: &ParallelExecutionMetrics) {
        if let Ok(mut history) = self.performance_history.lock() {
            history.push(metrics.clone());

            // Keep only recent history
            if history.len() > 1000 {
                history.remove(0);
            }
        }
    }

    /// Get performance statistics
    pub fn get_performance_statistics(&self) -> PerformanceStatistics {
        if let Ok(history) = self.performance_history.lock() {
            let total_operations = history.len();
            let avg_speedup = if !history.is_empty() {
                history.iter().map(|m| m.speedup).sum::<f64>() / history.len() as f64
            } else {
                0.0
            };

            let avg_efficiency = if !history.is_empty() {
                history.iter().map(|m| m.parallel_efficiency).sum::<f64>() / history.len() as f64
            } else {
                0.0
            };

            PerformanceStatistics {
                total_operations,
                average_speedup: avg_speedup,
                best_strategies: Vec::new(), // Would implement strategy tracking
                hardware_utilization: HardwareUtilization {
                    simd_utilization: 0.7, // Estimate
                    memory_bandwidth_utilization: 0.6,
                    cache_efficiency: avg_efficiency,
                    energy_efficiency: None,
                },
            }
        } else {
            PerformanceStatistics {
                total_operations: 0,
                average_speedup: 0.0,
                best_strategies: Vec::new(),
                hardware_utilization: HardwareUtilization {
                    simd_utilization: 0.0,
                    memory_bandwidth_utilization: 0.0,
                    cache_efficiency: 0.0,
                    energy_efficiency: None,
                },
            }
        }
    }
}

/// Convenience functions for advanced-parallel operations
#[allow(dead_code)]
pub fn create_advanced_parallel_processor() -> StatsResult<AdvancedParallelStatsProcessor> {
    AdvancedParallelStatsProcessor::default()
}

#[allow(dead_code)]
pub fn mean_advanced_parallel<F>(data: ArrayView1<F>) -> StatsResult<AdvancedParallelResult<F>>
where
    F: Float + NumCast + Send + Sync + Zero + std::iter::Sum + std::fmt::Display,
{
    let processor = AdvancedParallelStatsProcessor::default()?;
    processor.mean_advanced_parallel(data)
}

#[allow(dead_code)]
pub fn variance_advanced_parallel<F>(
    data: ArrayView1<F>,
    ddof: usize,
) -> StatsResult<AdvancedParallelResult<F>>
where
    F: Float + NumCast + Send + Sync + Zero + std::iter::Sum + std::fmt::Display,
{
    let processor = AdvancedParallelStatsProcessor::default()?;
    processor.variance_advanced_parallel(data, ddof)
}

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

    #[test]
    fn test_advanced_parallel_config() {
        let config = AdvancedParallelConfig::default();
        assert!(config.adaptive_work_distribution);
        assert!(config.numa_aware);
        assert!(config.performance_monitoring);
    }

    #[test]
    fn test_processor_creation() {
        let processor = AdvancedParallelStatsProcessor::default().expect("Operation failed");
        assert!(!processor.execution_contexts.is_empty());
    }

    #[test]
    fn test_optimization_strategy_selection() {
        let processor = AdvancedParallelStatsProcessor::default().expect("Operation failed");
        let strategy = processor
            .select_optimization_strategy("mean", 10000)
            .expect("Operation failed");

        assert!(!strategy.name.is_empty());
        assert!(strategy.thread_count > 0);
        assert!(strategy.expected_performance > 0.0);
    }

    #[test]
    fn test_work_unit_creation() {
        let processor = AdvancedParallelStatsProcessor::default().expect("Operation failed");
        let data = array![1.0f64, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
        let strategy = OptimizationStrategy {
            name: "test".to_string(),
            thread_count: 2,
            work_distribution: WorkDistributionMethod::EqualChunks,
            memory_layout: MemoryLayoutStrategy::Contiguous,
            expected_performance: 2.0,
        };

        let work_units = processor
            .create_work_units(&data.view(), &strategy)
            .expect("Operation failed");
        assert_eq!(work_units.len(), 2);
        assert!(!work_units[0].data.is_empty());
        assert!(!work_units[1].data.is_empty());
    }

    #[test]
    fn test_correlation_computation() {
        let processor = AdvancedParallelStatsProcessor::default().expect("Operation failed");
        let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let y = vec![2.0, 4.0, 6.0, 8.0, 10.0];

        let correlation = processor.compute_correlation(&x, &y).expect("Operation failed");
        assert!((correlation - 1.0).abs() < 1e-10); // Perfect positive correlation
    }

    #[test]
    fn test_performance_metrics_calculation() {
        let processor = AdvancedParallelStatsProcessor::default().expect("Operation failed");
        let work_units = vec![
            WorkUnit {
                id: 0,
                data: vec![1.0, 2.0],
                cost: 100.0,
                dependencies: Vec::new(),
                priority: WorkPriority::Normal,
                numa_node: None,
            },
            WorkUnit {
                id: 1,
                data: vec![3.0, 4.0],
                cost: 120.0,
                dependencies: Vec::new(),
                priority: WorkPriority::Normal,
                numa_node: None,
            },
        ];

        let metrics = processor
            .calculate_execution_metrics(Duration::from_millis(100), &work_units)
            .expect("Operation failed");

        assert_eq!(metrics.threads_used, 2);
        assert!(metrics.load_balance_efficiency > 0.0);
        assert!(metrics.speedup > 0.0);
    }
}