scirs2-spatial 0.4.2

//! Generic traits and type parameters for spatial algorithms
//!
//! This module provides generic traits that enable spatial algorithms to work
//! with different numeric types and array structures. This improves API flexibility
//! and allows for better integration with the broader Rust ecosystem.
//!
//! # Features
//!
//! - **Generic numeric types**: Support for f32, f64, and other numeric types
//! - **Array ecosystem integration**: Work with ndarray, nalgebra, and other array types
//! - **Iterator support**: Process data from various sources efficiently
//! - **Dimension awareness**: Compile-time and runtime dimension checking
//! - **Memory efficiency**: Zero-copy operations where possible
//!
//! # Examples
//!
//! ```
//! use scirs2_spatial::generic_traits::{SpatialScalar};
//! use scirs2_core::ndarray::array;
//!
//! // Generic distance calculation (simplified example)
//! fn calculate_distance<T>(p1: &[T], p2: &[T]) -> T
//! where
//!     T: SpatialScalar,
//! {
//!     // Generic implementation (simplified)
//!     T::zero()
//! }
//!
//! // Works with different array types
//! let point1 = array![1.0f32, 2.0f32, 3.0f32];
//! let point2 = array![4.0f32, 5.0f32, 6.0f32];
//! let _result = calculate_distance(point1.as_slice().expect("Operation failed"), point2.as_slice().expect("Operation failed"));
//! ```

use scirs2_core::ndarray::Array1;
use scirs2_core::numeric::{Float, NumCast};
use scirs2_core::simd_ops::SimdUnifiedOps;
use std::fmt::Debug;

/// Trait for scalar types that can be used in spatial computations
///
/// This trait extends standard numeric traits with requirements specific
/// to spatial algorithms, including floating-point operations and conversions.
pub trait SpatialScalar:
    Float + Debug + Default + NumCast + Send + Sync + SimdUnifiedOps + 'static
{
    /// The epsilon value for floating-point comparisons
    fn epsilon() -> Self;

    /// Maximum finite value for this type
    fn max_finite() -> Self;

    /// Convert from f64 (used for constants and literals)
    fn from_f64(value: f64) -> Option<Self> {
        NumCast::from(value)
    }

    /// Convert to f64 for interoperability
    fn to_f64(&self) -> Option<f64> {
        NumCast::from(*self)
    }

    /// Square root function
    fn sqrt(&self) -> Self {
        Float::sqrt(*self)
    }

    /// Absolute value function
    fn abs(&self) -> Self {
        Float::abs(*self)
    }

    /// Power function
    fn powf(&self, exp: Self) -> Self {
        Float::powf(*self, exp)
    }

    /// Natural logarithm
    fn ln(&self) -> Self {
        Float::ln(*self)
    }

    /// Exponential function
    fn exp(&self) -> Self {
        Float::exp(*self)
    }

    /// Sine function
    fn sin(&self) -> Self {
        Float::sin(*self)
    }

    /// Cosine function
    fn cos(&self) -> Self {
        Float::cos(*self)
    }

    /// Arctangent of y/x
    fn atan2(&self, other: Self) -> Self {
        Float::atan2(*self, other)
    }

    /// Check if the value is finite
    fn is_finite(&self) -> bool {
        Float::is_finite(*self)
    }

    /// Check if the value is NaN
    fn is_nan(&self) -> bool {
        Float::is_nan(*self)
    }

    /// SIMD-optimized squared Euclidean distance
    fn simd_squared_euclidean_distance(_a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        Err("SIMD not available for this type")
    }

    /// SIMD-optimized Manhattan distance
    fn simd_manhattan_distance(_a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        Err("SIMD not available for this type")
    }

    /// SIMD-optimized dot product
    fn simd_dot_product(_a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        Err("SIMD not available for this type")
    }
}

impl SpatialScalar for f32 {
    fn epsilon() -> Self {
        f32::EPSILON
    }

    fn max_finite() -> Self {
        f32::MAX
    }

