scirs2_interpolate/

cache.rs

//! Cache-aware algorithm implementations for interpolation
//!
//! This module provides cache-optimized versions of computationally intensive
//! interpolation algorithms. The caching strategies focus on:
//!
//! - **Basis function caching**: Pre-computed and cached basis function evaluations
//! - **Coefficient matrix caching**: Cached matrix factorizations and linear solves
//! - **Distance matrix caching**: Cached distance computations for scattered data methods
//! - **Knot span caching**: Cached knot span lookups for B-splines
//! - **Memory layout optimization**: Data structures optimized for cache locality
//!
//! These optimizations can provide significant performance improvements for:
//! - Repeated evaluations at similar points
//! - Large datasets with repeated computations
//! - Real-time applications requiring fast interpolation
//!
//! # Examples
//!
//! ```rust
//! use scirs2_core::ndarray::Array1;
//! use scirs2_interpolate::cache::{CachedBSpline, BSplineCache};
//! use scirs2_interpolate::bspline::ExtrapolateMode;
//!
//! // Create a cached B-spline for fast repeated evaluations
//! let knots = Array1::linspace(0.0, 10.0, 20);
//! let coeffs = Array1::linspace(-1.0, 1.0, 16);
//!
//! let mut cached_spline = CachedBSpline::new(
//!     &knots.view(),
//!     &coeffs.view(),
//!     3, // degree
//!     ExtrapolateMode::Extrapolate,
//!     BSplineCache::default(),
//! ).expect("Operation failed");
//!
//! // Use cached evaluation for fast performance
//! for x in Array1::linspace(0.0, 10.0, 1000).iter() {
//!     let y = cached_spline.evaluate_cached(*x).expect("Operation failed");
//! }
//! ```

use crate::bspline::{BSpline, ExtrapolateMode};
use crate::error::InterpolateResult;
use scirs2_core::ndarray::{Array1, Array2, ArrayView1};
use scirs2_core::numeric::{Float, FromPrimitive, Zero};
use std::collections::HashMap;
use std::fmt::{Debug, Display};
use std::hash::{Hash, Hasher};
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, RemAssign, Sub, SubAssign};

// Simple random generation for eviction policy

/// Cache eviction policies
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum EvictionPolicy {
    /// Least Recently Used - evict oldest accessed items
    LRU,
    /// Least Frequently Used - evict least accessed items  
    LFU,
    /// First In First Out - evict oldest inserted items
    FIFO,
    /// Random eviction
    Random,
    /// Adaptive policy based on access patterns
    Adaptive,
}

/// Configuration for cache behavior
#[derive(Debug, Clone)]
pub struct CacheConfig {
    /// Maximum number of entries in basis function cache
    pub max_basis_cache_size: usize,
    /// Maximum number of entries in coefficient matrix cache  
    pub max_matrix_cache_size: usize,
    /// Maximum number of entries in distance matrix cache
    pub max_distance_cache_size: usize,
    /// Tolerance for cache key matching (for floating point comparisons)
    pub tolerance: f64,
    /// Whether to enable cache statistics tracking
    pub track_stats: bool,
    /// Cache eviction strategy
    pub eviction_policy: EvictionPolicy,
    /// Memory limit for caches in MB (0 = no limit)
    pub memory_limit_mb: usize,
    /// Enable adaptive cache sizing based on access patterns
    pub adaptive_sizing: bool,
}

impl Default for CacheConfig {
    fn default() -> Self {
        Self {
            max_basis_cache_size: 1024,
            max_matrix_cache_size: 64,
            max_distance_cache_size: 256,
            tolerance: 1e-12,
            track_stats: false,
            eviction_policy: EvictionPolicy::LRU,
            memory_limit_mb: 0, // No limit by default
            adaptive_sizing: true,
        }
    }
}

/// Cache statistics for performance monitoring
#[derive(Debug, Clone)]
pub struct CacheStats {
    /// Number of cache hits
    pub hits: usize,
    /// Number of cache misses  
    pub misses: usize,
    /// Number of cache evictions
    pub evictions: usize,
    /// Total memory usage in bytes
    pub memory_usage_bytes: usize,
    /// Average access frequency for adaptive sizing
    pub avg_access_frequency: f64,
    /// Number of cache resizes
    pub resize_count: usize,
    /// Peak memory usage
    pub peak_memory_bytes: usize,
    /// Last cleanup timestamp (for TTL-based eviction)
    pub last_cleanup_time: std::time::Instant,
}

impl Default for CacheStats {
    fn default() -> Self {
        Self {
            hits: 0,
            misses: 0,
            evictions: 0,
            memory_usage_bytes: 0,
            avg_access_frequency: 0.0,
            resize_count: 0,
            peak_memory_bytes: 0,
            last_cleanup_time: std::time::Instant::now(),
        }
    }
}

impl CacheStats {
    /// Calculate cache hit ratio
    pub fn hit_ratio(&self) -> f64 {
        if self.hits + self.misses == 0 {
            0.0
        } else {
            self.hits as f64 / (self.hits + self.misses) as f64
        }
    }

