Mercurial > repos > alg_tools / file revision

/*!
Array containers that support vector space operations on floats.
For working with small vectors in $ℝ^2$ or $ℝ^3$.
*/

use std::ops::{Add,Sub,AddAssign,SubAssign,Mul,Div,MulAssign,DivAssign,Neg,Index,IndexMut};
use std::slice::{Iter,IterMut};
use crate::types::{Float,Num,SignedNum};
use crate::maputil::{FixedLength,FixedLengthMut,map1,map2,map1_mut,map2_mut};
use crate::euclidean::*;
use crate::norms::*;
use crate::linops::AXPY;
use serde::ser::{Serialize, Serializer, SerializeSeq};

/// A container type for (short) `N`-dimensional vectors of element type `F`.
///
/// Supports basic operations of an [`Euclidean`] space, several [`Norm`]s, and
/// fused [`AXPY`] operations, among others.
#[derive(Copy,Clone,Debug,PartialEq,Eq)]
pub struct Loc<F, const N : usize>(
    /// An array of the elements of the vector
    pub [F; N]
);

// Need to manually implement as [F; N] serialisation is provided only for some N.
impl<F, const N : usize> Serialize for Loc<F, N>
where
    F: Serialize,
{
    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
    where
        S: Serializer,
    {
        let mut seq = serializer.serialize_seq(Some(N))?;
        for e in self.iter() {
            seq.serialize_element(e)?;
        }
        seq.end()
    }
}

impl<F, const N : usize> Loc<F, N> {
    /// Creates a new `Loc` vector from an array.
    #[inline]
    pub fn new(arr : [F; N]) -> Self {
        Loc(arr)
    }

    /// Returns an iterator over the elements of the vector
    #[inline]
    pub fn iter(&self) -> Iter<'_, F> {
        self.0.iter()
    }

    /// Returns an iterator over mutable references to the elements of the vector
    #[inline]
    pub fn iter_mut(&mut self) -> IterMut<'_, F> {
        self.0.iter_mut()
    }
}

impl<F : Copy, const N : usize> Loc<F, N> {
    /// Maps `g` over the elements of the vector, returning a new [`Loc`] vector
    #[inline]
    pub fn map<H>(&self, g : impl Fn(F) -> H) -> Loc<H, N> {
        Loc::new(map1(self, |u| g(*u)))
    }

    /// Maps `g` over pairs of elements of two vectors, retuning a new one.
    #[inline]
    pub fn map2<H>(&self, other : &Self, g : impl Fn(F, F) -> H) -> Loc<H, N> {
        Loc::new(map2(self, other, |u, v| g(*u, *v)))
    }

    /// Maps `g` over mutable references to elements of the vector.
    #[inline]
    pub fn map_mut(&mut self, g : impl Fn(&mut F)) {
        map1_mut(self, g)
    }

    /// Maps `g` over pairs of mutable references to elements of `self, and elements
    /// of `other` vector.
    #[inline]
    pub fn map2_mut(&mut self, other : &Self, g : impl Fn(&mut F, F)) {
        map2_mut(self, other, |u, v| g(u, *v))
    }

    /// Maps `g` over the elements of `self` and returns the product of the results.
    #[inline]
    pub fn product_map<A : Num>(&self, g : impl Fn(F) -> A) -> A {
        match N {
            1 => g(unsafe { *self.0.get_unchecked(0) }),
            2 => g(unsafe { *self.0.get_unchecked(0) }) *
                 g(unsafe { *self.0.get_unchecked(1) }),
            3 => g(unsafe { *self.0.get_unchecked(0) }) *
                 g(unsafe { *self.0.get_unchecked(1) }) *
                 g(unsafe { *self.0.get_unchecked(2) }),
            _ => self.iter().fold(A::ONE, |m, &x| m * g(x))
        }
    }
}

/// Construct a [`Loc`].
///
/// Use as
/// ```
/// # use alg_tools::loc::Loc;
/// # use alg_tools::loc;
/// let x = loc![1.0, 2.0];
/// ```
#[macro_export]
macro_rules! loc {
    ($($x:expr),+ $(,)?) => { Loc::new([$($x),+]) }
}


impl<F, const N : usize> From<[F; N]> for Loc<F, N> {
    #[inline]
    fn from(other: [F; N]) -> Loc<F, N> {
        Loc(other)
    }
}

