scirs2_cluster/gpu/

acceleration.rs

//! Advanced GPU Acceleration Module for v0.2.0
//!
//! This module provides comprehensive GPU acceleration for clustering algorithms with:
//! - Multiple backend support (CUDA, OpenCL, ROCm, Metal, OneAPI)
//! - Advanced memory management strategies
//! - Tensor core and mixed precision support
//! - Automatic CPU fallback
//! - GPU-accelerated K-means clustering

use crate::error::{ClusteringError, Result};
use scirs2_core::ndarray::{Array1, Array2, ArrayView2, Axis};
use scirs2_core::numeric::{Float, FromPrimitive};
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::fmt;
use std::time::{Duration, Instant};

use super::core::{DeviceSelection, GpuBackend, GpuConfig, GpuContext, GpuDevice};
use super::memory::{GpuMemoryBlock, GpuMemoryManager, MemoryStats, MemoryStrategy};

// ============================================================================
// Advanced Memory Management
// ============================================================================

/// Advanced memory management strategy for GPU operations
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
pub enum AdvancedMemoryStrategy {
    /// Conservative: Minimize GPU memory usage, more host-device transfers
    Conservative,
    /// Aggressive: Maximize GPU memory usage for speed
    Aggressive,
    /// Adaptive: Dynamically adjust based on available memory and workload
    Adaptive,
    /// Streaming: Process data in chunks for datasets larger than GPU memory
    Streaming {
        /// Chunk size in bytes
        chunk_size: usize,
    },
    /// Unified: Use unified memory where available (CUDA managed memory)
    Unified,
    /// Pool: Use memory pool for fast allocations/deallocations
    Pool {
        /// Pool size in bytes
        pool_size: usize,
    },
}

impl Default for AdvancedMemoryStrategy {
    fn default() -> Self {
        AdvancedMemoryStrategy::Adaptive
    }
}

impl fmt::Display for AdvancedMemoryStrategy {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            AdvancedMemoryStrategy::Conservative => write!(f, "Conservative"),
            AdvancedMemoryStrategy::Aggressive => write!(f, "Aggressive"),
            AdvancedMemoryStrategy::Adaptive => write!(f, "Adaptive"),
            AdvancedMemoryStrategy::Streaming { chunk_size } => {
                write!(f, "Streaming({}MB)", chunk_size / (1024 * 1024))
            }
            AdvancedMemoryStrategy::Unified => write!(f, "Unified"),
            AdvancedMemoryStrategy::Pool { pool_size } => {
                write!(f, "Pool({}MB)", pool_size / (1024 * 1024))
            }
        }
    }
}

/// Advanced GPU memory manager with multiple strategies
#[derive(Debug)]
pub struct AdvancedGpuMemoryManager {
    /// Base memory manager
    base_manager: GpuMemoryManager,
    /// Current memory strategy
    strategy: AdvancedMemoryStrategy,
    /// Available GPU memory in bytes
    available_memory: usize,
    /// Memory allocation history for adaptive strategy
    allocation_history: Vec<AllocationRecord>,
    /// Memory pressure threshold (0.0 to 1.0)
    pressure_threshold: f64,
    /// Enable memory defragmentation
    enable_defrag: bool,
    /// Memory usage statistics
    usage_stats: MemoryUsageStats,
}

/// Record of memory allocation for adaptive management
#[derive(Debug, Clone)]
pub struct AllocationRecord {
    /// Size of allocation
    pub size: usize,
    /// Timestamp
    pub timestamp: Instant,
    /// Duration of use
    pub duration: Option<Duration>,
    /// Was allocation successful
    pub success: bool,
}

/// Memory usage statistics
#[derive(Debug, Clone, Default, Serialize, Deserialize)]
pub struct MemoryUsageStats {
    /// Total allocations
    pub total_allocations: usize,
    /// Successful allocations
    pub successful_allocations: usize,
    /// Failed allocations (out of memory)
    pub failed_allocations: usize,
    /// Total bytes allocated over time
    pub total_bytes_allocated: usize,
    /// Current bytes in use
    pub current_bytes_in_use: usize,
    /// Peak bytes in use
    pub peak_bytes_in_use: usize,
    /// Average allocation size
    pub avg_allocation_size: f64,
    /// Memory efficiency (successful/total)
    pub efficiency: f64,
}

impl AdvancedGpuMemoryManager {
    /// Create a new advanced memory manager
    pub fn new(strategy: AdvancedMemoryStrategy, available_memory: usize) -> Self {
        let alignment = 256; // 256-byte alignment for GPU
        let max_pool_size = match strategy {
            AdvancedMemoryStrategy::Pool { pool_size } => pool_size / (1024 * 1024),
            _ => 100, // Default pool size
        };

        Self {
            base_manager: GpuMemoryManager::new(alignment, max_pool_size),
            strategy,
            available_memory,
            allocation_history: Vec::new(),
            pressure_threshold: 0.85,
            enable_defrag: true,
            usage_stats: MemoryUsageStats::default(),
        }
    }

    /// Allocate memory with strategy-aware logic
    pub fn allocate(&mut self, size: usize) -> Result<GpuMemoryBlock> {
        self.usage_stats.total_allocations += 1;

        // Check memory pressure
        let memory_pressure = self.calculate_memory_pressure();

        // Handle based on strategy
        let result = match self.strategy {
            AdvancedMemoryStrategy::Conservative => self.allocate_conservative(size),
            AdvancedMemoryStrategy::Aggressive => self.allocate_aggressive(size),
            AdvancedMemoryStrategy::Adaptive => self.allocate_adaptive(size, memory_pressure),
            AdvancedMemoryStrategy::Streaming { chunk_size } => {
                self.allocate_streaming(size, chunk_size)
            }
            AdvancedMemoryStrategy::Unified => self.allocate_unified(size),
            AdvancedMemoryStrategy::Pool { .. } => self.base_manager.allocate(size),
        };

        // Record allocation
        let success = result.is_ok();
        self.allocation_history.push(AllocationRecord {
            size,
            timestamp: Instant::now(),
            duration: None,
            success,
        });

        if success {
            self.usage_stats.successful_allocations += 1;
            self.usage_stats.total_bytes_allocated += size;
            self.usage_stats.current_bytes_in_use += size;
            self.usage_stats.peak_bytes_in_use = self
                .usage_stats
                .peak_bytes_in_use
                .max(self.usage_stats.current_bytes_in_use);
        } else {
            self.usage_stats.failed_allocations += 1;
        }

        self.update_efficiency();
        result
    }

