# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1756953073 18000
# Node ID 53ab61a41d7092e552208085a743d35fe90f7d4c
# Parent  102421d462d1158a054535cfcbfabfb3a31cb37d
nalgebra matrix distances

diff -r 102421d462d1 -r 53ab61a41d70 src/euclidean/wrap.rs
--- a/src/euclidean/wrap.rs	Wed Sep 03 21:03:47 2025 -0500
+++ b/src/euclidean/wrap.rs	Wed Sep 03 21:31:13 2025 -0500
@@ -135,6 +135,8 @@
         $crate::wrap!(impl_scalar_mut $F, $type, std::ops::MulAssign, mul_assign where $($qual)*);
         $crate::wrap!(impl_scalar_mut $F, $type, std::ops::DivAssign, div_assign where $($qual)*);
 
+        $crate::self_ownable!($type where $($qual)*);
+
         impl<$($qual)*> $crate::euclidean::Euclidean<$F> for $type
         // where
         //     Self: $crate::euclidean::wrap::Wrapped<WrappedField = $F>
@@ -153,6 +155,8 @@
         //         + for<'b> std::ops::SubAssign<&'b Self>
         //         + std::ops::Neg<Output = <Self as $crate::linops::AXPY>::Owned>,
         {
+            type PrincipalE = Self;
+
             fn dot<I: $crate::instance::Instance<Self>>(&self, other: I) -> $F {
                 other.eval_decompose(|x| self.get_view().dot(&x.get_view()))
             }
@@ -166,6 +170,22 @@
             }
         }
 
+        impl<$($qual)*> $crate::linops::VectorSpace for $type
+        // where
+        //     Self : $crate::euclidean::wrap::Wrapped<WrappedField = $F>,
+        //     Self::Unwrapped : $crate::linops::AXPY<Field = F>,
+        //     Self: std::ops::MulAssign<F> + std::ops::DivAssign<F>,
+        //     Self::Unwrapped: std::ops::MulAssign<F> + std::ops::DivAssign<F>,
+        {
+            type Field = $F;
+            type PrincipalV = Self;
+
+            /// Return a similar zero as `self`.
+            fn similar_origin(&self) -> Self::PrincipalV {
+                Self::wrap(self.get_view().similar_origin())
+            }
+        }
+
         impl<$($qual)*> $crate::linops::AXPY for $type
         // where
         //     Self : $crate::euclidean::wrap::Wrapped<WrappedField = $F>,
@@ -173,10 +193,7 @@
         //     Self: std::ops::MulAssign<F> + std::ops::DivAssign<F>,
         //     Self::Unwrapped: std::ops::MulAssign<F> + std::ops::DivAssign<F>,
         {
-            type Field = $F;
-            type Owned = Self;
-
-            fn axpy<I: $crate::instance::Instance<Self>>(&mut self, α: $F, x: I, β: $F) {
+             fn axpy<I: $crate::instance::Instance<Self>>(&mut self, α: $F, x: I, β: $F) {
                 x.eval_decompose(|v| {
                     self.get_view_mut().axpy(α, v.get_view(), β)
                 })
@@ -194,11 +211,6 @@
                 })
             }
 
-            /// Return a similar zero as `self`.
-            fn similar_origin(&self) -> Self::Owned {
-                Self::wrap(self.get_view().similar_origin())
-            }
-
             /// Set self to zero.
             fn set_zero(&mut self) {
                 self.get_mut_view().set_zero()
@@ -207,6 +219,7 @@
 
         impl<$($qual)*> $crate::instance::Space for $type {
             type Decomp = <<Self as $crate::euclidean::wrap::Wrapped>::Unwrapped as $crate::instance::Space>::Decomp;
+            type Principal = Self;
         }
     };
 }
diff -r 102421d462d1 -r 53ab61a41d70 src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Wed Sep 03 21:03:47 2025 -0500
+++ b/src/nalgebra_support.rs	Wed Sep 03 21:31:13 2025 -0500
@@ -389,7 +389,7 @@
 
 impl<SM, M, E> Projection<E, Linfinity> for Vector<E, M, SM>
 where