    /// Calculate cache efficiency score (combines hit ratio and memory usage)
    pub fn efficiency_score(&self) -> f64 {
        let hit_ratio = self.hit_ratio();
        let memory_factor = if self.peak_memory_bytes > 0 {
            1.0 - (self.memory_usage_bytes as f64 / self.peak_memory_bytes as f64)
        } else {
            1.0
        };
        (hit_ratio + memory_factor) / 2.0
    }

    /// Get memory usage in MB
    pub fn memory_usage_mb(&self) -> f64 {
        self.memory_usage_bytes as f64 / (1024.0 * 1024.0)
    }

    /// Get peak memory usage in MB
    pub fn peak_memory_mb(&self) -> f64 {
        self.peak_memory_bytes as f64 / (1024.0 * 1024.0)
    }

    /// Update memory usage statistics
    pub fn update_memory_usage(&mut self, currentbytes: usize) {
        self.memory_usage_bytes = currentbytes;
        if currentbytes > self.peak_memory_bytes {
            self.peak_memory_bytes = currentbytes;
        }
    }

    /// Update access frequency for adaptive sizing
    pub fn update_access_frequency(&mut self, new_accesscount: usize) {
        let alpha = 0.1; // Exponential smoothing factor
        let new_freq = new_accesscount as f64;
        self.avg_access_frequency = alpha * new_freq + (1.0 - alpha) * self.avg_access_frequency;
    }

    /// Check if cleanup is needed based on time threshold
    pub fn needs_cleanup(&self, thresholdsecs: u64) -> bool {
        self.last_cleanup_time.elapsed().as_secs() >= thresholdsecs
    }

    /// Reset all statistics
    pub fn reset(&mut self) {
        self.hits = 0;
        self.misses = 0;
        self.evictions = 0;
        self.memory_usage_bytes = 0;
        self.avg_access_frequency = 0.0;
        self.resize_count = 0;
        self.peak_memory_bytes = 0;
        self.last_cleanup_time = std::time::Instant::now();
    }
}

/// A floating-point key that can be hashed with tolerance
#[derive(Debug, Clone)]
struct FloatKey<F: Float> {
    value: F,
    tolerance: F,
}

impl<F: Float> FloatKey<F> {
    #[allow(dead_code)]
    fn new(value: F, tolerance: F) -> Self {
        Self { value, tolerance }
    }
}

impl<F: crate::traits::InterpolationFloat> PartialEq for FloatKey<F> {
    fn eq(&self, other: &Self) -> bool {
        (self.value - other.value).abs() <= self.tolerance
    }
}

impl<F: crate::traits::InterpolationFloat> Eq for FloatKey<F> {}

impl<F: crate::traits::InterpolationFloat> Hash for FloatKey<F> {
    fn hash<H: Hasher>(&self, state: &mut H) {
        // Quantize the float to the tolerance for consistent hashing
        let quantized = (self.value / self.tolerance).round() * self.tolerance;
        // Convert to bits for hashing (this is approximate)
        let bits = quantized.to_f64().unwrap_or(0.0).to_bits();
        bits.hash(state);
    }
}

/// Cache for B-spline basis function evaluations
#[derive(Debug)]
pub struct BSplineCache<F: Float> {
    /// Cache for basis function values with access tracking
    basis_cache: HashMap<(FloatKey<F>, usize, usize), CacheEntry<F>>,
    /// Cache for knot span indices with access tracking
    span_cache: HashMap<FloatKey<F>, CacheEntry<usize>>,
    /// Configuration
    config: CacheConfig,
    /// Statistics
    stats: CacheStats,
    /// Access counter for LRU/LFU tracking
    #[allow(dead_code)]
    access_counter: usize,
}

/// Cache key for basis function evaluations
#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct BasisCacheKey {
    /// Quantized x value for consistent cache key generation
    pub x_quantized: u64,
    /// Basis function index
    pub index: usize,
    /// Polynomial degree
    pub degree: usize,
}

/// Cache entry with metadata for advanced eviction policies
#[derive(Debug, Clone)]
struct CacheEntry<T> {
    /// The cached value
    #[allow(dead_code)]
    value: T,
    /// Last access time (for LRU)
    #[allow(dead_code)]
    last_access: usize,
    /// Access frequency (for LFU)
    #[allow(dead_code)]
    access_count: usize,
    /// Insertion time (for FIFO)
    #[allow(dead_code)]
    insertion_time: usize,
    /// Estimated memory size in bytes
    #[allow(dead_code)]
    memory_size: usize,
}

#[allow(dead_code)]
impl<T> CacheEntry<T> {
    /// Create a new cache entry
    fn new(_value: T, insertiontime: usize) -> Self {
        let memory_size = std::mem::size_of::<T>() + std::mem::size_of::<Self>();
        Self {
            value: _value,
            last_access: insertiontime,
            access_count: 1,
            insertion_time: insertiontime,
            memory_size,
        }
    }