    /// Deallocate memory
    pub fn deallocate(&mut self, block: GpuMemoryBlock) -> Result<()> {
        let size = block.size;
        self.base_manager.deallocate(block)?;
        self.usage_stats.current_bytes_in_use =
            self.usage_stats.current_bytes_in_use.saturating_sub(size);
        Ok(())
    }

    /// Conservative allocation strategy
    fn allocate_conservative(&mut self, size: usize) -> Result<GpuMemoryBlock> {
        // Check if we have enough memory before allocating
        if self.usage_stats.current_bytes_in_use + size > self.available_memory {
            // Try to free unused memory first
            self.compact_memory()?;
        }
        self.base_manager.allocate(size)
    }

    /// Aggressive allocation strategy
    fn allocate_aggressive(&mut self, size: usize) -> Result<GpuMemoryBlock> {
        // Allocate without checking, rely on GPU driver
        self.base_manager.allocate(size)
    }

    /// Adaptive allocation strategy
    fn allocate_adaptive(&mut self, size: usize, memory_pressure: f64) -> Result<GpuMemoryBlock> {
        if memory_pressure > self.pressure_threshold {
            // High pressure: use conservative approach
            self.allocate_conservative(size)
        } else {
            // Low pressure: use aggressive approach
            self.allocate_aggressive(size)
        }
    }

    /// Streaming allocation for large datasets
    fn allocate_streaming(&mut self, size: usize, chunk_size: usize) -> Result<GpuMemoryBlock> {
        let actual_size = size.min(chunk_size);
        self.base_manager.allocate(actual_size)
    }

    /// Unified memory allocation
    fn allocate_unified(&mut self, size: usize) -> Result<GpuMemoryBlock> {
        // In real implementation, would use CUDA managed memory
        self.base_manager.allocate(size)
    }

    /// Calculate current memory pressure (0.0 to 1.0)
    fn calculate_memory_pressure(&self) -> f64 {
        if self.available_memory == 0 {
            return 1.0;
        }
        self.usage_stats.current_bytes_in_use as f64 / self.available_memory as f64
    }

    /// Compact memory by freeing unused allocations
    fn compact_memory(&mut self) -> Result<()> {
        self.base_manager.clear_pools()
    }

    /// Update efficiency statistics
    fn update_efficiency(&mut self) {
        if self.usage_stats.total_allocations > 0 {
            self.usage_stats.efficiency = self.usage_stats.successful_allocations as f64
                / self.usage_stats.total_allocations as f64;
            self.usage_stats.avg_allocation_size = self.usage_stats.total_bytes_allocated as f64
                / self.usage_stats.successful_allocations.max(1) as f64;
        }
    }

    /// Get memory statistics
    pub fn get_stats(&self) -> &MemoryUsageStats {
        &self.usage_stats
    }

    /// Get current memory strategy
    pub fn strategy(&self) -> AdvancedMemoryStrategy {
        self.strategy
    }

    /// Set memory strategy
    pub fn set_strategy(&mut self, strategy: AdvancedMemoryStrategy) {
        self.strategy = strategy;
    }

    /// Get memory pressure threshold
    pub fn pressure_threshold(&self) -> f64 {
        self.pressure_threshold
    }

    /// Set memory pressure threshold
    pub fn set_pressure_threshold(&mut self, threshold: f64) {
        self.pressure_threshold = threshold.clamp(0.0, 1.0);
    }

    /// Check if defragmentation is enabled
    pub fn is_defrag_enabled(&self) -> bool {
        self.enable_defrag
    }

    /// Enable or disable defragmentation
    pub fn set_defrag_enabled(&mut self, enabled: bool) {
        self.enable_defrag = enabled;
    }
}

// ============================================================================
// Tensor Core and Mixed Precision Support
// ============================================================================

/// Precision mode for GPU computations
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
pub enum PrecisionMode {
    /// Full precision (f64)
    Full,
    /// Single precision (f32)
    Single,
    /// Half precision (f16)
    Half,
    /// Mixed precision (f16 compute, f32 accumulator)
    Mixed,
    /// Brain floating point (bf16)
    BFloat16,
    /// Tensor float 32 (TF32) for NVIDIA Ampere+
    TensorFloat32,
    /// Automatic selection based on hardware
    Auto,
}

impl Default for PrecisionMode {
    fn default() -> Self {
        PrecisionMode::Auto
    }
}

impl fmt::Display for PrecisionMode {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            PrecisionMode::Full => write!(f, "Full (f64)"),
            PrecisionMode::Single => write!(f, "Single (f32)"),
            PrecisionMode::Half => write!(f, "Half (f16)"),
            PrecisionMode::Mixed => write!(f, "Mixed (f16/f32)"),
            PrecisionMode::BFloat16 => write!(f, "BFloat16"),
            PrecisionMode::TensorFloat32 => write!(f, "TF32"),
            PrecisionMode::Auto => write!(f, "Auto"),
        }
    }
}

/// Tensor core configuration
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct TensorCoreConfig {
    /// Enable tensor cores if available
    pub enabled: bool,
    /// Precision mode
    pub precision: PrecisionMode,
    /// Tile size for tensor core operations (M, N, K)
    pub tile_size: (usize, usize, usize),
    /// Use structured sparsity if available (NVIDIA Ampere+)
    pub use_sparsity: bool,
    /// Sparsity ratio (e.g., 0.5 for 2:4 sparsity)
    pub sparsity_ratio: f64,
    /// Enable automatic precision scaling
    pub auto_scale: bool,
    /// Loss scaling factor for mixed precision training
    pub loss_scale: f64,
}