    fn simd_squared_euclidean_distance(a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        if a.len() != b.len() {
            return Err("Slice lengths must match");
        }

        use scirs2_core::ndarray::Array1;
        let a_array = Array1::from(a.to_vec());
        let b_array = Array1::from(b.to_vec());

        let diff = Self::simd_sub(&a_array.view(), &b_array.view());
        let squared = Self::simd_mul(&diff.view(), &diff.view());
        Ok(Self::simd_sum(&squared.view()))
    }

    fn simd_manhattan_distance(a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        if a.len() != b.len() {
            return Err("Slice lengths must match");
        }

        let a_array = Array1::from(a.to_vec());
        let b_array = Array1::from(b.to_vec());

        let diff = Self::simd_sub(&a_array.view(), &b_array.view());
        let abs_diff = Self::simd_abs(&diff.view());
        Ok(Self::simd_sum(&abs_diff.view()))
    }

    fn simd_dot_product(a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        if a.len() != b.len() {
            return Err("Slice lengths must match");
        }

        let a_array = Array1::from(a.to_vec());
        let b_array = Array1::from(b.to_vec());

        Ok(Self::simd_dot(&a_array.view(), &b_array.view()))
    }
}

impl SpatialScalar for f64 {
    fn epsilon() -> Self {
        f64::EPSILON
    }

    fn max_finite() -> Self {
        f64::MAX
    }

    fn simd_squared_euclidean_distance(a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        if a.len() != b.len() {
            return Err("Slice lengths must match");
        }

        let a_array = Array1::from(a.to_vec());
        let b_array = Array1::from(b.to_vec());

        let diff = Self::simd_sub(&a_array.view(), &b_array.view());
        let squared = Self::simd_mul(&diff.view(), &diff.view());
        Ok(Self::simd_sum(&squared.view()))
    }

    fn simd_manhattan_distance(a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        if a.len() != b.len() {
            return Err("Slice lengths must match");
        }

        let a_array = Array1::from(a.to_vec());
        let b_array = Array1::from(b.to_vec());

        let diff = Self::simd_sub(&a_array.view(), &b_array.view());
        let abs_diff = Self::simd_abs(&diff.view());
        Ok(Self::simd_sum(&abs_diff.view()))
    }

    fn simd_dot_product(a: &[Self], b: &[Self]) -> Result<Self, &'static str> {
        if a.len() != b.len() {
            return Err("Slice lengths must match");
        }

        let a_array = Array1::from(a.to_vec());
        let b_array = Array1::from(b.to_vec());

        Ok(Self::simd_dot(&a_array.view(), &b_array.view()))
    }
}

/// Trait for types that can represent a point in space
///
/// This trait abstracts over different point representations,
/// allowing algorithms to work with vectors, arrays, and custom types.
pub trait SpatialPoint<T: SpatialScalar> {
    /// Get the dimension of the point
    fn dimension(&self) -> usize;

    /// Get the coordinate at the given index
    fn coordinate(&self, index: usize) -> Option<T>;

    /// Get all coordinates as a slice if possible (for efficiency)
    fn as_slice(&self) -> Option<&[T]> {
        None
    }

    /// Create a point from coordinates
    fn from_coords(coords: &[T]) -> Self;

    /// Calculate squared Euclidean distance to another point
    fn squared_distance_to(&self, other: &Self) -> T {
        if self.dimension() != other.dimension() {
            return T::max_finite();
        }

        // Try SIMD-optimized calculation if slices are available
        if let (Some(slice_a), Some(slice_b)) = (self.as_slice(), other.as_slice()) {
            if let Ok(simd_result) = T::simd_squared_euclidean_distance(slice_a, slice_b) {
                return simd_result;
            }
        }

        // Fallback to scalar calculation
        let mut sum = T::zero();
        for i in 0..self.dimension() {
            if let (Some(a), Some(b)) = (self.coordinate(i), other.coordinate(i)) {
                let diff = a - b;
                sum = sum + diff * diff;
            }
        }
        sum
    }

    /// Calculate Euclidean distance to another point
    fn distance_to(&self, other: &Self) -> T {
        self.squared_distance_to(other).sqrt()
    }