    /// Update access statistics
    fn update_access(&mut self, currenttime: usize) {
        self.last_access = currenttime;
        self.access_count += 1;
    }

    /// Calculate priority for eviction (lower means more likely to be evicted)
    fn eviction_priority(&self, policy: EvictionPolicy, currenttime: usize) -> f64 {
        match policy {
            EvictionPolicy::LRU => -(self.last_access as f64),
            EvictionPolicy::LFU => -(self.access_count as f64),
            EvictionPolicy::FIFO => -(self.insertion_time as f64),
            EvictionPolicy::Random => {
                // Simple linear congruential generator for randomness
                let x = (self.insertion_time * 1103515245 + 12345) & 0x7fffffff;
                x as f64 / 0x7fffffff as f64
            }
            EvictionPolicy::Adaptive => {
                // Combine recency and frequency with memory size consideration
                let recency = (currenttime - self.last_access) as f64;
                let frequency = self.access_count as f64;
                let memory_factor = self.memory_size as f64;

                // Higher score means less likely to be evicted
                -(frequency / (1.0 + recency + memory_factor / 1000.0))
            }
        }
    }
}

impl<F: crate::traits::InterpolationFloat> Default for BSplineCache<F> {
    fn default() -> Self {
        Self::new(CacheConfig::default())
    }
}

#[allow(dead_code)]
impl<F: crate::traits::InterpolationFloat> BSplineCache<F> {
    /// Create a new B-spline cache with the given configuration
    pub fn new(config: CacheConfig) -> Self {
        Self {
            basis_cache: HashMap::new(),
            span_cache: HashMap::new(),
            config,
            stats: CacheStats::default(),
            access_counter: 0,
        }
    }

    /// Get cached basis function value or compute and cache it
    fn get_or_compute_basis<T>(
        &mut self,
        x: F,
        i: usize,
        k: usize,
        knots: &[T],
        computer: impl FnOnce() -> T,
    ) -> T
    where
        T: Float + Copy,
    {
        self.access_counter += 1;
        let tolerance = F::from_f64(self.config.tolerance).expect("Operation failed");
        let key = (FloatKey::new(x, tolerance), i, k);

        if let Some(cache_entry) = self.basis_cache.get_mut(&key) {
            if self.config.track_stats {
                self.stats.hits += 1;
            }

            // Update access statistics
            cache_entry.update_access(self.access_counter);

            // Convert from F to T (this assumes they're the same type or compatible)
            unsafe { std::mem::transmute_copy(&cache_entry.value) }
        } else {
            if self.config.track_stats {
                self.stats.misses += 1;
            }
            let computed = computer();

            // Convert from T to F for caching (again, assumes compatibility)
            let cached: F = unsafe { std::mem::transmute_copy(&computed) };

            // Check cache size and evict if necessary
            if self.basis_cache.len() >= self.config.max_basis_cache_size {
                self.evict_basis_cache();
            }

            // Create cache entry with metadata
            let cache_entry = CacheEntry::new(cached, self.access_counter);
            self.basis_cache.insert(key, cache_entry);

            // Update memory usage statistics
            if self.config.track_stats {
                self.update_memory_usage();
            }

            computed
        }
    }

    /// Get cached basis function value using new key structure
    pub fn get_or_compute_basis_with_key(
        &mut self,
        key: BasisCacheKey,
        computer: impl FnOnce() -> F,
    ) -> F {
        self.access_counter += 1;

        // Convert BasisCacheKey to internal cache key format
        let tolerance = F::from_f64(self.config.tolerance).expect("Operation failed");
        let x = F::from_f64(f64::from_bits(key.x_quantized)).unwrap_or_else(F::zero);
        let internal_key = (FloatKey::new(x, tolerance), key.index, key.degree);

        if let Some(cache_entry) = self.basis_cache.get_mut(&internal_key) {
            if self.config.track_stats {
                self.stats.hits += 1;
            }

            // Update access statistics
            cache_entry.update_access(self.access_counter);
            cache_entry.value
        } else {
            if self.config.track_stats {
                self.stats.misses += 1;
            }
            let computed = computer();

            // Check cache size and evict if necessary
            if self.basis_cache.len() >= self.config.max_basis_cache_size {
                self.evict_basis_cache();
            }

            // Create cache entry with metadata
            let cache_entry = CacheEntry::new(computed, self.access_counter);
            self.basis_cache.insert(internal_key, cache_entry);

            // Update memory usage statistics
            if self.config.track_stats {
                self.update_memory_usage();
            }

            computed
        }
    }