impl Default for TensorCoreConfig {
    fn default() -> Self {
        Self {
            enabled: true,
            precision: PrecisionMode::Auto,
            tile_size: (16, 16, 16),
            use_sparsity: false,
            sparsity_ratio: 0.5,
            auto_scale: true,
            loss_scale: 1.0,
        }
    }
}

/// Tensor core capabilities detection
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct TensorCoreCapabilities {
    /// Tensor cores available
    pub available: bool,
    /// Supported precision modes
    pub supported_precisions: Vec<PrecisionMode>,
    /// Supported tile sizes
    pub supported_tile_sizes: Vec<(usize, usize, usize)>,
    /// Supports structured sparsity
    pub supports_sparsity: bool,
    /// Peak TOPS (Tera Operations Per Second)
    pub peak_tops: Option<f64>,
    /// Architecture name
    pub architecture: String,
}

impl Default for TensorCoreCapabilities {
    fn default() -> Self {
        Self {
            available: false,
            supported_precisions: vec![PrecisionMode::Single],
            supported_tile_sizes: vec![(16, 16, 16)],
            supports_sparsity: false,
            peak_tops: None,
            architecture: "Unknown".to_string(),
        }
    }
}

/// Detect tensor core capabilities for a GPU device
pub fn detect_tensor_core_capabilities(device: &GpuDevice) -> TensorCoreCapabilities {
    match device.backend {
        GpuBackend::Cuda => {
            // NVIDIA Tensor Core detection
            TensorCoreCapabilities {
                available: true,
                supported_precisions: vec![
                    PrecisionMode::Half,
                    PrecisionMode::Mixed,
                    PrecisionMode::BFloat16,
                    PrecisionMode::TensorFloat32,
                ],
                supported_tile_sizes: vec![(16, 16, 16), (32, 8, 16), (8, 32, 16)],
                supports_sparsity: true, // Ampere+ supports 2:4 sparsity
                peak_tops: Some(312.0),  // Example for A100
                architecture: "NVIDIA Tensor Cores".to_string(),
            }
        }
        GpuBackend::Rocm => {
            // AMD Matrix Cores
            TensorCoreCapabilities {
                available: true,
                supported_precisions: vec![
                    PrecisionMode::Half,
                    PrecisionMode::Mixed,
                    PrecisionMode::BFloat16,
                ],
                supported_tile_sizes: vec![(32, 32, 8), (16, 16, 16)],
                supports_sparsity: false,
                peak_tops: Some(383.0), // Example for MI250X
                architecture: "AMD Matrix Cores".to_string(),
            }
        }
        GpuBackend::Metal => {
            // Apple Neural Engine / GPU
            TensorCoreCapabilities {
                available: true,
                supported_precisions: vec![PrecisionMode::Half, PrecisionMode::Single],
                supported_tile_sizes: vec![(16, 16, 16)],
                supports_sparsity: false,
                peak_tops: Some(15.8), // Example for M1
                architecture: "Apple Neural Engine".to_string(),
            }
        }
        _ => TensorCoreCapabilities::default(),
    }
}

// ============================================================================
// Device Selection Strategies
// ============================================================================

/// Advanced device selection strategy
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum AdvancedDeviceSelection {
    /// Use the first available device
    First,
    /// Use the device with most available memory
    MostMemory,
    /// Use the device with highest compute capability
    HighestCompute,
    /// Use a specific device by ID
    Specific(u32),
    /// Automatic selection based on workload
    Auto,
    /// Use fastest device for current workload (benchmarked)
    Fastest,
    /// Use device with best power efficiency
    MostEfficient,
    /// Round-robin across available devices
    RoundRobin,
    /// Load-balanced across devices based on utilization
    LoadBalanced,
    /// Multi-GPU: use all available GPUs
    MultiGpu {
        /// Maximum number of GPUs to use
        max_gpus: usize,
    },
}

impl Default for AdvancedDeviceSelection {
    fn default() -> Self {
        AdvancedDeviceSelection::Auto
    }
}

impl From<AdvancedDeviceSelection> for DeviceSelection {
    fn from(adv: AdvancedDeviceSelection) -> Self {
        match adv {
            AdvancedDeviceSelection::First => DeviceSelection::First,
            AdvancedDeviceSelection::MostMemory => DeviceSelection::MostMemory,
            AdvancedDeviceSelection::HighestCompute => DeviceSelection::HighestCompute,
            AdvancedDeviceSelection::Specific(id) => DeviceSelection::Specific(id),
            AdvancedDeviceSelection::Auto => DeviceSelection::Auto,
            AdvancedDeviceSelection::Fastest => DeviceSelection::Fastest,
            _ => DeviceSelection::Auto, // Default for advanced strategies
        }
    }
}

/// Device selector for multi-GPU operations
#[derive(Debug)]
pub struct DeviceSelector {
    /// Available devices
    devices: Vec<GpuDevice>,
    /// Selection strategy
    strategy: AdvancedDeviceSelection,
    /// Device utilization tracking
    utilization: HashMap<u32, f64>,
    /// Round-robin counter
    round_robin_idx: usize,
    /// Device benchmark results
    benchmarks: HashMap<u32, DeviceBenchmark>,
}

/// Device benchmark result
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct DeviceBenchmark {
    /// Device ID
    pub device_id: u32,
    /// Distance computation throughput (GFLOPS)
    pub distance_throughput: f64,
    /// K-means iteration time (ms)
    pub kmeans_time_ms: f64,
    /// Memory bandwidth (GB/s)
    pub memory_bandwidth: f64,
    /// Power consumption (W)
    pub power_consumption: Option<f64>,
    /// Benchmark timestamp
    pub timestamp: std::time::SystemTime,
}

impl DeviceSelector {
    /// Create a new device selector
    pub fn new(strategy: AdvancedDeviceSelection) -> Self {
        Self {
            devices: Vec::new(),
            strategy,
            utilization: HashMap::new(),
            round_robin_idx: 0,
            benchmarks: HashMap::new(),
        }
    }

    /// Add a device to the selector
    pub fn add_device(&mut self, device: GpuDevice) {
        self.utilization.insert(device.device_id, 0.0);
        self.devices.push(device);
    }