    /// Calculate Manhattan distance to another point
    fn manhattan_distance_to(&self, other: &Self) -> T {
        if self.dimension() != other.dimension() {
            return T::max_finite();
        }

        // Try SIMD-optimized calculation if slices are available
        if let (Some(slice_a), Some(slice_b)) = (self.as_slice(), other.as_slice()) {
            if let Ok(simd_result) = T::simd_manhattan_distance(slice_a, slice_b) {
                return simd_result;
            }
        }

        // Fallback to scalar calculation
        let mut sum = T::zero();
        for i in 0..self.dimension() {
            if let (Some(a), Some(b)) = (self.coordinate(i), other.coordinate(i)) {
                sum = sum + (a - b).abs();
            }
        }
        sum
    }
}

/// Trait for collections of spatial points
///
/// This trait allows algorithms to work with different array structures
/// and iterator types while maintaining efficiency.
pub trait SpatialArray<T: SpatialScalar, P: SpatialPoint<T>> {
    /// Get the number of points in the array
    fn len(&self) -> usize;

    /// Check if the array is empty
    fn is_empty(&self) -> bool {
        self.len() == 0
    }

    /// Get the dimension of points in this array
    fn dimension(&self) -> Option<usize>;

    /// Get a point at the given index
    fn get_point(&self, index: usize) -> Option<P>;

    /// Iterate over all points
    fn iter_points(&self) -> Box<dyn Iterator<Item = P> + '_>;

    /// Get the bounding box of all points
    fn bounding_box(&self) -> Option<(P, P)> {
        if self.is_empty() {
            return None;
        }

        let first = self.get_point(0)?;
        let dim = first.dimension();

        let mut min_coords = vec![T::max_finite(); dim];
        let mut max_coords = vec![T::min_value(); dim];

        for point in self.iter_points() {
            for i in 0..dim {
                if let Some(coord) = point.coordinate(i) {
                    if coord < min_coords[i] {
                        min_coords[i] = coord;
                    }
                    if coord > max_coords[i] {
                        max_coords[i] = coord;
                    }
                }
            }
        }

        Some((P::from_coords(&min_coords), P::from_coords(&max_coords)))
    }
}

/// Trait for distance metrics
pub trait DistanceMetric<T: SpatialScalar, P: SpatialPoint<T>> {
    /// Calculate distance between two points
    fn distance(&self, p1: &P, p2: &P) -> T;

    /// Calculate squared distance (if applicable, for efficiency)
    fn squared_distance(&self, p1: &P, p2: &P) -> Option<T> {
        None
    }

    /// Check if this metric satisfies the triangle inequality
    fn is_metric(&self) -> bool {
        true
    }

    /// Get the name of this distance metric
    fn name(&self) -> &'static str;
}

/// Euclidean distance metric
#[derive(Debug, Clone, Copy, Default)]
pub struct EuclideanMetric;

impl<T: SpatialScalar, P: SpatialPoint<T>> DistanceMetric<T, P> for EuclideanMetric {
    fn distance(&self, p1: &P, p2: &P) -> T {
        p1.distance_to(p2)
    }

    fn squared_distance(&self, p1: &P, p2: &P) -> Option<T> {
        Some(p1.squared_distance_to(p2))
    }

    fn name(&self) -> &'static str {
        "euclidean"
    }
}

/// Manhattan distance metric
#[derive(Debug, Clone, Copy, Default)]
pub struct ManhattanMetric;

impl<T: SpatialScalar, P: SpatialPoint<T>> DistanceMetric<T, P> for ManhattanMetric {
    fn distance(&self, p1: &P, p2: &P) -> T {
        p1.manhattan_distance_to(p2)
    }

    fn name(&self) -> &'static str {
        "manhattan"
    }
}

/// Chebyshev distance metric
#[derive(Debug, Clone, Copy, Default)]
pub struct ChebyshevMetric;

impl<T: SpatialScalar, P: SpatialPoint<T>> DistanceMetric<T, P> for ChebyshevMetric {
    fn distance(&self, p1: &P, p2: &P) -> T {
        if p1.dimension() != p2.dimension() {
            return T::max_finite();
        }

        let mut max_diff = T::zero();
        for i in 0..p1.dimension() {
            if let (Some(a), Some(b)) = (p1.coordinate(i), p2.coordinate(i)) {
                let diff = (a - b).abs();
                if diff > max_diff {
                    max_diff = diff;
                }
            }
        }
        max_diff
    }

