zipora 3.1.4

//! Adaptive Strategy Selection for Rank/Select Operations
//!
//! This module implements adaptive selection of optimal rank/select implementations
//! based on data density analysis, inspired by advanced adaptive systems.
//!
//! The adaptive system analyzes data characteristics and automatically selects
//! the most appropriate implementation:
//!
//! - **Sparse data (< 5% density)**: Uses `RankSelectFew` for space efficiency
//! - **Dense data (> 90% density)**: Uses optimized dense implementations  
//! - **Balanced data**: Uses general-purpose implementations with optimal block sizes
//! - **Multi-dimensional**: Uses mixed implementations for related bit vectors
//!
//! # Examples
//!
//! ```rust
//! use zipora::succinct::{BitVector, rank_select::{AdaptiveRankSelect, RankSelectOps}};
//!
//! // Create sparse bit vector (every 100th bit set)
//! let mut sparse_bv = BitVector::new();
//! for i in 0..10000 {
//!     sparse_bv.push(i % 100 == 0).unwrap();
//! }
//!
//! // Adaptive selection automatically chooses RankSelectFew
//! let adaptive_rs = AdaptiveRankSelect::new(sparse_bv).unwrap();
//! println!("Selected: {}", adaptive_rs.implementation_name());
//! 
//! let rank = adaptive_rs.rank1(5000);
//! let pos = adaptive_rs.select1(25).unwrap();
//! ```

use crate::error::{Result, ZiporaError};

/// Runtime BMI2 capabilities detection and caching
#[derive(Debug, Clone)]
pub struct AdaptiveCpuCapabilities {
    /// BMI1 instruction set available (TZCNT, LZCNT, ANDN)
    pub has_bmi1: bool,
    /// BMI2 instruction set available (PDEP, PEXT, BZHI)
    pub has_bmi2: bool,
    /// POPCNT instruction available
    pub has_popcnt: bool,
    /// AVX2 instruction set available
    pub has_avx2: bool,
    /// Optimization tier (0=none, 1=basic, 2=BMI1, 3=BMI2, 4=BMI2+AVX2)
    pub optimization_tier: u8,
    /// Recommended chunk size for bulk operations
    pub chunk_size: usize,
}

impl AdaptiveCpuCapabilities {
    /// Detect available CPU features and determine optimization tier
    pub fn detect() -> Self {
        let has_bmi1 = Self::detect_bmi1();
        let has_bmi2 = Self::detect_bmi2();
        let has_popcnt = Self::detect_popcnt();
        let has_avx2 = Self::detect_avx2();
        
        let optimization_tier = match (has_popcnt, has_bmi1, has_bmi2, has_avx2) {
            (true, true, true, true) => 4,  // Full BMI2 + AVX2
            (true, true, true, false) => 3, // BMI2 without AVX2
            (true, true, false, _) => 2,    // BMI1 only
            (true, false, false, _) => 1,   // POPCNT only
            _ => 0,                          // No acceleration
        };
        
        let chunk_size = match optimization_tier {
            4 => 1024,  // Large chunks for AVX2
            3 => 512,   // Medium chunks for BMI2
            2 => 256,   // Smaller chunks for BMI1
            1 => 128,   // Small chunks for POPCNT
            _ => 64,    // Minimal chunks for scalar
        };
        
        Self {
            has_bmi1,
            has_bmi2,
            has_popcnt,
            has_avx2,
            optimization_tier,
            chunk_size,
        }
    }
    
    /// Get cached capabilities (thread-safe singleton)
    #[inline]
    pub fn get() -> &'static Self {
        static CAPABILITIES: std::sync::OnceLock<AdaptiveCpuCapabilities> = std::sync::OnceLock::new();
        CAPABILITIES.get_or_init(Self::detect)
    }
    
    #[cfg(target_arch = "x86_64")]
    fn detect_bmi1() -> bool {
        std::is_x86_feature_detected!("bmi1")
    }
    
    #[cfg(not(target_arch = "x86_64"))]
    fn detect_bmi1() -> bool {
        false
    }
    
    #[cfg(target_arch = "x86_64")]
    fn detect_bmi2() -> bool {
        std::is_x86_feature_detected!("bmi2")
    }
    
    #[cfg(not(target_arch = "x86_64"))]
    fn detect_bmi2() -> bool {
        false
    }
    
    #[cfg(target_arch = "x86_64")]
    fn detect_popcnt() -> bool {
        std::is_x86_feature_detected!("popcnt")
    }
    
    #[cfg(not(target_arch = "x86_64"))]
    fn detect_popcnt() -> bool {
        false
    }
    
    #[cfg(target_arch = "x86_64")]
    fn detect_avx2() -> bool {
        std::is_x86_feature_detected!("avx2")
    }
    
    #[cfg(not(target_arch = "x86_64"))]
    fn detect_avx2() -> bool {
        false
    }
}
use crate::succinct::BitVector;
use super::{
    RankSelectOps, RankSelectInterleaved256,
};
use std::fmt;