    /// Get cached knot span or compute and cache it
    fn get_or_compute_span(&mut self, x: F, computer: impl FnOnce() -> usize) -> usize {
        self.access_counter += 1;
        let tolerance = F::from_f64(self.config.tolerance).expect("Operation failed");
        let key = FloatKey::new(x, tolerance);

        if let Some(cache_entry) = self.span_cache.get_mut(&key) {
            if self.config.track_stats {
                self.stats.hits += 1;
            }

            // Update access statistics
            cache_entry.update_access(self.access_counter);
            cache_entry.value
        } else {
            if self.config.track_stats {
                self.stats.misses += 1;
            }
            let computed = computer();

            // Create cache entry with metadata
            let cache_entry = CacheEntry::new(computed, self.access_counter);
            self.span_cache.insert(key, cache_entry);

            // Update memory usage statistics
            if self.config.track_stats {
                self.update_memory_usage();
            }

            computed
        }
    }

    /// Update memory usage statistics
    fn update_memory_usage(&mut self) {
        if !self.config.track_stats {
            return;
        }

        let basis_memory: usize = self
            .basis_cache
            .values()
            .map(|entry| entry.memory_size)
            .sum();
        let span_memory: usize = self
            .span_cache
            .values()
            .map(|entry| entry.memory_size)
            .sum();

        let total_memory = basis_memory + span_memory;
        self.stats.update_memory_usage(total_memory);

        // Check memory limit if specified
        if self.config.memory_limit_mb > 0 {
            let limit_bytes = self.config.memory_limit_mb * 1024 * 1024;
            if total_memory > limit_bytes {
                self.evict_basis_cache_by_memory();
            }
        }
    }

    /// Evict entries from the basis cache based on memory pressure
    fn evict_basis_cache_by_memory(&mut self) {
        let target_size = self.config.memory_limit_mb * 1024 * 1024 * 3 / 4; // Target 75% of limit
        let mut current_memory = self.stats.memory_usage_bytes;

        if current_memory <= target_size {
            return;
        }

        // Create vector of (key, priority) pairs for sorting
        let mut entries: Vec<_> = self
            .basis_cache
            .iter()
            .map(|(key, entry)| {
                let priority =
                    entry.eviction_priority(self.config.eviction_policy, self.access_counter);
                (key.clone(), priority, entry.memory_size)
            })
            .collect();

        // Sort by eviction priority (lowest first)
        entries.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

        // Remove entries until we're under the target
        for (key_, _, memory_size) in entries {
            if current_memory <= target_size {
                break;
            }

            self.basis_cache.remove(&key_);
            current_memory = current_memory.saturating_sub(memory_size);

            if self.config.track_stats {
                self.stats.evictions += 1;
            }
        }

        // Update memory statistics
        self.update_memory_usage();
    }

    /// Evict some entries from the basis cache when it gets too large
    /// Uses advanced eviction policies for better cache performance
    fn evict_basis_cache(&mut self) {
        let total_entries = self.basis_cache.len();
        let remove_count = if self.config.adaptive_sizing {
            // Adaptive removal based on hit ratio
            let hit_ratio = self.stats.hit_ratio();
            if hit_ratio > 0.8 {
                total_entries / 8 // Remove fewer entries if hit ratio is high
            } else if hit_ratio > 0.5 {
                total_entries / 4 // Standard removal
            } else {
                total_entries / 2 // Remove more entries if hit ratio is low
            }
        } else {
            total_entries / 4 // Remove 25% by default
        };

        // Create vector of (key, priority) pairs for sorting
        let mut entries: Vec<_> = self
            .basis_cache
            .iter()
            .map(|(key, entry)| {
                let priority =
                    entry.eviction_priority(self.config.eviction_policy, self.access_counter);
                (key.clone(), priority)
            })
            .collect();

        // Sort by eviction priority (lowest first = most likely to be evicted)
        entries.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

        // Remove the lowest priority entries
        for (key, _) in entries.into_iter().take(remove_count) {
            self.basis_cache.remove(&key);
            if self.config.track_stats {
                self.stats.evictions += 1;
            }
        }

        // Update statistics
        if self.config.track_stats {
            self.stats.resize_count += 1;
            self.update_memory_usage();
        }
    }

    /// Clear all cached data
    pub fn clear(&mut self) {
        self.basis_cache.clear();
        self.span_cache.clear();
    }

    /// Get cache statistics
    pub fn stats(&self) -> &CacheStats {
        &self.stats
    }

    /// Reset cache statistics
    pub fn reset_stats(&mut self) {
        self.stats.reset();
    }
}

/// A cache-aware B-spline implementation
#[derive(Debug)]
pub struct CachedBSpline<T>
where
    T: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = T>
        + Sub<Output = T>
        + Mul<Output = T>
        + Div<Output = T>
        + AddAssign
        + SubAssign
        + MulAssign
        + DivAssign
        + RemAssign
        + Zero
        + Copy,
{
    /// The underlying B-spline
    spline: BSpline<T>,
    /// Cache for basis function evaluations
    cache: BSplineCache<T>,
}