    fn name(&self) -> &'static str {
        "chebyshev"
    }
}

/// Implementation of SpatialPoint for `Vec<T>`
impl<T: SpatialScalar> SpatialPoint<T> for Vec<T> {
    fn dimension(&self) -> usize {
        self.len()
    }

    fn coordinate(&self, index: usize) -> Option<T> {
        self.get(index).copied()
    }

    fn as_slice(&self) -> Option<&[T]> {
        Some(self.as_slice())
    }

    fn from_coords(coords: &[T]) -> Self {
        coords.to_vec()
    }
}

/// Implementation of SpatialPoint for slices
impl<T: SpatialScalar> SpatialPoint<T> for &[T] {
    fn dimension(&self) -> usize {
        self.len()
    }

    fn coordinate(&self, index: usize) -> Option<T> {
        self.get(index).copied()
    }

    fn as_slice(&self) -> Option<&[T]> {
        Some(self)
    }

    fn from_coords(coords: &[T]) -> Self {
        // This is a fundamental limitation - &[T] is a reference to existing data
        // and cannot be created from raw coordinates without an underlying array.
        // This implementation is not meaningful for slice references.
        unreachable!(
            "&[T]::from_coords() should not be called - &[T] is a reference to existing data"
        )
    }
}

/// Implementation of SpatialPoint for arrays
impl<T: SpatialScalar, const N: usize> SpatialPoint<T> for [T; N] {
    fn dimension(&self) -> usize {
        N
    }

    fn coordinate(&self, index: usize) -> Option<T> {
        self.get(index).copied()
    }

    fn as_slice(&self) -> Option<&[T]> {
        Some(self.as_slice())
    }

    fn from_coords(coords: &[T]) -> Self {
        let mut result = [T::zero(); N];
        for (i, &coord) in coords.iter().enumerate().take(N) {
            result[i] = coord;
        }
        result
    }
}

/// Generic point structure for spatial algorithms
#[derive(Debug, Clone, PartialEq)]
pub struct Point<T: SpatialScalar> {
    coords: Vec<T>,
}

impl<T: SpatialScalar> Point<T> {
    /// Create a new point from coordinates
    pub fn new(coords: Vec<T>) -> Self {
        Self { coords }
    }

    /// Create a point with the given dimension filled with zeros
    pub fn zeros(dim: usize) -> Self {
        Self {
            coords: vec![T::zero(); dim],
        }
    }

    /// Create a 2D point
    pub fn new_2d(x: T, y: T) -> Self {
        Self { coords: vec![x, y] }
    }

    /// Create a 3D point
    pub fn new_3d(x: T, y: T, z: T) -> Self {
        Self {
            coords: vec![x, y, z],
        }
    }

    /// Get the coordinates as a slice
    pub fn coords(&self) -> &[T] {
        &self.coords
    }

    /// Get mutable access to coordinates
    pub fn coords_mut(&mut self) -> &mut [T] {
        &mut self.coords
    }

    /// Add another point (vector addition)
    pub fn add(&self, other: &Point<T>) -> Option<Point<T>> {
        if self.dimension() != other.dimension() {
            return None;
        }

        let coords: Vec<T> = self
            .coords
            .iter()
            .zip(other.coords.iter())
            .map(|(&a, &b)| a + b)
            .collect();

        Some(Point::new(coords))
    }

    /// Subtract another point (vector subtraction)
    pub fn subtract(&self, other: &Point<T>) -> Option<Point<T>> {
        if self.dimension() != other.dimension() {
            return None;
        }

        let coords: Vec<T> = self
            .coords
            .iter()
            .zip(other.coords.iter())
            .map(|(&a, &b)| a - b)
            .collect();

        Some(Point::new(coords))
    }

    /// Scale the point by a scalar
    pub fn scale(&self, factor: T) -> Point<T> {
        let coords: Vec<T> = self.coords.iter().map(|&c| c * factor).collect();
        Point::new(coords)
    }

