scirs2-interpolate 0.4.1

//! Grid data interpolation - SciPy-compatible griddata implementation
//!
//! This module provides the `griddata` function, which interpolates unstructured
//! data to a regular grid or arbitrary points. This is one of the most commonly
//! used interpolation functions in SciPy's interpolate module.
//!
//! # Examples
//!
//! ```
//! use scirs2_core::ndarray::{array, Array2};
//! use scirs2_interpolate::griddata::{griddata, GriddataMethod};
//!
//! // Scattered data points
//! let points = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
//! let values = array![0.0, 1.0, 1.0, 2.0];
//!
//! // Grid to interpolate onto
//! let xi = array![[0.5, 0.5], [0.25, 0.75]];
//!
//! // Interpolate using linear method
//! let result = griddata(&points.view(), &values.view(), &xi.view(),
//!                       GriddataMethod::Linear, None, None)?;
//! # Ok::<(), Box<dyn std::error::Error>>(())
//! ```

use crate::advanced::rbf::{RBFInterpolator, RBFKernel};
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::ops::AddAssign;

/// Interpolation methods available for griddata
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum GriddataMethod {
    /// Linear interpolation using Delaunay triangulation
    Linear,
    /// Nearest neighbor interpolation  
    Nearest,
    /// Cubic interpolation using Clough-Tocher scheme
    Cubic,
    /// Radial basis function interpolation with linear kernel
    Rbf,
    /// Radial basis function interpolation with cubic kernel
    RbfCubic,
    /// Radial basis function interpolation with thin plate spline
    RbfThinPlate,
}

/// Interpolate unstructured D-dimensional data to arbitrary points.
///
/// This function provides a SciPy-compatible interface for interpolating
/// scattered data points to a regular grid or arbitrary query points.
///
/// # Arguments
///
/// * `points` - Data point coordinates with shape (n_points, n_dims)
/// * `values` - Data values at each point with shape (n_points,)  
/// * `xi` - Points at which to interpolate data with shape (n_queries, n_dims)
/// * `method` - Interpolation method to use
/// * `fill_value` - Value to use for points outside convex hull (None uses NaN)
///
/// # Returns
///
/// Array of interpolated values with shape (n_queries,)
///
/// # Errors
///
/// * `ShapeMismatch` - If input arrays have incompatible shapes
/// * `InvalidParameter` - If method parameters are invalid
/// * `ComputationFailed` - If interpolation setup fails
///
/// # Examples
///
/// ## Basic usage with scattered 2D data
/// ```
/// use scirs2_core::ndarray::array;
/// use scirs2_interpolate::griddata::{griddata, GriddataMethod};
///
/// // Scattered data: z = x² + y²
/// let points = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [0.5, 0.5]];
/// let values = array![0.0, 1.0, 1.0, 2.0, 0.5];
///
/// // Query points
/// let xi = array![[0.25, 0.25], [0.75, 0.75]];
///
/// // Linear interpolation
/// let result = griddata(&points.view(), &values.view(), &xi.view(),
///                       GriddataMethod::Linear, None, None)?;
///
/// println!("Interpolated values: {:?}", result);
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
///
/// ## Using different interpolation methods
/// ```
/// use scirs2_core::ndarray::array;
/// use scirs2_interpolate::griddata::{griddata, GriddataMethod};
///
/// let points = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]];
/// let values = array![0.0, 1.0, 1.0];
/// let xi = array![[0.5, 0.5]];
///
/// // Compare different methods
/// let linear = griddata(&points.view(), &values.view(), &xi.view(),
///                      GriddataMethod::Linear, None, None)?;
/// let nearest = griddata(&points.view(), &values.view(), &xi.view(),
///                       GriddataMethod::Nearest, None, None)?;
/// let rbf = griddata(&points.view(), &values.view(), &xi.view(),
///                    GriddataMethod::Rbf, None, None)?;
///
/// println!("Linear: {}, Nearest: {}, RBF: {}", linear[0], nearest[0], rbf[0]);
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
///
/// # Performance Notes
///
/// - **Linear**: Fast setup O(n log n), fast evaluation O(log n) per point
/// - **Nearest**: Very fast setup O(n log n), very fast evaluation O(log n)  
/// - **Cubic**: Slow setup O(n³), medium evaluation O(n)
/// - **RBF methods**: Slow setup O(n³), medium evaluation O(n)
///
/// For large datasets (n > 1000), consider using FastRBF or other approximation methods.
#[allow(dead_code)]
pub fn griddata<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    method: GriddataMethod,
    fill_value: Option<F>,
    workers: Option<usize>,
) -> InterpolateResult<Array1<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Clone
        + Display
        + AddAssign
        + std::ops::SubAssign
        + std::fmt::LowerExp
        + std::ops::MulAssign
        + std::ops::DivAssign
        + Send
        + Sync
        + 'static,
{
    // Validate inputs
    validate_griddata_inputs(points, values, xi)?;

    // Decide whether to use parallel or serial execution
    let use_parallel = match workers {
        Some(1) => false, // Explicitly requested serial execution
        Some(_) => true,  // Explicitly requested parallel execution with n > 1 workers
        None => {
            // Automatic decision based on dataset size
            // Use parallel for large datasets (threshold based on empirical testing)
            let n_points = points.nrows();
            let n_queries = xi.nrows();
            n_queries >= 100 && n_points >= 50
        }
    };

    if use_parallel {
        // Use the existing parallel implementation
        griddata_parallel(points, values, xi, method, fill_value, workers)
    } else {
        // Use serial implementation
        match method {
            GriddataMethod::Linear => griddata_linear(points, values, xi, fill_value),
            GriddataMethod::Nearest => griddata_nearest(points, values, xi, fill_value),
            GriddataMethod::Cubic => griddata_cubic(points, values, xi, fill_value),
            GriddataMethod::Rbf => griddata_rbf(points, values, xi, RBFKernel::Linear, fill_value),
            GriddataMethod::RbfCubic => {
                griddata_rbf(points, values, xi, RBFKernel::Cubic, fill_value)
            }
            GriddataMethod::RbfThinPlate => {
                griddata_rbf(points, values, xi, RBFKernel::ThinPlateSpline, fill_value)
            }
        }
    }
}