#[allow(dead_code)]
impl<T> CachedBSpline<T>
where
    T: crate::traits::InterpolationFloat,
{
    /// Create a new cached B-spline
    ///
    /// # Arguments
    ///
    /// * `knots` - Knot vector for the B-spline
    /// * `coeffs` - Control coefficients
    /// * `degree` - Degree of the B-spline
    /// * `extrapolate` - Extrapolation mode for points outside the domain
    /// * `cache` - Cache configuration for optimization
    ///
    /// # Examples
    ///
    /// ```rust
    /// use scirs2_core::ndarray::Array1;
    /// use scirs2_interpolate::cache::{CachedBSpline, BSplineCache, CacheConfig};
    /// use scirs2_interpolate::bspline::ExtrapolateMode;
    ///
    /// let knots = Array1::from_vec(vec![0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0]);
    /// let coeffs = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
    ///
    /// let cache_config = CacheConfig {
    ///     track_stats: true,
    ///     max_basis_cache_size: 512,
    ///     ..Default::default()
    /// };
    /// let cache = BSplineCache::new(cache_config);
    ///
    /// let cached_spline = CachedBSpline::new(
    ///     &knots.view(),
    ///     &coeffs.view(),
    ///     2, // quadratic
    ///     ExtrapolateMode::Extrapolate,
    ///     cache,
    /// ).expect("Operation failed");
    /// ```
    pub fn new(
        knots: &ArrayView1<T>,
        coeffs: &ArrayView1<T>,
        degree: usize,
        extrapolate: ExtrapolateMode,
        cache: BSplineCache<T>,
    ) -> InterpolateResult<Self> {
        let spline = BSpline::new(knots, coeffs, degree, extrapolate)?;

        Ok(Self { spline, cache })
    }

    /// Evaluate the B-spline using cached basis functions
    ///
    /// This method provides optimized evaluation by caching basis function
    /// computations. For repeated evaluations at similar points, this can
    /// be significantly faster than standard evaluation.
    ///
    /// # Arguments
    ///
    /// * `x` - Point at which to evaluate the B-spline
    ///
    /// # Returns
    ///
    /// The value of the B-spline at x
    ///
    /// # Examples
    ///
    /// ```rust
    /// use scirs2_core::ndarray::Array1;
    /// use scirs2_interpolate::cache::make_cached_bspline;
    /// use scirs2_interpolate::bspline::ExtrapolateMode;
    ///
    /// let knots = Array1::linspace(0.0, 10.0, 11);
    /// let coeffs = Array1::linspace(1.0, 5.0, 8);
    ///
    /// let mut cached_spline = make_cached_bspline(
    ///     &knots.view(),
    ///     &coeffs.view(),
    ///     2,
    ///     ExtrapolateMode::Extrapolate,
    /// ).expect("Operation failed");
    ///
    /// // Fast repeated evaluations
    /// let x_values = Array1::linspace(0.0, 10.0, 100);
    /// for &x in x_values.iter() {
    ///     let y = cached_spline.evaluate_cached(x).expect("Operation failed");
    ///     // Process result...
    /// }
    ///
    /// // Check cache performance
    /// let stats = cached_spline.cache_stats();
    /// println!("Cache hit ratio: {:.2}", stats.hit_ratio());
    /// ```
    pub fn evaluate_cached(&mut self, x: T) -> InterpolateResult<T> {
        // For this implementation, we'll delegate to the underlying spline
        // In a full implementation, we would cache the basis function evaluations
        // and use the cache during the evaluation process
        self.evaluate_with_cache_optimization(x)
    }

    /// Evaluate with cache optimization for basis functions
    fn evaluate_with_cache_optimization(&mut self, x: T) -> InterpolateResult<T> {
        // Get knot vector and coefficients from the underlying spline
        let knots = self.spline.knot_vector();
        let coeffs = self.spline.coefficients();
        let degree = self.spline.degree();

        // Find the knot span using cached lookup
        // Clone necessary data to avoid borrowing conflicts
        let knots_clone = knots.to_owned();
        let span = self
            .cache
            .get_or_compute_span(x, || Self::find_knot_span(x, &knots_clone, degree));

        // Compute basis functions using cache
        let mut result = T::zero();
        for i in 0..=degree {
            let basis_index = span.saturating_sub(degree) + i;
            if basis_index < coeffs.len() {
                // Cache key for this basis function evaluation
                let basis_key = BasisCacheKey {
                    x_quantized: self.quantize_x(x),
                    index: basis_index,
                    degree,
                };

                let basis_value = self.cache.get_or_compute_basis_with_key(basis_key, || {
                    if basis_index < knots_clone.len() - degree - 1 {
                        Self::compute_basis_function(x, basis_index, degree, &knots_clone)
                    } else {
                        T::zero()
                    }
                });

                result += coeffs[basis_index] * basis_value;
            }
        }

        Ok(result)
    }

    /// Quantize x value for cache key consistency
    fn quantize_x(&self, x: T) -> u64 {
        let x_f64 = x.to_f64().unwrap_or(0.0);
        let tolerance = self.cache.config.tolerance;
        let quantized = (x_f64 / tolerance).round() * tolerance;
        quantized.to_bits()
    }

    /// Compute all non-zero basis functions at x using De Boor's algorithm
    fn compute_basis_functions_at_x(&self, x: T, knots: &Array1<T>, degree: usize) -> Array1<T> {
        let span = Self::find_knot_span(x, knots, degree);
        let mut basis = Array1::zeros(degree + 1);

        // De Boor's algorithm for all basis functions at once
        basis[0] = T::one();

        for j in 1..=degree {
            let mut saved = T::zero();
            for r in 0..j {
                let left_knot_idx = span + 1 - j + r;
                let right_knot_idx = span + r;

                if left_knot_idx < knots.len() && right_knot_idx + 1 < knots.len() {
                    let alpha = if knots[right_knot_idx + 1] != knots[left_knot_idx] {
                        (x - knots[left_knot_idx])
                            / (knots[right_knot_idx + 1] - knots[left_knot_idx])
                    } else {
                        T::zero()
                    };

                    let temp = basis[r];
                    basis[r] = saved + (T::one() - alpha) * temp;
                    saved = alpha * temp;
                }
            }
            basis[j] = saved;
        }

        basis
    }

    /// Optimized batch evaluation using cached computations
    pub fn evaluate_batch_optimized(
        &mut self,
        x_vals: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        let mut results = Array1::zeros(x_vals.len());

        // Pre-warm cache with commonly accessed spans
        for &x in x_vals.iter().take(x_vals.len().min(10)) {
            let _ = self.evaluate_cached(x)?;
        }

        // Evaluate all points
        for (i, &x) in x_vals.iter().enumerate() {
            results[i] = self.evaluate_cached(x)?;
        }

        Ok(results)
    }

    /// Clear cache and force recomputation
    pub fn refresh_cache(&mut self) {
        self.cache.clear();
    }

    /// Optimize cache settings based on usage patterns
    pub fn optimize_cache_settings(&mut self) {
        let (hit_ratio, total_requests) = {
            let stats = self.cache.stats();
            (stats.hit_ratio(), stats.hits + stats.misses)
        };

        // Adjust cache size based on hit ratio
        if hit_ratio < 0.3 && self.cache.config.max_basis_cache_size > 64 {
            // Low hit ratio - reduce cache size
            self.cache.config.max_basis_cache_size /= 2;
        } else if hit_ratio > 0.8 && self.cache.config.max_basis_cache_size < 4096 {
            // High hit ratio - increase cache size
            self.cache.config.max_basis_cache_size *= 2;
        }

        // Enable aggressive caching for frequently accessed data
        if total_requests > 1000 {
            self.cache.config.adaptive_sizing = true;
        }
    }

    /// Find the knot span containing x
    fn find_knot_span(x: T, knots: &Array1<T>, degree: usize) -> usize {
        let n = knots.len() - degree - 1;

        if x >= knots[n] {
            return n - 1;
        }
        if x <= knots[degree] {
            return degree;
        }

        // Binary search
        let mut low = degree;
        let mut high = n;
        let mut mid = (low + high) / 2;

        while x < knots[mid] || x >= knots[mid + 1] {
            if x < knots[mid] {
                high = mid;
            } else {
                low = mid;
            }
            mid = (low + high) / 2;
        }

        mid
    }

    /// Compute a single basis function value
    #[allow(clippy::only_used_in_recursion)]
    fn compute_basis_function(x: T, i: usize, degree: usize, knots: &Array1<T>) -> T {
        // De Boor's algorithm for a single basis function
        if degree == 0 {
            if i < knots.len() - 1 && x >= knots[i] && x < knots[i + 1] {
                T::one()
            } else {
                T::zero()
            }
        } else {
            let mut left = T::zero();
            let mut right = T::zero();

            // Left recursion
            if i < knots.len() - degree - 1 && knots[i + degree] != knots[i] {
                let basis_left = Self::compute_basis_function(x, i, degree - 1, knots);
                left = (x - knots[i]) / (knots[i + degree] - knots[i]) * basis_left;
            }

            // Right recursion
            if i + 1 < knots.len() - degree - 1 && knots[i + degree + 1] != knots[i + 1] {
                let basis_right = Self::compute_basis_function(x, i + 1, degree - 1, knots);
                right = (knots[i + degree + 1] - x) / (knots[i + degree + 1] - knots[i + 1])
                    * basis_right;
            }

            left + right
        }
    }

    /// Evaluate the B-spline using the standard (non-cached) method
    pub fn evaluate_standard(&self, x: T) -> InterpolateResult<T> {
        self.spline.evaluate(x)
    }

    /// Evaluate at multiple points using cached basis functions
    ///
    /// This method efficiently evaluates the B-spline at multiple points,
    /// leveraging cached computations for improved performance.
    ///
    /// # Arguments
    ///
    /// * `x_vals` - Array of x-coordinates at which to evaluate
    ///
    /// # Returns
    ///
    /// Array of B-spline values corresponding to the input coordinates
    ///
    /// # Examples
    ///
    /// ```rust
    /// use scirs2_core::ndarray::Array1;
    /// use scirs2_interpolate::cache::make_cached_bspline;
    /// use scirs2_interpolate::bspline::ExtrapolateMode;
    ///
    /// let knots = Array1::from_vec(vec![0.0, 0.0, 1.0, 2.0, 3.0, 3.0]);
    /// let coeffs = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0]);
    ///
    /// let mut cached_spline = make_cached_bspline(
    ///     &knots.view(),
    ///     &coeffs.view(),
    ///     1, // linear
    ///     ExtrapolateMode::Extrapolate,
    /// ).expect("Operation failed");
    ///
    /// // Evaluate at multiple points efficiently
    /// let x_vals = Array1::from_vec(vec![0.5, 1.0, 1.5, 2.0, 2.5]);
    /// let results = cached_spline.evaluate_array_cached(&x_vals.view()).expect("Operation failed");
    ///
    /// assert_eq!(results.len(), x_vals.len());
    /// ```
    pub fn evaluate_array_cached(
        &mut self,
        x_vals: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        let mut results = Array1::zeros(x_vals.len());
        for (i, &x) in x_vals.iter().enumerate() {
            results[i] = self.evaluate_cached(x)?;
        }
        Ok(results)
    }

    /// Get access to the cache statistics
    pub fn cache_stats(&self) -> &CacheStats {
        self.cache.stats()
    }

    /// Reset cache statistics
    pub fn reset_cache_stats(&mut self) {
        self.cache.reset_stats();
    }

    /// Clear all cached data
    pub fn clear_cache(&mut self) {
        self.cache.clear();
    }

    /// Get the underlying B-spline
    pub fn spline(&self) -> &BSpline<T> {
        &self.spline
    }
}