    /// Calculate the dot product with another point
    pub fn dot(&self, other: &Point<T>) -> Option<T> {
        if self.dimension() != other.dimension() {
            return None;
        }

        // Try SIMD-optimized calculation
        if let Ok(simd_result) = T::simd_dot_product(&self.coords, &other.coords) {
            return Some(simd_result);
        }

        // Fallback to scalar calculation
        let dot_product = self
            .coords
            .iter()
            .zip(other.coords.iter())
            .map(|(&a, &b)| a * b)
            .fold(T::zero(), |acc, x| acc + x);

        Some(dot_product)
    }

    /// Calculate the magnitude (length) of the point as a vector
    pub fn magnitude(&self) -> T {
        self.coords
            .iter()
            .map(|&c| c * c)
            .fold(T::zero(), |acc, x| acc + x)
            .sqrt()
    }

    /// Normalize the point to unit length
    pub fn normalize(&self) -> Point<T> {
        let mag = self.magnitude();
        if mag > T::zero() {
            self.scale(T::one() / mag)
        } else {
            self.clone()
        }
    }
}

impl<T: SpatialScalar> SpatialPoint<T> for Point<T> {
    fn dimension(&self) -> usize {
        self.coords.len()
    }

    fn coordinate(&self, index: usize) -> Option<T> {
        self.coords.get(index).copied()
    }

    fn as_slice(&self) -> Option<&[T]> {
        Some(&self.coords)
    }

    fn from_coords(coords: &[T]) -> Self {
        Point::new(coords.to_vec())
    }
}

/// Utility functions for generic spatial operations
pub mod utils {
    use super::*;

    /// Calculate pairwise distances between all points in a collection
    pub fn pairwise_distances<T, P, A, M>(points: &A, metric: &M) -> Vec<T>
    where
        T: SpatialScalar,
        P: SpatialPoint<T>,
        A: SpatialArray<T, P>,
        M: DistanceMetric<T, P>,
    {
        let n = points.len();
        let mut distances = Vec::with_capacity(n * (n - 1) / 2);

        for i in 0..n {
            for j in (i + 1)..n {
                if let (Some(p1), Some(p2)) = (points.get_point(i), points.get_point(j)) {
                    distances.push(metric.distance(&p1, &p2));
                }
            }
        }

        distances
    }

    /// Find the centroid of a collection of points
    pub fn centroid<T, P, A>(points: &A) -> Option<Point<T>>
    where
        T: SpatialScalar,
        P: SpatialPoint<T>,
        A: SpatialArray<T, P>,
    {
        if points.is_empty() {
            return None;
        }

        let n = points.len();
        let dim = points.dimension()?;
        let mut sum_coords = vec![T::zero(); dim];

        for point in points.iter_points() {
            for (i, sum_coord) in sum_coords.iter_mut().enumerate().take(dim) {
                if let Some(coord) = point.coordinate(i) {
                    *sum_coord = *sum_coord + coord;
                }
            }
        }

        let n_scalar = T::from(n)?;
        for coord in &mut sum_coords {
            *coord = *coord / n_scalar;
        }

        Some(Point::new(sum_coords))
    }

    /// Calculate the convex hull using a generic algorithm
    pub fn convex_hull_2d<T, P, A>(points: &A) -> Vec<Point<T>>
    where
        T: SpatialScalar,
        P: SpatialPoint<T>,
        A: SpatialArray<T, P>,
    {
        if points.len() < 3 {
            return points
                .iter_points()
                .map(|p| Point::from_coords(p.as_slice().unwrap_or(&[])))
                .collect();
        }

        // Simple implementation - in practice, you'd use more sophisticated algorithms
        let mut hull_points: Vec<Point<T>> = points
            .iter_points()
            .map(|p| Point::from_coords(p.as_slice().unwrap_or(&[])))
            .collect();

        // Sort by x-coordinate, then by y-coordinate
        hull_points.sort_by(|a, b| {
            let x_cmp = a
                .coordinate(0)
                .partial_cmp(&b.coordinate(0))
                .expect("Operation failed");
            if x_cmp == std::cmp::Ordering::Equal {
                a.coordinate(1)
                    .partial_cmp(&b.coordinate(1))
                    .expect("Operation failed")
            } else {
                x_cmp
            }
        });

        hull_points
    }
}