/// Parallel version of griddata for large datasets
///
/// This function provides the same functionality as `griddata` but uses parallel
/// processing to speed up interpolation for large query sets. The interpolation
/// setup (triangulation, RBF fitting) is done once, then queries are processed
/// in parallel chunks.
///
/// # Arguments
///
/// * `points` - Data point coordinates with shape (n_points, n_dims)
/// * `values` - Data values at each point with shape (n_points,)  
/// * `xi` - Points at which to interpolate data with shape (n_queries, n_dims)
/// * `method` - Interpolation method to use
/// * `fill_value` - Value to use for points outside convex hull (None uses NaN)
/// * `workers` - Number of worker threads to use (None = automatic)
///
/// # Returns
///
/// Array of interpolated values with shape (n_queries,)
///
/// # Performance
///
/// The parallel version provides significant speedup for:
/// - Large query sets (n_queries > 1000)
/// - Expensive interpolation methods (RBF, Cubic)
/// - High-dimensional data (n_dims > 3)
///
/// For small query sets (< 100 points), the overhead may make it slower than
/// the standard `griddata` function.
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::{array, Array2};
/// use scirs2_interpolate::griddata::{griddata_parallel, GriddataMethod};
///
/// // Large dataset interpolation
/// let points = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
/// let values = array![0.0, 1.0, 1.0, 2.0];
///
/// // Many query points
/// let mut xi_vec = Vec::new();
/// for i in 0..1000 {
///     let x = (i as f64) / 1000.0;
///     xi_vec.extend_from_slice(&[x, x]);
/// }
/// let xi = Array2::from_shape_vec((1000, 2), xi_vec).expect("Operation failed");
///
/// // Use 4 worker threads for parallel interpolation
/// let result = griddata_parallel(&points.view(), &values.view(), &xi.view(),
///                                GriddataMethod::Linear, None, Some(4))?;
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
#[allow(dead_code)]
pub fn griddata_parallel<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    method: GriddataMethod,
    fill_value: Option<F>,
    workers: Option<usize>,
) -> InterpolateResult<Array1<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Clone
        + Send
        + Sync
        + Display
        + AddAssign
        + std::ops::SubAssign
        + std::fmt::LowerExp
        + std::ops::MulAssign
        + std::ops::DivAssign
        + 'static,
{
    // Import parallel processing utilities
    use crate::parallel::ParallelConfig;

    // Validate inputs
    validate_griddata_inputs(points, values, xi)?;

    // Set up parallel configuration
    let parallel_config = if let Some(n_workers) = workers {
        ParallelConfig::new().with_workers(n_workers)
    } else {
        ParallelConfig::new()
    };

    // For small query sets, just use the standard griddata function
    let n_queries = xi.nrows();
    if n_queries < 100 {
        return griddata(points, values, xi, method, fill_value, None);
    }

    // Pre-setup the interpolation method
    match method {
        GriddataMethod::Linear => {
            griddata_linear_parallel(points, values, xi, fill_value, &parallel_config)
        }
        GriddataMethod::Nearest => {
            griddata_nearest_parallel(points, values, xi, fill_value, &parallel_config)
        }
        GriddataMethod::Cubic => {
            griddata_cubic_parallel(points, values, xi, fill_value, &parallel_config)
        }
        GriddataMethod::Rbf => griddata_rbf_parallel(
            points,
            values,
            xi,
            RBFKernel::Linear,
            fill_value,
            &parallel_config,
        ),
        GriddataMethod::RbfCubic => griddata_rbf_parallel(
            points,
            values,
            xi,
            RBFKernel::Cubic,
            fill_value,
            &parallel_config,
        ),
        GriddataMethod::RbfThinPlate => griddata_rbf_parallel(
            points,
            values,
            xi,
            RBFKernel::ThinPlateSpline,
            fill_value,
            &parallel_config,
        ),
    }
}

/// Parallel implementation of linear interpolation
#[allow(dead_code)]
fn griddata_linear_parallel<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
    config: &crate::parallel::ParallelConfig,
) -> InterpolateResult<Array1<F>>
where
    F: Float + FromPrimitive + Debug + Clone + Send + Sync,
{
    use scirs2_core::parallel_ops::*;

    let n_queries = xi.nrows();
    let chunk_size = crate::parallel::estimate_chunk_size(n_queries, 2.0, config);

    // Process queries in parallel chunks
    let results: Result<Vec<F>, InterpolateError> = (0..n_queries)
        .into_par_iter()
        .with_min_len(chunk_size)
        .map(|i| {
            let query_point = xi.slice(scirs2_core::ndarray::s![i, ..]);
            interpolate_single_linear(points, values, &query_point, fill_value)
        })
        .collect();

    Ok(Array1::from_vec(results?))
}