/// Cache for distance matrices used in scattered data interpolation
#[derive(Debug)]
pub struct DistanceMatrixCache<F: Float> {
    /// Cache for computed distance matrices
    matrix_cache: HashMap<u64, Array2<F>>,
    /// Configuration
    config: CacheConfig,
    /// Statistics
    stats: CacheStats,
}

impl<F: crate::traits::InterpolationFloat> DistanceMatrixCache<F> {
    /// Create a new distance matrix cache
    pub fn new(config: CacheConfig) -> Self {
        Self {
            matrix_cache: HashMap::new(),
            config,
            stats: CacheStats::default(),
        }
    }

    /// Get a cached distance matrix or compute and cache it
    pub fn get_or_compute_distance_matrix<T>(
        &mut self,
        points: &Array2<T>,
        computer: impl FnOnce(&Array2<T>) -> Array2<F>,
    ) -> Array2<F>
    where
        T: Float,
    {
        // Create a hash of the points array for cache key
        let key = self.hash_points_safe(points);

        if let Some(cached_matrix) = self.matrix_cache.get(&key) {
            if self.config.track_stats {
                self.stats.hits += 1;
            }
            cached_matrix.clone()
        } else {
            if self.config.track_stats {
                self.stats.misses += 1;
            }

            let computed = computer(points);

            // Check cache size and evict if necessary
            if self.matrix_cache.len() >= self.config.max_distance_cache_size {
                self.evict_matrix_cache();
            }

            self.matrix_cache.insert(key, computed.clone());
            computed
        }
    }