-    SM: StorageMut<E, M> + Clone,
+    SM: StorageMut<E, M>,
     M: Dim,
     E: Scalar + Zero + One + Float + RealField,
     DefaultAllocator: Allocator<M>,
@@ -405,7 +405,7 @@
 
 impl<SM, M, E> ProjectionMut<E, Linfinity> for Vector<E, M, SM>
 where
-    SM: StorageMut<E, M> + Clone,
+    SM: StorageMut<E, M>,
     M: Dim,
     E: Scalar + Zero + One + Copy + Float + RealField,
     DefaultAllocator: Allocator<M>,
@@ -464,26 +464,25 @@
     })
 }
 
-// TODO: should allow different input storages in `Euclidean`.
-
-impl<E, M, S> Euclidean<E> for Vector<E, M, S>
+impl<E, M, N, S> Euclidean<E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
-    type PrincipalE = OVector<E, M>;
+    type PrincipalE = OMatrix<E, M, N>;
 
     #[inline]
     fn dot<I: Instance<Self>>(&self, other: I) -> E {
-        other.eval_ref(|ref r| Vector::<E, M, S>::dot(self, r))
+        other.eval_ref(|ref r| Matrix::<E, M, N, S>::dot(self, r))
     }
 
     #[inline]
     fn norm2_squared(&self) -> E {
-        Vector::<E, M, S>::norm_squared(self)
+        Matrix::<E, M, N, S>::norm_squared(self)
     }
 
     #[inline]
@@ -507,13 +506,14 @@
 }
 
 /// The default norm for `Vector` is [`L2`].
-impl<E, M, S> Normed<E> for Vector<E, M, S>
+impl<E, M, N, S> Normed<E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     type NormExp = L2;
 
@@ -524,22 +524,24 @@
 
     #[inline]
     fn is_zero(&self) -> bool {
-        Vector::<E, M, S>::norm_squared(self) == E::ZERO
+        Matrix::<E, M, N, S>::norm_squared(self) == E::ZERO
     }
 }
 
-impl<E, M, S> HasDual<E> for Vector<E, M, S>
+impl<E, M, N, S> HasDual<E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
-    type DualSpace = OVector<E, M>;
+    type DualSpace = OMatrix<E, M, N>;
 
-    fn dual_origin(&self) -> OVector<E, M> {
-        OVector::zeros_generic(M::from_usize(self.len()), Const)
+    fn dual_origin(&self) -> OMatrix<E, M, N> {
+        let (m, n) = self.shape_generic();
+        OMatrix::zeros_generic(m, n)
     }
 }
 
@@ -560,7 +562,7 @@
 impl<E, M, S> Dist<L1, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    S: Storage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
     ShapeConstraint: StridesOk<E, M, U1, S>,
@@ -571,13 +573,14 @@
     }
 }
 
-impl<E, M, S> Norm<L2, E> for Vector<E, M, S>
+impl<E, M, N, S> Norm<L2, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn norm(&self, _: L2) -> E {
@@ -585,13 +588,14 @@
     }
 }
 
-impl<E, M, S> Dist<L2, E> for Vector<E, M, S>
+impl<E, M, N, S> Dist<L2, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn dist<I: Instance<Self>>(&self, other: I, _: L2) -> E {
@@ -599,13 +603,14 @@
     }
 }
 
-impl<E, M, S> Norm<Linfinity, E> for Vector<E, M, S>
+impl<E, M, N, S> Norm<Linfinity, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn norm(&self, _: Linfinity) -> E {
@@ -613,13 +618,14 @@
     }
 }
 
-impl<E, M, S> Dist<Linfinity, E> for Vector<E, M, S>
+impl<E, M, N, S> Dist<Linfinity, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn dist<I: Instance<Self>>(&self, other: I, _: Linfinity) -> E {