    /// Select the best device based on current strategy
    pub fn select_device(&mut self) -> Option<&GpuDevice> {
        if self.devices.is_empty() {
            return None;
        }

        match &self.strategy {
            AdvancedDeviceSelection::First => self.devices.first(),
            AdvancedDeviceSelection::MostMemory => {
                self.devices.iter().max_by_key(|d| d.available_memory)
            }
            AdvancedDeviceSelection::HighestCompute => {
                self.devices.iter().max_by_key(|d| d.compute_units)
            }
            AdvancedDeviceSelection::Specific(id) => {
                self.devices.iter().find(|d| d.device_id == *id)
            }
            AdvancedDeviceSelection::Auto => {
                // Score-based selection
                self.devices.iter().max_by(|a, b| {
                    a.get_device_score()
                        .partial_cmp(&b.get_device_score())
                        .unwrap_or(std::cmp::Ordering::Equal)
                })
            }
            AdvancedDeviceSelection::Fastest => {
                // Use benchmark results if available
                if self.benchmarks.is_empty() {
                    self.devices.first()
                } else {
                    let fastest_id = self
                        .benchmarks
                        .iter()
                        .min_by(|a, b| {
                            a.1.kmeans_time_ms
                                .partial_cmp(&b.1.kmeans_time_ms)
                                .unwrap_or(std::cmp::Ordering::Equal)
                        })
                        .map(|(id, _)| *id);

                    fastest_id.and_then(|id| self.devices.iter().find(|d| d.device_id == id))
                }
            }
            AdvancedDeviceSelection::MostEfficient => {
                // Use power efficiency if available
                if self.benchmarks.is_empty() {
                    self.devices.first()
                } else {
                    let most_efficient_id = self
                        .benchmarks
                        .iter()
                        .filter_map(|(id, bench)| {
                            bench.power_consumption.map(|power| {
                                let efficiency = bench.distance_throughput / power;
                                (*id, efficiency)
                            })
                        })
                        .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))
                        .map(|(id, _)| id);

                    most_efficient_id.and_then(|id| self.devices.iter().find(|d| d.device_id == id))
                }
            }
            AdvancedDeviceSelection::RoundRobin => {
                let idx = self.round_robin_idx % self.devices.len();
                self.round_robin_idx += 1;
                self.devices.get(idx)
            }
            AdvancedDeviceSelection::LoadBalanced => {
                // Select least utilized device
                let least_utilized = self
                    .utilization
                    .iter()
                    .min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
                    .map(|(id, _)| *id);

                least_utilized.and_then(|id| self.devices.iter().find(|d| d.device_id == id))
            }
            AdvancedDeviceSelection::MultiGpu { max_gpus } => {
                // For multi-GPU, return the first device (caller handles multi-GPU logic)
                self.devices.iter().take(*max_gpus).next()
            }
        }
    }

    /// Update device utilization
    pub fn update_utilization(&mut self, device_id: u32, utilization: f64) {
        self.utilization
            .insert(device_id, utilization.clamp(0.0, 1.0));
    }

    /// Add benchmark result
    pub fn add_benchmark(&mut self, benchmark: DeviceBenchmark) {
        self.benchmarks.insert(benchmark.device_id, benchmark);
    }

    /// Get all devices
    pub fn devices(&self) -> &[GpuDevice] {
        &self.devices
    }

    /// Get current strategy
    pub fn strategy(&self) -> &AdvancedDeviceSelection {
        &self.strategy
    }

    /// Set selection strategy
    pub fn set_strategy(&mut self, strategy: AdvancedDeviceSelection) {
        self.strategy = strategy;
    }
}

// ============================================================================
// GPU Acceleration Configuration
// ============================================================================

/// Comprehensive GPU acceleration configuration
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct GpuAccelerationConfig {
    /// Enable GPU acceleration
    pub enabled: bool,
    /// Preferred backend
    pub backend: GpuBackend,
    /// Device selection strategy
    pub device_selection: AdvancedDeviceSelection,
    /// Memory management strategy
    pub memory_strategy: AdvancedMemoryStrategy,
    /// Tensor core configuration
    pub tensor_cores: TensorCoreConfig,
    /// Enable automatic CPU fallback
    pub auto_fallback: bool,
    /// Minimum problem size for GPU acceleration
    pub min_problem_size: usize,
    /// Tile size for blocked algorithms
    pub tile_size: usize,
    /// Enable asynchronous execution
    pub async_execution: bool,
    /// Number of CUDA streams / OpenCL queues
    pub num_streams: usize,
    /// Enable profiling
    pub enable_profiling: bool,
    /// Custom kernel optimizations
    pub kernel_optimizations: KernelOptimizations,
}

impl Default for GpuAccelerationConfig {
    fn default() -> Self {
        Self {
            enabled: true,
            backend: GpuBackend::CpuFallback,
            device_selection: AdvancedDeviceSelection::Auto,
            memory_strategy: AdvancedMemoryStrategy::Adaptive,
            tensor_cores: TensorCoreConfig::default(),
            auto_fallback: true,
            min_problem_size: 1000,
            tile_size: 256,
            async_execution: true,
            num_streams: 4,
            enable_profiling: false,
            kernel_optimizations: KernelOptimizations::default(),
        }
    }
}

impl GpuAccelerationConfig {
    /// Create CUDA configuration
    pub fn cuda() -> Self {
        Self {
            backend: GpuBackend::Cuda,
            tensor_cores: TensorCoreConfig {
                enabled: true,
                precision: PrecisionMode::Mixed,
                ..Default::default()
            },
            ..Default::default()
        }
    }

    /// Create OpenCL configuration
    pub fn opencl() -> Self {
        Self {
            backend: GpuBackend::OpenCl,
            tensor_cores: TensorCoreConfig {
                enabled: false,
                precision: PrecisionMode::Single,
                ..Default::default()
            },
            ..Default::default()
        }
    }

    /// Create ROCm configuration
    pub fn rocm() -> Self {
        Self {
            backend: GpuBackend::Rocm,
            tensor_cores: TensorCoreConfig {
                enabled: true,
                precision: PrecisionMode::Mixed,
                ..Default::default()
            },
            ..Default::default()
        }
    }