    /// Create a hash of the points array using a safe approach for floating point
    fn hash_points_safe<T: Float>(&self, points: &Array2<T>) -> u64 {
        use std::collections::hash_map::DefaultHasher;
        let mut hasher = DefaultHasher::new();

        // Hash the shape
        points.shape()[0].hash(&mut hasher);
        points.shape()[1].hash(&mut hasher);

        // Hash a subset of points for efficiency (or all if small)
        let hash_stride = if points.len() > 1000 {
            points.len() / 100
        } else {
            1
        };

        for (i, &val) in points.iter().enumerate() {
            if i % hash_stride == 0 {
                // Convert to f64 and then to bits for consistent hashing
                let val_f64 = val.to_f64().unwrap_or(0.0);
                // Quantize to tolerance for consistent hashing of similar values
                let quantized = (val_f64 / self.config.tolerance).round() * self.config.tolerance;
                let bits = quantized.to_bits();
                bits.hash(&mut hasher);
            }
        }

        hasher.finish()
    }

    /// Evict some entries from the matrix cache
    fn evict_matrix_cache(&mut self) {
        let remove_count = self.matrix_cache.len() / 4; // Remove 25%
        let keys_to_remove: Vec<_> = self
            .matrix_cache
            .keys()
            .take(remove_count)
            .cloned()
            .collect();

        for key in keys_to_remove {
            self.matrix_cache.remove(&key);
            if self.config.track_stats {
                self.stats.evictions += 1;
            }
        }
    }

    /// Clear all cached data
    pub fn clear(&mut self) {
        self.matrix_cache.clear();
    }

