Mercurial > repos > alg_tools / changeset

--- a/src/instance.rs	Tue Sep 02 12:37:53 2025 -0500
+++ b/src/instance.rs	Tue Sep 02 15:11:35 2025 -0500
@@ -155,7 +155,7 @@
     /// Type for a lightweight intermediate conversion that does not own the original variable.
     /// Usually this is just a reference, but may also be a lightweight structure that
     /// contains references; see the implementation for [`crate::direct_product::Pair`].
-    type Reference<'b>: Instance<X, Self> + Copy
+    type Reference<'b>: Instance<X, Self>
     where
         X: 'b;
 

--- a/src/nalgebra_support.rs	Tue Sep 02 12:37:53 2025 -0500
+++ b/src/nalgebra_support.rs	Tue Sep 02 15:11:35 2025 -0500
@@ -19,7 +19,7 @@
 use nalgebra::{
     ClosedAddAssign, ClosedMulAssign, DefaultAllocator, Dim, LpNorm, Matrix, MatrixView, OMatrix,
     OVector, RawStorage, RealField, Scalar, SimdComplexField, Storage, StorageMut, UniformNorm,
-    Vector,
+    Vector, U1,
 };
 use num_traits::identities::{One, Zero};
 use std::ops::Mul;
@@ -45,6 +45,39 @@
     }
 }
 
+trait StridesOk<E, N, M = U1, S = <DefaultAllocator as Allocator<N, M>>::Buffer<E>>:
+    DimEq<Dyn, S::RStride>
+    + DimEq<Dyn, S::CStride>
+    + DimEq<Dyn, <<DefaultAllocator as Allocator<N, M>>::Buffer<E> as RawStorage<E, N, M>>::RStride>
+    + DimEq<Dyn, <<DefaultAllocator as Allocator<N, M>>::Buffer<E> as RawStorage<E, N, M>>::CStride>
+where
+    S: RawStorage<E, N, M>,
+    E: Scalar,
+    N: Dim,
+    M: Dim,
+    DefaultAllocator: Allocator<N, M>,
+{
+}
+
+impl<S, E, N, M> StridesOk<E, N, M, S> for ShapeConstraint
+where
+    ShapeConstraint: DimEq<Dyn, S::RStride>
+        + DimEq<Dyn, S::CStride>
+        + DimEq<
+            Dyn,
+            <<DefaultAllocator as Allocator<N, M>>::Buffer<E> as RawStorage<E, N, M>>::RStride,
+        > + DimEq<
+            Dyn,
+            <<DefaultAllocator as Allocator<N, M>>::Buffer<E> as RawStorage<E, N, M>>::CStride,
+        >,
+    S: Storage<E, N, M>,
+    E: Scalar,
+    N: Dim,
+    M: Dim,
+    DefaultAllocator: Allocator<N, M>,
+{
+}
+
 impl<SM, N, M, E> Space for Matrix<E, N, M, SM>
 where
     SM: Storage<E, N, M>,
@@ -52,6 +85,7 @@
     M: Dim,
     E: Scalar + Zero + One,
     DefaultAllocator: Allocator<N, M>,