/// Data characteristics profile for adaptive selection
#[derive(Debug, Clone)]
pub struct DataProfile {
    /// Total number of bits
    pub total_bits: usize,
    /// Number of set bits (ones)
    pub ones_count: usize,
    /// Bit density (ratio of ones to total bits)
    pub density: f64,
    /// Estimated access pattern preference
    pub access_pattern: AccessPattern,
    /// Data size category
    pub size_category: SizeCategory,
    /// Pattern complexity score (0.0 = very regular, 1.0 = random)
    pub pattern_complexity: f64,
    /// Clustering coefficient for spatial locality
    pub clustering_coefficient: f64,
    /// Entropy of bit distribution (bits)
    pub entropy: f64,
    /// Run length statistics
    pub run_length_stats: RunLengthStats,
    /// Hardware acceleration tier available
    pub hardware_tier: u8,
}

/// Expected access pattern for optimization
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum AccessPattern {
    /// Mixed rank and select operations
    Mixed,
    /// Primarily rank operations
    RankHeavy,
    /// Primarily select operations  
    SelectHeavy,
    /// Sequential access pattern
    Sequential,
    /// Random access pattern
    Random,
}

/// Data size categories for implementation selection
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum SizeCategory {
    /// Small datasets (< 10K bits)
    Small,
    /// Medium datasets (10K - 1M bits)
    Medium,
    /// Large datasets (1M - 100M bits)
    Large,
    /// Very large datasets (> 100M bits)
    VeryLarge,
}

/// Run length statistics for pattern analysis
#[derive(Debug, Clone)]
pub struct RunLengthStats {
    /// Average run length of ones
    pub avg_ones_run: f64,
    /// Average run length of zeros
    pub avg_zeros_run: f64,
    /// Maximum run length of ones
    pub max_ones_run: usize,
    /// Maximum run length of zeros
    pub max_zeros_run: usize,
    /// Number of alternations (pattern changes)
    pub alternations: usize,
}

/// Strategy selection criteria based on optimization patterns
#[derive(Debug, Clone)]
pub struct SelectionCriteria {
    /// Sparsity threshold for RankSelectFew (adaptive: 0.02-0.15 based on complexity)
    pub sparse_threshold: f64,
    /// Dense threshold for dense optimizations (adaptive: 0.85-0.95 based on patterns)
    pub dense_threshold: f64,
    /// Small dataset threshold in bits (default: 10,000)
    pub small_dataset_threshold: usize,
    /// Large dataset threshold in bits (default: 1,000,000)
    pub large_dataset_threshold: usize,
    /// Very large dataset threshold in bits (default: 100,000,000)
    pub very_large_dataset_threshold: usize,
    /// Preferred access pattern for optimization
    pub access_pattern: AccessPattern,
    /// Enable select cache optimization
    pub enable_select_cache: bool,
    /// Prefer space efficiency over speed
    pub prefer_space: bool,
    /// Enable automatic threshold tuning based on data patterns
    pub enable_adaptive_thresholds: bool,
    /// Hardware acceleration preference (0=none, 1=basic, 2=BMI1, 3=BMI2, 4=BMI2+AVX2)
    pub min_hardware_tier: u8,
    /// Pattern complexity weight for selection (0.0-1.0)
    pub pattern_complexity_weight: f64,
    /// Clustering weight for spatial locality (0.0-1.0)
    pub clustering_weight: f64,
}

impl Default for SelectionCriteria {
    fn default() -> Self {
        Self {
            sparse_threshold: 0.05,
            dense_threshold: 0.90,
            small_dataset_threshold: 10_000,
            large_dataset_threshold: 1_000_000,
            very_large_dataset_threshold: 100_000_000,
            access_pattern: AccessPattern::Mixed,
            enable_select_cache: true,
            prefer_space: false,
            enable_adaptive_thresholds: true,
            min_hardware_tier: 0,
            pattern_complexity_weight: 0.3,
            clustering_weight: 0.2,
        }
    }
}

impl Default for RunLengthStats {
    fn default() -> Self {
        Self {
            avg_ones_run: 1.0,
            avg_zeros_run: 1.0,
            max_ones_run: 1,
            max_zeros_run: 1,
            alternations: 0,
        }
    }
}

/// Adaptive rank/select implementation that automatically selects optimal strategy.
///
/// Currently always uses `RankSelectInterleaved256` (verified as best performer
/// at 121-302 Mops/s, 50-150x over other implementations). Direct storage
/// eliminates vtable indirection on hot-path rank/select operations.
pub struct AdaptiveRankSelect {
    /// The actual implementation — always RankSelectInterleaved256
    implementation: RankSelectInterleaved256,
    /// Name of the selected implementation
    implementation_name: String,
    /// Data profile used for selection
    profile: DataProfile,
    /// Selection criteria used
    criteria: SelectionCriteria,
}

impl AdaptiveRankSelect {
    /// Create new adaptive rank/select with automatic strategy selection
    pub fn new(bit_vector: BitVector) -> Result<Self> {
        Self::with_criteria(bit_vector, SelectionCriteria::default())
    }

    /// Create new adaptive rank/select with custom selection criteria
    pub fn with_criteria(bit_vector: BitVector, criteria: SelectionCriteria) -> Result<Self> {
        let profile = Self::analyze_data(&bit_vector, &criteria);
        let (implementation, name) = Self::select_implementation(bit_vector, &profile, &criteria)?;
        
        Ok(Self {
            implementation,
            implementation_name: name,
            profile,
            criteria,
        })
    }

