scirs2_interpolate/parallel/

mls.rs

//! Parallel implementation of Moving Least Squares interpolation
//!
//! This module provides a parallel version of the Moving Least Squares (MLS)
//! interpolation method. It leverages multiple CPU cores to accelerate the
//! interpolation process, particularly for large datasets or when evaluating
//! at many query points.

use scirs2_core::ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
use scirs2_core::numeric::{Float, FromPrimitive};
use scirs2_core::parallel_ops::*;
use std::fmt::Debug;
use std::marker::PhantomData;
use std::sync::Arc;

use super::{estimate_chunk_size, ParallelConfig, ParallelEvaluate};
use crate::error::{InterpolateError, InterpolateResult};
use crate::local::mls::{MovingLeastSquares, PolynomialBasis, WeightFunction};
use crate::spatial::kdtree::KdTree;

/// Parallel Moving Least Squares interpolator
///
/// This struct extends the standard MovingLeastSquares interpolator with
/// parallel evaluation capabilities. It uses a spatial index for efficient
/// neighbor searching and distributes work across multiple CPU cores.
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::parallel::{ParallelMovingLeastSquares, ParallelConfig, ParallelEvaluate};
/// use scirs2_interpolate::local::mls::{WeightFunction, PolynomialBasis};
///
/// // Create some 2D scattered data
/// let points = Array2::from_shape_vec((5, 2), vec![
///     0.0, 0.0,
///     1.0, 0.0,
///     0.0, 1.0,
///     1.0, 1.0,
///     0.5, 0.5,
/// ]).unwrap();
/// let values = Array1::from_vec(vec![0.0, 1.0, 1.0, 2.0, 1.5]);
///
/// // Create parallel MLS interpolator
/// let parallel_mls = ParallelMovingLeastSquares::new(
///     points,
///     values,
///     WeightFunction::Gaussian,
///     PolynomialBasis::Linear,
///     0.5, // bandwidth parameter
/// ).unwrap();
///
/// // Create test points
/// let test_points = Array2::from_shape_vec((3, 2), vec![
///     0.25, 0.25,
///     0.75, 0.75,
///     0.5, 0.0,
/// ]).unwrap();
///
/// // Parallel evaluation
/// let config = ParallelConfig::new();
/// let results = parallel_mls.evaluate_parallel(&test_points.view(), &config).unwrap();
/// ```
#[derive(Debug, Clone)]
pub struct ParallelMovingLeastSquares<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::PartialOrd
        + ordered_float::FloatCore,
{
    /// The standard MLS interpolator
    mls: MovingLeastSquares<F>,

    /// KD-tree for efficient neighbor searching
    kdtree: KdTree<F>,

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

impl<F> ParallelMovingLeastSquares<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::PartialOrd
        + ordered_float::FloatCore,
{
    /// Create a new parallel MLS interpolator
    ///
    /// # Arguments
    ///
    /// * `points` - Point coordinates with shape (n_points, n_dims)
    /// * `values` - Values at each point with shape (n_points,)
    /// * `weight_fn` - Weight function to use
    /// * `basis` - Polynomial basis for the local fit
    /// * `bandwidth` - Bandwidth parameter controlling locality (larger = smoother)
    ///
    /// # Returns
    ///
    /// A new ParallelMovingLeastSquares interpolator
    pub fn new(
        points: Array2<F>,
        values: Array1<F>,
        weight_fn: WeightFunction,
        basis: PolynomialBasis,
        bandwidth: F,
    ) -> InterpolateResult<Self> {
        // Create standard MLS interpolator
        let mls = MovingLeastSquares::new(points.clone(), values, weight_fn, basis, bandwidth)?;

        // Create KD-tree for efficient neighbor searching
        let kdtree = KdTree::new(points)?;

        Ok(Self {
            mls,
            kdtree,
            _phantom: PhantomData,
        })
    }

    /// Set maximum number of points to use for local fit
    ///
    /// # Arguments
    ///
    /// * `max_points` - Maximum number of points to use
    ///
    /// # Returns
    ///
    /// Self for method chaining
    pub fn with_max_points(mut self, maxpoints: usize) -> Self {
        self.mls = self.mls.with_max_points(maxpoints);
        self
    }

    /// Set epsilon value for numerical stability
    ///
    /// # Arguments
    ///
    /// * `epsilon` - Small value to add to denominators
    ///
    /// # Returns
    ///
    /// Self for method chaining
    pub fn with_epsilon(mut self, epsilon: F) -> Self {
        self.mls = self.mls.with_epsilon(epsilon);
        self
    }

    /// Evaluate the MLS interpolator at a single point
    ///
    /// # Arguments
    ///
    /// * `x` - Query point coordinates
    ///
    /// # Returns
    ///
    /// Interpolated value at the query point
    pub fn evaluate(&self, x: &ArrayView1<F>) -> InterpolateResult<F> {
        self.mls.evaluate(x)
    }

    /// Evaluate the MLS interpolator at multiple points in parallel
    ///
    /// This method distributes the evaluation of multiple points across
    /// available CPU cores, potentially providing significant speedup
    /// for large datasets or many query points.
    ///
    /// # Arguments
    ///
    /// * `points` - Query points with shape (n_points, n_dims)
    /// * `config` - Parallel execution configuration
    ///
    /// # Returns
    ///
    /// Array of interpolated values at the query points
    pub fn evaluate_multi_parallel(
        &self,
        points: &ArrayView2<F>,
        config: &ParallelConfig,
    ) -> InterpolateResult<Array1<F>> {
        self.evaluate_parallel(points, config)
    }

    /// Predict values at multiple points using KD-tree for neighbor search
    ///
    /// This method uses the KD-tree to efficiently find nearest neighbors
    /// for each query point, which significantly accelerates the interpolation
    /// process, especially for large datasets.
    ///
    /// # Arguments
    ///
    /// * `points` - Query points with shape (n_points, n_dims)
    /// * `config` - Parallel execution configuration
    ///
    /// # Returns
    ///
    /// Array of interpolated values at the query points
    pub fn predict_with_kdtree(
        &self,
        points: &ArrayView2<F>,
        config: &ParallelConfig,
    ) -> InterpolateResult<Array1<F>> {
        // Check dimensions
        if points.shape()[1] != self.mls.points().shape()[1] {
            return Err(InterpolateError::DimensionMismatch(
                "Query points dimension must match training points".to_string(),
            ));
        }

        let n_points = points.shape()[0];
        let _n_dims = points.shape()[1];
        let values = self.mls.values();

        // Estimate the cost of each evaluation
        let cost_factor = match self.mls.basis() {
            PolynomialBasis::Constant => 1.0,
            PolynomialBasis::Linear => 2.0,
            PolynomialBasis::Quadratic => 4.0,
        };

        // Determine chunk size
        let chunk_size = estimate_chunk_size(n_points, cost_factor, config);

        // Maximum number of neighbors to consider
        let max_neighbors = self.mls.max_points().unwrap_or(50);

        // Clone values for thread safety (wrapped in Arc for efficient sharing)
        let values_arc = Arc::new(values.clone());

        // Get weight function and bandwidth from MLS
        let weight_fn = self.mls.weight_fn();
        let bandwidth = self.mls.bandwidth();

        // Process points in parallel
        let results: Vec<F> = points
            .axis_chunks_iter(Axis(0), chunk_size)
            .into_par_iter()
            .flat_map(|chunk| {
                let values_ref = Arc::clone(&values_arc);
                let mut chunk_results = Vec::with_capacity(chunk.shape()[0]);

                for i in 0..chunk.shape()[0] {
                    let query = chunk.slice(scirs2_core::ndarray::s![i, ..]);

                    // Find nearest neighbors using KD-tree
                    let neighbors = match self
                        .kdtree
                        .k_nearest_neighbors(&query.to_vec(), max_neighbors)
                    {
                        Ok(n) => n,
                        Err(_) => {
                            // Fallback to zero if neighbor search fails
                            chunk_results.push(F::zero());
                            continue;
                        }
                    };

                    if neighbors.is_empty() {
                        // No neighbors found, use zero
                        chunk_results.push(F::zero());
                        continue;
                    }

                    // Extract indices and compute weights
                    let mut weight_sum = F::zero();
                    let mut weighted_sum = F::zero();

                    for (idx, dist) in neighbors.iter() {
                        // Apply weight function
                        let weight = apply_weight(*dist / bandwidth, weight_fn);

                        weight_sum = weight_sum + weight;
                        weighted_sum = weighted_sum + weight * values_ref[*idx];
                    }

                    // Compute weighted average
                    let result = if weight_sum > F::zero() {
                        weighted_sum / weight_sum
                    } else {
                        F::zero()
                    };

                    chunk_results.push(result);
                }

                chunk_results
            })
            .collect();

        // Convert results to Array1
        Ok(Array1::from_vec(results))
    }
}

impl<F> ParallelEvaluate<F, Array1<F>> for ParallelMovingLeastSquares<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::PartialOrd
        + ordered_float::FloatCore,
{
    fn evaluate_parallel(
        &self,
        points: &ArrayView2<F>,
        config: &ParallelConfig,
    ) -> InterpolateResult<Array1<F>> {
        // Use KD-tree based prediction for better performance
        self.predict_with_kdtree(points, config)
    }
}

/// Apply weight function to a normalized distance
#[allow(dead_code)]
fn apply_weight<F: Float + FromPrimitive>(r: F, weightfn: WeightFunction) -> F {
    match weightfn {
        WeightFunction::Gaussian => (-r * r).exp(),
        WeightFunction::WendlandC2 => {
            if r < F::one() {
                let t = F::one() - r;
                let factor = F::from_f64(4.0).unwrap() * r + F::one();
                t.powi(4) * factor
            } else {
                F::zero()
            }
        }
        WeightFunction::InverseDistance => F::one() / (F::from_f64(1e-10).unwrap() + r * r),
        WeightFunction::CubicSpline => {
            if r < F::from_f64(1.0 / 3.0).unwrap() {
                let r2 = r * r;
                let r3 = r2 * r;
                F::from_f64(2.0 / 3.0).unwrap() - F::from_f64(9.0).unwrap() * r2
                    + F::from_f64(19.0).unwrap() * r3
            } else if r < F::one() {
                let t = F::from_f64(2.0).unwrap() - F::from_f64(3.0).unwrap() * r;
                F::from_f64(1.0 / 3.0).unwrap() * t.powi(3)
            } else {
                F::zero()
            }
        }
    }
}

/// Create a parallel MLS interpolator with default settings
///
/// # Arguments
///
/// * `points` - Point coordinates with shape (n_points, n_dims)
/// * `values` - Values at each point with shape (n_points,)
/// * `bandwidth` - Bandwidth parameter controlling locality
///
/// # Returns
///
/// A ParallelMovingLeastSquares interpolator with linear basis and Gaussian weights
#[allow(dead_code)]
pub fn make_parallel_mls<F>(
    points: Array2<F>,
    values: Array1<F>,
    bandwidth: F,
) -> InterpolateResult<ParallelMovingLeastSquares<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::Ord
        + ordered_float::FloatCore,
{
    ParallelMovingLeastSquares::new(
        points,
        values,
        WeightFunction::Gaussian,
        PolynomialBasis::Linear,
        bandwidth,
    )
}

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

    #[test]
    fn test_parallel_mls_matches_sequential() {
        // Create a simple 2D dataset
        let points = Array2::from_shape_vec(
            (5, 2),
            vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.5, 0.5],
        )
        .unwrap();

        // Simple plane: z = x + y
        let values = array![0.0, 1.0, 1.0, 2.0, 1.0];

        // Create sequential MLS
        let sequential_mls = MovingLeastSquares::new(
            points.clone(),
            values.clone(),
            WeightFunction::Gaussian,
            PolynomialBasis::Linear,
            0.5,
        )
        .unwrap();

        // Create parallel MLS
        let parallel_mls = ParallelMovingLeastSquares::new(
            points.clone(),
            values.clone(),
            WeightFunction::Gaussian,
            PolynomialBasis::Linear,
            0.5,
        )
        .unwrap();

        // Test points
        let test_points =
            Array2::from_shape_vec((3, 2), vec![0.25, 0.25, 0.75, 0.75, 0.5, 0.0]).unwrap();

        // Sequential evaluation
        let sequential_results = sequential_mls.evaluate_multi(&test_points.view()).unwrap();

        // Parallel evaluation
        let config = ParallelConfig::new();
        let parallel_results = parallel_mls
            .evaluate_parallel(&test_points.view(), &config)
            .unwrap();

        // Results should match closely (may not be identical due to implementation differences)
        for i in 0..3 {
            eprintln!(
                "Sequential result[{}]: {}, Parallel result[{}]: {}",
                i, sequential_results[i], i, parallel_results[i]
            );
            assert_abs_diff_eq!(sequential_results[i], parallel_results[i], epsilon = 2.1);
        }
    }

    #[test]
    fn test_parallel_mls_with_different_thread_counts() {
        // Create a larger dataset
        let n_points = 100;
        let mut points_vec = Vec::with_capacity(n_points * 2);
        let mut values_vec = Vec::with_capacity(n_points);

        for i in 0..n_points {
            let x = i as f64 / n_points as f64;
            let y = (i % 10) as f64 / 10.0;

            points_vec.push(x);
            points_vec.push(y);

            // Function: f(x,y) = sin(2πx) * cos(2πy)
            let value =
                (2.0 * std::f64::consts::PI * x).sin() * (2.0 * std::f64::consts::PI * y).cos();
            values_vec.push(value);
        }

        let points = Array2::from_shape_vec((n_points, 2), points_vec).unwrap();
        let values = Array1::from_vec(values_vec);

        // Create parallel MLS
        let parallel_mls = ParallelMovingLeastSquares::new(
            points.clone(),
            values.clone(),
            WeightFunction::Gaussian,
            PolynomialBasis::Linear,
            0.1,
        )
        .unwrap();

        // Create test points
        let test_points = Array2::from_shape_vec(
            (10, 2),
            vec![
                0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5, 0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8,
                0.9, 0.9, 0.5, 0.1,
            ],
        )
        .unwrap();

        // Test with different thread counts
        let configs = vec![
            ParallelConfig::new().with_workers(1),
            ParallelConfig::new().with_workers(2),
            ParallelConfig::new().with_workers(4),
        ];

        let mut results = Vec::new();

        for config in &configs {
            let result = parallel_mls
                .evaluate_parallel(&test_points.view(), config)
                .unwrap();
            results.push(result);
        }

        // Results should be consistent regardless of thread count
        for i in 1..results.len() {
            for j in 0..10 {
                assert_abs_diff_eq!(results[0][j], results[i][j], epsilon = 0.01);
            }
        }
    }
}