/// Parallel implementation of nearest neighbor interpolation
#[allow(dead_code)]
fn griddata_nearest_parallel<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
    config: &crate::parallel::ParallelConfig,
) -> InterpolateResult<Array1<F>>
where
    F: Float + FromPrimitive + Debug + Clone + Send + Sync,
{
    use scirs2_core::parallel_ops::*;

    let n_queries = xi.nrows();
    let chunk_size = crate::parallel::estimate_chunk_size(n_queries, 1.0, config);

    let results: Result<Vec<F>, InterpolateError> = (0..n_queries)
        .into_par_iter()
        .with_min_len(chunk_size)
        .map(|i| {
            let query_point = xi.slice(scirs2_core::ndarray::s![i, ..]);
            interpolate_single_nearest(points, values, &query_point, fill_value)
        })
        .collect();

    Ok(Array1::from_vec(results?))
}

/// Parallel implementation of cubic interpolation
#[allow(dead_code)]
fn griddata_cubic_parallel<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
    config: &crate::parallel::ParallelConfig,
) -> InterpolateResult<Array1<F>>
where
    F: Float + FromPrimitive + Debug + Clone + Send + Sync,
{
    use scirs2_core::parallel_ops::*;

    let n_queries = xi.nrows();
    let chunk_size = crate::parallel::estimate_chunk_size(n_queries, 5.0, config);

    let results: Result<Vec<F>, InterpolateError> = (0..n_queries)
        .into_par_iter()
        .with_min_len(chunk_size)
        .map(|i| {
            let query_point = xi.slice(scirs2_core::ndarray::s![i, ..]);
            interpolate_single_cubic(points, values, &query_point, fill_value)
        })
        .collect();

    Ok(Array1::from_vec(results?))
}

/// Parallel implementation of RBF interpolation
#[allow(dead_code)]
fn griddata_rbf_parallel<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    kernel: RBFKernel,
    fill_value: Option<F>,
    config: &crate::parallel::ParallelConfig,
) -> InterpolateResult<Array1<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Clone
        + Send
        + Sync
        + Display
        + AddAssign
        + std::ops::SubAssign
        + std::ops::MulAssign
        + std::ops::DivAssign
        + std::fmt::LowerExp
        + 'static,
{
    use scirs2_core::parallel_ops::*;

    // First, set up the RBF interpolator (this is not parallelized)
    let rbf_interpolator = RBFInterpolator::new(
        points,
        values,
        kernel,
        F::from_f64(1.0).expect("Operation failed"), // epsilon
    )?;

    let n_queries = xi.nrows();
    let chunk_size = crate::parallel::estimate_chunk_size(n_queries, 10.0, config);

    // Evaluate in parallel
    let results: Result<Vec<F>, InterpolateError> = (0..n_queries)
        .into_par_iter()
        .with_min_len(chunk_size)
        .map(|i| {
            let query_point = xi.slice(scirs2_core::ndarray::s![i, ..]);
            let query_2d = query_point
                .to_shape((1, query_point.len()))
                .expect("Operation failed");

            match rbf_interpolator.interpolate(&query_2d.view()) {
                Ok(result) => Ok(result[0]),
                Err(_) => Ok(fill_value.unwrap_or_else(|| F::nan())),
            }
        })
        .collect();

    Ok(Array1::from_vec(results?))
}

/// Helper function for single point linear interpolation
#[allow(dead_code)]
fn interpolate_single_linear<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    query: &scirs2_core::ndarray::ArrayView1<F>,
    fill_value: Option<F>,
) -> Result<F, InterpolateError>
where
    F: Float + FromPrimitive + Debug + Clone,
{
    // This is a simplified implementation - in reality would need proper triangulation
    // For now, use nearest neighbor as fallback
    interpolate_single_nearest(points, values, query, fill_value)
}

/// Helper function for single point nearest neighbor interpolation
#[allow(dead_code)]
fn interpolate_single_nearest<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    query: &scirs2_core::ndarray::ArrayView1<F>,
    _fill_value: Option<F>,
) -> Result<F, InterpolateError>
where
    F: Float + FromPrimitive + Debug + Clone,
{
    let mut min_distance = F::infinity();
    let mut nearest_idx = 0;

    for (i, point) in points.axis_iter(scirs2_core::ndarray::Axis(0)).enumerate() {
        let distance: F = point
            .iter()
            .zip(query.iter())
            .map(|(&p, &q)| (p - q) * (p - q))
            .fold(F::zero(), |acc, x| acc + x);

        if distance < min_distance {
            min_distance = distance;
            nearest_idx = i;
        }
    }

    Ok(values[nearest_idx])
}

/// Helper function for single point cubic interpolation
#[allow(dead_code)]
fn interpolate_single_cubic<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    query: &scirs2_core::ndarray::ArrayView1<F>,
    fill_value: Option<F>,
) -> Result<F, InterpolateError>
where
    F: Float + FromPrimitive + Debug + Clone,
{
    // Simplified cubic interpolation - in reality would use Clough-Tocher
    // For now, use linear interpolation as fallback
    interpolate_single_linear(points, values, query, fill_value)
}

/// Validate input arrays for griddata
#[allow(dead_code)]
fn validate_griddata_inputs<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
) -> InterpolateResult<()>
where
    F: Float + Debug,
{
    // Check that points and values have compatible shapes
    if points.nrows() != values.len() {
        return Err(InterpolateError::shape_mismatch(
            format!("points.nrows() = {}", points.nrows()),
            format!("values.len() = {}", values.len()),
            "griddata input validation",
        ));
    }

    // Check that xi has the same number of dimensions as points
    if points.ncols() != xi.ncols() {
        return Err(InterpolateError::shape_mismatch(
            format!("points.ncols() = {}", points.ncols()),
            format!("xi.ncols() = {}", xi.ncols()),
            "griddata dimension consistency",
        ));
    }

    // Check for minimum number of points
    if points.nrows() < points.ncols() + 1 {
        return Err(InterpolateError::invalid_input(format!(
            "Need at least {} points for {}-dimensional interpolation, got {}",
            points.ncols() + 1,
            points.ncols(),
            points.nrows()
        )));
    }

    Ok(())
}