    /// Analyze bit vector characteristics for strategy selection
    fn analyze_data(bit_vector: &BitVector, criteria: &SelectionCriteria) -> DataProfile {
        let total_bits = bit_vector.len();
        let ones_count = bit_vector.count_ones();
        let density = if total_bits == 0 { 0.0 } else { ones_count as f64 / total_bits as f64 };

        let size_category = match total_bits {
            0..=9999 => SizeCategory::Small,
            10000..=999999 => SizeCategory::Medium,
            1000000..=99999999 => SizeCategory::Large,
            _ => SizeCategory::VeryLarge,
        };

        // Advanced pattern analysis
        let run_length_stats = Self::analyze_run_lengths(bit_vector);
        let pattern_complexity = Self::calculate_pattern_complexity(bit_vector, &run_length_stats);
        let clustering_coefficient = Self::calculate_clustering_coefficient(bit_vector);
        let entropy = Self::calculate_entropy(bit_vector);
        
        // Detect hardware acceleration capabilities
        let hardware_tier = AdaptiveCpuCapabilities::get().optimization_tier;
        
        // Use provided access pattern or infer from data characteristics
        let access_pattern = Self::infer_access_pattern(criteria.access_pattern, &run_length_stats, pattern_complexity);

        DataProfile {
            total_bits,
            ones_count,
            density,
            access_pattern,
            size_category,
            pattern_complexity,
            clustering_coefficient,
            entropy,
            run_length_stats,
            hardware_tier,
        }
    }
    
    /// Analyze run lengths to understand bit patterns
    fn analyze_run_lengths(bit_vector: &BitVector) -> RunLengthStats {
        if bit_vector.len() == 0 {
            return RunLengthStats::default();
        }
        
        let mut ones_runs = Vec::new();
        let mut zeros_runs = Vec::new();
        let mut current_run = 1;
        let mut current_bit = bit_vector.get(0).unwrap_or(false);
        let mut alternations = 0;
        
        for i in 1..bit_vector.len() {
            let bit = bit_vector.get(i).unwrap_or(false);
            if bit == current_bit {
                current_run += 1;
            } else {
                // Run ended, record it
                if current_bit {
                    ones_runs.push(current_run);
                } else {
                    zeros_runs.push(current_run);
                }
                current_run = 1;
                current_bit = bit;
                alternations += 1;
            }
        }
        
        // Record final run
        if current_bit {
            ones_runs.push(current_run);
        } else {
            zeros_runs.push(current_run);
        }
        
        let avg_ones_run = if ones_runs.is_empty() { 0.0 } else {
            ones_runs.iter().sum::<usize>() as f64 / ones_runs.len() as f64
        };
        
        let avg_zeros_run = if zeros_runs.is_empty() { 0.0 } else {
            zeros_runs.iter().sum::<usize>() as f64 / zeros_runs.len() as f64
        };
        
        RunLengthStats {
            avg_ones_run,
            avg_zeros_run,
            max_ones_run: ones_runs.iter().max().copied().unwrap_or(0),
            max_zeros_run: zeros_runs.iter().max().copied().unwrap_or(0),
            alternations,
        }
    }
    
    /// Calculate pattern complexity score (0.0 = very regular, 1.0 = random)
    fn calculate_pattern_complexity(bit_vector: &BitVector, run_stats: &RunLengthStats) -> f64 {
        if bit_vector.len() == 0 {
            return 0.0;
        }
        
        // Factors that increase complexity:
        // 1. High number of alternations relative to size
        // 2. Variance in run lengths
        // 3. Lack of clear patterns
        
        let total_bits = bit_vector.len() as f64;
        let alternation_density = run_stats.alternations as f64 / total_bits;
        
        // Calculate run length variance (high variance = more complex)
        let avg_run_length = (run_stats.avg_ones_run + run_stats.avg_zeros_run) / 2.0;
        let max_run_length = run_stats.max_ones_run.max(run_stats.max_zeros_run) as f64;
        
        let run_variance = if avg_run_length > 0.0 {
            (max_run_length - avg_run_length) / avg_run_length
        } else {
            0.0
        };
        
        // Combine factors (normalized to 0-1)
        let complexity = (alternation_density * 2.0 + run_variance.min(1.0)) / 3.0;
        complexity.min(1.0)
    }
    
    /// Calculate clustering coefficient for spatial locality
    fn calculate_clustering_coefficient(bit_vector: &BitVector) -> f64 {
        if bit_vector.len() < 3 {
            return 1.0; // Perfect clustering for small data
        }
        
        let mut clusters = 0;
        let window_size = 64; // Analyze in 64-bit windows
        
        for start in (0..bit_vector.len()).step_by(window_size) {
            let end = (start + window_size).min(bit_vector.len());
            let mut local_transitions = 0;
            
            for i in start..end-1 {
                if bit_vector.get(i) != bit_vector.get(i + 1) {
                    local_transitions += 1;
                }
            }
            
            if local_transitions < window_size / 4 {
                clusters += 1;
            }
        }
        
        let num_windows = (bit_vector.len() + window_size - 1) / window_size;
        if num_windows == 0 { 1.0 } else { clusters as f64 / num_windows as f64 }
    }
    