/*impl<F : Copy, const N : usize> From<&[F; N]> for Loc<F, N> {
    #[inline]
    fn from(other: &[F; N]) -> Loc<F, N> {
        Loc(*other)
    }
}*/

impl<F> From<F> for Loc<F, 1> {
    #[inline]
    fn from(other: F) -> Loc<F, 1> {
        Loc([other])
    }
}

impl<F, const N : usize> From<Loc<F, N>> for [F; N] {
    #[inline]
    fn from(other : Loc<F, N>) -> [F; N] {
        other.0
    }
}

/*impl<F : Copy, const N : usize> From<&Loc<F, N>> for [F; N] {
    #[inline]
    fn from(other : &Loc<F, N>) -> [F; N] {
        other.0
    }
}*/


impl<F, const N : usize> IntoIterator for Loc<F, N> {
    type Item = <[F; N] as IntoIterator>::Item;
    type IntoIter = <[F; N] as IntoIterator>::IntoIter;

    #[inline]
    fn into_iter(self) -> Self::IntoIter {
        self.0.into_iter()
    }
}

// Indexing

impl<F, Ix, const N : usize> Index<Ix> for Loc<F,N>
where [F; N] : Index<Ix> {
    type Output = <[F; N] as Index<Ix>>::Output;

    #[inline]
    fn index(&self, ix : Ix) -> &Self::Output {
        self.0.index(ix)
    }
}

impl<F, Ix, const N : usize> IndexMut<Ix> for Loc<F,N>
where [F; N] : IndexMut<Ix> {
    #[inline]
    fn index_mut(&mut self, ix : Ix) -> &mut Self::Output {
        self.0.index_mut(ix)
    }
}

// Arithmetic

macro_rules! make_binop {
    ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
        impl<F : Num, const N : usize> $trait<Loc<F,N>> for Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(mut self, other : Loc<F, N>) -> Self::Output {
                self.$fn_assign(other);
                self
            }
        }

        impl<'a, F : Num, const N : usize> $trait<&'a Loc<F,N>> for Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(mut self, other : &'a Loc<F, N>) -> Self::Output {
                self.$fn_assign(other);
                self
            }
        }

        impl<'b, F : Num, const N : usize> $trait<Loc<F,N>> for &'b Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(self, other : Loc<F, N>) -> Self::Output {
                self.map2(&other, |a, b| a.$fn(b))
            }
        }

        impl<'a, 'b, F : Num, const N : usize> $trait<&'a Loc<F,N>> for &'b Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(self, other : &'a Loc<F, N>) -> Self::Output {
                self.map2(other, |a, b| a.$fn(b))
            }
        }

       impl<F : Num, const N : usize> $trait_assign<Loc<F,N>> for Loc<F, N> {
            #[inline]
            fn $fn_assign(&mut self, other : Loc<F, N>) {
                self.map2_mut(&other, |a, b| a.$fn_assign(b))
            }
        }

        impl<'a, F : Num, const N : usize> $trait_assign<&'a Loc<F,N>> for Loc<F, N> {
            #[inline]
            fn $fn_assign(&mut self, other : &'a Loc<F, N>) {
                self.map2_mut(other, |a, b| a.$fn_assign(b))
            }
        }
    }
}

make_binop!(Add, add, AddAssign, add_assign);
make_binop!(Sub, sub, SubAssign, sub_assign);

macro_rules! make_scalarop_rhs {
    ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
        impl<F : Num, const N : usize> $trait<F> for Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(self, b : F) -> Self::Output {
                self.map(|a| a.$fn(b))
            }
        }

        impl<'a, F : Num, const N : usize> $trait<&'a F> for Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(self, b : &'a F) -> Self::Output {
                self.map(|a| a.$fn(*b))
            }
        }

        impl<'b, F : Num, const N : usize> $trait<F> for &'b Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(self, b : F) -> Self::Output {
                self.map(|a| a.$fn(b))
            }
        }

        impl<'a, 'b, F : Float, const N : usize> $trait<&'a F> for &'b Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(self, b : &'a F) -> Self::Output {
                self.map(|a| a.$fn(*b))
            }
        }

        impl<F : Num, const N : usize> $trait_assign<F> for Loc<F, N> {
            #[inline]
            fn $fn_assign(&mut self, b : F) {
                self.map_mut(|a| a.$fn_assign(b));
            }
        }

        impl<'a, F : Num, const N : usize> $trait_assign<&'a F> for Loc<F, N> {
            #[inline]
            fn $fn_assign(&mut self, b : &'a F) {
                self.map_mut(|a| a.$fn_assign(*b));
            }
        }
    }
}