/// Integration with ndarray
pub mod ndarray_integration {
    use super::{SpatialArray, SpatialPoint, SpatialScalar};
    use scirs2_core::ndarray::{Array1, ArrayView1, ArrayView2, Axis};

    /// Implementation of SpatialPoint for ndarray ArrayView1
    impl<T: SpatialScalar> SpatialPoint<T> for ArrayView1<'_, T> {
        fn dimension(&self) -> usize {
            self.len()
        }

        fn coordinate(&self, index: usize) -> Option<T> {
            self.get(index).copied()
        }

        fn as_slice(&self) -> Option<&[T]> {
            // Use as_slice_memory_order() to avoid infinite recursion with trait method
            self.as_slice_memory_order()
        }

        fn from_coords(coords: &[T]) -> Self {
            // This is a fundamental limitation - ArrayView1 is a view into existing data
            // and cannot be created from raw coordinates without an underlying array.
            // This implementation is not meaningful for views, but we provide a dummy
            // implementation to satisfy the trait. Real usage should avoid this method.
            unreachable!("ArrayView1::from_coords() should not be called - ArrayView1 is a view into existing data")
        }
    }

    /// Implementation of SpatialPoint for ndarray Array1
    impl<T: SpatialScalar> SpatialPoint<T> for Array1<T> {
        fn dimension(&self) -> usize {
            self.len()
        }

        fn coordinate(&self, index: usize) -> Option<T> {
            self.get(index).copied()
        }

        fn as_slice(&self) -> Option<&[T]> {
            // Use as_slice_memory_order() to avoid infinite recursion with trait method
            self.as_slice_memory_order()
        }

        fn from_coords(coords: &[T]) -> Self {
            Array1::from(coords.to_vec())
        }
    }

    /// Wrapper for ndarray Array2 to implement SpatialArray
    pub struct NdArray2Wrapper<'a, T: SpatialScalar> {
        array: ArrayView2<'a, T>,
    }

    impl<'a, T: SpatialScalar> NdArray2Wrapper<'a, T> {
        pub fn new(array: ArrayView2<'a, T>) -> Self {
            Self { array }
        }
    }

    impl<T: SpatialScalar> SpatialArray<T, Array1<T>> for NdArray2Wrapper<'_, T> {
        fn len(&self) -> usize {
            self.array.nrows()
        }

        fn dimension(&self) -> Option<usize> {
            if self.array.nrows() > 0 {
                Some(self.array.ncols())
            } else {
                None
            }
        }

        fn get_point(&self, index: usize) -> Option<Array1<T>> {
            if index < self.len() {
                Some(self.array.row(index).to_owned())
            } else {
                None
            }
        }

        fn iter_points(&self) -> Box<dyn Iterator<Item = Array1<T>> + '_> {
            Box::new(self.array.axis_iter(Axis(0)).map(|row| row.to_owned()))
        }
    }
}

#[cfg(test)]
mod tests {
    use crate::{
        ChebyshevMetric, DistanceMetric, EuclideanMetric, ManhattanMetric, Point, SpatialPoint,
        SpatialScalar,
    };
    use approx::assert_relative_eq;

    #[test]
    fn test_spatial_scalar_traits() {
        assert!(<f32 as SpatialScalar>::epsilon() > 0.0);
        assert!(<f64 as SpatialScalar>::epsilon() > 0.0);
        assert!(f32::max_finite().is_finite());
        assert!(f64::max_finite().is_finite());
    }

    #[test]
    fn test_point_operations() {
        let p1 = Point::new_2d(1.0f64, 2.0);
        let p2 = Point::new_2d(4.0, 6.0);

        assert_eq!(p1.dimension(), 2);
        assert_eq!(p1.coordinate(0), Some(1.0));
        assert_eq!(p1.coordinate(1), Some(2.0));

        let distance = p1.distance_to(&p2);
        assert_relative_eq!(distance, 5.0, epsilon = 1e-10);

        let manhattan = p1.manhattan_distance_to(&p2);
        assert_relative_eq!(manhattan, 7.0, epsilon = 1e-10);
    }