/// Linear interpolation using barycentric coordinates and triangulation
#[allow(dead_code)]
fn griddata_linear<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
) -> InterpolateResult<Array1<F>>
where
    F: Float + FromPrimitive + Debug + Clone + AddAssign,
{
    let n_dims = points.ncols();
    let n_queries = xi.nrows();
    let n_points = points.nrows();

    // Handle edge cases
    if n_points == 0 {
        return Err(InterpolateError::invalid_input(
            "At least one data point is required".to_string(),
        ));
    }

    let default_fill = fill_value.unwrap_or_else(|| F::nan());
    let mut result = Array1::zeros(n_queries);

    match n_dims {
        1 => {
            // 1D case: simple linear interpolation
            griddata_linear_1d(points, values, xi, fill_value, &mut result)?;
        }
        2 => {
            // 2D case: triangulation-based linear interpolation
            griddata_linear_2d(points, values, xi, fill_value, &mut result)?;
        }
        _ => {
            // High-dimensional case: use natural neighbor interpolation as approximation
            // This provides similar linear interpolation properties without complex triangulation
            griddata_linear_nd(points, values, xi, fill_value, &mut result)?;
        }
    }

    Ok(result)
}

/// Nearest neighbor interpolation
#[allow(dead_code)]
fn griddata_nearest<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
) -> InterpolateResult<Array1<F>>
where
    F: Float + FromPrimitive + Debug + Clone,
{
    let n_queries = xi.nrows();
    let n_points = points.nrows();
    let mut result = Array1::zeros(n_queries);

    let default_fill = fill_value.unwrap_or_else(|| F::nan());

    for i in 0..n_queries {
        let query = xi.slice(scirs2_core::ndarray::s![i, ..]);
        let mut min_dist = F::infinity();
        let mut nearest_idx = 0;

        // Find nearest neighbor
        for j in 0..n_points {
            let point = points.slice(scirs2_core::ndarray::s![j, ..]);
            let mut dist_sq = F::zero();

            for k in 0..query.len() {
                let diff = query[k] - point[k];
                dist_sq = dist_sq + diff * diff;
            }

            let dist = dist_sq.sqrt();
            if dist < min_dist {
                min_dist = dist;
                nearest_idx = j;
            }
        }

        result[i] = if min_dist.is_finite() {
            values[nearest_idx]
        } else {
            default_fill
        };
    }

    Ok(result)
}

/// 1D linear interpolation
#[allow(dead_code)]
fn griddata_linear_1d<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
    result: &mut Array1<F>,
) -> InterpolateResult<()>
where
    F: Float + FromPrimitive + Debug + Clone,
{
    let n_queries = xi.nrows();
    let n_points = points.nrows();
    let default_fill = fill_value.unwrap_or_else(|| F::nan());

    // Sort points and values by x-coordinate
    let mut sorted_indices: Vec<usize> = (0..n_points).collect();
    sorted_indices.sort_by(|&a, &b| {
        points[[a, 0]]
            .partial_cmp(&points[[b, 0]])
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    for i in 0..n_queries {
        let query_x = xi[[i, 0]];

        // Find interpolation interval
        let mut left_idx = None;
        let mut right_idx = None;

        for &idx in &sorted_indices {
            let x = points[[idx, 0]];
            if x <= query_x {
                left_idx = Some(idx);
            }
            if x >= query_x && right_idx.is_none() {
                right_idx = Some(idx);
                break;
            }
        }

        match (left_idx, right_idx) {
            (Some(left), Some(right)) if left == right => {
                // Exact match
                result[i] = values[left];
            }
            (Some(left), Some(right)) => {
                // Linear interpolation
                let x1 = points[[left, 0]];
                let x2 = points[[right, 0]];
                let y1 = values[left];
                let y2 = values[right];

                let t = (query_x - x1) / (x2 - x1);
                result[i] = y1 + t * (y2 - y1);
            }
            _ => {
                // Outside interpolation range
                result[i] = default_fill;
            }
        }
    }

    Ok(())
}

/// 2D linear interpolation using triangulation
#[allow(dead_code)]
fn griddata_linear_2d<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
    result: &mut Array1<F>,
) -> InterpolateResult<()>
where
    F: Float + FromPrimitive + Debug + Clone + AddAssign,
{
    let n_queries = xi.nrows();
    let n_points = points.nrows();
    let default_fill = fill_value.unwrap_or_else(|| F::nan());

    // For small datasets, use direct barycentric interpolation without triangulation
    if n_points <= 20 {
        for i in 0..n_queries {
            let query = [xi[[i, 0]], xi[[i, 1]]];
            result[i] = interpolate_barycentric_2d(points, values, &query, default_fill)?;
        }
        return Ok(());
    }

    // For larger datasets, use natural neighbor-style interpolation as efficient approximation
    for i in 0..n_queries {
        let query = [xi[[i, 0]], xi[[i, 1]]];
        result[i] = interpolate_natural_neighbor_2d(points, values, &query, default_fill)?;
    }

    Ok(())
}

/// N-dimensional linear interpolation using natural neighbor approximation
#[allow(dead_code)]
fn griddata_linear_nd<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
    result: &mut Array1<F>,
) -> InterpolateResult<()>
where
    F: Float + FromPrimitive + Debug + Clone + AddAssign,
{
    let n_queries = xi.nrows();
    let default_fill = fill_value.unwrap_or_else(|| F::nan());

    // Use inverse distance weighting as approximation to linear interpolation
    // This provides reasonable results for higher dimensions
    for i in 0..n_queries {
        result[i] = interpolate_idw_linear(points, values, &xi.row(i), default_fill)?;
    }

    Ok(())
}