+    ShapeConstraint: StridesOk<E, N, M, SM> + StridesOk<E, N, M>,
 {
     type OwnedSpace = OMatrix<E, N, M>;
     type Decomp = MatrixDecomposition;
@@ -60,13 +94,60 @@
 #[derive(Copy, Clone, Debug)]
 pub struct MatrixDecomposition;
 
-impl<E, M, K, S> Decomposition<Matrix<E, M, K, S>> for MatrixDecomposition
+/// nalgebra [`Storage`] that's convertible to dynamic strides (practically all, but we need this
+/// to work around the type system).
+trait ViewableStorage<E, M, K = U1>: Storage<E, M, K>
+where
+    M: Dim,
+    K: Dim,
+    E: Scalar,
+    ShapeConstraint: StridesOk<E, M, K, Self>,
+    DefaultAllocator: Allocator<M, K>,
+{
+}
+
+impl<M, K, E, S> ViewableStorage<E, M, K> for S
 where
     S: Storage<E, M, K>,
     M: Dim,
     K: Dim,
+    E: Scalar,
+    ShapeConstraint: StridesOk<E, M, K, Self>,
+    DefaultAllocator: Allocator<M, K>,
+{
+}
+
+/// nalgebra [`StorageMut`] that's convertible to dynamic strides (practically all, but we need this
+/// to work around the type system).
+trait ViewableStorageMut<E, M, K = U1>: StorageMut<E, M, K>
+where
+    M: Dim,
+    K: Dim,
+    E: Scalar,
+    ShapeConstraint: StridesOk<E, M, K, Self>,
+    DefaultAllocator: Allocator<M, K>,
+{
+}
+
+impl<M, K, E, S> ViewableStorageMut<E, M, K> for S
+where
+    S: StorageMut<E, M, K>,
+    M: Dim,
+    K: Dim,
+    E: Scalar,
+    ShapeConstraint: StridesOk<E, M, K, Self>,
+    DefaultAllocator: Allocator<M, K>,
+{
+}
+
+impl<E, M, K, S> Decomposition<Matrix<E, M, K, S>> for MatrixDecomposition
+where
+    S: ViewableStorage<E, M, K>,
+    M: Dim,
+    K: Dim,
     E: Scalar + Zero + One,
     DefaultAllocator: Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K, S> + StridesOk<E, M, K>,
 {
     type OwnedInstance = OMatrix<E, M, K>;
 
@@ -88,96 +169,85 @@
     }
 }
 
-macro_rules! impl_instances {
-    ($rstride:ty, $cstride:ty where $($qual:tt)*) => {
-        impl<$($qual)* S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition> for Matrix<E, M, K, S2>
-        where
-            S1: Storage<E, M, K>,
-            S2: RawStorage<E, M, K, RStride=$rstride, CStride=$cstride>,
-            //ShapeConstraint: DimEq<Dyn, <S2 as RawStorage<E, M, K>>::RStride>
-            //    + DimEq<Dyn, <S2 as RawStorage<E, M, K>>::CStride>,
-            M: Dim,
-            K: Dim,
-            E: Scalar + Zero + One,
-            DefaultAllocator: Allocator<M, K>,
-        {
-            fn eval_decompose<'b, R>(self, f: impl FnOnce(OMatrix<E, M, K>) -> R) -> R
-            where
-                Self: 'b,
-            {
-                f(self.into_owned())
-            }
-
-            fn eval_ref_decompose<'b, R>(
-                &'b self,
-                f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
-            ) -> R
-            where
-                Self: 'b,
-                Matrix<E, M, K, S1>: 'b,
-            {
-                f(self.as_view::<M, K, Dyn, Dyn>())
-            }
+impl<S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition> for Matrix<E, M, K, S2>
+where
+    S1: ViewableStorage<E, M, K>,
+    S2: ViewableStorage<E, M, K>,
+    M: Dim,
+    K: Dim,
+    E: Scalar + Zero + One,
+    DefaultAllocator: Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K, S1> + StridesOk<E, M, K, S2>,
+{
+    fn eval_decompose<'b, R>(self, f: impl FnOnce(OMatrix<E, M, K>) -> R) -> R
+    where
+        Self: 'b,
+    {
+        f(self.into_owned())
+    }
 
-            #[inline]
-            fn own(self) -> OMatrix<E, M, K> {
-                self.into_owned()
-            }
-        }
+    fn eval_ref_decompose<'b, R>(
+        &'b self,
+        f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
+    ) -> R
+    where
+        Self: 'b,
+        Matrix<E, M, K, S1>: 'b,
+    {
+        f(self.as_view::<M, K, Dyn, Dyn>())
+    }
 
-        impl<'a, $($qual)* S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition>
-            for &'a Matrix<E, M, K, S2>
-        where
-            S1: Storage<E, M, K>,
-            S2: RawStorage<E, M, K, RStride=$rstride, CStride=$cstride>,
-            M: Dim,
-            K: Dim,
-            E: Scalar + Zero + One,
-            DefaultAllocator: Allocator<M, K>,
-        {
-            fn eval_decompose<'b, R>(self, f: impl FnOnce(OMatrix<E, M, K>) -> R) -> R
-            where
-                Self: 'b,
-            {
-                f(self.into_owned())
-            }
-
-            fn eval_ref_decompose<'b, R>(
-                &'b self,
-                f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
-            ) -> R
-            where
-                Self: 'b,
-                Matrix<E, M, K, S1>: 'b,
-            {
-                f((*self).as_view::<M, K, Dyn, Dyn>())
-            }
-
-            #[inline]
-            fn own(self) -> OMatrix<E, M, K> {
-                self.into_owned()
-            }
-        }
-    };
+    #[inline]
+    fn own(self) -> OMatrix<E, M, K> {
+        self.into_owned()
+    }
 }
 