    /// Calculate Shannon entropy of bit distribution
    fn calculate_entropy(bit_vector: &BitVector) -> f64 {
        if bit_vector.len() == 0 {
            return 0.0;
        }
        
        let ones = bit_vector.count_ones() as f64;
        let zeros = (bit_vector.len() - bit_vector.count_ones()) as f64;
        let total = bit_vector.len() as f64;
        
        if ones == 0.0 || zeros == 0.0 {
            return 0.0; // No entropy for all-same bits
        }
        
        let p1 = ones / total;
        let p0 = zeros / total;
        
        -(p1 * p1.log2() + p0 * p0.log2())
    }
    
    /// Infer optimal access pattern from data characteristics
    fn infer_access_pattern(
        provided: AccessPattern,
        run_stats: &RunLengthStats,
        complexity: f64,
    ) -> AccessPattern {
        match provided {
            AccessPattern::Mixed => {
                // Infer pattern from data characteristics
                if run_stats.avg_ones_run > 32.0 || run_stats.avg_zeros_run > 32.0 {
                    AccessPattern::Sequential
                } else if complexity < 0.3 {
                    AccessPattern::Sequential
                } else if complexity > 0.7 {
                    AccessPattern::Random
                } else {
                    AccessPattern::Mixed
                }
            }
            _ => provided, // Use explicitly provided pattern
        }
    }

    /// Select optimal implementation based on data profile (simplified to use best performer)
    fn select_implementation(
        bit_vector: BitVector,
        profile: &DataProfile,
        criteria: &SelectionCriteria,
    ) -> Result<(RankSelectInterleaved256, String)> {
        let description = format!("RankSelectInterleaved256 ({})",
            Self::get_description(profile, criteria));

        Ok((
            RankSelectInterleaved256::new(bit_vector)?,
            description,
        ))
    }

    /// Get descriptive text for the selected configuration
    fn get_description(profile: &DataProfile, criteria: &SelectionCriteria) -> String {
        let (adaptive_sparse_threshold, adaptive_dense_threshold) = if criteria.enable_adaptive_thresholds {
            Self::tune_thresholds(profile, criteria)
        } else {
            (criteria.sparse_threshold, criteria.dense_threshold)
        };

        if profile.ones_count == 0 {
            "all zeros"
        } else if profile.ones_count == profile.total_bits {
            "all ones"
        } else if profile.density < adaptive_sparse_threshold {
            if profile.ones_count * 2 < profile.total_bits {
                "sparse ones optimized"
            } else {
                "sparse zeros optimized"
            }
        } else if profile.density > adaptive_dense_threshold {
            "dense data optimized"
        } else {
            match (profile.size_category, profile.access_pattern) {
                (SizeCategory::Small, _) => "small dataset optimized",
                (SizeCategory::Medium, AccessPattern::Sequential) => "medium sequential optimized",
                (SizeCategory::Medium, AccessPattern::Mixed) => "medium mixed optimized",
                (SizeCategory::Medium, _) => "medium dataset optimized",
                (SizeCategory::Large, AccessPattern::Sequential) => "large sequential optimized",
                (SizeCategory::Large, _) => "large dataset optimized",
                (SizeCategory::VeryLarge, _) => "very large dataset optimized",
            }
        }.to_string()
    }
    
    /// Tune thresholds based on data pattern analysis
    fn tune_thresholds(profile: &DataProfile, criteria: &SelectionCriteria) -> (f64, f64) {
        let mut sparse_threshold = criteria.sparse_threshold;
        let mut dense_threshold = criteria.dense_threshold;

        // CRITICAL OPTIMIZATION: Apply referenced project's cache-aware threshold strategy
        // referenced project uses Q=1,4 parameters for different cache optimization levels
        let cache_optimization_level = Self::determine_cache_level(profile);

        match cache_optimization_level {
            1 => {
                // Q=1: Linear search optimization for small ranges (≤32)
                // Better branch prediction, prefer implementations with linear search paths
                sparse_threshold *= 1.5; // More aggressive sparse selection
                dense_threshold *= 0.7;   // Earlier switch to dense for cache efficiency
            },
            4 => {
                // Q=4: Binary search optimization for larger ranges
                // Better cache locality, prefer implementations with binary search
                sparse_threshold *= 0.8;  // Less aggressive sparse selection
                dense_threshold *= 1.2;   // Later switch to dense to leverage binary search
            },
            _ => {
                // Default: balanced approach
                sparse_threshold *= 1.0;
                dense_threshold *= 1.0;
            }
        }

        // Adjust thresholds based on pattern complexity
        let complexity_factor = profile.pattern_complexity * criteria.pattern_complexity_weight;
        let clustering_factor = profile.clustering_coefficient * criteria.clustering_weight;

        // For highly clustered data, increase sparse threshold (more likely to benefit from sparse representation)
        if profile.clustering_coefficient > 0.7 {
            sparse_threshold *= 1.0 + clustering_factor;
        }

        // For complex patterns, be more conservative about sparse representation
        if profile.pattern_complexity > 0.6 {
            sparse_threshold *= 1.0 - complexity_factor;
        }

        // Adjust based on run length characteristics - referenced project optimized
        if profile.run_length_stats.avg_ones_run > 32.0 || profile.run_length_stats.avg_zeros_run > 32.0 {
            // Use referenced project's 32-element threshold for linear vs binary search
            dense_threshold *= 0.85; // More aggressive dense for longer runs
        }

        // Hardware acceleration can handle denser data more efficiently
        if profile.hardware_tier >= 3 { // BMI2 available
            sparse_threshold *= 1.3; // Can handle denser sparse data with BMI2
            dense_threshold *= 0.8;   // More aggressive dense optimization with BMI2
        }
        
        // Clamp thresholds to reasonable ranges
        sparse_threshold = sparse_threshold.clamp(0.01, 0.25);
        dense_threshold = dense_threshold.clamp(0.75, 0.98);
        
        (sparse_threshold, dense_threshold)
    }

