--- a/src/nalgebra_support.rs Tue Feb 20 12:33:16 2024 -0500
+++ b/src/nalgebra_support.rs Mon Feb 03 19:22:16 2025 -0500
@@ -1,7 +1,7 @@
/*!
Integration with nalgebra.
-This module mainly implements [`Euclidean`], [`Norm`], [`Dot`], [`Linear`], etc. for [`nalgebra`]
+This module mainly implements [`Euclidean`], [`Norm`], [`Linear`], etc. for [`nalgebra`]
matrices and vectors.
It also provides [`ToNalgebraRealField`] as a vomit-inducingly ugly workaround to nalgebra
force-feeding its own versions of the same basic mathematical methods on `f32` and `f64` as
@@ -10,10 +10,9 @@
use nalgebra::{
Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar,
- ClosedMul, ClosedAdd, SimdComplexField, Vector, OVector, RealField,
+ ClosedAddAssign, ClosedMulAssign, SimdComplexField, Vector, OVector, RealField,
LpNorm, UniformNorm
};
-use nalgebra::Norm as NalgebraNorm;
use nalgebra::base::constraint::{
ShapeConstraint, SameNumberOfRows, SameNumberOfColumns
};
@@ -23,102 +22,127 @@
use num_traits::identities::{Zero, One};
use crate::linops::*;
use crate::euclidean::*;
+use crate::mapping::{Space, BasicDecomposition};
use crate::types::Float;
use crate::norms::*;
+use crate::instance::Instance;
-impl<SM,SV,N,M,K,E> Apply<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV: Storage<E,M,K>,
- N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One,
- DefaultAllocator : Allocator<E,N,K>,
- DefaultAllocator : Allocator<E,M,K>,
- DefaultAllocator : Allocator<E,N,M>,
- DefaultAllocator : Allocator<E,M,N> {
- type Output = OMatrix<E,N,K>;
+impl<SM,N,M,E> Space for Matrix<E,N,M,SM>
+where
+ SM: Storage<E,N,M> + Clone,
+ N : Dim, M : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
+ DefaultAllocator : Allocator<N,M>,
+{
+ type Decomp = BasicDecomposition;
+}
+
+impl<SM,SV,N,M,K,E> Mapping<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
+where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone,
+ N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
+ DefaultAllocator : Allocator<N,K>,
+ DefaultAllocator : Allocator<M,K>,
+ DefaultAllocator : Allocator<N,M>,
+ DefaultAllocator : Allocator<M,N> {
+ type Codomain = OMatrix<E,N,K>;
#[inline]
- fn apply(&self, x : Matrix<E,M,K,SV>) -> Self::Output {
- self.mul(x)
+ fn apply<I : Instance<Matrix<E,M,K,SV>>>(
+ &self, x : I
+ ) -> Self::Codomain {
+ x.either(|owned| self.mul(owned), |refr| self.mul(refr))
}
}
-impl<'a, SM,SV,N,M,K,E> Apply<&'a Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV: Storage<E,M,K>,
- N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One,
- DefaultAllocator : Allocator<E,N,K>,
- DefaultAllocator : Allocator<E,M,K>,
- DefaultAllocator : Allocator<E,N,M>,
- DefaultAllocator : Allocator<E,M,N> {
- type Output = OMatrix<E,N,K>;
-
- #[inline]
- fn apply(&self, x : &'a Matrix<E,M,K,SV>) -> Self::Output {
- self.mul(x)
- }
-}
impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV: Storage<E,M,K>,
- N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One,
- DefaultAllocator : Allocator<E,N,K>,
- DefaultAllocator : Allocator<E,M,K>,
- DefaultAllocator : Allocator<E,N,M>,
- DefaultAllocator : Allocator<E,M,N> {
- type Codomain = OMatrix<E,N,K>;
+where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone,
+ N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
+ DefaultAllocator : Allocator<N,K>,
+ DefaultAllocator : Allocator<M,K>,
+ DefaultAllocator : Allocator<N,M>,
+ DefaultAllocator : Allocator<M,N> {
}
impl<SM,SV1,SV2,N,M,K,E> GEMV<E, Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: StorageMut<E,N,K>,
- N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float,
- DefaultAllocator : Allocator<E,N,K>,
- DefaultAllocator : Allocator<E,M,K>,
- DefaultAllocator : Allocator<E,N,M>,
- DefaultAllocator : Allocator<E,M,N> {
+where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: StorageMut<E,N,K>,
+ N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float,
+ DefaultAllocator : Allocator<N,K>,
+ DefaultAllocator : Allocator<M,K>,
+ DefaultAllocator : Allocator<N,M>,
+ DefaultAllocator : Allocator<M,N> {
#[inline]
- fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : &Matrix<E,M,K,SV1>, β : E) {
- Matrix::gemm(y, α, self, x, β)
+ fn gemv<I : Instance<Matrix<E,M,K,SV1>>>(
+ &self, y : &mut Matrix<E,N,K,SV2>, α : E, x : I, β : E
+ ) {
+ x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β))
}
#[inline]
- fn apply_mut<'a>(&self, y : &mut Matrix<E,N,K,SV2>, x : &Matrix<E,M,K,SV1>) {
- self.mul_to(x, y)
+ fn apply_mut<'a, I : Instance<Matrix<E,M,K,SV1>>>(&self, y : &mut Matrix<E,N,K,SV2>, x : I) {
+ x.eval(|x̃| self.mul_to(x̃, y))
}
}
impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM>
-where SM: StorageMut<E,M>, SV1: Storage<E,M>,
- M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float,
- DefaultAllocator : Allocator<E,M> {
+where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone,
+ M : Dim, E : Scalar + Zero + One + Float,
+ DefaultAllocator : Allocator<M> {
+ type Owned = OVector<E, M>;
#[inline]
- fn axpy(&mut self, α : E, x : &Vector<E,M,SV1>, β : E) {
- Matrix::axpy(self, α, x, β)
+ fn axpy<I : Instance<Vector<E,M,SV1>>>(&mut self, α : E, x : I, β : E) {
+ x.eval(|x̃| Matrix::axpy(self, α, x̃, β))
+ }
+
+ #[inline]
+ fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) {
+ y.eval(|ỹ| Matrix::copy_from(self, ỹ))
+ }
+
+ #[inline]
+ fn set_zero(&mut self) {
+ self.iter_mut().for_each(|e| *e = E::ZERO);
}
#[inline]
- fn copy_from(&mut self, y : &Vector<E,M,SV1>) {
- Matrix::copy_from(self, y)
+ fn similar_origin(&self) -> Self::Owned {
+ OVector::zeros_generic(M::from_usize(self.len()), Const)
}
}
+/* Implemented automatically as Euclidean.
+impl<SM,M,E> Projection<E, L2> for Vector<E,M,SM>
+where SM: StorageMut<E,M> + Clone,
+ M : Dim, E : Scalar + Zero + One + Float + RealField,
+ DefaultAllocator : Allocator<M> {
+ #[inline]
+ fn proj_ball_mut(&mut self, ρ : E, _ : L2) {
+ let n = self.norm(L2);
+ if n > ρ {
+ self.iter_mut().for_each(|v| *v *= ρ/n)
+ }
+ }
+}*/
+
impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM>
-where SM: StorageMut<E,M>,
- M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float + RealField,
- DefaultAllocator : Allocator<E,M> {
+where SM: StorageMut<E,M> + Clone,
+ M : Dim, E : Scalar + Zero + One + Float + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) {
self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ))
}
}
-impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>,Matrix<E,N,K,SV2>>
+impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>>
for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: Storage<E,N,K>,
- N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + SimdComplexField,
- DefaultAllocator : Allocator<E,N,K>,
- DefaultAllocator : Allocator<E,M,K>,
- DefaultAllocator : Allocator<E,N,M>,
- DefaultAllocator : Allocator<E,M,N> {
+where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: Storage<E,N,K> + Clone,
+ N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + SimdComplexField,
+ DefaultAllocator : Allocator<N,K>,
+ DefaultAllocator : Allocator<M,K>,
+ DefaultAllocator : Allocator<N,M>,
+ DefaultAllocator : Allocator<M,N> {
type AdjointCodomain = OMatrix<E,M,K>;
type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a;
@@ -128,20 +152,6 @@
}
}
-impl<E,M,S,Si> Dot<Vector<E,M,Si>,E>
-for Vector<E,M,S>
-where M : Dim,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One,
- S : Storage<E,M>,
- Si : Storage<E,M>,
- DefaultAllocator : Allocator<E,M> {
-
- #[inline]
- fn dot(&self, other : &Vector<E,M,Si>) -> E {
- Vector::<E,M,S>::dot(self, other)
- }
-}
-
/// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`.
#[inline]
fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>(
@@ -170,15 +180,15 @@
impl<E,M,S> Euclidean<E>
for Vector<E,M,S>
where M : Dim,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : StorageMut<E,M> + Clone,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
type Output = OVector<E, M>;
-
+
#[inline]
- fn similar_origin(&self) -> OVector<E, M> {
- OVector::zeros_generic(M::from_usize(self.len()), Const)
+ fn dot<I : Instance<Self>>(&self, other : I) -> E {
+ Vector::<E,M,S>::dot(self, other.ref_instance())
}
#[inline]
@@ -187,17 +197,17 @@
}
#[inline]
- fn dist2_squared(&self, other : &Self) -> E {
- metric_distance_squared(self, other)
+ fn dist2_squared<I : Instance<Self>>(&self, other : I) -> E {
+ metric_distance_squared(self, other.ref_instance())
}
}
impl<E,M,S> StaticEuclidean<E>
for Vector<E,M,S>
where M : DimName,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : StorageMut<E,M> + Clone,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
fn origin() -> OVector<E, M> {
@@ -205,78 +215,109 @@
}
}
+/// The default norm for `Vector` is [`L2`].
+impl<E,M,S> Normed<E>
+for Vector<E,M,S>
+where M : Dim,
+ S : Storage<E,M> + Clone,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
+
+ type NormExp = L2;
+
+ #[inline]
+ fn norm_exponent(&self) -> Self::NormExp {
+ L2
+ }
+
+ #[inline]
+ fn is_zero(&self) -> bool {
+ Vector::<E,M,S>::norm_squared(self) == E::ZERO
+ }
+}
+
+impl<E,M,S> HasDual<E>
+for Vector<E,M,S>
+where M : Dim,
+ S : Storage<E,M> + Clone,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
+ // TODO: Doesn't work with different storage formats.
+ type DualSpace = Self;
+}
+
impl<E,M,S> Norm<E, L1>
for Vector<E,M,S>
where M : Dim,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : Storage<E,M>,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
fn norm(&self, _ : L1) -> E {
- LpNorm(1).norm(self)
+ nalgebra::Norm::norm(&LpNorm(1), self)
}
}
impl<E,M,S> Dist<E, L1>
for Vector<E,M,S>
where M : Dim,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : Storage<E,M> + Clone,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
- fn dist(&self, other : &Self, _ : L1) -> E {
- LpNorm(1).metric_distance(self, other)
+ fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> E {
+ nalgebra::Norm::metric_distance(&LpNorm(1), self, other.ref_instance())
}
}
impl<E,M,S> Norm<E, L2>
for Vector<E,M,S>
where M : Dim,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : Storage<E,M>,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
fn norm(&self, _ : L2) -> E {
- LpNorm(2).norm(self)
+ nalgebra::Norm::norm(&LpNorm(2), self)
}
}
impl<E,M,S> Dist<E, L2>
for Vector<E,M,S>
where M : Dim,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : Storage<E,M> + Clone,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
- fn dist(&self, other : &Self, _ : L2) -> E {
- LpNorm(2).metric_distance(self, other)
+ fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> E {
+ nalgebra::Norm::metric_distance(&LpNorm(2), self, other.ref_instance())
}
}
impl<E,M,S> Norm<E, Linfinity>
for Vector<E,M,S>
where M : Dim,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : Storage<E,M>,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
fn norm(&self, _ : Linfinity) -> E {
- UniformNorm.norm(self)
+ nalgebra::Norm::norm(&UniformNorm, self)
}
}
impl<E,M,S> Dist<E, Linfinity>
for Vector<E,M,S>
where M : Dim,
- S : StorageMut<E,M>,
- E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField,
- DefaultAllocator : Allocator<E,M> {
+ S : Storage<E,M> + Clone,
+ E : Float + Scalar + Zero + One + RealField,
+ DefaultAllocator : Allocator<M> {
#[inline]
- fn dist(&self, other : &Self, _ : Linfinity) -> E {
- UniformNorm.metric_distance(self, other)
+ fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> E {
+ nalgebra::Norm::metric_distance(&UniformNorm, self, other.ref_instance())
}
}