-impl_instances!(Dyn, Dyn where );
-impl_instances!(Const<RS>, Dyn where const RS : usize,);
-impl_instances!(Dyn, Const<CS> where const CS : usize,);
-impl_instances!(Const<RS>, Const<CS> where const RS : usize, const CS : usize,);
+impl<'a, S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition>
+    for &'a Matrix<E, M, K, S2>
+where
+    S1: ViewableStorage<E, M, K>,
+    S2: ViewableStorage<E, M, K>,
+    M: Dim,
+    K: Dim,
+    E: Scalar + Zero + One,
+    DefaultAllocator: Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K, S1> + StridesOk<E, M, K, S2>,
+{
+    fn eval_decompose<'b, R>(self, f: impl FnOnce(OMatrix<E, M, K>) -> R) -> R
+    where
+        Self: 'b,
+    {
+        f(self.into_owned())
+    }
+
+    fn eval_ref_decompose<'b, R>(
+        &'b self,
+        f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
+    ) -> R
+    where
+        Self: 'b,
+        Matrix<E, M, K, S1>: 'b,
+    {
+        f((*self).as_view::<M, K, Dyn, Dyn>())
+    }
+
+    #[inline]
+    fn own(self) -> OMatrix<E, M, K> {
+        self.into_owned()
+    }
+}
 
 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,