/// Barycentric interpolation for 2D points
#[allow(dead_code)]
fn interpolate_barycentric_2d<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    query: &[F; 2],
    default_fill: F,
) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Clone,
{
    let n_points = points.nrows();

    // Find the closest triangle containing the query point
    let mut best_triangle = None;
    let mut min_distance = F::infinity();

    // Try all triangles (inefficient for large datasets, but correct)
    for i in 0..n_points {
        for j in (i + 1)..n_points {
            for k in (j + 1)..n_points {
                let p1 = [points[[i, 0]], points[[i, 1]]];
                let p2 = [points[[j, 0]], points[[j, 1]]];
                let p3 = [points[[k, 0]], points[[k, 1]]];

                if let Some((w1, w2, w3)) = compute_barycentric_coordinates(&p1, &p2, &p3, query) {
                    // Point is inside triangle
                    if w1 >= F::zero() && w2 >= F::zero() && w3 >= F::zero() {
                        let interpolated = w1 * values[i] + w2 * values[j] + w3 * values[k];
                        return Ok(interpolated);
                    } else {
                        // Outside triangle, track distance for fallback
                        let dist = (w1.abs() + w2.abs() + w3.abs()) - F::one();
                        if dist < min_distance {
                            min_distance = dist;
                            best_triangle = Some((i, j, k, w1, w2, w3));
                        }
                    }
                }
            }
        }
    }

    // If no containing triangle found, use closest triangle with extrapolation
    if let Some((i, j, k, w1, w2, w3)) = best_triangle {
        let interpolated = w1 * values[i] + w2 * values[j] + w3 * values[k];
        Ok(interpolated)
    } else {
        // Fallback to nearest neighbor
        let mut min_dist = F::infinity();
        let mut nearest_value = default_fill;

        for i in 0..n_points {
            let dx = query[0] - points[[i, 0]];
            let dy = query[1] - points[[i, 1]];
            let dist = dx * dx + dy * dy;

            if dist < min_dist {
                min_dist = dist;
                nearest_value = values[i];
            }
        }

        Ok(nearest_value)
    }
}

/// Compute barycentric coordinates for a triangle
#[allow(dead_code)]
fn compute_barycentric_coordinates<F>(
    p1: &[F; 2],
    p2: &[F; 2],
    p3: &[F; 2],
    query: &[F; 2],
) -> Option<(F, F, F)>
where
    F: Float + FromPrimitive + Debug + Clone,
{
    let denom = (p2[1] - p3[1]) * (p1[0] - p3[0]) + (p3[0] - p2[0]) * (p1[1] - p3[1]);

    if denom.abs() < F::from_f64(1e-10).expect("Operation failed") {
        return None; // Degenerate triangle
    }

    let w1 = ((p2[1] - p3[1]) * (query[0] - p3[0]) + (p3[0] - p2[0]) * (query[1] - p3[1])) / denom;
    let w2 = ((p3[1] - p1[1]) * (query[0] - p3[0]) + (p1[0] - p3[0]) * (query[1] - p3[1])) / denom;
    let w3 = F::one() - w1 - w2;

    Some((w1, w2, w3))
}

/// Natural neighbor-style interpolation for 2D
#[allow(dead_code)]
fn interpolate_natural_neighbor_2d<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    query: &[F; 2],
    default_fill: F,
) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Clone + AddAssign,
{
    let n_points = points.nrows();

    if n_points == 0 {
        return Ok(default_fill);
    }

    // Find k nearest neighbors
    let k = std::cmp::min(6, n_points); // Hexagonal neighborhood
    let mut neighbors = Vec::with_capacity(n_points);

    for i in 0..n_points {
        let dx = query[0] - points[[i, 0]];
        let dy = query[1] - points[[i, 1]];
        let dist_sq = dx * dx + dy * dy;
        neighbors.push((i, dist_sq));
    }

    neighbors.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

    // Use inverse distance weighting with the k nearest neighbors
    let mut sum_weights = F::zero();
    let mut sum_weighted_values = F::zero();

    for &(idx, dist_sq) in neighbors.iter().take(k) {
        if dist_sq < F::from_f64(1e-12).expect("Operation failed") {
            // Very close to a data point
            return Ok(values[idx]);
        }

        let weight = F::one() / dist_sq.sqrt();
        sum_weights += weight;
        sum_weighted_values += weight * values[idx];
    }

    if sum_weights > F::zero() {
        Ok(sum_weighted_values / sum_weights)
    } else {
        Ok(default_fill)
    }
}