    /// Get cache statistics
    pub fn stats(&self) -> &CacheStats {
        &self.stats
    }

    /// Reset cache statistics
    pub fn reset_stats(&mut self) {
        self.stats.reset();
    }
}

/// Create a cached B-spline with default cache settings
#[allow(dead_code)]
pub fn make_cached_bspline<T>(
    knots: &ArrayView1<T>,
    coeffs: &ArrayView1<T>,
    degree: usize,
    extrapolate: ExtrapolateMode,
) -> InterpolateResult<CachedBSpline<T>>
where
    T: crate::traits::InterpolationFloat,
{
    let cache = BSplineCache::default();
    CachedBSpline::new(knots, coeffs, degree, extrapolate, cache)
}

/// Create a cached B-spline with custom cache configuration
#[allow(dead_code)]
pub fn make_cached_bspline_with_config<T>(
    knots: &ArrayView1<T>,
    coeffs: &ArrayView1<T>,
    degree: usize,
    extrapolate: ExtrapolateMode,
    cache_config: CacheConfig,
) -> InterpolateResult<CachedBSpline<T>>
where
    T: crate::traits::InterpolationFloat,
{
    let cache = BSplineCache::new(cache_config);
    CachedBSpline::new(knots, coeffs, degree, extrapolate, cache)
}

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

    #[test]
    fn test_cached_bspline_evaluation() {
        // Create a simple B-spline
        let knots = array![0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0];
        let coeffs = array![1.0, 2.0, 3.0, 4.0, 5.0];

        let mut cached_spline = make_cached_bspline(
            &knots.view(),
            &coeffs.view(),
            2, // quadratic
            ExtrapolateMode::Extrapolate,
        )
        .expect("Operation failed");

        // Test evaluation at a few points
        let test_points = array![0.5, 1.0, 1.5, 2.0, 2.5];

        for &x in test_points.iter() {
            let cached_result = cached_spline.evaluate_cached(x).expect("Operation failed");
            let standard_result = cached_spline
                .evaluate_standard(x)
                .expect("Operation failed");

            // Results should be very close
            assert_relative_eq!(cached_result, standard_result, epsilon = 1e-10);
        }
    }

    #[test]
    fn test_cache_statistics() {
        let knots = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
        let coeffs = array![1.0, 2.0, 3.0, 4.0];

        let mut cached_spline = make_cached_bspline_with_config(
            &knots.view(),
            &coeffs.view(),
            2,
            ExtrapolateMode::Extrapolate,
            CacheConfig {
                track_stats: true,
                ..Default::default()
            },
        )
        .expect("Operation failed");

        // First evaluation should result in cache misses
        let _ = cached_spline
            .evaluate_cached(1.5)
            .expect("Operation failed");
        let stats_after_first = cached_spline.cache_stats();
        assert!(stats_after_first.misses > 0);

        // Second evaluation at the same point should result in cache hits
        let _ = cached_spline
            .evaluate_cached(1.5)
            .expect("Operation failed");
        let stats_after_second = cached_spline.cache_stats();
        assert!(stats_after_second.hits > 0);
    }

    #[test]
    fn test_distance_matrix_cache() {
        let config = CacheConfig {
            track_stats: true,
            max_distance_cache_size: 10,
            ..Default::default()
        };
        let mut cache = DistanceMatrixCache::<f64>::new(config);

        // Create test points
        let points = Array2::from_shape_vec((3, 2), vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0])
            .expect("Operation failed");

        // First computation should be a cache miss
        let result1 = cache.get_or_compute_distance_matrix(&points, |pts| {
            let n = pts.nrows();
            let mut distances = Array2::zeros((n, n));
            for i in 0..n {
                for j in 0..n {
                    let diff = &pts.slice(scirs2_core::ndarray::s![i, ..])
                        - &pts.slice(scirs2_core::ndarray::s![j, ..]);
                    distances[[i, j]] = diff.iter().map(|&x| x * x).sum::<f64>().sqrt();
                }
            }
            distances
        });

        assert_eq!(result1.shape(), &[3, 3]);
        assert_eq!(cache.stats().misses, 1);
        assert_eq!(cache.stats().hits, 0);

        // Second computation with same points should be a cache hit
        let result2 = cache.get_or_compute_distance_matrix(&points, |_| {
            panic!("Should not be called on cache hit");
        });

        assert_eq!(result1, result2);
        assert_eq!(cache.stats().misses, 1);
        assert_eq!(cache.stats().hits, 1);

        // Different points should be a cache miss
        let different_points =
            Array2::from_shape_vec((2, 2), vec![2.0, 2.0, 3.0, 3.0]).expect("Operation failed");
        let _result3 = cache.get_or_compute_distance_matrix(&different_points, |pts| {
            let n = pts.nrows();
            Array2::zeros((n, n))
        });

        assert_eq!(cache.stats().misses, 2);
        assert_eq!(cache.stats().hits, 1);
    }
}