+    SM: ViewableStorage<E, N, M>,
+    SV: ViewableStorage<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>,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<N, M> + Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, N, M, SM> + StridesOk<E, M, K, SV> + StridesOk<E, N, K>,
 {
     type Codomain = OMatrix<E, N, K>;
 
@@ -189,33 +259,29 @@
 
 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> + Clone,
+    SM: ViewableStorage<E, N, M>,
+    SV: ViewableStorage<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>,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<N, M> + Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, N, M, SM> + StridesOk<E, M, K, SV> + StridesOk<E, N, K>,
 {
 }
 
 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> + Clone,
-    SV2: StorageMut<E, N, K>,
+    SM: ViewableStorage<E, N, M>,
+    SV1: ViewableStorage<E, M, K> + Clone,
+    SV2: ViewableStorageMut<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>,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<N, M> + Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, N, M, SM> + StridesOk<E, M, K, SV1> + StridesOk<E, N, K, SV2>,
 {
     #[inline]
     fn gemv<I: Instance<Matrix<E, M, K, SV1>>>(
@@ -236,11 +302,12 @@
 
 impl<S, M, N, E> VectorSpace for Matrix<E, M, N, S>
 where
-    S: Storage<E, M, N>,
+    S: ViewableStorage<E, M, N>,
     M: Dim,
     N: Dim,
     E: Scalar + Zero + One + Float,
     DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     type Field = E;
     type Owned = OMatrix<E, M, N>;
@@ -254,12 +321,13 @@
 
 impl<SM, SV1, M, N, E> AXPY<Matrix<E, M, N, SV1>> for Matrix<E, M, N, SM>
 where
-    SM: StorageMut<E, M, N>,
-    SV1: Storage<E, M, N>,
+    SM: ViewableStorageMut<E, M, N>,
+    SV1: ViewableStorage<E, M, N>,
     M: Dim,
     N: Dim,
     E: Scalar + Zero + One + Float,
     DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, SM> + StridesOk<E, M, N, SV1>,
 {
     #[inline]
     fn axpy<I: Instance<Matrix<E, M, N, SV1>>>(&mut self, α: E, x: I, β: E) {
@@ -300,10 +368,11 @@
 
 impl<SM, M, E> Projection<E, Linfinity> for Vector<E, M, SM>
 where
-    SM: Storage<E, M> + Clone,
+    SM: ViewableStorageMut<E, M> + Clone,
     M: Dim,
     E: Scalar + Zero + One + Float + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, SM>,
 {
     #[inline]
     fn proj_ball(self, ρ: E, exp: Linfinity) -> <Self as Space>::OwnedSpace {
@@ -315,10 +384,11 @@
 
 impl<SM, M, E> ProjectionMut<E, Linfinity> for Vector<E, M, SM>
 where
-    SM: StorageMut<E, M> + Clone,
+    SM: ViewableStorageMut<E, M> + Clone,
     M: Dim,
     E: Scalar + Zero + One + Float + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, SM>,
 {
     #[inline]
     fn proj_ball_mut(&mut self, ρ: E, _: Linfinity) {
@@ -330,17 +400,18 @@
 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> + Clone,
-    SV2: Storage<E, N, K> + Clone,
+    SM: ViewableStorage<E, N, M>,
+    SV1: ViewableStorage<E, M, K> + Clone,
+    SV2: ViewableStorage<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>,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<N, M> + Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, N, M, SM>
+        + StridesOk<E, M, K, SV1>
+        + StridesOk<E, N, K, SV2>
+        + StridesOk<E, M, N>,
 {
     type AdjointCodomain = OMatrix<E, M, K>;
     type Adjoint<'a>
@@ -382,9 +453,10 @@
 impl<E, M, S> Euclidean<E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    S: ViewableStorage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     type OwnedEuclidean = OVector<E, M>;
 
@@ -407,9 +479,10 @@
 impl<E, M, S> StaticEuclidean<E> for Vector<E, M, S>
 where
     M: DimName,
-    S: Storage<E, M>,
+    S: ViewableStorage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     #[inline]
     fn origin() -> OVector<E, M> {
@@ -421,9 +494,10 @@
 impl<E, M, S> Normed<E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    S: ViewableStorage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     type NormExp = L2;
 
@@ -441,9 +515,10 @@
 impl<E, M, S> HasDual<E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    S: ViewableStorage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     type DualSpace = OVector<E, M>;
 
@@ -455,9 +530,10 @@
 impl<E, M, S> Norm<L1, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    S: ViewableStorage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     #[inline]
     fn norm(&self, _: L1) -> E {
@@ -468,9 +544,10 @@
 impl<E, M, S> Dist<L1, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    S: ViewableStorage<E, M> + Clone,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     #[inline]
     fn dist<I: Instance<Self>>(&self, other: I, _: L1) -> E {
@@ -481,9 +558,10 @@
 impl<E, M, S> Norm<L2, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    S: ViewableStorage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     #[inline]
     fn norm(&self, _: L2) -> E {
@@ -494,9 +572,10 @@
 impl<E, M, S> Dist<L2, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    S: ViewableStorage<E, M> + Clone,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     #[inline]
     fn dist<I: Instance<Self>>(&self, other: I, _: L2) -> E {
@@ -507,9 +586,10 @@
 impl<E, M, S> Norm<Linfinity, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    S: ViewableStorage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     #[inline]
     fn norm(&self, _: Linfinity) -> E {
@@ -520,9 +600,10 @@
 impl<E, M, S> Dist<Linfinity, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    S: ViewableStorage<E, M> + Clone,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M, U1, S>,
 {
     #[inline]
     fn dist<I: Instance<Self>>(&self, other: I, _: Linfinity) -> E {