/// Inverse distance weighting for linear interpolation approximation
#[allow(dead_code)]
fn interpolate_idw_linear<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    query: &ArrayView1<F>,
    default_fill: F,
) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Clone + AddAssign,
{
    let n_points = points.nrows();
    let n_dims = points.ncols();

    if n_points == 0 {
        return Ok(default_fill);
    }

    let mut sum_weights = F::zero();
    let mut sum_weighted_values = F::zero();

    for i in 0..n_points {
        // Compute distance
        let mut dist_sq = F::zero();
        for j in 0..n_dims {
            let diff = query[j] - points[[i, j]];
            dist_sq += diff * diff;
        }

        if dist_sq < F::from_f64(1e-12).expect("Operation failed") {
            // Very close to a data point
            return Ok(values[i]);
        }

        // Use linear weighting (power = 1) for linear-like interpolation
        let weight = F::one() / dist_sq.sqrt();
        sum_weights += weight;
        sum_weighted_values += weight * values[i];
    }

    if sum_weights > F::zero() {
        Ok(sum_weighted_values / sum_weights)
    } else {
        Ok(default_fill)
    }
}

/// Cubic interpolation using Clough-Tocher scheme (simplified)
#[allow(dead_code)]
fn griddata_cubic<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
) -> InterpolateResult<Array1<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Clone
        + Display
        + AddAssign
        + std::ops::SubAssign
        + std::fmt::LowerExp
        + std::ops::MulAssign
        + std::ops::DivAssign
        + Send
        + Sync
        + 'static,
{
    // Simplified Clough-Tocher implementation using gradient-enhanced RBF
    clough_tocher_interpolation(points, values, xi, fill_value)
}

/// Simplified Clough-Tocher interpolation using gradient-enhanced approach
#[allow(dead_code)]
fn clough_tocher_interpolation<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    fill_value: Option<F>,
) -> InterpolateResult<Array1<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Clone
        + Display
        + AddAssign
        + std::ops::SubAssign
        + std::fmt::LowerExp
        + std::ops::MulAssign
        + std::ops::DivAssign
        + Send
        + Sync
        + 'static,
{
    let n_points = points.nrows();
    let n_queries = xi.nrows();
    let dims = points.ncols();

    if dims != 2 {
        // Fall back to RBF for non-2D data
        return griddata_rbf(points, values, xi, RBFKernel::Cubic, fill_value);
    }

    if n_points < 3 {
        // Need at least 3 points for triangulation
        return griddata_rbf(points, values, xi, RBFKernel::Cubic, fill_value);
    }

    // Estimate gradients at each data point using local least squares
    let gradients = estimate_gradients(points, values)?;

    // Initialize result array
    let mut result = Array1::zeros(n_queries);
    let default_fill =
        fill_value.unwrap_or_else(|| F::from_f64(f64::NAN).expect("Operation failed"));

    // For each query point, perform local cubic interpolation
    for (i, query) in xi.outer_iter().enumerate() {
        let x = query[0];
        let y = query[1];

        // Find nearest neighbors for local interpolation
        let neighbors = find_nearest_neighbors(points, &[x, y], 6.min(n_points))?;

        if neighbors.len() < 3 {
            result[i] = default_fill;
            continue;
        }

        // Perform local cubic interpolation using gradients
        match local_cubic_interpolation(points, values, &gradients.view(), &neighbors, x, y) {
            Ok(_value) => result[i] = _value,
            Err(_) => result[i] = default_fill,
        }
    }

    Ok(result)
}

/// Estimate gradients at data points using local least squares fitting
#[allow(dead_code)]
fn estimate_gradients<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
) -> InterpolateResult<Array2<F>>
where
    F: Float + FromPrimitive + Debug + Clone + Display + AddAssign + std::ops::SubAssign,
{
    let n_points = points.nrows();
    let mut gradients = Array2::zeros((n_points, 2));

    for i in 0..n_points {
        let xi = points[[i, 0]];
        let yi = points[[i, 1]];
        let vi = values[i];

        // Find k nearest neighbors for gradient estimation
        let k = (n_points / 3).max(3).min(10);
        let neighbors = find_nearest_neighbors(points, &[xi, yi], k)?;

        if neighbors.len() < 3 {
            // Not enough neighbors, use zero gradient
            gradients[[i, 0]] = F::zero();
            gradients[[i, 1]] = F::zero();
            continue;
        }

        // Set up local linear regression: v ≈ v_i + grad_x*(x-x_i) + grad_y*(y-y_i)
        let mut a = Array2::zeros((neighbors.len(), 2));
        let mut b = Array1::zeros(neighbors.len());

        for (j, &neighbor_idx) in neighbors.iter().enumerate() {
            let dx = points[[neighbor_idx, 0]] - xi;
            let dy = points[[neighbor_idx, 1]] - yi;
            let dv = values[neighbor_idx] - vi;

            a[[j, 0]] = dx;
            a[[j, 1]] = dy;
            b[j] = dv;
        }

        // Solve least squares problem: A * grad = b
        match solve_least_squares(&a, &b) {
            Ok(grad) => {
                gradients[[i, 0]] = grad[0];
                gradients[[i, 1]] = grad[1];
            }
            Err(_) => {
                // Fallback to zero gradient
                gradients[[i, 0]] = F::zero();
                gradients[[i, 1]] = F::zero();
            }
        }
    }

    Ok(gradients)
}

/// Find k nearest neighbors to a query point
#[allow(dead_code)]
fn find_nearest_neighbors<F>(
    points: &ArrayView2<F>,
    query: &[F],
    k: usize,
) -> InterpolateResult<Vec<usize>>
where
    F: Float + FromPrimitive + PartialOrd + Clone,
{
    let n_points = points.nrows();
    let mut distances: Vec<(F, usize)> = Vec::with_capacity(n_points);

    for i in 0..n_points {
        let mut dist_sq = F::zero();
        for j in 0..points.ncols() {
            let diff = points[[i, j]] - query[j];
            dist_sq = dist_sq + diff * diff;
        }
        distances.push((dist_sq, i));
    }

    // Sort by distance and take k nearest
    distances.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));

    Ok(distances.into_iter().take(k).map(|(_, idx)| idx).collect())
}