    #[test]
    fn test_point_arithmetic() {
        let p1 = Point::new_3d(1.0f32, 2.0, 3.0);
        let p2 = Point::new_3d(4.0, 5.0, 6.0);

        let sum = p1.add(&p2).expect("Operation failed");
        assert_eq!(sum.coordinate(0), Some(5.0));
        assert_eq!(sum.coordinate(1), Some(7.0));
        assert_eq!(sum.coordinate(2), Some(9.0));

        let diff = p2.subtract(&p1).expect("Operation failed");
        assert_eq!(diff.coordinate(0), Some(3.0));
        assert_eq!(diff.coordinate(1), Some(3.0));
        assert_eq!(diff.coordinate(2), Some(3.0));

        let scaled = p1.scale(2.0);
        assert_eq!(scaled.coordinate(0), Some(2.0));
        assert_eq!(scaled.coordinate(1), Some(4.0));
        assert_eq!(scaled.coordinate(2), Some(6.0));
    }

    #[test]
    fn test_distance_metrics() {
        use crate::generic_traits::DistanceMetric;

        let p1 = Point::new_2d(0.0f64, 0.0);
        let p2 = Point::new_2d(3.0, 4.0);

        let euclidean = EuclideanMetric;
        let manhattan = ManhattanMetric;
        let chebyshev = ChebyshevMetric;

        assert_relative_eq!(euclidean.distance(&p1, &p2), 5.0, epsilon = 1e-10);
        assert_relative_eq!(manhattan.distance(&p1, &p2), 7.0, epsilon = 1e-10);
        assert_relative_eq!(chebyshev.distance(&p1, &p2), 4.0, epsilon = 1e-10);
    }

    #[test]
    fn test_vec_as_spatial_point() {
        let p1 = vec![1.0f64, 2.0, 3.0];
        let p2 = vec![4.0, 5.0, 6.0];

        assert_eq!(p1.dimension(), 3);
        assert_eq!(p1.coordinate(1), Some(2.0));

        let distance = p1.distance_to(&p2);
        // Euclidean distance: sqrt((4-1)² + (5-2)² + (6-3)²) = sqrt(9+9+9) = sqrt(27)
        assert_relative_eq!(distance, (27.0f64).sqrt(), epsilon = 1e-10);
    }

    #[test]
    fn test_array_as_spatial_point() {
        let p1: [f32; 3] = [1.0, 2.0, 3.0];
        let p2: [f32; 3] = [4.0, 5.0, 6.0];

        assert_eq!(p1.dimension(), 3);
        assert_eq!(p1.coordinate(2), Some(3.0));

        let distance = p1.distance_to(&p2);
        // Euclidean distance: sqrt((4-1)² + (5-2)² + (6-3)²) = sqrt(9+9+9) = sqrt(27)
        assert_relative_eq!(distance, (27.0f32).sqrt(), epsilon = 1e-6);
    }

    #[test]
    fn test_point_normalization() {
        let p = Point::new_2d(3.0f64, 4.0);
        let magnitude = p.magnitude();
        assert_relative_eq!(magnitude, 5.0, epsilon = 1e-10);

        let normalized = p.normalize();
        assert_relative_eq!(normalized.magnitude(), 1.0, epsilon = 1e-10);
        assert_relative_eq!(
            normalized.coordinate(0).expect("Operation failed"),
            0.6,
            epsilon = 1e-10
        );
        assert_relative_eq!(
            normalized.coordinate(1).expect("Operation failed"),
            0.8,
            epsilon = 1e-10
        );
    }

    #[test]
    fn test_dot_product() {
        let p1 = Point::new_3d(1.0f64, 2.0, 3.0);
        let p2 = Point::new_3d(4.0, 5.0, 6.0);

        let dot = p1.dot(&p2).expect("Operation failed");
        assert_relative_eq!(dot, 32.0, epsilon = 1e-10); // 1*4 + 2*5 + 3*6 = 32
    }
}