make_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
make_scalarop_rhs!(Div, div, DivAssign, div_assign);

macro_rules! make_unaryop {
    ($trait:ident, $fn:ident) => {
        impl<F : SignedNum, const N : usize> $trait for Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(mut self) -> Self::Output {
                self.map_mut(|a| *a = (*a).$fn());
                self
            }
        }

        impl<'a, F : SignedNum, const N : usize> $trait for &'a Loc<F, N> {
            type Output = Loc<F, N>;
            #[inline]
            fn $fn(self) -> Self::Output {
                self.map(|a| a.$fn())
            }
        }
    }
}

make_unaryop!(Neg, neg);

macro_rules! make_scalarop_lhs {
    ($trait:ident, $fn:ident; $($f:ident)+) => { $(
        impl<const N : usize> $trait<Loc<$f,N>> for $f {
            type Output = Loc<$f, N>;
            #[inline]
            fn $fn(self, v : Loc<$f,N>) -> Self::Output {
                v.map(|b| self.$fn(b))
            }
        }

        impl<'a, const N : usize> $trait<&'a Loc<$f,N>> for $f {
            type Output = Loc<$f, N>;
            #[inline]
            fn $fn(self, v : &'a Loc<$f,N>) -> Self::Output {
                v.map(|b| self.$fn(b))
            }
        }

        impl<'b, const N : usize> $trait<Loc<$f,N>> for &'b $f {
            type Output = Loc<$f, N>;
            #[inline]
            fn $fn(self, v : Loc<$f,N>) -> Self::Output {
                v.map(|b| self.$fn(b))
            }
        }

        impl<'a, 'b, const N : usize> $trait<&'a Loc<$f,N>> for &'b $f {
            type Output = Loc<$f, N>;
            #[inline]
            fn $fn(self, v : &'a Loc<$f, N>) -> Self::Output {
                v.map(|b| self.$fn(b))
            }
        }
    )+ }
}

make_scalarop_lhs!(Mul, mul; f32 f64 i8 i16 i32 i64 isize u8 u16 u32 u64 usize);
make_scalarop_lhs!(Div, div; f32 f64 i8 i16 i32 i64 isize u8 u16 u32 u64 usize);

// Norms

macro_rules! domination {
    ($norm:ident, $dominates:ident) => {
        impl<F : Float, const N : usize> Dominated<F, $dominates, Loc<F, N>> for $norm {
            #[inline]
            fn norm_factor(&self, _p : $dominates) -> F {
                F::ONE
            }
            #[inline]
            fn from_norm(&self, p_norm : F, _p : $dominates) -> F {
                p_norm
            }
        }
    };
    ($norm:ident, $dominates:ident, $fn:path) => {
        impl<F : Float, const N : usize> Dominated<F, $dominates, Loc<F, N>> for $norm {
            #[inline]
            fn norm_factor(&self, _p : $dominates) -> F {
                $fn(F::cast_from(N))
            }
        }
    };
}

domination!(L1, L1);
domination!(L2, L2);
domination!(Linfinity, Linfinity);

domination!(L1, L2, F::sqrt);
domination!(L2, Linfinity, F::sqrt);
domination!(L1, Linfinity, std::convert::identity);

domination!(Linfinity, L1);
domination!(Linfinity, L2);
domination!(L2, L1);

impl<F : Num,const N : usize> Dot<Loc<F, N>,F> for Loc<F, N> {
    /// This implementation is not stabilised as it's meant to be used for very small vectors.
    /// Use [`nalgebra`] for larger vectors.
    #[inline]
    fn dot(&self, other : &Loc<F, N>) -> F {
        self.0.iter()
              .zip(other.0.iter())
              .fold(F::ZERO, |m, (&v, &w)| m + v * w)
    }
}

impl<F : Float,const N : usize> Euclidean<F> for Loc<F, N> {
    type Output = Self;

    #[inline]
    fn similar_origin(&self) -> Self {
        Self::ORIGIN
    }

    /// This implementation is not stabilised as it's meant to be used for very small vectors.
    /// Use [`nalgebra`] for larger vectors.
    #[inline]
    fn norm2_squared(&self) -> F {
        self.iter().fold(F::ZERO, |m, &v| m + v * v)
    }

