scirs2-interpolate 0.4.1

//! Fast Approximate Kriging for Large Datasets
//!
//! This module provides computationally efficient kriging algorithms for large spatial datasets.
//! Standard kriging requires O(n³) operations for fitting and O(n²) for prediction,
//! which becomes prohibitively expensive for large datasets. When fully implemented, this module will provide:
//!
//! 1. **Local Kriging**: Uses only nearby points for each prediction (O(k³) complexity per prediction)
//! 2. **Fixed Rank Kriging**: Low-rank approximation of the covariance matrix (O(nr²) fitting and O(r²) prediction)
//! 3. **Tapering**: Covariance tapering for sparse approximations (reduces complexity through sparsity)
//! 4. **HODLR**: Hierarchical Off-Diagonal Low-Rank approximation (O(n log² n) complexity)
//!
//! These methods trade some accuracy for substantial performance improvements,
//! making kriging feasible for datasets with thousands to millions of points.

// IMPLEMENTATION STATUS:
//
// This module is currently undergoing significant refactoring to improve maintainability,
// performance, and integration with the rest of the library. The API is defined and stable,
// but implementations are temporarily placeholders returning NotImplemented errors.
//
// Future versions will complete this implementation with optimized algorithms leveraging:
// - Efficient matrix-free algorithms
// - Sparse matrix operations
// - Parallel execution for larger datasets
// - GPU acceleration through CUDA integration (as a feature flag)
//
// For now, users should use the standard kriging implementation for smaller datasets,
// or external libraries for very large datasets requiring approximate methods.
//
// The modular refactoring is in progress in the fast_kriging/ subdirectory.

use crate::advanced::enhanced_kriging::{AnisotropicCovariance, TrendFunction};
use crate::advanced::kriging::CovarianceFunction;
use crate::error::{InterpolateError, InterpolateResult};
use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2};
use scirs2_core::numeric::{Float, FromPrimitive};
use std::fmt::{Debug, Display};
use std::marker::PhantomData;
use std::ops::{Add, Div, Mul, Sub};

/// Type alias for sparse matrix representation
type SparseComponents<F> = (Vec<(usize, usize)>, Vec<F>);

/// Maximum number of neighbors to consider in local kriging
const DEFAULT_MAX_NEIGHBORS: usize = 50;

/// Default radius multiplier for local neighborhood search
const DEFAULT_RADIUS_MULTIPLIER: f64 = 3.0;

/// Fast kriging approximation methods for large datasets
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum FastKrigingMethod {
    /// Local kriging using only nearby points for prediction
    /// Provides O(k³) complexity per prediction where k is the neighborhood size
    Local,

    /// Fixed Rank Kriging with low-rank approximation
    /// Provides O(nr²) fitting and O(r²) prediction where r is the rank
    FixedRank(usize), // Rank parameter

    /// Tapering approach that zeros out small covariance values
    /// Creates sparse matrices for efficient computation
    Tapering(f64), // Taper range parameter

    /// Hierarchical off-diagonal low-rank approximation
    /// Balances accuracy and computational efficiency
    HODLR(usize), // Maximum leaf size
}