    /// Determine cache optimization level based on referenced project's Q parameter strategy
    fn determine_cache_level(profile: &DataProfile) -> u8 {
        // referenced project uses Q=1 for linear search optimization (≤32 elements)
        // and Q=4 for binary search optimization (larger ranges)

        // Estimate typical access range based on data characteristics
        let estimated_range = if profile.clustering_coefficient > 0.7 {
            // Highly clustered data: smaller effective ranges
            (profile.ones_count as f64 / profile.clustering_coefficient).min(128.0) as usize
        } else {
            // Scattered data: larger ranges
            (profile.total_bits as f64 / (profile.ones_count as f64 + 1.0)).min(512.0) as usize
        };

        // Apply referenced project's threshold strategy
        if estimated_range <= 32 {
            // Q=1: Linear search optimization for small ranges
            // Better branch prediction, prefer linear algorithms
            1
        } else if estimated_range <= 128 {
            // Q=2: Hybrid approach for medium ranges
            2
        } else {
            // Q=4: Binary search optimization for large ranges
            // Better cache locality, prefer binary search algorithms
            4
        }
    }

    /// Get the name of the selected implementation
    pub fn implementation_name(&self) -> &str {
        &self.implementation_name
    }

    /// Get the data profile used for selection
    pub fn data_profile(&self) -> &DataProfile {
        &self.profile
    }

    /// Get the selection criteria used
    pub fn selection_criteria(&self) -> &SelectionCriteria {
        &self.criteria
    }

    /// Get optimization statistics
    pub fn optimization_stats(&self) -> OptimizationStats {
        OptimizationStats {
            total_bits: self.profile.total_bits,
            ones_count: self.profile.ones_count,
            density: self.profile.density,
            implementation: self.implementation_name.clone(),
            space_overhead_percent: self.implementation.space_overhead_percent(),
            estimated_performance_tier: self.estimate_performance_tier(),
        }
    }

    /// Estimate performance tier of selected implementation
    fn estimate_performance_tier(&self) -> PerformanceTier {
        // Since we always use RankSelectInterleaved256 (the best performer),
        // always return HighPerformance tier
        PerformanceTier::HighPerformance
    }
}

/// Performance tier classification
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum PerformanceTier {
    /// Optimized for space efficiency
    SpaceOptimized,
    /// High-performance for mixed operations
    HighPerformance,
    /// Optimized for sequential access
    Sequential,
    /// Balanced space and performance
    Balanced,
}

/// Optimization statistics for adaptive selection
#[derive(Debug, Clone)]
pub struct OptimizationStats {
    /// Total number of bits
    pub total_bits: usize,
    /// Number of set bits
    pub ones_count: usize,
    /// Bit density
    pub density: f64,
    /// Selected implementation name
    pub implementation: String,
    /// Space overhead percentage
    pub space_overhead_percent: f64,
    /// Estimated performance tier
    pub estimated_performance_tier: PerformanceTier,
}

impl fmt::Display for OptimizationStats {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "AdaptiveRankSelect Stats:\n\
             Total bits: {}\n\
             Ones count: {} ({:.2}% density)\n\
             Implementation: {}\n\
             Space overhead: {:.2}%\n\
             Performance tier: {:?}",
            self.total_bits,
            self.ones_count,
            self.density * 100.0,
            self.implementation,
            self.space_overhead_percent,
            self.estimated_performance_tier
        )
    }
}

// Implement RankSelectOps trait by delegating to selected implementation
impl RankSelectOps for AdaptiveRankSelect {
    #[inline]
    fn rank1(&self, pos: usize) -> usize {
        self.implementation.rank1(pos)
    }

    #[inline]
    fn rank0(&self, pos: usize) -> usize {
        self.implementation.rank0(pos)
    }

    #[inline]
    fn select1(&self, k: usize) -> Result<usize> {
        self.implementation.select1(k)
    }

    #[inline]
    fn select0(&self, k: usize) -> Result<usize> {
        self.implementation.select0(k)
    }

    fn len(&self) -> usize {
        self.implementation.len()
    }

    fn count_ones(&self) -> usize {
        self.implementation.count_ones()
    }

    fn get(&self, index: usize) -> Option<bool> {
        self.implementation.get(index)
    }

    fn space_overhead_percent(&self) -> f64 {
        self.implementation.space_overhead_percent()
    }
}

/// Multi-dimensional adaptive selection for related bit vectors.
///
/// Direct storage eliminates vtable indirection.
pub struct AdaptiveMultiDimensional {
    /// The selected implementation — always RankSelectInterleaved256
    implementation: RankSelectInterleaved256,
    /// Implementation name
    implementation_name: String,
    /// Number of dimensions
    dimensions: usize,
}