    fn dist2_squared(&self, other : &Self) -> F {
        self.iter()
            .zip(other.iter())
            .fold(F::ZERO, |m, (&v, &w)| { let d = v - w; m + d * d })
    }

    #[inline]
    fn norm2(&self) -> F {
        // Optimisation for N==1 that avoids squaring and square rooting.
        if N==1 {
            unsafe { self.0.get_unchecked(0) }.abs()
        } else {
            self.norm2_squared().sqrt()
        }
    }

    #[inline]
    fn dist2(&self, other : &Self) -> F {
        // Optimisation for N==1 that avoids squaring and square rooting.
        if N==1 {
            unsafe { *self.0.get_unchecked(0) - *other.0.get_unchecked(0) }.abs()
        } else {
            self.dist2_squared(other).sqrt()
        }
    }
}

impl<F : Num, const N : usize> Loc<F, N> {
    pub const ORIGIN : Self = Loc([F::ZERO; N]);
}

impl<F : Float,const N : usize> StaticEuclidean<F> for Loc<F, N> {
    #[inline]
    fn origin() -> Self {
        Self::ORIGIN
    }
}

impl<F : Float, const N : usize> Norm<F, L2> for Loc<F, N> {
    #[inline]
    fn norm(&self, _ : L2) -> F { self.norm2() }
}

impl<F : Float, const N : usize> Dist<F, L2> for Loc<F, N> {
    #[inline]
    fn dist(&self, other : &Self, _ : L2) -> F { self.dist2(other) }
}

impl<F : Float, const N : usize> Norm<F, L1> for Loc<F, N> {
    /// This implementation is not stabilised as it's meant to be used for very small vectors.
    /// Use [`nalgebra`] for larger vectors.
    #[inline]
    fn norm(&self, _ : L1) -> F {
        self.iter().fold(F::ZERO, |m, v| m + v.abs())
    }
}

impl<F : Float, const N : usize> Dist<F, L1> for Loc<F, N> {
    #[inline]
    fn dist(&self, other : &Self, _ : L1) -> F {
        self.iter()
            .zip(other.iter())
            .fold(F::ZERO, |m, (&v, &w)| m + (v-w).abs() )
    }
}

impl<F : Float, const N : usize> Projection<F, Linfinity> for Loc<F, N> {
    #[inline]
    fn proj_ball_mut(&mut self, ρ : F, _ : Linfinity) {
        self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ))
    }
}

impl<F : Float, const N : usize> Norm<F, Linfinity> for Loc<F, N> {
    /// This implementation is not stabilised as it's meant to be used for very small vectors.
    /// Use [`nalgebra`] for larger vectors.
    #[inline]
    fn norm(&self, _ : Linfinity) -> F {
        self.iter().fold(F::ZERO, |m, v| m.max(v.abs()))
    }
}

impl<F : Float, const N : usize> Dist<F, Linfinity> for Loc<F, N> {
    #[inline]
    fn dist(&self, other : &Self, _ : Linfinity) -> F {
        self.iter()
            .zip(other.iter())
            .fold(F::ZERO, |m, (&v, &w)| m.max((v-w).abs()))
    }
}


// Misc.

impl<A, const N : usize> FixedLength<N> for Loc<A,N> {
    type Iter = std::array::IntoIter<A, N>;
    type Elem = A;
    #[inline]
    fn fl_iter(self) -> Self::Iter {
        self.into_iter()
    }
}

impl<A, const N : usize> FixedLengthMut<N> for Loc<A,N> {
    type IterMut<'a> = std::slice::IterMut<'a, A> where A : 'a;
    #[inline]
    fn fl_iter_mut(&mut self) -> Self::IterMut<'_> {
        self.iter_mut()
    }
}

impl<'a, A, const N : usize> FixedLength<N> for &'a Loc<A,N> {
    type Iter = std::slice::Iter<'a, A>;
    type Elem = &'a A;
    #[inline]
    fn fl_iter(self) -> Self::Iter {
        self.iter()
    }
}

impl<F : Num, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> {

    #[inline]
    fn axpy(&mut self, α : F, x : &Loc<F, N>, β : F) {
        if β == F::ZERO {
            map2_mut(self, x, |yi, xi| { *yi = α * (*xi) })
        } else {
            map2_mut(self, x, |yi, xi| { *yi = β * (*yi) + α * (*xi) })
        }
    }

    #[inline]
    fn copy_from(&mut self, x : &Loc<F, N>) {
        map2_mut(self, x, |yi, xi| *yi = *xi )
    }
}