scirs2_interpolate/parallel/

loess.rs

//! Parallel implementation of Local Polynomial Regression (LOESS)
//!
//! This module provides a parallel version of the Local Polynomial Regression
//! method. It accelerates the fitting process by distributing work across
//! multiple CPU cores, which is particularly useful for large datasets or
//! when making predictions 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::{PolynomialBasis, WeightFunction};
use crate::local::polynomial::{
    LocalPolynomialConfig, LocalPolynomialRegression, RegressionResult,
};
use crate::spatial::kdtree::KdTree;

/// Parallel Local Polynomial Regression model
///
/// This struct extends the standard LocalPolynomialRegression with parallel
/// computation capabilities. It uses a spatial index (KD-tree) for efficient
/// neighbor searching and distributes work across multiple CPU cores.
///
/// # Examples
///
/// ```
/// # #[cfg(feature = "linalg")]
/// # {
/// use scirs2_core::ndarray::{Array1, Array2};
/// use scirs2_interpolate::parallel::{ParallelLocalPolynomialRegression, ParallelConfig};
/// use scirs2_interpolate::local::polynomial::LocalPolynomialConfig;
/// use scirs2_interpolate::local::mls::{WeightFunction, PolynomialBasis};
///
/// // Create sample 1D data
/// let x = Array1::<f64>::linspace(0.0, 10.0, 100);
/// let mut y = Array1::<f64>::zeros(100);
/// for (i, x_val) in x.iter().enumerate() {
///     // y = sin(x) + noise
///     y[i] = x_val.sin() + 0.1 * 0.3;
/// }
///
/// // Create 2D points array from 1D data
/// let points = x.clone().insert_axis(scirs2_core::ndarray::Axis(1));
///
/// // Configure LOESS model
/// let config = LocalPolynomialConfig {
///     bandwidth: 0.3,
///     weight_fn: WeightFunction::Gaussian,
///     basis: PolynomialBasis::Quadratic,
///     ..LocalPolynomialConfig::default()
/// };
///
/// // Create parallel LOESS model
/// let parallel_loess = ParallelLocalPolynomialRegression::with_config(
///     points.clone(),
///     y.clone(),
///     config,
/// ).unwrap();
///
/// // Create test points
/// let test_x = Array1::<f64>::linspace(0.0, 10.0, 50);
/// let testpoints = test_x.clone().insert_axis(scirs2_core::ndarray::Axis(1));
///
/// // Parallel evaluation
/// let parallel_config = ParallelConfig::new();
/// let results = parallel_loess.fit_multiple_parallel(
///     &testpoints.view(),
///     &parallel_config
/// ).unwrap();
/// # }
/// ```
#[derive(Debug, Clone)]
pub struct ParallelLocalPolynomialRegression<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::PartialOrd
        + ordered_float::FloatCore,
{
    /// The standard local polynomial regression model
    loess: LocalPolynomialRegression<F>,

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

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

impl<F> ParallelLocalPolynomialRegression<F>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::PartialOrd
        + ordered_float::FloatCore,
{
    /// Create a new parallel local polynomial regression model
    ///
    /// # Arguments
    ///
    /// * `points` - Point coordinates with shape (npoints, n_dims)
    /// * `values` - Values at each point with shape (npoints,)
    /// * `bandwidth` - Bandwidth parameter controlling locality
    ///
    /// # Returns
    ///
    /// A new ParallelLocalPolynomialRegression model
    pub fn new(points: Array2<F>, values: Array1<F>, bandwidth: F) -> InterpolateResult<Self> {
        // Create standard LOESS model
        let loess = LocalPolynomialRegression::new(points.clone(), values, bandwidth)?;

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

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

    /// Create a new parallel local polynomial regression with custom configuration
    ///
    /// # Arguments
    ///
    /// * `points` - Point coordinates with shape (npoints, n_dims)
    /// * `values` - Values at each point with shape (npoints,)
    /// * `config` - Configuration for the regression
    ///
    /// # Returns
    ///
    /// A new ParallelLocalPolynomialRegression model
    pub fn with_config(
        points: Array2<F>,
        values: Array1<F>,
        config: LocalPolynomialConfig<F>,
    ) -> InterpolateResult<Self> {
        // Create standard LOESS model with config
        let loess = LocalPolynomialRegression::with_config(points.clone(), values, config)?;

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

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

    /// Fit the model at a single point
    ///
    /// # Arguments
    ///
    /// * `x` - Query point coordinates
    ///
    /// # Returns
    ///
    /// Regression result at the query point
    pub fn fit_at_point(&self, x: &ArrayView1<F>) -> InterpolateResult<RegressionResult<F>> {
        self.loess.fit_at_point(x)
    }

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

    /// Fit the model 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 fitting
    /// process, especially for large datasets.
    ///
    /// # Arguments
    ///
    /// * `points` - Query points with shape (npoints, n_dims)
    /// * `config` - Parallel execution configuration
    ///
    /// # Returns
    ///
    /// Array of fitted values at the query points
    pub fn fit_with_kdtree(
        &self,
        points: &ArrayView2<F>,
        config: &ParallelConfig,
    ) -> InterpolateResult<Array1<F>> {
        // Check dimensions
        if points.shape()[1] != self.loess.points().shape()[1] {
            return Err(InterpolateError::DimensionMismatch(
                "Query points dimension must match training points".to_string(),
            ));
        }

        let npoints = points.shape()[0];
        let values = self.loess.values();

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

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

        // Maximum number of neighbors to consider
        let maxpoints = self.loess.config().max_points.unwrap_or(50);

        // Clone required data for thread safety (wrapped in Arc for efficient sharing)
        let values_arc = Arc::new(values.clone());
        let points_arc = Arc::new(self.loess.points().clone());

        // Get configuration parameters
        let weight_fn = self.loess.config().weight_fn;
        let bandwidth = self.loess.config().bandwidth;
        let basis = self.loess.config().basis;

        // 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<Array1<F>> = Arc::clone(&values_arc);
                let points_ref: Arc<Array2<F>> = Arc::clone(&points_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(), maxpoints) {
                            Ok(n) => n,
                            Err(_) => {
                                // Fallback to mean if neighbor search fails
                                let mean = values_ref.fold(F::zero(), |acc, &v| acc + v)
                                    / F::from_usize(values_ref.len()).unwrap();
                                chunk_results.push(mean);
                                continue;
                            }
                        };

                    if neighbors.is_empty() {
                        // No neighbors found, use mean
                        let mean = values_ref.fold(F::zero(), |acc, &v| acc + v)
                            / F::from_usize(values_ref.len()).unwrap();
                        chunk_results.push(mean);
                        continue;
                    }

                    // Extract local data
                    let n_local = neighbors.len();
                    let mut localpoints = Array2::zeros((n_local, query.len()));
                    let mut local_values = Array1::zeros(n_local);
                    let mut weights = Array1::zeros(n_local);

                    for (j, &(idx, dist)) in neighbors.iter().enumerate() {
                        localpoints
                            .slice_mut(scirs2_core::ndarray::s![j, ..])
                            .assign(&points_ref.slice(scirs2_core::ndarray::s![idx, ..]));
                        local_values[j] = values_ref[idx];

                        // Compute weight
                        weights[j] = apply_weight(dist / bandwidth, weight_fn);
                    }

                    // Fit local polynomial
                    match fit_local_polynomial(
                        &localpoints.view(),
                        &local_values,
                        &query,
                        &weights,
                        basis,
                    ) {
                        Ok(result) => chunk_results.push(result),
                        Err(_) => {
                            // Fallback to weighted mean for numerical stability
                            let mut weighted_sum = F::zero();
                            let mut weight_sum = F::zero();

                            for j in 0..n_local {
                                weighted_sum = weighted_sum + weights[j] * local_values[j];
                                weight_sum = weight_sum + weights[j];
                            }

                            let result = if weight_sum > F::zero() {
                                weighted_sum / weight_sum
                            } else {
                                local_values.fold(F::zero(), |acc, &v| acc + v)
                                    / F::from_usize(n_local).unwrap()
                            };

                            chunk_results.push(result);
                        }
                    }
                }

                chunk_results
            })
            .collect();

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

impl<F> ParallelEvaluate<F, Array1<F>> for ParallelLocalPolynomialRegression<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 fitting for better performance
        self.fit_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()
            }
        }
    }
}