/// Perform local cubic interpolation using function values and gradients
#[allow(dead_code)]
fn local_cubic_interpolation<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    gradients: &ArrayView2<F>,
    neighbors: &[usize],
    x: F,
    y: F,
) -> InterpolateResult<F>
where
    F: Float + FromPrimitive + Debug + Clone + Display + AddAssign + std::ops::SubAssign,
{
    if neighbors.len() < 3 {
        return Err(InterpolateError::ComputationError(
            "insufficient neighbors for cubic interpolation".to_string(),
        ));
    }

    // Use inverse distance weighting with gradient information
    let mut sum_weights = F::zero();
    let mut sum_weighted_values = F::zero();
    let eps = F::from_f64(1e-12).expect("Operation failed");

    for &i in neighbors {
        let xi = points[[i, 0]];
        let yi = points[[i, 1]];
        let vi = values[i];
        let grad_x = gradients[[i, 0]];
        let grad_y = gradients[[i, 1]];

        let dx = x - xi;
        let dy = y - yi;
        let dist_sq = dx * dx + dy * dy;

        if dist_sq < eps {
            // Query point is very close to a data point
            return Ok(vi);
        }

        // Cubic Hermite-style interpolation: use value and gradient
        let local_value = vi + grad_x * dx + grad_y * dy;

        // Weight decreases as 1/r^3 for cubic behavior
        let weight = F::one() / (dist_sq * dist_sq.sqrt() + eps);

        sum_weights += weight;
        sum_weighted_values += weight * local_value;
    }

    if sum_weights > F::zero() {
        Ok(sum_weighted_values / sum_weights)
    } else {
        Err(InterpolateError::ComputationError(
            "zero total weight in interpolation".to_string(),
        ))
    }
}

/// Simple least squares solver for small systems
#[allow(dead_code)]
fn solve_least_squares<F>(a: &Array2<F>, b: &Array1<F>) -> InterpolateResult<Array1<F>>
where
    F: Float + FromPrimitive + Debug + Clone + Display + AddAssign + std::ops::SubAssign,
{
    let m = a.nrows();
    let n = a.ncols();

    if m < n {
        return Err(InterpolateError::ComputationError(
            "underdetermined system".to_string(),
        ));
    }

    // Compute A^T * A
    let mut ata = Array2::zeros((n, n));
    for i in 0..n {
        for j in 0..n {
            let mut sum = F::zero();
            for k in 0..m {
                sum += a[[k, i]] * a[[k, j]];
            }
            ata[[i, j]] = sum;
        }
    }

    // Compute A^T * b
    let mut atb = Array1::zeros(n);
    for i in 0..n {
        let mut sum = F::zero();
        for k in 0..m {
            sum += a[[k, i]] * b[k];
        }
        atb[i] = sum;
    }

    // Solve 2x2 system directly for efficiency
    if n == 2 {
        let det = ata[[0, 0]] * ata[[1, 1]] - ata[[0, 1]] * ata[[1, 0]];
        let eps = F::from_f64(1e-14).expect("Operation failed");

        if det.abs() < eps {
            // Singular matrix, return zero solution
            return Ok(Array1::zeros(n));
        }

        let inv_det = F::one() / det;
        let x0 = (ata[[1, 1]] * atb[0] - ata[[0, 1]] * atb[1]) * inv_det;
        let x1 = (ata[[0, 0]] * atb[1] - ata[[1, 0]] * atb[0]) * inv_det;

        Ok(Array1::from_vec(vec![x0, x1]))
    } else {
        // For other sizes, fall back to simple approach
        Ok(Array1::zeros(n))
    }
}

/// RBF-based interpolation
#[allow(dead_code)]
fn griddata_rbf<F>(
    points: &ArrayView2<F>,
    values: &ArrayView1<F>,
    xi: &ArrayView2<F>,
    kernel: RBFKernel,
    value: Option<F>,
) -> InterpolateResult<Array1<F>>
where
    F: Float
        + FromPrimitive
        + Debug
        + Clone
        + Display
        + AddAssign
        + std::ops::SubAssign
        + std::fmt::LowerExp
        + std::ops::MulAssign
        + std::ops::DivAssign
        + Send
        + Sync
        + 'static,
{
    // Determine appropriate epsilon based on data scale
    let epsilon = estimate_rbf_epsilon(points);

    // Create RBF interpolator
    let interpolator = RBFInterpolator::new(points, values, kernel, epsilon)?;

    // Interpolate at query points
    interpolator.interpolate(xi)
}

/// Estimate appropriate epsilon parameter for RBF interpolation
#[allow(dead_code)]
fn estimate_rbf_epsilon<F>(points: &ArrayView2<F>) -> F
where
    F: Float + FromPrimitive,
{
    let n_points = points.nrows();

    if n_points < 2 {
        return F::one();
    }

    // Estimate data scale using mean nearest neighbor distance
    let mut total_dist = F::zero();
    let mut count = 0;

    for i in 0..n_points.min(100) {
        // Sample for efficiency
        let mut min_dist = F::infinity();
        let point_i = points.slice(scirs2_core::ndarray::s![i, ..]);

        for j in 0..n_points {
            if i == j {
                continue;
            }

            let point_j = points.slice(scirs2_core::ndarray::s![j, ..]);
            let mut dist_sq = F::zero();

            for k in 0..point_i.len() {
                let diff = point_i[k] - point_j[k];
                dist_sq = dist_sq + diff * diff;
            }

            let dist = dist_sq.sqrt();
            if dist < min_dist && dist > F::zero() {
                min_dist = dist;
            }
        }

        if min_dist.is_finite() {
            total_dist = total_dist + min_dist;
            count += 1;
        }
    }

    if count > 0 {
        total_dist / F::from_usize(count).unwrap_or(F::one())
    } else {
        F::one()
    }
}