/// Fast approximate kriging interpolator for large datasets
///
/// This struct provides methods for efficient kriging interpolation
/// when dealing with large spatial datasets (thousands to millions of points).
/// It trades some accuracy for significant computational savings.
///
/// # Examples
///
/// ```no_run
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::advanced::fast_kriging::{
///     FastKriging, FastKrigingMethod, FastKrigingBuilder
/// };
/// use scirs2_interpolate::advanced::kriging::CovarianceFunction;
///
/// // Create sample data
/// let n_points = 100; // Reduced for testing
/// let points = Array2::<f64>::zeros((n_points, 2));
/// let values = Array1::<f64>::zeros(n_points);
///
/// // Create a fast kriging model using local approximation
/// let local_kriging = FastKrigingBuilder::<f64>::new()
///     .points(points.clone())
///     .values(values.clone())
///     .covariance_function(CovarianceFunction::Matern52)
///     .approximation_method(FastKrigingMethod::Local)
///     .max_neighbors(50)
///     .build()
///     .expect("doc example: should succeed");
///
/// // Predict at new points
/// let query_points = Array2::<f64>::zeros((10, 2));
/// let predictions = local_kriging.predict(&query_points.view()).expect("Operation failed");
///
/// // Create a model using fixed rank approximation
/// let low_rank_kriging = FastKrigingBuilder::<f64>::new()
///     .points(points.clone())
///     .values(values.clone())
///     .covariance_function(CovarianceFunction::Exponential)
///     .approximation_method(FastKrigingMethod::FixedRank(10))
///     .build()
///     .expect("doc example: should succeed");
/// # }
/// ```
#[derive(Debug, Clone)]
#[allow(dead_code)]
pub struct FastKriging<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Div<Output = F>
        + Mul<Output = F>
        + Sub<Output = F>
        + Add<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign,
{
    /// Points coordinates (input locations)
    points: Array2<F>,

    /// Values at points (observations)
    values: Array1<F>,

    /// Anisotropic covariance specification
    anisotropic_cov: AnisotropicCovariance<F>,

    /// Trend function type for Universal Kriging
    trend_fn: TrendFunction,

    /// Approximation method to use for fast computation
    approx_method: FastKrigingMethod,

    /// Maximum number of neighbors for local kriging
    max_neighbors: usize,

    /// Search radius multiplier for local kriging
    radius_multiplier: F,

    /// Pre-computed low-rank approximation components
    /// For FixedRank: [U, S, V] where K ≈ U * S * V^T
    low_rank_components: Option<(Array2<F>, Array1<F>, Array2<F>)>,

    /// Sparse representation for Tapering method
    /// Indices and values for sparse covariance matrix
    sparse_components: Option<SparseComponents<F>>,

    /// Weights for kriging predictions
    /// Different format based on approximation method
    weights: Array1<F>,

    /// Pre-computed basis functions for trend model
    basis_functions: Option<Array2<F>>,

    /// Trend coefficients for Universal Kriging
    trend_coeffs: Option<Array1<F>>,

    /// Flag indicating whether to optimize parameters
    optimize_parameters: bool,

    /// Flag indicating whether exact variance computation is needed
    compute_exact_variance: bool,

    /// Marker for generic type parameter
    _phantom: PhantomData<F>,
}

/// Builder for constructing FastKriging models
///
/// This builder provides a method-chaining interface for configuring and
/// creating FastKriging models for large datasets.
///
/// # Examples
///
/// ```no_run
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::advanced::fast_kriging::{
///     FastKrigingBuilder, FastKrigingMethod
/// };
/// use scirs2_interpolate::advanced::kriging::CovarianceFunction;
///
/// // Create sample data
/// let points = Array2::<f64>::zeros((100, 2));
/// let values = Array1::<f64>::zeros(100);
///
/// // Build a fast kriging model with local approximation
/// let kriging = FastKrigingBuilder::<f64>::new()
///     .points(points.clone())
///     .values(values.clone())
///     .covariance_function(CovarianceFunction::Matern52)
///     .approximation_method(FastKrigingMethod::Local)
///     .max_neighbors(30)
///     .radius_multiplier(2.5)
///     .build()
///     .expect("doc example: should succeed");
/// # }
/// ```
#[derive(Debug, Clone)]
pub struct FastKrigingBuilder<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Div<Output = F>
        + Mul<Output = F>
        + Sub<Output = F>
        + Add<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign,
{
    /// Points coordinates
    points: Option<Array2<F>>,

    /// Values at points
    values: Option<Array1<F>>,

    /// Covariance function
    cov_fn: CovarianceFunction,

    /// Directional length scales
    length_scales: Option<Array1<F>>,

    /// Signal variance parameter
    sigma_sq: F,

    /// Nugget parameter
    nugget: F,

    /// Trend function type
    trend_fn: TrendFunction,

    /// Approximation method
    approx_method: FastKrigingMethod,

    /// Maximum number of neighbors
    max_neighbors: usize,

    /// Search radius multiplier
    radius_multiplier: F,

    /// Marker for generic type parameter
    _phantom: PhantomData<F>,
}

/// Result type for FastKriging predictions
#[derive(Debug, Clone)]
pub struct FastPredictionResult<F: Float> {
    /// Predicted values
    pub value: Array1<F>,

    /// Approximate prediction variances
    pub variance: Array1<F>,

    /// Method used for computation
    pub method: FastKrigingMethod,

    /// Computation time in milliseconds (if available)
    pub computation_time_ms: Option<f64>,
}

impl<F> Default for FastKrigingBuilder<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = F>
        + Sub<Output = F>
        + Mul<Output = F>
        + Div<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + 'static,
{
    fn default() -> Self {
        Self::new()
    }
}

impl<F> FastKrigingBuilder<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = F>
        + Sub<Output = F>
        + Mul<Output = F>
        + Div<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + 'static,
{
    /// Create a new builder for FastKriging
    pub fn new() -> Self {
        Self {
            points: None,
            values: None,
            cov_fn: CovarianceFunction::Matern52,
            length_scales: None,
            sigma_sq: F::from_f64(1.0).expect("Operation failed"),
            nugget: F::from_f64(1e-6).expect("Operation failed"),
            trend_fn: TrendFunction::Constant,
            approx_method: FastKrigingMethod::Local,
            max_neighbors: DEFAULT_MAX_NEIGHBORS,
            radius_multiplier: F::from_f64(DEFAULT_RADIUS_MULTIPLIER).expect("Operation failed"),
            _phantom: PhantomData,
        }
    }

    /// Set points for interpolation
    pub fn points(mut self, points: Array2<F>) -> Self {
        self.points = Some(points);
        self
    }

    /// Set values for interpolation
    pub fn values(mut self, values: Array1<F>) -> Self {
        self.values = Some(values);
        self
    }

    /// Set covariance function
    pub fn covariance_function(mut self, covfn: CovarianceFunction) -> Self {
        self.cov_fn = covfn;
        self
    }

    /// Set length scales (one per dimension)
    pub fn length_scales(mut self, lengthscales: Array1<F>) -> Self {
        self.length_scales = Some(lengthscales);
        self
    }

    /// Set a single isotropic length scale
    pub fn length_scale(mut self, lengthscale: F) -> Self {
        self.sigma_sq = lengthscale;
        self
    }

    /// Set signal variance
    pub fn sigma_sq(mut self, sigmasq: F) -> Self {
        self.sigma_sq = sigmasq;
        self
    }

    /// Set nugget parameter
    pub fn nugget(mut self, nugget: F) -> Self {
        self.nugget = nugget;
        self
    }

    /// Set trend function type
    pub fn trend_function(mut self, trendfn: TrendFunction) -> Self {
        self.trend_fn = trendfn;
        self
    }

    /// Set approximation method
    pub fn approximation_method(mut self, method: FastKrigingMethod) -> Self {
        self.approx_method = method;
        self
    }

    /// Set maximum number of neighbors for local kriging
    pub fn max_neighbors(mut self, maxneighbors: usize) -> Self {
        self.max_neighbors = maxneighbors;
        self
    }

    /// Set radius multiplier for neighborhood search
    pub fn radius_multiplier(mut self, multiplier: F) -> Self {
        self.radius_multiplier = multiplier;
        self
    }

    /// Build the FastKriging model
    pub fn build(self) -> InterpolateResult<FastKriging<F>> {
        // Implementation placeholder - will be completed in future update
        Err(InterpolateError::NotImplemented(
            "FastKriging implementation is currently being refactored for improved performance and maintainability. Please use standard kriging instead for now.".to_string(),
        ))
    }
}

impl<F> FastKriging<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = F>
        + Sub<Output = F>
        + Mul<Output = F>
        + Div<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + 'static,
{
    /// Create a new builder for FastKriging
    pub fn builder() -> FastKrigingBuilder<F> {
        FastKrigingBuilder::new()
    }

    /// Predict values at new points using fast approximation
    pub fn predict(
        &self,
        _query_points: &ArrayView2<F>,
    ) -> InterpolateResult<FastPredictionResult<F>> {
        // Implementation placeholder - will be completed in future update
        Err(InterpolateError::NotImplemented(
            "FastKriging implementation is currently being refactored for improved performance and maintainability. Please use standard kriging instead for now.".to_string(),
        ))
    }
}

/// Creates a new FastKriging model with local approximation
///
/// This convenience function provides a simpler interface for creating
/// a FastKriging model using the local approximation method, which is
/// suitable for most large datasets.
///
/// # Arguments
///
/// * `points` - Coordinates of sample points (n_points × n_dimensions)
/// * `values` - Values at sample points (n_points)
/// * `cov_fn` - Covariance function to use
/// * `length_scale` - Isotropic length scale parameter
/// * `max_neighbors` - Maximum number of neighbors to use
///
/// # Returns
///
/// A FastKriging interpolator with local approximation method
///
/// # Example
///
/// ```no_run
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::advanced::fast_kriging::make_local_kriging;
/// use scirs2_interpolate::advanced::kriging::CovarianceFunction;
///
/// // Create sample data
/// let points = Array2::<f64>::zeros((100, 2));
/// let values = Array1::<f64>::zeros(100);
///
/// // Create a local kriging model
/// let kriging = make_local_kriging(
///     &points.view(),
///     &values.view(),
///     CovarianceFunction::Matern52,
///     1.0,  // length_scale
///     50    // max_neighbors
/// ).expect("Operation failed");
///
/// // Make a prediction
/// let query_point = Array2::<f64>::zeros((1, 2));
/// let pred = kriging.predict(&query_point.view()).expect("Operation failed");
/// # }
/// ```
#[allow(dead_code)]
pub fn make_local_kriging<
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = F>
        + Sub<Output = F>
        + Mul<Output = F>
        + Div<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + 'static,
>(
    _points: &ArrayView2<F>,
    _values: &ArrayView1<F>,
    _cov_fn: CovarianceFunction,
    _scale: F,
    neighbors: usize,
) -> InterpolateResult<FastKriging<F>> {
    // Implementation placeholder - will be completed in future update
    Err(InterpolateError::NotImplemented(
        "FastKriging implementation is currently being refactored for improved performance and maintainability. Please use standard kriging instead for now.".to_string(),
    ))
}

/// Creates a new FastKriging model with fixed rank approximation
///
/// This function creates a kriging model that uses a low-rank approximation
/// of the covariance matrix to enable fast computation with large datasets.
///
/// # Arguments
///
/// * `points` - Coordinates of sample points (n_points × n_dimensions)
/// * `values` - Values at sample points (n_points)
/// * `cov_fn` - Covariance function to use
/// * `length_scale` - Isotropic length scale parameter
/// * `rank` - Rank of the approximation (smaller rank = faster but less accurate)
///
/// # Returns
///
/// A FastKriging interpolator with fixed rank approximation method
///
/// # Example
///
/// ```no_run
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::advanced::fast_kriging::make_fixed_rank_kriging;
/// use scirs2_interpolate::advanced::kriging::CovarianceFunction;
///
/// // Create sample data
/// let points = Array2::<f64>::zeros((100, 2));
/// let values = Array1::<f64>::zeros(100);
///
/// // Create a fixed rank kriging model
/// let kriging = make_fixed_rank_kriging(
///     &points.view(),
///     &values.view(),
///     CovarianceFunction::Matern52,
///     1.0,  // length_scale
///     10    // rank
/// ).expect("Operation failed");
///
/// // Make a prediction
/// let query_point = Array2::<f64>::zeros((1, 2));
/// let pred = kriging.predict(&query_point.view()).expect("Operation failed");
/// # }
/// ```
#[allow(dead_code)]
pub fn make_fixed_rank_kriging<
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = F>
        + Sub<Output = F>
        + Mul<Output = F>
        + Div<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + 'static,
>(
    _points: &ArrayView2<F>,
    _values: &ArrayView1<F>,
    _cov_fn: CovarianceFunction,
    _scale: F,
    rank: usize,
) -> InterpolateResult<FastKriging<F>> {
    // Implementation placeholder - will be completed in future update
    Err(InterpolateError::NotImplemented(
        "FastKriging implementation is currently being refactored for improved performance and maintainability. Please use standard kriging instead for now.".to_string(),
    ))
}

/// Creates a new FastKriging model with tapering approximation
///
/// This function creates a kriging model that uses covariance tapering
/// to create a sparse representation, enabling fast computation with large datasets.
///
/// # Arguments
///
/// * `points` - Coordinates of sample points (n_points × n_dimensions)
/// * `values` - Values at sample points (n_points)
/// * `cov_fn` - Covariance function to use
/// * `length_scale` - Isotropic length scale parameter
/// * `taper_range` - Range beyond which covariance is set to zero
///
/// # Returns
///
/// A FastKriging interpolator with tapering approximation method
///
/// # Example
///
/// ```no_run
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::advanced::fast_kriging::make_tapered_kriging;
/// use scirs2_interpolate::advanced::kriging::CovarianceFunction;
///
/// // Create sample data
/// let points = Array2::<f64>::zeros((100, 2));
/// let values = Array1::<f64>::zeros(100);
///
/// // Create a tapered kriging model
/// let kriging = make_tapered_kriging(
///     &points.view(),
///     &values.view(),
///     CovarianceFunction::Matern52,
///     1.0,  // length_scale
///     3.0   // taper_range
/// ).expect("Operation failed");
///
/// // Make a prediction
/// let query_point = Array2::<f64>::zeros((1, 2));
/// let pred = kriging.predict(&query_point.view()).expect("Operation failed");
/// # }
/// ```
#[allow(dead_code)]
pub fn make_tapered_kriging<
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = F>
        + Sub<Output = F>
        + Mul<Output = F>
        + Div<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + 'static,
>(
    _points: &ArrayView2<F>,
    _values: &ArrayView1<F>,
    _cov_fn: CovarianceFunction,
    _scale: F,
    range: F,
) -> InterpolateResult<FastKriging<F>> {
    // Implementation placeholder - will be completed in future update
    Err(InterpolateError::NotImplemented(
        "FastKriging implementation is currently being refactored for improved performance and maintainability. Please use standard kriging instead for now.".to_string(),
    ))
}

/// Creates a new FastKriging model with HODLR approximation
///
/// This function creates a kriging model that uses a hierarchical off-diagonal
/// low-rank approximation for fast computation with large datasets.
///
/// # Arguments
///
/// * `points` - Coordinates of sample points (n_points × n_dimensions)
/// * `values` - Values at sample points (n_points)
/// * `cov_fn` - Covariance function to use
/// * `length_scale` - Isotropic length scale parameter
/// * `leaf_size` - Size of leaf nodes in the hierarchical decomposition
///
/// # Returns
///
/// A FastKriging interpolator with HODLR approximation method
///
/// # Example
///
/// ```no_run
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::advanced::fast_kriging::make_hodlr_kriging;
/// use scirs2_interpolate::advanced::kriging::CovarianceFunction;
///
/// // Create sample data
/// let points = Array2::<f64>::zeros((100, 2));
/// let values = Array1::<f64>::zeros(100);
///
/// // Create a HODLR kriging model
/// let kriging = make_hodlr_kriging(
///     &points.view(),
///     &values.view(),
///     CovarianceFunction::Matern52,
///     1.0,  // length_scale
///     32    // leaf_size
/// ).expect("Operation failed");
///
/// // Make a prediction
/// let query_point = Array2::<f64>::zeros((1, 2));
/// let pred = kriging.predict(&query_point.view()).expect("Operation failed");
/// # }
/// ```
#[allow(dead_code)]
pub fn make_hodlr_kriging<
    F: Float
        + FromPrimitive
        + Debug
        + Display
        + Add<Output = F>
        + Sub<Output = F>
        + Mul<Output = F>
        + Div<Output = F>
        + std::ops::AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::ops::RemAssign
        + 'static,
>(
    _points: &ArrayView2<F>,
    _values: &ArrayView1<F>,
    _cov_fn: CovarianceFunction,
    _scale: F,
    size: usize,
) -> InterpolateResult<FastKriging<F>> {
    // Implementation placeholder - will be completed in future update
    Err(InterpolateError::NotImplemented(
        "FastKriging implementation is currently being refactored for improved performance and maintainability. Please use standard kriging instead for now.".to_string(),
    ))
}

/// Automatically choose the best approximation method based on dataset size
///
/// This function selects an appropriate approximation method based on
/// the size of the dataset, balancing accuracy and computational efficiency.
///
/// # Returns
///
/// A FastKrigingMethod appropriate for the given dataset size
///
/// # Example
///
/// ```
/// use scirs2_interpolate::advanced::fast_kriging::select_approximation_method;
///
/// // Get recommended method for a dataset with 10,000 points
/// let method = select_approximation_method(10_000);
/// ```
#[allow(dead_code)]
pub fn select_approximation_method(_npoints: usize) -> FastKrigingMethod {
    if _npoints < 500 {
        // For small datasets, local kriging is accurate and fast enough
        FastKrigingMethod::Local
    } else if _npoints < 5_000 {
        // For medium datasets, use fixed rank with moderate rank
        FastKrigingMethod::FixedRank(50)
    } else if _npoints < 50_000 {
        // For large datasets, use tapering with moderate range
        FastKrigingMethod::Tapering(3.0)
    } else {
        // For very large datasets, use HODLR with appropriate leaf size
        FastKrigingMethod::HODLR(64)
    }
}