impl AdaptiveMultiDimensional {
    /// Create adaptive multi-dimensional rank/select for two related bit vectors
    pub fn new_dual(bit_vector1: BitVector, bit_vector2: BitVector) -> Result<Self> {
        if bit_vector1.len() != bit_vector2.len() {
            return Err(ZiporaError::invalid_data(
                "Bit vectors must have the same length".to_string()
            ));
        }

        // Since we only have RankSelectInterleaved256 available and it's the best performer,
        // we use it for the first bit vector. Multi-dimensional support can be enhanced later.
        let implementation = RankSelectInterleaved256::new(bit_vector1)?;

        Ok(Self {
            implementation,
            implementation_name: "RankSelectInterleaved256 (dual dimension)".to_string(),
            dimensions: 2,
        })
    }

    /// Get the name of the selected implementation
    pub fn implementation_name(&self) -> &str {
        &self.implementation_name
    }

    /// Get the number of dimensions
    pub fn dimensions(&self) -> usize {
        self.dimensions
    }
}

impl RankSelectOps for AdaptiveMultiDimensional {
    #[inline]
    fn rank1(&self, pos: usize) -> usize {
        self.implementation.rank1(pos)
    }

    #[inline]
    fn rank0(&self, pos: usize) -> usize {
        self.implementation.rank0(pos)
    }

    #[inline]
    fn select1(&self, k: usize) -> Result<usize> {
        self.implementation.select1(k)
    }

    #[inline]
    fn select0(&self, k: usize) -> Result<usize> {
        self.implementation.select0(k)
    }

    fn len(&self) -> usize {
        self.implementation.len()
    }

    fn count_ones(&self) -> usize {
        self.implementation.count_ones()
    }

    fn get(&self, index: usize) -> Option<bool> {
        self.implementation.get(index)
    }