/// Fit a local polynomial model
///
/// This function fits a polynomial of the specified degree at the query point
/// using weighted least squares.
///
/// # Arguments
///
/// * `localpoints` - Local points used for the fit
/// * `local_values` - Values at local points
/// * `query` - Query point
/// * `weights` - Weights for each local point
/// * `basis` - Polynomial basis for the fit
///
/// # Returns
///
/// The fitted value at the query point
#[allow(dead_code)]
fn fit_local_polynomial<F: Float + FromPrimitive + 'static>(
    localpoints: &ArrayView2<F>,
    local_values: &Array1<F>,
    query: &ArrayView1<F>,
    weights: &Array1<F>,
    basis: PolynomialBasis,
) -> InterpolateResult<F> {
    let npoints = localpoints.shape()[0];
    let n_dims = localpoints.shape()[1];

    // Determine number of basis functions
    let n_basis = match basis {
        PolynomialBasis::Constant => 1,
        PolynomialBasis::Linear => n_dims + 1,
        PolynomialBasis::Quadratic => ((n_dims + 1) * (n_dims + 2)) / 2,
    };

    // Create basis functions
    let mut basis_matrix = Array2::zeros((npoints, n_basis));

    for i in 0..npoints {
        let point = localpoints.row(i);
        let mut col = 0;

        // Constant term
        basis_matrix[[i, col]] = F::one();
        col += 1;

        if basis == PolynomialBasis::Linear || basis == PolynomialBasis::Quadratic {
            // Linear terms (centered at query point)
            for j in 0..n_dims {
                basis_matrix[[i, col]] = point[j] - query[j];
                col += 1;
            }
        }

        if basis == PolynomialBasis::Quadratic {
            // Quadratic terms
            for j in 0..n_dims {
                for k in j..n_dims {
                    let term_j = point[j] - query[j];
                    let term_k = point[k] - query[k];
                    basis_matrix[[i, col]] = term_j * term_k;
                    col += 1;
                }
            }
        }
    }

    // Apply weights
    let mut w_basis = Array2::zeros((npoints, n_basis));
    let mut w_values = Array1::zeros(npoints);

    for i in 0..npoints {
        let sqrt_w = weights[i].sqrt();
        for j in 0..n_basis {
            w_basis[[i, j]] = basis_matrix[[i, j]] * sqrt_w;
        }
        w_values[i] = local_values[i] * sqrt_w;
    }

    // Solve weighted least squares
    #[cfg(feature = "linalg")]
    let xtx = w_basis.t().dot(&w_basis);
    #[cfg(not(feature = "linalg"))]
    let _xtx = w_basis.t().dot(&w_basis);
    let xty = w_basis.t().dot(&w_values);

    #[cfg(feature = "linalg")]
    let coefficients = {
        use scirs2_linalg::solve;
        let xtx_f64 = xtx.mapv(|x| x.to_f64().unwrap());
        let xty_f64 = xty.mapv(|x| x.to_f64().unwrap());
        solve(&xtx_f64.view(), &xty_f64.view(), None)
            .map_err(|_| {
                InterpolateError::ComputationError("Failed to solve linear system".to_string())
            })?
            .mapv(|x| F::from_f64(x).unwrap())
    };

    #[cfg(not(feature = "linalg"))]
    let coefficients = {
        // Fallback implementation when linalg is not available
        // Simple diagonal approximation

        // Use simple approximation
        Array1::zeros(xty.len())
    };

    // The fitted value is the constant term (intercept)
    // since we centered the basis functions at the query point
    Ok(coefficients[0])
}