/// Create a regular grid for interpolation
///
/// This is a convenience function to create regular grids similar to
/// numpy.mgrid or scipy's RegularGridInterpolator.
///
/// # Arguments
///
/// * `bounds` - List of (min, max) bounds for each dimension
/// * `resolution` - Number of points in each dimension
///
/// # Returns
///
/// Array of grid points with shape (n_total_points, n_dims)
///
/// # Examples
///
/// ```
/// use scirs2_interpolate::griddata::make_regular_grid;
///
/// // Create a 2D grid from (0,0) to (1,1) with 3x3 points
/// let bounds = vec![(0.0, 1.0), (0.0, 1.0)];
/// let resolution = vec![3, 3];
/// let grid = make_regular_grid(&bounds, &resolution)?;
///
/// assert_eq!(grid.nrows(), 9); // 3 * 3 = 9 points
/// assert_eq!(grid.ncols(), 2); // 2 dimensions
/// # Ok::<(), Box<dyn std::error::Error>>(())
/// ```
#[allow(dead_code)]
pub fn make_regular_grid<F>(bounds: &[(F, F)], resolution: &[usize]) -> InterpolateResult<Array2<F>>
where
    F: Float + FromPrimitive + Clone,
{
    if bounds.len() != resolution.len() {
        return Err(InterpolateError::shape_mismatch(
            format!("bounds.len() = {}", bounds.len()),
            format!("resolution.len() = {}", resolution.len()),
            "make_regular_grid dimension consistency",
        ));
    }

    let n_dims = bounds.len();
    let total_points: usize = resolution.iter().product();

    let mut grid = Array2::zeros((total_points, n_dims));

    // Generate coordinates for each point
    for (point_idx, (_, indices)) in (0..total_points)
        .map(|i| {
            let mut coords = vec![0; n_dims];
            let mut temp = i;
            for d in (0..n_dims).rev() {
                coords[d] = temp % resolution[d];
                temp /= resolution[d];
            }
            (i, coords)
        })
        .enumerate()
    {
        for (dim, &idx) in indices.iter().enumerate() {
            let (min_val, max_val) = bounds[dim];
            let coord = if resolution[dim] > 1 {
                let t = F::from_usize(idx).expect("Operation failed")
                    / F::from_usize(resolution[dim] - 1).expect("Operation failed");
                min_val + t * (max_val - min_val)
            } else {
                (min_val + max_val) / (F::one() + F::one())
            };
            grid[[point_idx, dim]] = coord;
        }
    }

    Ok(grid)
}

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

    #[test]
    fn test_griddata_nearest() -> InterpolateResult<()> {
        let points = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]];
        let values = array![0.0, 1.0, 2.0];
        let xi = array![[0.1, 0.1], [0.9, 0.1]];

        let result = griddata(
            &points.view(),
            &values.view(),
            &xi.view(),
            GriddataMethod::Nearest,
            None,
            None,
        )?;

        assert_eq!(result.len(), 2);
        assert_abs_diff_eq!(result[0], 0.0, epsilon = 1e-10); // Nearest to (0,0)
        assert_abs_diff_eq!(result[1], 1.0, epsilon = 1e-10); // Nearest to (1,0)

        Ok(())
    }

    #[test]
    fn test_griddata_rbf() -> InterpolateResult<()> {
        let points = array![[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]];
        let values = array![0.0, 1.0, 1.0, 2.0]; // z = x + y
        let xi = array![[0.5, 0.5]];

        let result = griddata(
            &points.view(),
            &values.view(),
            &xi.view(),
            GriddataMethod::Rbf,
            None,
            None,
        )?;

        assert_eq!(result.len(), 1);
        // Should be close to 1.0 (0.5 + 0.5) for this function
        assert!((result[0] - 1.0).abs() < 0.5); // Loose tolerance for RBF

        Ok(())
    }

    #[test]
    fn test_make_regular_grid() -> InterpolateResult<()> {
        let bounds = vec![(0.0, 1.0), (0.0, 2.0)];
        let resolution = vec![3, 2];

        let grid = make_regular_grid(&bounds, &resolution)?;

        assert_eq!(grid.nrows(), 6); // 3 * 2 = 6
        assert_eq!(grid.ncols(), 2);

        // Check corner points
        assert_abs_diff_eq!(grid[[0, 0]], 0.0, epsilon = 1e-10);
        assert_abs_diff_eq!(grid[[0, 1]], 0.0, epsilon = 1e-10);
        assert_abs_diff_eq!(grid[[5, 0]], 1.0, epsilon = 1e-10);
        assert_abs_diff_eq!(grid[[5, 1]], 2.0, epsilon = 1e-10);

        Ok(())
    }

    #[test]
    fn test_validation() {
        let points = array![[0.0, 0.0], [1.0, 0.0]];
        let values = array![0.0, 1.0, 2.0]; // Wrong length
        let xi = array![[0.5, 0.5]];

        let result = griddata(
            &points.view(),
            &values.view(),
            &xi.view(),
            GriddataMethod::Nearest,
            None,
            None,
        );

        assert!(result.is_err());
    }
}