    fn space_overhead_percent(&self) -> f64 {
        self.implementation.space_overhead_percent()
    }
}

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

    fn create_sparse_bitvector(size: usize, sparsity: usize) -> BitVector {
        let mut bv = BitVector::new();
        for i in 0..size {
            bv.push(i % sparsity == 0).unwrap();
        }
        bv
    }

    fn create_dense_bitvector(size: usize) -> BitVector {
        let mut bv = BitVector::new();
        for i in 0..size {
            bv.push(i % 10 != 0).unwrap(); // 90% density
        }
        bv
    }

    fn create_balanced_bitvector(size: usize) -> BitVector {
        let mut bv = BitVector::new();
        for i in 0..size {
            bv.push(i % 2 == 0).unwrap(); // 50% density
        }
        bv
    }

    #[test]
    fn test_sparse_data_selection() {
        // Very sparse data (1% density)
        let sparse_bv = create_sparse_bitvector(10000, 100);
        let adaptive = AdaptiveRankSelect::new(sparse_bv).unwrap();
        
        assert!(adaptive.implementation_name().contains("RankSelectInterleaved256"));
        assert!(adaptive.data_profile().density < 0.05);
        
        // Test basic operations
        let rank = adaptive.rank1(1000);
        assert_eq!(rank, 10); // Every 100th bit, so 10 bits set in first 1000
        
        let pos = adaptive.select1(5).unwrap();
        assert_eq!(pos, 500); // 5th set bit is at position 500
    }

    #[test]
    fn test_dense_data_selection() {
        // Dense data (90% density)
        let dense_bv = create_dense_bitvector(1000);
        let adaptive = AdaptiveRankSelect::new(dense_bv).unwrap();
        
        println!("Dense test - Selected: {}, Density: {:.3}, Size category: {:?}", 
                 adaptive.implementation_name(), 
                 adaptive.data_profile().density,
                 adaptive.data_profile().size_category);
        
        // Small dense dataset should select either cache-efficient or separated implementation  
        assert!(adaptive.implementation_name().contains("RankSelectInterleaved256"));
        assert!(adaptive.data_profile().density > 0.8);
        
        // Test basic operations
        let rank = adaptive.rank1(100);
        assert_eq!(rank, 90); // 90% of 100 bits
    }

    #[test]
    fn test_balanced_data_selection() {
        // Balanced data (50% density)
        let balanced_bv = create_balanced_bitvector(50000);
        let adaptive = AdaptiveRankSelect::new(balanced_bv).unwrap();
        
        // Should select based on size (medium dataset)
        assert!(adaptive.implementation_name().contains("RankSelect"));
        assert!((adaptive.data_profile().density - 0.5).abs() < 0.1);
        
        // Test basic operations
        let rank = adaptive.rank1(1000);
        assert_eq!(rank, 500); // 50% of 1000 bits
    }

    #[test]
    fn test_small_dataset_selection() {
        // Small dataset should prefer cache-efficient implementation
        let small_bv = create_balanced_bitvector(5000);
        let adaptive = AdaptiveRankSelect::new(small_bv).unwrap();
        
        assert!(adaptive.implementation_name().contains("RankSelectInterleaved256"));
        assert_eq!(adaptive.data_profile().size_category, SizeCategory::Small);
    }

    #[test]
    fn test_large_dataset_selection() {
        // Large dataset should prefer separated implementation
        let large_bv = create_balanced_bitvector(5_000_000);
        let adaptive = AdaptiveRankSelect::new(large_bv).unwrap();
        
        assert!(adaptive.implementation_name().contains("RankSelectInterleaved256"));
        assert_eq!(adaptive.data_profile().size_category, SizeCategory::Large);
    }

    #[test]
    fn test_custom_criteria() {
        let bv = create_sparse_bitvector(1000, 50); // 2% density
        
        // Custom criteria with higher sparse threshold
        let criteria = SelectionCriteria {
            sparse_threshold: 0.025, // 2.5% - should trigger sparse optimization for 2% data
            ..Default::default()
        };
        
        let adaptive = AdaptiveRankSelect::with_criteria(bv.clone(), criteria).unwrap();
        
        println!("Custom criteria test - Selected: {}, Density: {:.3}", 
                 adaptive.implementation_name(), 
                 adaptive.data_profile().density);
        
        assert!(adaptive.implementation_name().contains("RankSelectInterleaved256"));
    }

    #[test]
    fn test_edge_cases() {
        // All zeros
        let all_zeros = BitVector::with_size(1000, false).unwrap();
        let adaptive = AdaptiveRankSelect::new(all_zeros).unwrap();
        assert!(adaptive.implementation_name().contains("RankSelectInterleaved256"));
        assert_eq!(adaptive.count_ones(), 0);
        
        // All ones
        let all_ones = BitVector::with_size(1000, true).unwrap();
        let adaptive = AdaptiveRankSelect::new(all_ones).unwrap();
        assert!(adaptive.implementation_name().contains("RankSelectInterleaved256"));
        assert_eq!(adaptive.count_ones(), 1000);
    }

    #[test]
    fn test_optimization_stats() {
        let sparse_bv = create_sparse_bitvector(10000, 100);
        let adaptive = AdaptiveRankSelect::new(sparse_bv).unwrap();
        
        let stats = adaptive.optimization_stats();
        assert_eq!(stats.total_bits, 10000);
        assert_eq!(stats.ones_count, 100);
        assert!((stats.density - 0.01).abs() < 0.001);
        assert!(stats.implementation.contains("RankSelectInterleaved256"));
        assert_eq!(stats.estimated_performance_tier, PerformanceTier::HighPerformance);
        
        // Test display
        let display_str = format!("{}", stats);
        assert!(display_str.contains("AdaptiveRankSelect Stats"));
        assert!(display_str.contains("RankSelectInterleaved256"));
    }

    #[test]
    fn test_multi_dimensional_adaptive() {
        let bv1 = create_balanced_bitvector(1000);
        let bv2 = create_sparse_bitvector(1000, 50);
        
        let multi = AdaptiveMultiDimensional::new_dual(bv1, bv2).unwrap();
        assert!(multi.implementation_name().contains("RankSelectInterleaved256"));
        assert_eq!(multi.dimensions(), 2);
        
        // Test basic operations
        assert_eq!(multi.len(), 1000);
        let rank = multi.rank1(100);
        assert!(rank > 0);
    }

    #[test]
    fn test_multi_dimensional_length_mismatch() {
        let bv1 = create_balanced_bitvector(1000);
        let bv2 = create_balanced_bitvector(500); // Different length
        
        let result = AdaptiveMultiDimensional::new_dual(bv1, bv2);
        assert!(result.is_err());
    }

    #[test]
    fn test_data_profile_analysis() {
        let bv = create_sparse_bitvector(50000, 200); // 0.5% density
        let criteria = SelectionCriteria::default();
        let profile = AdaptiveRankSelect::analyze_data(&bv, &criteria);
        
        assert_eq!(profile.total_bits, 50000);
        assert_eq!(profile.ones_count, 250);
        assert!((profile.density - 0.005).abs() < 0.001);
        assert_eq!(profile.size_category, SizeCategory::Medium);
        // Access pattern may be inferred based on data characteristics, so check it's valid
        assert!(matches!(profile.access_pattern, AccessPattern::Mixed | AccessPattern::Sequential | AccessPattern::Random));
    }

    #[test]
    fn test_performance_consistency() {
        // Ensure adaptive selection provides consistent results
        let bv = create_balanced_bitvector(10000);
        let adaptive = AdaptiveRankSelect::new(bv.clone()).unwrap();
        let reference = RankSelectInterleaved256::new(bv).unwrap();
        
        // Test multiple positions
        for pos in [0, 1000, 5000, 9999] {
            assert_eq!(adaptive.rank1(pos), reference.rank1(pos));
            assert_eq!(adaptive.rank0(pos), reference.rank0(pos));
        }
        
        // Test select operations
        let ones_count = adaptive.count_ones();
        for k in [0, ones_count / 4, ones_count / 2, ones_count - 1] {
            if k < ones_count {
                assert_eq!(adaptive.select1(k).unwrap(), reference.select1(k).unwrap());
            }
        }
    }
    
    #[test]
    fn test_pattern_analysis() {
        // Test run length analysis
        let mut regular_bv = BitVector::new();
        for i in 0..1000 {
            regular_bv.push(i % 8 < 4).unwrap(); // Regular pattern: 4 ones, 4 zeros
        }
        
        let criteria = SelectionCriteria::default();
        let profile = AdaptiveRankSelect::analyze_data(&regular_bv, &criteria);
        
        // Should detect regular pattern with low complexity
        assert!(profile.pattern_complexity < 0.5);
        assert!(profile.run_length_stats.avg_ones_run > 3.0);
        assert!(profile.run_length_stats.avg_zeros_run > 3.0);
        assert_eq!(profile.access_pattern, AccessPattern::Sequential);
        
        // Test with random pattern
        let mut random_bv = BitVector::new();
        for i in 0..1000 {
            random_bv.push((i * 31 + 17) % 3 == 0).unwrap(); // Pseudo-random pattern
        }
        
        let random_profile = AdaptiveRankSelect::analyze_data(&random_bv, &criteria);
        assert!(random_profile.pattern_complexity > profile.pattern_complexity);
    }
    
    #[test]
    fn test_adaptive_thresholds() {
        // Test adaptive threshold tuning
        let mut clustered_bv = BitVector::new();
        for i in 0..1000 {
            // Create clustered pattern: groups of 50 ones followed by groups of 50 zeros
            clustered_bv.push((i / 50) % 2 == 0).unwrap();
        }
        
        let criteria = SelectionCriteria {
            enable_adaptive_thresholds: true,
            ..Default::default()
        };
        
        let profile = AdaptiveRankSelect::analyze_data(&clustered_bv, &criteria);
        let (adaptive_sparse, adaptive_dense) = AdaptiveRankSelect::tune_thresholds(&profile, &criteria);
        
        // Adaptive thresholds should be different from defaults for clustered data
        assert!(profile.clustering_coefficient > 0.5);
        
        // Test hardware-aware threshold adjustment
        if profile.hardware_tier >= 3 {
            // With BMI2, sparse threshold should be higher (can handle denser sparse data)
            assert!(adaptive_sparse >= criteria.sparse_threshold);
        }
        
        println!("Clustering: {:.3}, Original sparse: {:.3}, Adaptive sparse: {:.3}", 
                 profile.clustering_coefficient, criteria.sparse_threshold, adaptive_sparse);
    }
    
    #[test]
    fn test_entropy_calculation() {
        // Test entropy calculation
        let all_ones = BitVector::with_size(1000, true).unwrap();
        let criteria = SelectionCriteria::default();
        let profile_ones = AdaptiveRankSelect::analyze_data(&all_ones, &criteria);
        assert_eq!(profile_ones.entropy, 0.0); // No entropy for uniform data
        
        let balanced = create_balanced_bitvector(1000);
        let profile_balanced = AdaptiveRankSelect::analyze_data(&balanced, &criteria);
        assert!(profile_balanced.entropy > 0.9); // High entropy for balanced data
        assert!(profile_balanced.entropy <= 1.0); // Max entropy is 1.0 for binary
    }
    
    #[test]
    fn test_hardware_aware_selection() {
        let bv = create_balanced_bitvector(50000);
        
        // Test with hardware acceleration preference
        let hardware_criteria = SelectionCriteria {
            min_hardware_tier: 3, // Require BMI2
            enable_adaptive_thresholds: true,
            ..Default::default()
        };
        
        let profile = AdaptiveRankSelect::analyze_data(&bv, &hardware_criteria);
        
        // Should detect available hardware
        assert!(profile.hardware_tier <= 4); // Max tier is 4 (BMI2+AVX2)
        
        println!("Detected hardware tier: {}", profile.hardware_tier);
        
        // Hardware-aware adaptive should work
        let adaptive = AdaptiveRankSelect::with_criteria(bv, hardware_criteria);
        if adaptive.is_ok() {
            let adaptive = adaptive.unwrap();
            println!("Hardware-aware selection: {}", adaptive.implementation_name());
            
            // Test basic operations
            assert_eq!(adaptive.rank1(1000), 500); // 50% density
        }
    }
    
    #[test]
    fn test_complex_pattern_detection() {
        // Create a complex alternating pattern
        let mut complex_bv = BitVector::new();
        for i in 0..10000 {
            // Complex pattern: alternating with some regularity
            let bit = match i % 7 {
                0 | 1 | 2 => true,
                3 | 4 => false,
                5 => (i / 7) % 3 == 0,
                _ => false,
            };
            complex_bv.push(bit).unwrap();
        }
        
        let criteria = SelectionCriteria {
            enable_adaptive_thresholds: true,
            pattern_complexity_weight: 0.5,
            clustering_weight: 0.3,
            ..Default::default()
        };
        
        let profile = AdaptiveRankSelect::analyze_data(&complex_bv, &criteria);
        
        // Should detect moderate complexity
        assert!(profile.pattern_complexity > 0.2);
        assert!(profile.pattern_complexity < 0.8);
        
        // Should have reasonable alternations
        assert!(profile.run_length_stats.alternations > 1000);
        
        println!("Complex pattern - Complexity: {:.3}, Alternations: {}, Avg run lengths: {:.1}/{:.1}",
                 profile.pattern_complexity, 
                 profile.run_length_stats.alternations,
                 profile.run_length_stats.avg_ones_run,
                 profile.run_length_stats.avg_zeros_run);
        
        // Adaptive selection should handle this appropriately
        let adaptive = AdaptiveRankSelect::with_criteria(complex_bv.clone(), criteria).unwrap();
        
        // Test correctness
        let interleaved = RankSelectInterleaved256::new(complex_bv).unwrap();
        for pos in [0, 1000, 5000, 9999] {
            assert_eq!(adaptive.rank1(pos), interleaved.rank1(pos));
        }
    }
}