    /// Create Metal configuration
    pub fn metal() -> Self {
        Self {
            backend: GpuBackend::Metal,
            tensor_cores: TensorCoreConfig {
                enabled: true,
                precision: PrecisionMode::Half,
                ..Default::default()
            },
            ..Default::default()
        }
    }

    /// Create CPU fallback configuration
    pub fn cpu() -> Self {
        Self {
            enabled: false,
            backend: GpuBackend::CpuFallback,
            ..Default::default()
        }
    }

    /// Convert to basic GpuConfig
    pub fn to_basic_config(&self) -> GpuConfig {
        GpuConfig {
            preferred_backend: self.backend,
            device_selection: self.device_selection.clone().into(),
            auto_fallback: self.auto_fallback,
            memory_pool_size: match self.memory_strategy {
                AdvancedMemoryStrategy::Pool { pool_size } => Some(pool_size),
                _ => None,
            },
            optimize_memory: true,
            backend_options: HashMap::new(),
        }
    }
}

/// Kernel optimization settings
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct KernelOptimizations {
    /// Use loop unrolling
    pub loop_unrolling: bool,
    /// Use shared memory tiling
    pub shared_memory_tiling: bool,
    /// Use register blocking
    pub register_blocking: bool,
    /// Use vectorized loads (e.g., float4)
    pub vectorized_loads: bool,
    /// Use texture memory for read-only data
    pub texture_memory: bool,
    /// Use constant memory for frequently accessed constants
    pub constant_memory: bool,
    /// Occupancy optimization level (0-3)
    pub occupancy_level: u8,
}

impl Default for KernelOptimizations {
    fn default() -> Self {
        Self {
            loop_unrolling: true,
            shared_memory_tiling: true,
            register_blocking: true,
            vectorized_loads: true,
            texture_memory: false,
            constant_memory: true,
            occupancy_level: 2,
        }
    }
}

// ============================================================================
// GPU Accelerated K-Means
// ============================================================================

/// GPU-accelerated K-means clustering
#[derive(Debug)]
pub struct GpuKMeans<F: Float> {
    /// Configuration
    config: GpuAccelerationConfig,
    /// GPU context (if available)
    context: Option<GpuContext>,
    /// Memory manager
    memory_manager: AdvancedGpuMemoryManager,
    /// Device selector
    device_selector: DeviceSelector,
    /// Tensor core capabilities
    tensor_caps: TensorCoreCapabilities,
    /// Is GPU actually available
    gpu_available: bool,
    /// Profiling data
    profiling_data: Vec<ProfilingRecord>,
    /// Phantom data for type
    _phantom: std::marker::PhantomData<F>,
}

/// Profiling record for GPU operations
#[derive(Debug, Clone)]
pub struct ProfilingRecord {
    /// Operation name
    pub operation: String,
    /// Duration in microseconds
    pub duration_us: u64,
    /// Memory transferred (bytes)
    pub memory_transferred: usize,
    /// Compute operations performed
    pub compute_ops: usize,
    /// Timestamp
    pub timestamp: Instant,
}

/// K-means result from GPU computation
#[derive(Debug, Clone)]
pub struct GpuKMeansResult<F: Float> {
    /// Final centroids
    pub centroids: Array2<F>,
    /// Cluster assignments for each point
    pub labels: Array1<usize>,
    /// Inertia (sum of squared distances to centroids)
    pub inertia: F,
    /// Number of iterations
    pub n_iterations: usize,
    /// Whether converged
    pub converged: bool,
    /// Computation metrics
    pub metrics: KMeansMetrics,
}

/// K-means computation metrics
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct KMeansMetrics {
    /// Total computation time (ms)
    pub total_time_ms: f64,
    /// Time for distance computation (ms)
    pub distance_time_ms: f64,
    /// Time for centroid update (ms)
    pub centroid_update_time_ms: f64,
    /// Time for label assignment (ms)
    pub label_assignment_time_ms: f64,
    /// Data transfer time (ms)
    pub transfer_time_ms: f64,
    /// Used GPU acceleration
    pub used_gpu: bool,
    /// Backend used
    pub backend: String,
    /// Memory used (bytes)
    pub memory_used: usize,
    /// Throughput (samples/second)
    pub throughput: f64,
}