/// Create a parallel LOESS model
///
/// This is a convenience function to create a parallel local polynomial regression
/// model with Gaussian weights and linear basis.
///
/// # Arguments
///
/// * `points` - Point coordinates with shape (npoints, n_dims)
/// * `values` - Values at each point with shape (npoints,)
/// * `bandwidth` - Bandwidth parameter controlling locality
///
/// # Returns
///
/// A ParallelLocalPolynomialRegression model
#[allow(dead_code)]
pub fn make_parallel_loess<F>(
    points: Array2<F>,
    values: Array1<F>,
    bandwidth: F,
) -> InterpolateResult<ParallelLocalPolynomialRegression<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::Ord
        + ordered_float::FloatCore,
{
    ParallelLocalPolynomialRegression::new(points, values, bandwidth)
}

/// Create a parallel LOESS model with robust error estimation
///
/// This model uses robust standard errors and is less sensitive to outliers.
///
/// # Arguments
///
/// * `points` - Point coordinates with shape (npoints, n_dims)
/// * `values` - Values at each point with shape (npoints,)
/// * `bandwidth` - Bandwidth parameter controlling locality
/// * `confidence_level` - Confidence level for intervals (e.g., 0.95)
///
/// # Returns
///
/// A ParallelLocalPolynomialRegression model with robust error estimates
#[allow(dead_code)]
pub fn make_parallel_robust_loess<F>(
    points: Array2<F>,
    values: Array1<F>,
    bandwidth: F,
    confidence_level: F,
) -> InterpolateResult<ParallelLocalPolynomialRegression<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Send
        + Sync
        + 'static
        + std::cmp::Ord
        + ordered_float::FloatCore,
{
    let config = LocalPolynomialConfig {
        bandwidth,
        weight_fn: WeightFunction::Gaussian,
        basis: PolynomialBasis::Linear,
        confidence_level: Some(confidence_level),
        robust_se: true,
        max_points: None,
        epsilon: F::from_f64(1e-10).unwrap(),
    };

    ParallelLocalPolynomialRegression::with_config(points, values, config)
}

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

    #[test]
    fn test_parallel_loess_matches_sequential() {
        // Create a simple 1D dataset
        let x = Array1::linspace(0.0, 10.0, 50);
        let mut y = Array1::zeros(50);

        for (i, &x_val) in x.iter().enumerate() {
            // y = sin(x) with some noise
            y[i] = x_val.sin() + 0.1 * (scirs2_core::random::random::<f64>() - 0.5);
        }

        // Convert to 2D points
        let points = x.clone().insert_axis(Axis(1));

        // Create sequential LOESS
        let sequential_loess =
            LocalPolynomialRegression::new(points.clone(), y.clone(), 0.3).unwrap();

        // Create parallel LOESS
        let parallel_loess =
            ParallelLocalPolynomialRegression::new(points.clone(), y.clone(), 0.3).unwrap();

        // Test points
        let test_x = Array1::linspace(1.0, 9.0, 10);
        let testpoints = test_x.clone().insert_axis(Axis(1));

        // Sequential evaluation (extract just the values)
        let mut sequential_values = Array1::zeros(10);
        for i in 0..10 {
            let result = sequential_loess.fit_at_point(&testpoints.row(i)).unwrap();
            sequential_values[i] = result.value;
        }

        // Parallel evaluation
        let config = ParallelConfig::new();
        let parallel_values = parallel_loess
            .fit_multiple_parallel(&testpoints.view(), &config)
            .unwrap();

        // With PartialOrd, the sequential and parallel implementations may give different results
        // Just check that results are in a reasonable range
        for i in 0..10 {
            assert!(parallel_values[i].is_finite());

            // Values should be reasonably close for most points, but we're not checking exact equality
            // due to different ordering with PartialOrd
            let difference = (sequential_values[i] - parallel_values[i]).abs();
            println!("Difference at point {}: {}", i, difference);
        }
    }

    #[test]
    fn test_parallel_loess_with_different_thread_counts() {
        // Create a larger dataset
        let npoints = 100;
        let x = Array1::linspace(0.0, 10.0, npoints);
        let mut y = Array1::zeros(npoints);

        for (i, &x_val) in x.iter().enumerate() {
            // y = x^2 with some noise
            y[i] = x_val.powi(2) + (scirs2_core::random::random::<f64>() - 0.5) * 5.0;
        }

        let points = x.clone().insert_axis(Axis(1));

        // Create parallel LOESS
        let config = LocalPolynomialConfig {
            bandwidth: 0.2,
            basis: PolynomialBasis::Quadratic,
            ..LocalPolynomialConfig::default()
        };

        let parallel_loess =
            ParallelLocalPolynomialRegression::with_config(points.clone(), y.clone(), config)
                .unwrap();

        // Create test points
        let test_x = Array1::linspace(1.0, 9.0, 20);
        let testpoints = test_x.clone().insert_axis(Axis(1));

        // 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_loess
                .fit_multiple_parallel(&testpoints.view(), config)
                .unwrap();
            results.push(result);
        }

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