impl<F: Float + FromPrimitive + Send + Sync + 'static> GpuKMeans<F> {
    /// Create new GPU-accelerated K-means
    pub fn new(config: GpuAccelerationConfig) -> Result<Self> {
        let device_selector = DeviceSelector::new(config.device_selection.clone());

        // Try to create GPU context
        let (context, gpu_available, tensor_caps) = Self::try_create_context(&config)?;

        let available_memory = context
            .as_ref()
            .map(|ctx| ctx.device.available_memory)
            .unwrap_or(1024 * 1024 * 1024); // 1GB default

        let memory_manager =
            AdvancedGpuMemoryManager::new(config.memory_strategy, available_memory);

        Ok(Self {
            config,
            context,
            memory_manager,
            device_selector,
            tensor_caps,
            gpu_available,
            profiling_data: Vec::new(),
            _phantom: std::marker::PhantomData,
        })
    }

    /// Try to create GPU context
    fn try_create_context(
        config: &GpuAccelerationConfig,
    ) -> Result<(Option<GpuContext>, bool, TensorCoreCapabilities)> {
        if !config.enabled || config.backend == GpuBackend::CpuFallback {
            return Ok((None, false, TensorCoreCapabilities::default()));
        }

        // Create device and context
        let device = GpuDevice::new(
            0,
            format!("{} Device", config.backend),
            8_000_000_000,
            6_000_000_000,
            "1.0".to_string(),
            1024,
            config.backend,
            true,
        );

        let tensor_caps = detect_tensor_core_capabilities(&device);
        let basic_config = config.to_basic_config();

        match GpuContext::new(device.clone(), basic_config) {
            Ok(ctx) => Ok((Some(ctx), true, tensor_caps)),
            Err(_) if config.auto_fallback => Ok((None, false, TensorCoreCapabilities::default())),
            Err(e) => Err(e),
        }
    }

    /// Fit K-means to data
    pub fn fit(
        &mut self,
        data: ArrayView2<F>,
        k: usize,
        max_iter: usize,
        tol: F,
    ) -> Result<GpuKMeansResult<F>> {
        let start_time = Instant::now();
        let n_samples = data.nrows();
        let n_features = data.ncols();

        // Validate inputs
        if k == 0 || k > n_samples {
            return Err(ClusteringError::InvalidInput(format!(
                "k must be between 1 and n_samples ({}), got {}",
                n_samples, k
            )));
        }

        // Decide whether to use GPU
        let use_gpu = self.should_use_gpu(n_samples, n_features);

        if use_gpu && self.gpu_available {
            self.fit_gpu(data, k, max_iter, tol, start_time)
        } else {
            self.fit_cpu(data, k, max_iter, tol, start_time)
        }
    }

    /// Decide whether to use GPU based on problem size
    fn should_use_gpu(&self, n_samples: usize, n_features: usize) -> bool {
        let problem_size = n_samples * n_features;
        problem_size >= self.config.min_problem_size && self.config.enabled
    }

    /// GPU implementation of K-means
    fn fit_gpu(
        &mut self,
        data: ArrayView2<F>,
        k: usize,
        max_iter: usize,
        tol: F,
        start_time: Instant,
    ) -> Result<GpuKMeansResult<F>> {
        let n_samples = data.nrows();
        let n_features = data.ncols();

        // Initialize centroids using K-means++
        let mut centroids = self.initialize_centroids_gpu(data, k)?;
        let mut labels = Array1::zeros(n_samples);
        let mut inertia = F::infinity();
        let mut converged = false;
        let mut n_iterations = 0;

        let mut distance_time = Duration::ZERO;
        let mut centroid_time = Duration::ZERO;
        let mut label_time = Duration::ZERO;

        // Main K-means loop
        for iter in 0..max_iter {
            n_iterations = iter + 1;

            // Step 1: Compute distances and assign labels
            let label_start = Instant::now();
            let (new_labels, distances) = self.compute_labels_gpu(data, centroids.view())?;
            labels = new_labels;
            label_time += label_start.elapsed();

            // Step 2: Compute new centroids
            let centroid_start = Instant::now();
            let new_centroids = self.compute_centroids_gpu(data, &labels, k)?;
            centroid_time += centroid_start.elapsed();

            // Step 3: Check convergence
            let new_inertia = self.compute_inertia(&distances);
            let centroid_shift =
                self.compute_centroid_shift(centroids.view(), new_centroids.view());

            centroids = new_centroids;

            if centroid_shift <= tol
                || (inertia - new_inertia).abs() < tol * F::from(0.01).unwrap_or(tol)
            {
                converged = true;
                inertia = new_inertia;
                break;
            }

            inertia = new_inertia;
        }

        let total_time = start_time.elapsed();

        let metrics = KMeansMetrics {
            total_time_ms: total_time.as_secs_f64() * 1000.0,
            distance_time_ms: distance_time.as_secs_f64() * 1000.0,
            centroid_update_time_ms: centroid_time.as_secs_f64() * 1000.0,
            label_assignment_time_ms: label_time.as_secs_f64() * 1000.0,
            transfer_time_ms: 0.0, // Would be populated with actual transfer times
            used_gpu: true,
            backend: format!("{}", self.config.backend),
            memory_used: self.memory_manager.get_stats().current_bytes_in_use,
            throughput: n_samples as f64 / total_time.as_secs_f64(),
        };

        Ok(GpuKMeansResult {
            centroids,
            labels,
            inertia,
            n_iterations,
            converged,
            metrics,
        })
    }

    /// CPU fallback implementation of K-means
    fn fit_cpu(
        &self,
        data: ArrayView2<F>,
        k: usize,
        max_iter: usize,
        tol: F,
        start_time: Instant,
    ) -> Result<GpuKMeansResult<F>> {
        let n_samples = data.nrows();
        let n_features = data.ncols();

        // Initialize centroids using K-means++
        let mut centroids = self.initialize_centroids_cpu(data, k)?;
        let mut labels = Array1::zeros(n_samples);
        let mut inertia = F::infinity();
        let mut converged = false;
        let mut n_iterations = 0;

        // Main K-means loop
        for iter in 0..max_iter {
            n_iterations = iter + 1;

            // Step 1: Assign labels
            let (new_labels, distances) = self.assign_labels_cpu(data, centroids.view())?;
            labels = new_labels;

            // Step 2: Update centroids
            let new_centroids = self.update_centroids_cpu(data, &labels, k, n_features)?;

            // Step 3: Check convergence
            let new_inertia = self.compute_inertia(&distances);
            let centroid_shift =
                self.compute_centroid_shift(centroids.view(), new_centroids.view());

            centroids = new_centroids;

            if centroid_shift <= tol {
                converged = true;
                inertia = new_inertia;
                break;
            }

            inertia = new_inertia;
        }

        let total_time = start_time.elapsed();

        let metrics = KMeansMetrics {
            total_time_ms: total_time.as_secs_f64() * 1000.0,
            distance_time_ms: 0.0,
            centroid_update_time_ms: 0.0,
            label_assignment_time_ms: 0.0,
            transfer_time_ms: 0.0,
            used_gpu: false,
            backend: "CPU".to_string(),
            memory_used: 0,
            throughput: n_samples as f64 / total_time.as_secs_f64(),
        };

        Ok(GpuKMeansResult {
            centroids,
            labels,
            inertia,
            n_iterations,
            converged,
            metrics,
        })
    }

    /// Initialize centroids using K-means++ on GPU
    fn initialize_centroids_gpu(&self, data: ArrayView2<F>, k: usize) -> Result<Array2<F>> {
        // For now, use CPU initialization (GPU K-means++ is complex)
        self.initialize_centroids_cpu(data, k)
    }

    /// Initialize centroids using K-means++ on CPU
    fn initialize_centroids_cpu(&self, data: ArrayView2<F>, k: usize) -> Result<Array2<F>> {
        let n_samples = data.nrows();
        let n_features = data.ncols();
        let mut centroids = Array2::zeros((k, n_features));
        let mut rng = scirs2_core::random::rng();

        // Choose first centroid randomly
        let first_idx = scirs2_core::random::RngExt::random_range(&mut rng, 0..n_samples);
        for j in 0..n_features {
            centroids[[0, j]] = data[[first_idx, j]];
        }

        if k == 1 {
            return Ok(centroids);
        }

        // Choose remaining centroids with K-means++
        let mut min_distances = Array1::from_elem(n_samples, F::infinity());

        for i in 1..k {
            // Update minimum distances
            for sample_idx in 0..n_samples {
                let dist =
                    self.euclidean_distance_squared(data.row(sample_idx), centroids.row(i - 1));
                if dist < min_distances[sample_idx] {
                    min_distances[sample_idx] = dist;
                }
            }

            // Compute probability distribution
            let sum_distances: F = min_distances.iter().copied().fold(F::zero(), |a, b| a + b);
            if sum_distances <= F::zero() {
                // All points are at centroids, pick random
                let idx = scirs2_core::random::RngExt::random_range(&mut rng, 0..n_samples);
                for j in 0..n_features {
                    centroids[[i, j]] = data[[idx, j]];
                }
                continue;
            }

            // Sample next centroid
            let threshold = F::from(scirs2_core::random::RngExt::random_range(
                &mut rng,
                0.0..1.0,
            ))
            .unwrap_or(F::zero())
                * sum_distances;
            let mut cumsum = F::zero();
            let mut next_idx = 0;

            for (idx, &dist) in min_distances.iter().enumerate() {
                cumsum = cumsum + dist;
                if cumsum >= threshold {
                    next_idx = idx;
                    break;
                }
            }

            for j in 0..n_features {
                centroids[[i, j]] = data[[next_idx, j]];
            }
        }

        Ok(centroids)
    }

    /// Compute labels using GPU acceleration
    fn compute_labels_gpu(
        &self,
        data: ArrayView2<F>,
        centroids: ArrayView2<F>,
    ) -> Result<(Array1<usize>, Array1<F>)> {
        // GPU-accelerated label assignment
        // For now, use CPU implementation with the structure for GPU
        self.assign_labels_cpu(data, centroids)
    }

    /// Compute centroids using GPU acceleration
    fn compute_centroids_gpu(
        &self,
        data: ArrayView2<F>,
        labels: &Array1<usize>,
        k: usize,
    ) -> Result<Array2<F>> {
        let n_features = data.ncols();
        self.update_centroids_cpu(data, labels, k, n_features)
    }

    /// Assign labels to each point (CPU implementation)
    fn assign_labels_cpu(
        &self,
        data: ArrayView2<F>,
        centroids: ArrayView2<F>,
    ) -> Result<(Array1<usize>, Array1<F>)> {
        let n_samples = data.nrows();
        let n_centroids = centroids.nrows();
        let mut labels = Array1::zeros(n_samples);
        let mut distances = Array1::zeros(n_samples);

        for i in 0..n_samples {
            let mut min_dist = F::infinity();
            let mut min_label = 0;

            for j in 0..n_centroids {
                let dist = self.euclidean_distance_squared(data.row(i), centroids.row(j));
                if dist < min_dist {
                    min_dist = dist;
                    min_label = j;
                }
            }

            labels[i] = min_label;
            distances[i] = min_dist;
        }

        Ok((labels, distances))
    }

    /// Update centroids (CPU implementation)
    fn update_centroids_cpu(
        &self,
        data: ArrayView2<F>,
        labels: &Array1<usize>,
        k: usize,
        n_features: usize,
    ) -> Result<Array2<F>> {
        let mut centroids = Array2::zeros((k, n_features));
        let mut counts = vec![0usize; k];

        // Sum points in each cluster
        for (i, &label) in labels.iter().enumerate() {
            if label < k {
                for j in 0..n_features {
                    centroids[[label, j]] = centroids[[label, j]] + data[[i, j]];
                }
                counts[label] += 1;
            }
        }

        // Divide by counts
        for i in 0..k {
            if counts[i] > 0 {
                let count = F::from(counts[i]).unwrap_or(F::one());
                for j in 0..n_features {
                    centroids[[i, j]] = centroids[[i, j]] / count;
                }
            }
        }

        Ok(centroids)
    }

    /// Compute squared Euclidean distance
    fn euclidean_distance_squared(
        &self,
        a: scirs2_core::ndarray::ArrayView1<F>,
        b: scirs2_core::ndarray::ArrayView1<F>,
    ) -> F {
        a.iter()
            .zip(b.iter())
            .map(|(&x, &y)| {
                let diff = x - y;
                diff * diff
            })
            .fold(F::zero(), |acc, x| acc + x)
    }

    /// Compute inertia from distances
    fn compute_inertia(&self, distances: &Array1<F>) -> F {
        distances.iter().copied().fold(F::zero(), |a, b| a + b)
    }

    /// Compute centroid shift
    fn compute_centroid_shift(&self, old: ArrayView2<F>, new: ArrayView2<F>) -> F {
        let mut max_shift = F::zero();
        for i in 0..old.nrows() {
            let shift = self
                .euclidean_distance_squared(old.row(i), new.row(i))
                .sqrt();
            if shift > max_shift {
                max_shift = shift;
            }
        }
        max_shift
    }

    /// Get configuration
    pub fn config(&self) -> &GpuAccelerationConfig {
        &self.config
    }

    /// Check if GPU is available
    pub fn is_gpu_available(&self) -> bool {
        self.gpu_available
    }

    /// Get tensor core capabilities
    pub fn tensor_core_capabilities(&self) -> &TensorCoreCapabilities {
        &self.tensor_caps
    }

    /// Get memory statistics
    pub fn memory_stats(&self) -> &MemoryUsageStats {
        self.memory_manager.get_stats()
    }

    /// Get profiling data
    pub fn profiling_data(&self) -> &[ProfilingRecord] {
        &self.profiling_data
    }
}

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

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

    #[test]
    fn test_advanced_memory_strategy_display() {
        assert_eq!(
            AdvancedMemoryStrategy::Conservative.to_string(),
            "Conservative"
        );
        assert_eq!(
            AdvancedMemoryStrategy::Streaming {
                chunk_size: 1024 * 1024
            }
            .to_string(),
            "Streaming(1MB)"
        );
    }

    #[test]
    fn test_advanced_memory_manager_creation() {
        let manager = AdvancedGpuMemoryManager::new(
            AdvancedMemoryStrategy::Adaptive,
            4 * 1024 * 1024 * 1024, // 4GB
        );
        assert_eq!(manager.strategy(), AdvancedMemoryStrategy::Adaptive);
    }

    #[test]
    fn test_advanced_memory_allocation() {
        let mut manager = AdvancedGpuMemoryManager::new(
            AdvancedMemoryStrategy::Conservative,
            1024 * 1024 * 1024, // 1GB
        );

        let result = manager.allocate(1024);
        assert!(result.is_ok());

        let stats = manager.get_stats();
        assert_eq!(stats.total_allocations, 1);
        assert_eq!(stats.successful_allocations, 1);
    }

    #[test]
    fn test_precision_mode_display() {
        assert_eq!(PrecisionMode::Mixed.to_string(), "Mixed (f16/f32)");
        assert_eq!(PrecisionMode::TensorFloat32.to_string(), "TF32");
    }

    #[test]
    fn test_tensor_core_config_default() {
        let config = TensorCoreConfig::default();
        assert!(config.enabled);
        assert_eq!(config.precision, PrecisionMode::Auto);
        assert!(config.auto_scale);
    }

    #[test]
    fn test_device_selector_creation() {
        let selector = DeviceSelector::new(AdvancedDeviceSelection::Auto);
        assert!(selector.devices().is_empty());
    }

    #[test]
    fn test_device_selector_add_device() {
        let mut selector = DeviceSelector::new(AdvancedDeviceSelection::MostMemory);

        let device = GpuDevice::new(
            0,
            "Test GPU".to_string(),
            8_000_000_000,
            6_000_000_000,
            "1.0".to_string(),
            1024,
            GpuBackend::Cuda,
            true,
        );

        selector.add_device(device);
        assert_eq!(selector.devices().len(), 1);
    }

    #[test]
    fn test_gpu_acceleration_config_default() {
        let config = GpuAccelerationConfig::default();
        assert!(config.enabled);
        assert!(config.auto_fallback);
    }

    #[test]
    fn test_gpu_acceleration_config_cuda() {
        let config = GpuAccelerationConfig::cuda();
        assert_eq!(config.backend, GpuBackend::Cuda);
        assert!(config.tensor_cores.enabled);
    }

    #[test]
    fn test_gpu_kmeans_creation() {
        let config = GpuAccelerationConfig::cpu();
        let kmeans = GpuKMeans::<f64>::new(config);
        assert!(kmeans.is_ok());
    }

    #[test]
    fn test_gpu_kmeans_fit_cpu_fallback() {
        let config = GpuAccelerationConfig::cpu();
        let mut kmeans = GpuKMeans::<f64>::new(config).expect("Failed to create GpuKMeans");

        // Create test data with two clear clusters
        let data = Array2::from_shape_vec(
            (6, 2),
            vec![1.0, 2.0, 1.2, 1.8, 0.8, 1.9, 4.0, 5.0, 4.2, 4.8, 3.9, 5.1],
        )
        .expect("Failed to create test data");

        let result = kmeans.fit(data.view(), 2, 100, 1e-4);
        assert!(result.is_ok());

        let result = result.expect("Failed to fit");
        assert_eq!(result.centroids.nrows(), 2);
        assert_eq!(result.labels.len(), 6);
        assert!(!result.metrics.used_gpu);
    }

    #[test]
    fn test_gpu_kmeans_convergence() {
        let config = GpuAccelerationConfig::cpu();
        let mut kmeans = GpuKMeans::<f64>::new(config).expect("Failed to create GpuKMeans");

        // Well-separated clusters should converge quickly
        let data = Array2::from_shape_vec(
            (8, 2),
            vec![
                0.0, 0.0, 0.1, 0.1, 0.0, 0.1, 0.1, 0.0, 10.0, 10.0, 10.1, 10.1, 10.0, 10.1, 10.1,
                10.0,
            ],
        )
        .expect("Failed to create test data");

        let result = kmeans.fit(data.view(), 2, 100, 1e-6);
        assert!(result.is_ok());

        let result = result.expect("Failed to fit");
        assert!(result.converged);
        assert!(result.n_iterations < 50);
    }

    #[test]
    fn test_memory_usage_stats() {
        let mut manager =
            AdvancedGpuMemoryManager::new(AdvancedMemoryStrategy::Aggressive, 1024 * 1024 * 1024);

        // Multiple allocations
        for _ in 0..5 {
            let _ = manager.allocate(1024);
        }

        let stats = manager.get_stats();
        assert_eq!(stats.total_allocations, 5);
        assert!(stats.efficiency > 0.0);
    }

    #[test]
    fn test_kernel_optimizations_default() {
        let opts = KernelOptimizations::default();
        assert!(opts.loop_unrolling);
        assert!(opts.shared_memory_tiling);
        assert_eq!(opts.occupancy_level, 2);
    }

    #[test]
    fn test_detect_tensor_core_capabilities() {
        let cuda_device = GpuDevice::new(
            0,
            "CUDA Device".to_string(),
            8_000_000_000,
            6_000_000_000,
            "8.0".to_string(),
            1024,
            GpuBackend::Cuda,
            true,
        );

        let caps = detect_tensor_core_capabilities(&cuda_device);
        assert!(caps.available);
        assert!(!caps.supported_precisions.is_empty());
    }

    #[test]
    fn test_profiling_record_creation() {
        let record = ProfilingRecord {
            operation: "distance_compute".to_string(),
            duration_us: 1000,
            memory_transferred: 1024 * 1024,
            compute_ops: 1000000,
            timestamp: Instant::now(),
        };

        assert_eq!(record.operation, "distance_compute");
        assert_eq!(record.duration_us, 1000);
    }
}