# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1735651802 18000
# Node ID 9226980e45a704d75ae1b3c8c808e318647792af
# Parent  1a38447a89fa5c573aef980d6cc92ffa34ab6dfd
Significantly simplify Mapping / Apply through Instance

diff -r 1a38447a89fa -r 9226980e45a7 Cargo.lock
--- a/Cargo.lock	Sat Dec 14 09:31:27 2024 -0500
+++ b/Cargo.lock	Tue Dec 31 08:30:02 2024 -0500
@@ -17,6 +17,7 @@
  "rayon",
  "serde",
  "serde_json",
+ "simba",
 ]
 
 [[package]]
diff -r 1a38447a89fa -r 9226980e45a7 Cargo.toml
--- a/Cargo.toml	Sat Dec 14 09:31:27 2024 -0500
+++ b/Cargo.toml	Tue Dec 31 08:30:02 2024 -0500
@@ -23,6 +23,7 @@
 cpu-time = "~1.0.0"
 serde_json = "~1.0.85"
 rayon = "1.5.3"
+simba = "0.9.0"
 
 [package.metadata.docs.rs]
 rustdoc-args = [ "--html-in-header", "katex-header.html" ]
@@ -32,7 +33,7 @@
 debug = true
 
 [features]
-default = []
+default = ["nightly"]
 use_custom_thread_pool = []
 nightly = [] # enable for higher-performance implementations of some things
 
diff -r 1a38447a89fa -r 9226980e45a7 src/bisection_tree/btfn.rs
--- a/src/bisection_tree/btfn.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/bisection_tree/btfn.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -4,7 +4,10 @@
 use std::marker::PhantomData;
 use std::sync::Arc;
 use crate::types::Float;
-use crate::mapping::{Apply, Mapping, Differentiable};
+use crate::mapping::{
+    Instance, Mapping, DifferentiableImpl, DifferentiableMapping, Space,
+    BasicDecomposition,
+};
 //use crate::linops::{Apply, Linear};
 use crate::sets::Set;
 use crate::sets::Cube;
@@ -40,10 +43,22 @@
 }
 
 impl<F : Float, G, BT, const N : usize>
+Space for BTFN<F, G, BT, N>
+where
+    G : SupportGenerator<F, N, Id=BT::Data>,
+    G::SupportType : LocalAnalysis<F, BT::Agg, N>,
+    BT : BTImpl<F, N>
+{
+    type Decomp = BasicDecomposition;
+}
+
+impl<F : Float, G, BT, const N : usize>
 BTFN<F, G, BT, N>
-where G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N>,
-      BT : BTImpl<F, N> {
+where
+    G : SupportGenerator<F, N, Id=BT::Data>,
+    G::SupportType : LocalAnalysis<F, BT::Agg, N>,
+    BT : BTImpl<F, N>
+{
 
     /// Create a new BTFN from a support generator and a pre-initialised bisection tree.
     ///
@@ -390,64 +405,40 @@
 // Apply, Mapping, Differentiate
 //
 
-impl<'a, F : Float, G, BT, V, const N : usize> Apply<&'a Loc<F, N>>
+impl<F : Float, G, BT, V, const N : usize> Mapping<Loc<F, N>>
 for BTFN<F, G, BT, N>
-where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<&'a Loc<F, N>, Output = V>,
-      V : Sum {
+where
+    BT : BTImpl<F, N>,
+    G : SupportGenerator<F, N, Id=BT::Data>,
+    G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F, N>, Codomain = V>,
+    V : Sum + Space,
+{
 
-    type Output = V;
+    type Codomain = V;
 
-    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
-        self.bt.iter_at(x)
-            .map(|&d| self.generator.support_for(d).apply(x)).sum()
+    fn apply<I : Instance<Loc<F,N>>>(&self, x : I) -> Self::Codomain {
+        let xc = x.cow();
+        self.bt.iter_at(&*xc)
+            .map(|&d| self.generator.support_for(d).apply(&*xc)).sum()
     }
 }
 
-impl<F : Float, G, BT, V, const N : usize> Apply<Loc<F, N>>
+impl<F : Float, G, BT, V, const N : usize> DifferentiableImpl<Loc<F, N>>
 for BTFN<F, G, BT, N>
-where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Apply<Loc<F, N>, Output = V>,
-      V : Sum {
-
-    type Output = V;
-
-    fn apply(&self, x : Loc<F, N>) -> Self::Output {
-        self.bt.iter_at(&x)
-            .map(|&d| self.generator.support_for(d).apply(x)).sum()
-    }
-}
-
-impl<'a, F : Float, G, BT, V, const N : usize> Differentiable<&'a Loc<F, N>>
-for BTFN<F, G, BT, N>
-where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Differentiable<&'a Loc<F, N>, Derivative = V>,
-      V : Sum {
+where
+    BT : BTImpl<F, N>,
+    G : SupportGenerator<F, N, Id=BT::Data>,
+    G::SupportType : LocalAnalysis<F, BT::Agg, N>
+        + DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>,
+    V : Sum + Space,
+{
 
     type Derivative = V;
 
-    fn differential(&self, x : &'a Loc<F, N>) -> Self::Derivative {
-        self.bt.iter_at(x)
-            .map(|&d| self.generator.support_for(d).differential(x))
-            .sum()
-    }
-}
-
-impl<F : Float, G, BT, V, const N : usize> Differentiable<Loc<F, N>>
-for BTFN<F, G, BT, N>
-where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> + Differentiable<Loc<F, N>, Derivative = V>,
-      V : Sum {
-
-    type Derivative = V;
-
-    fn differential(&self, x : Loc<F, N>) -> Self::Derivative {
-        self.bt.iter_at(&x)
-            .map(|&d| self.generator.support_for(d).differential(x))
+    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x :I) -> Self::Derivative {
+        let xc = x.cow();
+        self.bt.iter_at(&*xc)
+            .map(|&d| self.generator.support_for(d).differential(&*xc))
             .sum()
     }
 }
@@ -840,7 +831,7 @@
 impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N>
 where BT : BTSearch<F, N, Agg=Bounds<F>>,
       G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : Mapping<Loc<F, N>,Codomain=F>
+      G::SupportType : Mapping<Loc<F, N>, Codomain=F>
                        + LocalAnalysis<F, Bounds<F>, N>,
       Cube<F, N> : P2Minimise<Loc<F, N>, F> {
 
diff -r 1a38447a89fa -r 9226980e45a7 src/bisection_tree/either.rs
--- a/src/bisection_tree/either.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/bisection_tree/either.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -3,7 +3,13 @@
 use std::sync::Arc;
 
 use crate::types::*;
-use crate::mapping::{Apply, Differentiable};
+use crate::mapping::{
+    Instance,
+    Mapping,
+    DifferentiableImpl,
+    DifferentiableMapping,
+    Space,
+};
 use crate::iter::{Mappable, MapF, MapZ};
 use crate::sets::Cube;
 use crate::loc::Loc;
@@ -177,12 +183,17 @@
     }
 }
 
-impl<F, S1, S2, X> Apply<X> for EitherSupport<S1, S2>
-where S1 : Apply<X, Output=F>,
-      S2 : Apply<X, Output=F> {
-    type Output = F;
+impl<F, S1, S2, X> Mapping<X> for EitherSupport<S1, S2>
+where
+    F : Space,
+    X : Space,
+    S1 : Mapping<X, Codomain=F>,
+    S2 : Mapping<X, Codomain=F>,
+{
+    type Codomain = F;
+
     #[inline]
-    fn apply(&self, x : X) -> F {
+    fn apply<I : Instance<X>>(&self, x : I) -> F {
         match self {
             EitherSupport::Left(ref a) => a.apply(x),
             EitherSupport::Right(ref b) => b.apply(x),
@@ -190,12 +201,17 @@
     }
 }
 
-impl<F, S1, S2, X> Differentiable<X> for EitherSupport<S1, S2>
-where S1 : Differentiable<X, Derivative=F>,
-      S2 : Differentiable<X, Derivative=F> {
-    type Derivative = F;
+impl<X, S1, S2, O> DifferentiableImpl<X> for EitherSupport<S1, S2>
+where
+    O : Space,
+    X : Space,
+    S1 : DifferentiableMapping<X, DerivativeDomain=O>,
+    S2 : DifferentiableMapping<X, DerivativeDomain=O>,
+{
+    type Derivative = O;
+
     #[inline]
-    fn differential(&self, x : X) -> F {
+    fn differential_impl<I : Instance<X>>(&self, x : I) -> O {
         match self {
             EitherSupport::Left(ref a) => a.differential(x),
             EitherSupport::Right(ref b) => b.differential(x),
diff -r 1a38447a89fa -r 9226980e45a7 src/bisection_tree/support.rs
--- a/src/bisection_tree/support.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/bisection_tree/support.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -6,7 +6,9 @@
 use std::ops::{MulAssign,DivAssign,Neg};
 use crate::types::{Float, Num};
 use crate::maputil::map2;
-use crate::mapping::{Apply, Differentiable};
+use crate::mapping::{
+    Instance, Mapping, DifferentiableImpl, DifferentiableMapping, Space
+};
 use crate::sets::Cube;
 use crate::loc::Loc;
 use super::aggregator::Bounds;
@@ -127,39 +129,23 @@
     base_fn : T,
 }
 
-impl<'a, T, V, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Shift<T,F,N>
-where T : Apply<Loc<F, N>, Output=V> {
-    type Output = V;
+impl<'a, T, V : Space, F : Float, const N : usize> Mapping<Loc<F, N>> for Shift<T,F,N>
+where T : Mapping<Loc<F, N>, Codomain=V> {
+    type Codomain = V;
+
     #[inline]
-    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
-        self.base_fn.apply(x - &self.shift)
+    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+        self.base_fn.apply(x.own() - &self.shift)
     }
 }
 
-impl<'a, T, V, F : Float, const N : usize> Apply<Loc<F, N>> for Shift<T,F,N>
-where T : Apply<Loc<F, N>, Output=V> {
-    type Output = V;
-    #[inline]
-    fn apply(&self, x : Loc<F, N>) -> Self::Output {
-        self.base_fn.apply(x - &self.shift)
-    }
-}
-
-impl<'a, T, V, F : Float, const N : usize> Differentiable<&'a Loc<F, N>> for Shift<T,F,N>
-where T : Differentiable<Loc<F, N>, Derivative=V> {
+impl<'a, T, V : Space, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>> for Shift<T,F,N>
+where T : DifferentiableMapping<Loc<F, N>, DerivativeDomain=V> {
     type Derivative = V;
+
     #[inline]
-    fn differential(&self, x : &'a Loc<F, N>) -> Self::Derivative {
-        self.base_fn.differential(x - &self.shift)
-    }
-}
-
-impl<'a, T, V, F : Float, const N : usize> Differentiable<Loc<F, N>> for Shift<T,F,N>
-where T : Differentiable<Loc<F, N>, Derivative=V> {
-    type Derivative = V;
-    #[inline]
-    fn differential(&self, x : Loc<F, N>) -> Self::Derivative {
-        self.base_fn.differential(x - &self.shift)
+    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
+        self.base_fn.differential(x.own() - &self.shift)
     }
 }
 
@@ -228,46 +214,26 @@
     pub base_fn : T,
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Apply<&'a Loc<F, N>> for Weighted<T, C>
-where T : for<'b> Apply<&'b Loc<F, N>, Output=V>,
-      V : std::ops::Mul<F,Output=V>,
+impl<'a, T, V, F : Float, C, const N : usize> Mapping<Loc<F, N>> for Weighted<T, C>
+where T : Mapping<Loc<F, N>, Codomain=V>,
+      V : Space + std::ops::Mul<F,Output=V>,
       C : Constant<Type=F> {
-    type Output = V;
+    type Codomain = V;
+
     #[inline]
-    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
+    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
         self.base_fn.apply(x) * self.weight.value()
     }
 }
 
-impl<'a, T, V, F : Float, C, const N : usize> Apply<Loc<F, N>> for Weighted<T, C>
-where T : Apply<Loc<F, N>, Output=V>,
-      V : std::ops::Mul<F,Output=V>,
-      C : Constant<Type=F> {
-    type Output = V;
-    #[inline]
-    fn apply(&self, x : Loc<F, N>) -> Self::Output {
-        self.base_fn.apply(x) * self.weight.value()
-    }
-}
-
-impl<'a, T, V, F : Float, C, const N : usize> Differentiable<&'a Loc<F, N>> for Weighted<T, C>
-where T : for<'b> Differentiable<&'b Loc<F, N>, Derivative=V>,
-      V : std::ops::Mul<F, Output=V>,
+impl<'a, T, V, F : Float, C, const N : usize> DifferentiableImpl<Loc<F, N>> for Weighted<T, C>
+where T : DifferentiableMapping<Loc<F, N>, DerivativeDomain=V>,
+      V : Space + std::ops::Mul<F, Output=V>,
       C : Constant<Type=F> {
     type Derivative = V;
+
     #[inline]
-    fn differential(&self, x : &'a Loc<F, N>) -> Self::Derivative {
-        self.base_fn.differential(x) * self.weight.value()
-    }
-}
-
-impl<'a, T, V, F : Float, C, const N : usize> Differentiable<Loc<F, N>> for Weighted<T, C>
-where T : Differentiable<Loc<F, N>, Derivative=V>,
-      V : std::ops::Mul<F, Output=V>,
-      C : Constant<Type=F> {
-    type Derivative = V;
-    #[inline]
-    fn differential(&self, x : Loc<F, N>) -> Self::Derivative {
+    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
         self.base_fn.differential(x) * self.weight.value()
     }
 }
@@ -380,21 +346,12 @@
     pub T
 );
 
-impl<'a, T, F : Float, const N : usize> Apply<&'a Loc<F, N>> for Normalised<T>
-where T : Norm<F, L1> + for<'b> Apply<&'b Loc<F, N>, Output=F> {
-    type Output = F;
+impl<'a, T, F : Float, const N : usize> Mapping<Loc<F, N>> for Normalised<T>
+where T : Norm<F, L1> + Mapping<Loc<F,N>, Codomain=F> {
+    type Codomain = F;
+
     #[inline]
-    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
-        let w = self.0.norm(L1);
-        if w == F::ZERO { F::ZERO } else { self.0.apply(x) / w }
-    }
-}
-
-impl<'a, T, F : Float, const N : usize> Apply<Loc<F, N>> for Normalised<T>
-where T : Norm<F, L1> + Apply<Loc<F,N>, Output=F> {
-    type Output = F;
-    #[inline]
-    fn apply(&self, x : Loc<F, N>) -> Self::Output {
+    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
         let w = self.0.norm(L1);
         if w == F::ZERO { F::ZERO } else { self.0.apply(x) / w }
     }
diff -r 1a38447a89fa -r 9226980e45a7 src/collection.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/collection.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -0,0 +1,54 @@
+/*!
+Abstract collections of objects.
+*/
+
+use crate::loc::Loc;
+
+/// An abstract collection of elements.
+pub trait Collection : IntoIterator<Item = Self::Element> {
+    /// Type of elements of the collection
+    type Element;
+    /// Iterator over references to elements of the collection
+    type RefsIter<'a> : Iterator<Item=&'a Self::Element> where Self : 'a;
+
+    /// Returns an iterator over references to elements of the collection.
+    fn iter_refs(&self) -> Self::RefsIter<'_>;
+}
+
+/// An abstract collection of mutable elements.
+pub trait CollectionMut : Collection {
+    /// Iterator over references to elements of the collection
+    type RefsIterMut<'a> : Iterator<Item=&'a mut Self::Element> where Self : 'a;
+
+    /// Returns an iterator over references to elements of the collection.
+    fn iter_refs_mut(&mut self) -> Self::RefsIterMut<'_>;
+}
+
+/// Helps implement Collection and CollectionMut for slice-like collections.
+#[macro_export]
+macro_rules! slice_like_collection {
+    ($type : ty where $($where:tt)*) => {
+        impl<$($where)*> Collection for $type {
+            type Element = E;
+            type RefsIter<'_a> = std::slice::Iter<'_a, E> where Self : '_a;
+
+            #[inline]
+            fn iter_refs(&self) -> Self::RefsIter<'_> {
+                self.iter()
+            }
+        }
+
+        impl<$($where)*> CollectionMut for $type {
+            type RefsIterMut<'_a> = std::slice::IterMut<'_a, E> where Self : '_a;
+
+            #[inline]
+            fn iter_refs_mut(&mut self) -> Self::RefsIterMut<'_> {
+                self.iter_mut()
+            }
+        }
+    }
+}
+
+slice_like_collection!(Vec<E> where E);
+slice_like_collection!([E; N] where E, const N : usize);
+slice_like_collection!(Loc<E, N> where E, const N : usize);
diff -r 1a38447a89fa -r 9226980e45a7 src/convex.rs
--- a/src/convex.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/convex.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -2,20 +2,24 @@
 Some convex analysis basics
 */
 
-use crate::mapping::{Apply, Mapping};
+use crate::types::*;
+use crate::mapping::{Mapping, Space};
+use crate::instance::{Instance, InstanceMut, DecompositionMut};
+use crate::norms::*;
 
 /// Trait for convex mappings. Has no features, just serves as a constraint
 ///
 /// TODO: should constrain `Mapping::Codomain` to implement a partial order,
 /// but this makes everything complicated with little benefit.
-pub trait ConvexMapping<Domain> : Mapping<Domain> {}
+pub trait ConvexMapping<Domain : Space> : Mapping<Domain>
+{}
 
 /// Trait for mappings with a Fenchel conjugate
 ///
 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
 /// not be convex.
-pub trait Conjugable<Domain> : Mapping<Domain> {
-    type DualDomain;
+pub trait Conjugable<Domain : Space> : Mapping<Domain> {
+    type DualDomain : Space;
     type Conjugate<'a> : ConvexMapping<Self::DualDomain> where Self : 'a;
 
     fn conjugate(&self) -> Self::Conjugate<'_>;
@@ -25,8 +29,8 @@
 ///
 /// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`],
 /// but a `Preconjugable` mapping has be convex.
-pub trait Preconjugable<Domain> : ConvexMapping<Domain> {
-    type PredualDomain;
+pub trait Preconjugable<Domain : Space> : ConvexMapping<Domain> {
+    type PredualDomain : Space;
     type Preconjugate<'a> : Mapping<Self::PredualDomain> where Self : 'a;
 
     fn preconjugate(&self) -> Self::Preconjugate<'_>;
@@ -36,15 +40,78 @@
 ///
 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
 /// not be convex.
-pub trait HasProx<Domain> : Mapping<Domain> {
+pub trait Prox<Domain : Space> : Mapping<Domain> {
     type Prox<'a> : Mapping<Domain, Codomain=Domain> where Self : 'a;
 
     /// Returns a proximal mapping with weight τ
     fn prox_mapping(&self, τ : Self::Codomain) -> Self::Prox<'_>;
 
     /// Calculate the proximal mapping with weight τ
-    fn prox(&self, z : Domain, τ : Self::Codomain) -> Domain {
+    fn prox<I : Instance<Domain>>(&self, τ : Self::Codomain, z : I) -> Domain {
         self.prox_mapping(τ).apply(z)
     }
+
+    /// Calculate the proximal mapping with weight τ in-place
+    fn prox_mut<'b>(&self, τ : Self::Codomain, y : &'b mut Domain)
+    where
+        &'b mut Domain : InstanceMut<Domain>,
+        Domain:: Decomp : DecompositionMut<Domain>,
+        for<'a> &'a Domain : Instance<Domain>,
+    {
+        *y = self.prox(τ, &*y);
+    }
 }
 
+
+pub struct NormConjugate<F : Float, E : NormExponent>(NormMapping<F, E>);
+
+impl<Domain, E, F> ConvexMapping<Domain> for NormMapping<F, E>
+where
+    Domain : Space,
+    E : NormExponent,
+    F : Float,
+    Self : Mapping<Domain, Codomain=F> {}
+
+
+impl<Domain, E, F> ConvexMapping<Domain> for NormConjugate<F, E>
+where
+    Domain : Space,
+    E : NormExponent,
+    F : Float,
+    Self : Mapping<Domain, Codomain=F> {}
+
+
+impl<F, E, Domain> Mapping<Domain> for NormConjugate<F, E>
+where
+    Domain : Space + Norm<F, E>,
+    F : Float,
+    E : NormExponent,
+{
+    type Codomain = F;
+
+    fn apply<I : Instance<Domain>>(&self, d : I) -> F {
+        if d.eval(|x| x.norm(self.0.exponent)) <= F::ONE {
+            F::ZERO
+        } else {
+            F::INFINITY
+        }
+    }
+}
+
+
+
+impl<E, F, Domain> Conjugable<Domain> for NormMapping<F, E>
+where
+    E : NormExponent + Clone,
+    F : Float,
+    Domain : Norm<F, E> + Space,
+{
+
+    type DualDomain = Domain;
+    type Conjugate<'a> = NormConjugate<F, E> where Self : 'a;
+
+    fn conjugate(&self) -> Self::Conjugate<'_> {
+        NormConjugate(self.clone())
+    }
+}
+
diff -r 1a38447a89fa -r 9226980e45a7 src/direct_product.rs
--- a/src/direct_product.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/direct_product.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -9,65 +9,81 @@
 use std::clone::Clone;
 use serde::{Serialize, Deserialize};
 use crate::types::{Num, Float};
-use crate::{maybe_lifetime, maybe_ref, impl_vectorspace_ops};
+use crate::{maybe_lifetime, maybe_ref};
 use crate::euclidean::{Dot, Euclidean};
+use crate::instance::{Instance, InstanceMut, Decomposition, DecompositionMut, MyCow};
+use crate::mapping::Space;
+use crate::linops::AXPY;
+use crate::loc::Loc;
+use crate::norms::{Norm, PairNorm, NormExponent};
 
-#[derive(Debug,Clone,PartialEq,Eq,Serialize,Deserialize)]
+#[derive(Debug,Clone,Copy,PartialEq,Eq,Serialize,Deserialize)]
 pub struct Pair<A, B> (pub A, pub B);
 
 impl<A, B> Pair<A,B> {
-    pub fn new(a : A, b : B) -> Pair<A,B> { Pair{ 0 : a, 1 : b } }
+    pub fn new(a : A, b : B) -> Pair<A,B> { Pair(a, b) }
 }
 
 impl<A, B> From<(A,B)> for Pair<A,B> {
     #[inline]
-    fn from((a, b) : (A, B)) -> Pair<A,B> { Pair{ 0 : a, 1 : b } }
+    fn from((a, b) : (A, B)) -> Pair<A,B> { Pair(a, b) }
+}
+
+impl<A, B> From<Pair<A,B>> for (A,B) {
+    #[inline]
+    fn from(Pair(a, b) : Pair<A, B>) -> (A,B) { (a, b) }
 }
 
 macro_rules! impl_binop {
-    ($trait : ident, $fn : ident, $refl:ident, $refr:ident) => {
-        impl_binop!(@doit: $trait, $fn;
-                           maybe_lifetime!($refl, &'l Pair<A,B>),
-                               (maybe_lifetime!($refl, &'l A), maybe_lifetime!($refl, &'l B));
+    (($a : ty, $b : ty), $trait : ident, $fn : ident, $refl:ident, $refr:ident) => {
+        impl_binop!(@doit: $a, $b, $trait, $fn;
+                           maybe_lifetime!($refl, &'l Pair<$a,$b>),
+                           (maybe_lifetime!($refl, &'l $a),
+                            maybe_lifetime!($refl, &'l $b));
                            maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
-                               (maybe_lifetime!($refr, &'r Ai), maybe_lifetime!($refr, &'r Bi));
+                           (maybe_lifetime!($refr, &'r Ai),
+                            maybe_lifetime!($refr, &'r Bi));
                            $refl, $refr);
     };
 
-    (@doit: $trait:ident, $fn:ident;
+    (@doit: $a:ty, $b:ty,
+            $trait:ident, $fn:ident;
             $self:ty, ($aself:ty, $bself:ty);
             $in:ty, ($ain:ty, $bin:ty);
             $refl:ident, $refr:ident) => {
-        impl<'l, 'r, A, B, Ai, Bi> $trait<$in>
+        impl<'l, 'r, Ai, Bi> $trait<$in>
         for $self
         where $aself: $trait<$ain>,
               $bself: $trait<$bin> {
             type Output = Pair<<$aself as $trait<$ain>>::Output,
-                                       <$bself as $trait<$bin>>::Output>;
+                               <$bself as $trait<$bin>>::Output>;
 
             #[inline]
             fn $fn(self, y : $in) -> Self::Output {
-                Pair { 0 : maybe_ref!($refl, self.0).$fn(maybe_ref!($refr, y.0)),
-                               1 : maybe_ref!($refl, self.1).$fn(maybe_ref!($refr, y.1)) }
+                Pair(maybe_ref!($refl, self.0).$fn(maybe_ref!($refr, y.0)),
+                     maybe_ref!($refl, self.1).$fn(maybe_ref!($refr, y.1)))
             }
         }
     };
 }
 
 macro_rules! impl_assignop {
-    ($trait : ident, $fn : ident, $refr:ident) => {
-        impl_assignop!(@doit: $trait, $fn;
+    (($a : ty, $b : ty), $trait : ident, $fn : ident, $refr:ident) => {
+        impl_assignop!(@doit: $a, $b,
+                              $trait, $fn;
                               maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
-                                  (maybe_lifetime!($refr, &'r Ai), maybe_lifetime!($refr, &'r Bi));
+                              (maybe_lifetime!($refr, &'r Ai),
+                               maybe_lifetime!($refr, &'r Bi));
                               $refr);
     };
-    (@doit: $trait:ident, $fn:ident;
+    (@doit: $a : ty, $b : ty,
+            $trait:ident, $fn:ident;
             $in:ty, ($ain:ty, $bin:ty);
             $refr:ident) => {
-        impl<'r, A, B, Ai, Bi> $trait<$in>
-        for Pair<A,B>
-        where A: $trait<$ain>,
-              B: $trait<$bin> {
+        impl<'r, Ai, Bi> $trait<$in>
+        for Pair<$a,$b>
+        where $a: $trait<$ain>,
+              $b: $trait<$bin> {
             #[inline]
             fn $fn(&mut self, y : $in) -> () {
                 self.0.$fn(maybe_ref!($refr, y.0));
@@ -78,26 +94,29 @@
 }
 
 macro_rules! impl_scalarop {
-    ($trait : ident, $fn : ident, $refl:ident) => {
-        impl_scalarop!(@doit: $trait, $fn;
-                              maybe_lifetime!($refl, &'l Pair<A,B>),
-                                  (maybe_lifetime!($refl, &'l A), maybe_lifetime!($refl, &'l B));
+    (($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident, $refl:ident) => {
+        impl_scalarop!(@doit: $field,
+                              $trait, $fn;
+                              maybe_lifetime!($refl, &'l Pair<$a,$b>),
+                              (maybe_lifetime!($refl, &'l $a),
+                               maybe_lifetime!($refl, &'l $b));
                               $refl);
     };
-    (@doit: $trait:ident, $fn:ident;
+    (@doit: $field : ty,
+            $trait:ident, $fn:ident;
             $self:ty, ($aself:ty, $bself:ty);
             $refl:ident) => {
         // Scalar as Rhs
-        impl<'l, F : Num, A, B> $trait<F>
+        impl<'l> $trait<$field>
         for $self
-        where $aself: $trait<F>,
-              $bself: $trait<F> {
-            type Output = Pair<<$aself as $trait<F>>::Output,
-                                       <$bself as $trait<F>>::Output>;
+        where $aself: $trait<$field>,
+              $bself: $trait<$field> {
+            type Output = Pair<<$aself as $trait<$field>>::Output,
+                               <$bself as $trait<$field>>::Output>;
             #[inline]
-            fn $fn(self, a : F) -> Self::Output {
-                Pair{ 0 : maybe_ref!($refl, self.0).$fn(a),
-                              1 : maybe_ref!($refl, self.1).$fn(a)}
+            fn $fn(self, a : $field) -> Self::Output {
+                Pair(maybe_ref!($refl, self.0).$fn(a),
+                     maybe_ref!($refl, self.1).$fn(a))
             }
         }
     }
@@ -106,37 +125,38 @@
 // Not used due to compiler overflow
 #[allow(unused_macros)]
 macro_rules! impl_scalarlhs_op {
-    ($trait:ident, $fn:ident, $refr:ident, $field:ty) => {
-        impl_scalarlhs_op!(@doit: $trait, $fn, 
-                                  maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
-                                  (maybe_lifetime!($refr, &'r Ai), maybe_lifetime!($refr, &'r Bi));
+    (($a : ty, $b : ty), $field : ty, $trait:ident, $fn:ident, $refr:ident) => {
+        impl_scalarlhs_op!(@doit: $trait, $fn,
+                                  maybe_lifetime!($refr, &'r Pair<$a,$b>),
+                                  (maybe_lifetime!($refr, &'r $a),
+                                   maybe_lifetime!($refr, &'r $b));
                                   $refr, $field);
     };
-    (@doit: $trait:ident, $fn:ident, 
+    (@doit: $trait:ident, $fn:ident,
             $in:ty, ($ain:ty, $bin:ty);
             $refr:ident, $field:ty) => {
-        impl<'r, Ai, Bi> $trait<$in>
+        impl<'r> $trait<$in>
         for $field
         where $field : $trait<$ain>
                      + $trait<$bin> {
             type Output = Pair<<$field as $trait<$ain>>::Output,
-                                       <$field as $trait<$bin>>::Output>;
+                               <$field as $trait<$bin>>::Output>;
             #[inline]
             fn $fn(self, x : $in) -> Self::Output {
-                Pair{ 0 : self.$fn(maybe_ref!($refr, x.0)),
-                              1 : self.$fn(maybe_ref!($refr, x.1))}
+                Pair(self.$fn(maybe_ref!($refr, x.0)),
+                     self.$fn(maybe_ref!($refr, x.1)))
             }
         }
     };
 }
 
 macro_rules! impl_scalar_assignop {
-    ($trait : ident, $fn : ident) => {
-        impl<'r, F : Num, A, B> $trait<F>
-        for Pair<A,B>
-        where A: $trait<F>, B: $trait<F> {
+    (($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => {
+        impl<'r> $trait<$field>
+        for Pair<$a, $b>
+        where $a: $trait<$field>, $b: $trait<$field> {
             #[inline]
-            fn $fn(&mut self, a : F) -> () {
+            fn $fn(&mut self, a : $field) -> () {
                 self.0.$fn(a);
                 self.1.$fn(a);
             }
@@ -145,34 +165,77 @@
 }
 
 macro_rules! impl_unaryop {
-    ($trait:ident, $fn:ident, $refl:ident) => {
+    (($a : ty, $b : ty), $trait:ident, $fn:ident, $refl:ident) => {
         impl_unaryop!(@doit: $trait, $fn;
-                             maybe_lifetime!($refl, &'l Pair<A,B>),
-                                 (maybe_lifetime!($refl, &'l A), maybe_lifetime!($refl, &'l B));
+                             maybe_lifetime!($refl, &'l Pair<$a,$b>),
+                             (maybe_lifetime!($refl, &'l $a),
+                              maybe_lifetime!($refl, &'l $b));
                              $refl);
     };
     (@doit: $trait:ident, $fn:ident;
             $self:ty, ($aself:ty, $bself:ty);
             $refl : ident) => {
-        impl<'l, A, B> $trait
+        impl<'l> $trait
         for $self
         where $aself: $trait,
               $bself: $trait {
             type Output = Pair<<$aself as $trait>::Output,
-                                       <$bself as $trait>::Output>;
+                               <$bself as $trait>::Output>;
             #[inline]
             fn $fn(self) -> Self::Output {
-                Pair{ 0 : maybe_ref!($refl, self.0).$fn(),
-                              1 : maybe_ref!($refl, self.1).$fn()}
+                Pair(maybe_ref!($refl, self.0).$fn(),
+                     maybe_ref!($refl, self.1).$fn())
             }
         }
     }
 }
 
+#[macro_export]
+macro_rules! impl_pair_vectorspace_ops {
+    (($a:ty, $b:ty), $field:ty) => {
+        impl_pair_vectorspace_ops!(@binary, ($a, $b), Add, add);
+        impl_pair_vectorspace_ops!(@binary, ($a, $b), Sub, sub);
+        impl_pair_vectorspace_ops!(@assign, ($a, $b), AddAssign, add_assign);
+        impl_pair_vectorspace_ops!(@assign, ($a, $b), SubAssign, sub_assign);
+        impl_pair_vectorspace_ops!(@scalar, ($a, $b), $field, Mul, mul);
+        impl_pair_vectorspace_ops!(@scalar, ($a, $b), $field, Div, div);
+        // Compiler overflow
+        // $(
+        //     impl_pair_vectorspace_ops!(@scalar_lhs, ($a, $b), $field, $impl_scalarlhs_op, Mul, mul);
+        // )*
+        impl_pair_vectorspace_ops!(@scalar_assign, ($a, $b), $field, MulAssign, mul_assign);
+        impl_pair_vectorspace_ops!(@scalar_assign, ($a, $b), $field, DivAssign, div_assign);
+        impl_pair_vectorspace_ops!(@unary, ($a, $b), Neg, neg);
+    };
+    (@binary, ($a : ty, $b : ty), $trait : ident, $fn : ident) => {
+        impl_binop!(($a, $b), $trait, $fn, ref, ref);
+        impl_binop!(($a, $b), $trait, $fn, ref, noref);
+        impl_binop!(($a, $b), $trait, $fn, noref, ref);
+        impl_binop!(($a, $b), $trait, $fn, noref, noref);
+    };
+    (@assign, ($a : ty, $b : ty), $trait : ident, $fn :ident) => {
+        impl_assignop!(($a, $b), $trait, $fn, ref);
+        impl_assignop!(($a, $b), $trait, $fn, noref);
+    };
+    (@scalar, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn :ident) => {
+        impl_scalarop!(($a, $b), $field, $trait, $fn, ref);
+        impl_scalarop!(($a, $b), $field, $trait, $fn, noref);
+    };
+    (@scalar_lhs, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => {
+        impl_scalarlhs_op!(($a, $b), $field, $trait, $fn, ref);
+        impl_scalarlhs_op!(($a, $b), $field, $trait, $fn, noref);
+    };
+    (@scalar_assign, ($a : ty, $b : ty), $field : ty, $trait : ident, $fn : ident) => {
+        impl_scalar_assignop!(($a, $b), $field, $trait, $fn);
+    };
+    (@unary, ($a : ty, $b : ty), $trait : ident, $fn :  ident) => {
+        impl_unaryop!(($a, $b), $trait, $fn, ref);
+        impl_unaryop!(($a, $b), $trait, $fn, noref);
+    };
+}
 
-impl_vectorspace_ops!(impl_binop, impl_assignop, impl_scalarop, impl_scalarlhs_op,
-                      impl_scalar_assignop, impl_unaryop);
-
+impl_pair_vectorspace_ops!((f32, f32), f32);
+impl_pair_vectorspace_ops!((f64, f64), f64);
 
 impl<A, B, U, V, F> Dot<Pair<U, V>, F> for Pair<A, B>
 where
@@ -186,15 +249,28 @@
     }
 }
 
-impl<A, B, F> Euclidean<F>  for Pair<A, B>
+type PairOutput<F, A, B> = Pair<<A as Euclidean<F>>::Output, <B as Euclidean<F>>::Output>;
+
+impl<A, B, F> Euclidean<F> for Pair<A, B>
 where
     A : Euclidean<F>,
     B : Euclidean<F>,
-    F : Float
+    F : Float,
+    PairOutput<F, A, B> : Euclidean<F>,
+    Self : Sized + Dot<Self,F>
+          + Mul<F, Output=PairOutput<F, A, B>> + MulAssign<F>
+          + Div<F, Output=PairOutput<F, A, B>> + DivAssign<F>
+          + Add<Self, Output=PairOutput<F, A, B>>
+          + Sub<Self, Output=PairOutput<F, A, B>>
+          + for<'b> Add<&'b Self, Output=PairOutput<F, A, B>>
+          + for<'b> Sub<&'b Self, Output=PairOutput<F, A, B>>
+          + AddAssign<Self> + for<'b> AddAssign<&'b Self>
+          + SubAssign<Self> + for<'b> SubAssign<&'b Self>
+          + Neg<Output=PairOutput<F, A, B>>
 {
-    type Output = Pair<<A as Euclidean<F>>::Output, <B as Euclidean<F>>::Output>;
+    type Output = PairOutput<F, A, B>;
 
-    fn similar_origin(&self) -> <Self as Euclidean<F>>::Output {
+    fn similar_origin(&self) -> PairOutput<F, A, B> {
         Pair(self.0.similar_origin(), self.1.similar_origin())
     }
 
@@ -203,26 +279,192 @@
     }
 }
 
-use crate::linops::AXPY;
-
 impl<F, A, B, U, V> AXPY<F, Pair<U, V>> for Pair<A, B>
 where
+    U : Space,
+    V : Space,
     A : AXPY<F, U>,
     B : AXPY<F, V>,
     F : Num
 {
-    fn axpy(&mut self, α : F, x : &Pair<U, V>, β : F) {
-        self.0.axpy(α, &x.0, β);
-        self.1.axpy(α, &x.1, β);
+
+    fn axpy<I : Instance<Pair<U,V>>>(&mut self, α : F, x : I, β : F) {
+        let Pair(u, v) = x.decompose();
+        self.0.axpy(α, u, β);
+        self.1.axpy(α, v, β);
+    }
+
+    fn copy_from<I : Instance<Pair<U,V>>>(&mut self, x : I) {
+        let Pair(u, v) = x.decompose();
+        self.0.copy_from(u);
+        self.1.copy_from(v);
+    }
+
+    fn scale_from<I : Instance<Pair<U,V>>>(&mut self, α : F, x : I) {
+        let Pair(u, v) = x.decompose();
+        self.0.scale_from(α, u);
+        self.1.scale_from(α, v);
+    }
+}
+
+/// [`Decomposition`] for working with [`Pair`]s.
+#[derive(Copy, Clone, Debug)]
+pub struct PairDecomposition<D, Q>(D, Q);
+
+impl<A : Space, B : Space> Space for Pair<A, B> {
+    type Decomp = PairDecomposition<A::Decomp, B::Decomp>;
+}
+
+impl<A, B, D, Q> Decomposition<Pair<A, B>> for PairDecomposition<D,Q>
+where
+    A : Space,
+    B : Space,
+    D : Decomposition<A>,
+    Q : Decomposition<B>,
+{
+    type Decomposition<'b> = Pair<D::Decomposition<'b>, Q::Decomposition<'b>> where Pair<A, B> : 'b;
+    type Reference<'b> = Pair<D::Reference<'b>, Q::Reference<'b>> where Pair<A, B> : 'b;
+
+    #[inline]
+    fn lift<'b>(Pair(u, v) : Self::Reference<'b>) -> Self::Decomposition<'b> {
+        Pair(D::lift(u), Q::lift(v))
+    }
+}
+
+impl<A, B, U, V, D, Q> Instance<Pair<A, B>, PairDecomposition<D, Q>> for Pair<U, V>
+where
+    A : Space,
+    B : Space,
+    D : Decomposition<A>,
+    Q : Decomposition<B>,
+    U : Instance<A, D>,
+    V : Instance<B, Q>,
+{
+    #[inline]
+    fn decompose<'b>(self)
+        -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Decomposition<'b>
+    where Self : 'b, Pair<A, B> : 'b
+    {
+        Pair(self.0.decompose(), self.1.decompose())
+    }
+
+    #[inline]
+    fn ref_instance(&self)
+        -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Reference<'_>
+    {
+        Pair(self.0.ref_instance(), self.1.ref_instance())
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, Pair<A, B>> where Self : 'b{
+        MyCow::Owned(Pair(self.0.own(), self.1.own()))
     }
 
-    fn copy_from(&mut self, x : &Pair<U, V>,) {
-        self.0.copy_from(&x.0);
-        self.1.copy_from(&x.1);
+    #[inline]
+    fn own(self) -> Pair<A, B> {
+        Pair(self.0.own(), self.1.own())
+    }
+}
+
+
+impl<'a, A, B, U, V, D, Q> Instance<Pair<A, B>, PairDecomposition<D, Q>> for &'a Pair<U, V>
+where
+    A : Space,
+    B : Space,
+    D : Decomposition<A>,
+    Q : Decomposition<B>,
+    U : Instance<A, D>,
+    V : Instance<B, Q>,
+    &'a U : Instance<A, D>,
+    &'a V : Instance<B, Q>,
+{
+    #[inline]
+    fn decompose<'b>(self)
+        -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Decomposition<'b>
+    where Self : 'b, Pair<A, B> : 'b
+    {
+        Pair(D::lift(self.0.ref_instance()), Q::lift(self.1.ref_instance()))
+    }
+
+    #[inline]
+    fn ref_instance(&self)
+        -> <PairDecomposition<D, Q> as Decomposition<Pair<A, B>>>::Reference<'_>
+    {
+        Pair(self.0.ref_instance(), self.1.ref_instance())
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, Pair<A, B>> where Self : 'b {
+        MyCow::Owned(self.own())
+    }
+
+    #[inline]
+    fn own(self) -> Pair<A, B> {
+        let Pair(ref u, ref v) = self;
+        Pair(u.own(), v.own())
     }
 
-    fn scale_from(&mut self, α : F, x : &Pair<U, V>) {
-        self.0.scale_from(α, &x.0);
-        self.1.scale_from(α, &x.1);
+}
+
+impl<A, B, D, Q> DecompositionMut<Pair<A, B>> for PairDecomposition<D,Q>
+where
+    A : Space,
+    B : Space,
+    D : DecompositionMut<A>,
+    Q : DecompositionMut<B>,
+{
+    type ReferenceMut<'b> = Pair<D::ReferenceMut<'b>, Q::ReferenceMut<'b>> where Pair<A, B> : 'b;
+}
+
+impl<A, B, U, V, D, Q> InstanceMut<Pair<A, B>, PairDecomposition<D, Q>> for Pair<U, V>
+where
+    A : Space,
+    B : Space,
+    D : DecompositionMut<A>,
+    Q : DecompositionMut<B>,
+    U : InstanceMut<A, D>,
+    V : InstanceMut<B, Q>,
+{
+    #[inline]
+    fn ref_instance_mut(&mut self)
+        -> <PairDecomposition<D, Q> as DecompositionMut<Pair<A, B>>>::ReferenceMut<'_>
+    {
+        Pair(self.0.ref_instance_mut(), self.1.ref_instance_mut())
     }
-}
\ No newline at end of file
+}
+
+impl<'a, A, B, U, V, D, Q> InstanceMut<Pair<A, B>, PairDecomposition<D, Q>> for &'a mut Pair<U, V>
+where
+    A : Space,
+    B : Space,
+    D : DecompositionMut<A>,
+    Q : DecompositionMut<B>,
+    U : InstanceMut<A, D>,
+    V : InstanceMut<B, Q>,
+{
+    #[inline]
+    fn ref_instance_mut(&mut self)
+        -> <PairDecomposition<D, Q> as DecompositionMut<Pair<A, B>>>::ReferenceMut<'_>
+    {
+        Pair(self.0.ref_instance_mut(), self.1.ref_instance_mut())
+    }
+}
+
+
+impl<F, A, B, ExpA, ExpB, ExpJ> Norm<F, PairNorm<ExpA, ExpB, ExpJ>>
+for Pair<A,B>
+where
+    F : Num,
+    ExpA : NormExponent,
+    ExpB : NormExponent,
+    ExpJ : NormExponent,
+    A : Norm<F, ExpA>,
+    B : Norm<F, ExpB>,
+    Loc<F, 2> : Norm<F, ExpJ>,
+{
+    fn norm(&self, PairNorm(expa, expb, expj) : PairNorm<ExpA, ExpB, ExpJ>) -> F {
+        Loc([self.0.norm(expa), self.1.norm(expb)]).norm(expj)
+    }
+}
+
+
diff -r 1a38447a89fa -r 9226980e45a7 src/euclidean.rs
--- a/src/euclidean.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/euclidean.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -2,8 +2,9 @@
 Euclidean spaces.
 */
 
+use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg};
 use crate::types::*;
-use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg};
+use crate::mapping::Space;
 
 /// Space (type) with a defined dot product.
 ///
@@ -18,7 +19,7 @@
 /// The type should implement vector space operations (addition, subtraction, scalar
 /// multiplication and scalar division) along with their assignment versions, as well
 /// as the [`Dot`] product with respect to `Self`.
-pub trait Euclidean<F : Float> : Sized + Dot<Self,F>
+pub trait Euclidean<F : Float> : Space + Dot<Self,F>
         + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F>
         + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F>
         + Add<Self, Output=<Self as Euclidean<F>>::Output>
diff -r 1a38447a89fa -r 9226980e45a7 src/instance.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/instance.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -0,0 +1,282 @@
+/*!
+Helper traits to work with references or owned values of types and their subsets.
+*/
+
+#[derive(Clone, Copy)]
+pub enum EitherDecomp<A, B> {
+    Owned(A),
+    Borrowed(B),
+}
+
+/// A very basic implementation of [`Cow`] without a [`Clone`] trait dependency.
+pub type MyCow<'b, X> = EitherDecomp<X, &'b X>;
+
+impl<'b, X> std::ops::Deref for MyCow<'b, X> {
+    type Target = X;
+
+    #[inline]
+    fn deref(&self) -> &Self::Target {
+        match self {
+            EitherDecomp::Owned(x) => &x,
+            EitherDecomp::Borrowed(x) => x,
+        }
+    }
+}
+
+impl<'b, X> MyCow<'b, X> {
+    #[inline]
+    pub fn into_owned(self) -> X where X : Clone {
+        match self {
+            EitherDecomp::Owned(x) => x,
+            EitherDecomp::Borrowed(x) => x.clone(),
+        }
+    }
+}
+
+/// Trait for abitrary mathematical spaces.
+pub trait Space : Instance<Self, Self::Decomp> {
+    /// Default decomposition for the space
+    type Decomp : Decomposition<Self>;
+}
+
+#[macro_export]
+macro_rules! impl_basic_space {
+    ($($type:ty)*) => { $(
+        impl $crate::instance::Space for $type {
+            type Decomp = $crate::instance::BasicDecomposition;
+        }
+    )* };
+    ($type:ty where $($where:tt)*) => {
+        impl<$($where)*> $crate::instance::Space for $type {
+            type Decomp = $crate::instance::BasicDecomposition;
+        }
+    };
+}
+
+impl_basic_space!(u8 u16 u32 u64 u128 usize
+                  i8 i16 i32 i64 i128 isize
+                  f32 f64);
+
+/// Marker type for decompositions to be used with [`Instance`].
+pub trait Decomposition<X : Space> : Sized {
+    /// Possibly owned form of the decomposition
+    type Decomposition<'b> : Instance<X, Self> where X : 'b;
+    /// Unlikely owned form of the decomposition.
+    /// 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 where X : 'b;
+
+    /// Left the lightweight reference type into a full decomposition type.
+    fn lift<'b>(r : Self::Reference<'b>) -> Self::Decomposition<'b>;
+}
+
+/// Most common [`Decomposition`] (into `Either<X, &'b X>`) that allows working with owned
+/// values and all sorts of references.
+#[derive(Copy, Clone, Debug)]
+pub struct BasicDecomposition;
+
+impl<X : Space + Clone> Decomposition<X> for BasicDecomposition {
+    type Decomposition<'b> = MyCow<'b, X> where X : 'b;
+    type Reference<'b> = &'b X where X : 'b;
+
+    #[inline]
+    fn lift<'b>(r : Self::Reference<'b>) -> Self::Decomposition<'b> {
+        MyCow::Borrowed(r)
+    }
+}
+
+/// Helper trait for functions to work with either owned values or references to either the
+/// “principal type” `X` or types some present a subset of `X`. In the latter sense, this
+/// generalises [`std::borrow::ToOwned`], [`std::borrow::Borrow`], and [`std::borrow::Cow`].
+/// This type also includes iteratation facilities when `X` is a [`Collection`], to avoid
+/// the possibly costly conversion of such subset types into `X`.
+pub trait Instance<X : Space, D = <X as Space>::Decomp> : Sized where D : Decomposition<X> {
+    /// Decomposes self according to `decomposer`.
+    fn decompose<'b>(self) -> D::Decomposition<'b>
+    where Self : 'b, X : 'b;
+  
+    /// Returns a lightweight instance of `self`.
+    fn ref_instance(&self) -> D::Reference<'_>;
+
+    /// Returns an owned instance of `X`, cloning or converting non-true instances when necessary.
+    fn own(self) -> X;
+
+    // ************** automatically implemented methods below from here **************
+
+    /// Returns an owned instance or reference to `X`, converting non-true instances when necessary.
+    ///
+    /// Default implementation uses [`Self::own`]. Consumes the input.
+    fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b {
+        MyCow::Owned(self.own())
+    }
+    
+    #[inline]
+    /// Evaluates `f` on a reference to self.
+    ///
+    /// Default implementation uses [`Self::cow`]. Consumes the input.
+    fn eval<'b, R>(self, f : impl FnOnce(&X) -> R) -> R
+    where X : 'b, Self : 'b
+    {
+        f(&*self.cow())
+    }
+
+    #[inline]
+    /// Evaluates `f` or `g` depending on whether a reference or owned value is available.
+    ///
+    /// Default implementation uses [`Self::cow`]. Consumes the input.
+    fn either<'b, R>(
+        self,
+        f : impl FnOnce(X) -> R,
+        g : impl FnOnce(&X) -> R
+    ) -> R
+    where Self : 'b
+    {
+        match self.cow() {
+            EitherDecomp::Owned(x) => f(x),
+            EitherDecomp::Borrowed(x) => g(x),
+        }
+    }
+}
+
+
+impl<X : Space + Clone> Instance<X, BasicDecomposition> for X {
+    #[inline]
+    fn decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b>
+    where Self : 'b, X : 'b
+    {
+        MyCow::Owned(self)
+    }
+
+    #[inline]
+    fn own(self) -> X {
+        self
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b {
+        MyCow::Owned(self)
+    }
+
+    #[inline]
+    fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> {
+        self
+    }
+}
+
+impl<'a, X : Space + Clone> Instance<X, BasicDecomposition> for &'a X {
+    #[inline]
+    fn decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b>
+    where Self : 'b, X : 'b
+    {
+        MyCow::Borrowed(self)
+    }
+
+    #[inline]
+    fn own(self) -> X {
+        self.clone()
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b {
+        MyCow::Borrowed(self)
+    }
+
+    #[inline]
+    fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> {
+        *self
+    }
+}
+
+impl<'a, X : Space + Clone> Instance<X, BasicDecomposition> for &'a mut X {
+    #[inline]
+    fn  decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b>
+    where Self : 'b, X : 'b
+    {
+        EitherDecomp::Borrowed(self)
+    }
+
+    #[inline]
+    fn own(self) -> X {
+        self.clone()
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b, X : Clone {
+        EitherDecomp::Borrowed(self)
+    }
+
+    #[inline]
+    fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> {
+        *self
+    }
+}
+
+impl<'a, X : Space + Clone> Instance<X, BasicDecomposition> for MyCow<'a, X> {
+
+    #[inline]
+    fn  decompose<'b>(self) -> <BasicDecomposition as Decomposition<X>>::Decomposition<'b>
+    where Self : 'b, X : 'b
+    {
+        self
+    }
+
+    #[inline]
+    fn own(self) -> X {
+        match self {
+            MyCow::Borrowed(a) => a.own(),
+            MyCow::Owned(b) => b.own()
+        }
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, X> where Self : 'b {
+        match self {
+            MyCow::Borrowed(a) => a.cow(),
+            MyCow::Owned(b) => b.cow()
+        }
+    }
+
+    #[inline]
+    fn ref_instance(&self) -> <BasicDecomposition as Decomposition<X>>::Reference<'_> {
+        match self {
+            MyCow::Borrowed(a) => a,
+            MyCow::Owned(b) => &b,
+        }
+    }
+}
+
+/// Marker type for mutable decompositions to be used with [`InstanceMut`].
+pub trait DecompositionMut<X : Space> : Sized {
+    type ReferenceMut<'b> : InstanceMut<X, Self> where X : 'b;
+}
+
+
+/// Helper trait for functions to work with mutable references.
+pub trait InstanceMut<X : Space , D = <X as Space>::Decomp> : Sized where D : DecompositionMut<X> {
+    /// Returns a mutable decomposition of self.
+    fn ref_instance_mut(&mut self) -> D::ReferenceMut<'_>;
+}
+
+impl<X : Space> DecompositionMut<X> for BasicDecomposition {
+    type ReferenceMut<'b> = &'b mut X where X : 'b;
+}
+
+/// This impl may seem pointless, but allows throwaway mutable scratch variables
+impl<'a, X : Space> InstanceMut<X, BasicDecomposition> for X {
+    #[inline]
+    fn ref_instance_mut(&mut self)
+        -> <BasicDecomposition as DecompositionMut<X>>::ReferenceMut<'_>
+    {
+        self
+    }
+}
+
+impl<'a, X : Space> InstanceMut<X, BasicDecomposition> for &'a mut X {
+    #[inline]
+    fn ref_instance_mut(&mut self)
+        -> <BasicDecomposition as DecompositionMut<X>>::ReferenceMut<'_>
+    {
+        self
+    }
+}
diff -r 1a38447a89fa -r 9226980e45a7 src/lib.rs
--- a/src/lib.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/lib.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -11,13 +11,12 @@
     feature(maybe_uninit_uninit_array,maybe_uninit_array_assume_init,maybe_uninit_slice),
     feature(float_minimum_maximum),
     feature(get_mut_unchecked),
+    feature(cow_is_borrowed),
 )]
 
-// They don't work:
-//#![feature(negative_impls)]
-//#![feature(specialization)]
-
 pub mod types;
+pub mod instance;
+pub mod collection;
 pub mod nanleast;
 pub mod error;
 pub mod parallelism;
diff -r 1a38447a89fa -r 9226980e45a7 src/linops.rs
--- a/src/linops.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/linops.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -4,59 +4,69 @@
 
 use numeric_literals::replace_float_literals;
 use std::marker::PhantomData;
+use serde::Serialize;
 use crate::types::*;
-use serde::Serialize;
-pub use crate::mapping::Apply;
+pub use crate::mapping::{Mapping, Space};
 use crate::direct_product::Pair;
+use crate::instance::Instance;
+use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity};
 
 /// Trait for linear operators on `X`.
-pub trait Linear<X> : Apply<X, Output=Self::Codomain>
-                      + for<'a> Apply<&'a X, Output=Self::Codomain> {
-    type Codomain;
-}
+pub trait Linear<X : Space> : Mapping<X>
+{ }
 
 /// Efficient in-place summation.
 #[replace_float_literals(F::cast_from(literal))]
-pub trait AXPY<F : Num, X = Self> {
+pub trait AXPY<F, X = Self> : Space
+where
+    F : Num,
+    X : Space,
+{
     /// Computes  `y = βy + αx`, where `y` is `Self`.
-    fn axpy(&mut self, α : F, x : &X, β : F);
+    fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F);
 
     /// Copies `x` to `self`.
-    fn copy_from(&mut self, x : &X) {
+    fn copy_from<I : Instance<X>>(&mut self, x : I) {
         self.axpy(1.0, x, 0.0)
     }
 
     /// Computes  `y = αx`, where `y` is `Self`.
-    fn scale_from(&mut self, α : F, x : &X) {
+    fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) {
         self.axpy(α, x, 0.0)
     }
 }
 
 /// Efficient in-place application for [`Linear`] operators.
 #[replace_float_literals(F::cast_from(literal))]
-pub trait GEMV<F : Num, X, Y = <Self as Linear<X>>::Codomain> : Linear<X> {
+pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> {
     /// Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv(&self, y : &mut Y, α : F, x : &X, β : F);
+    fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F);
 
     /// Computes `y = Ax`, where `A` is `Self`
-    fn apply_mut(&self, y : &mut Y, x : &X){
+    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
         self.gemv(y, 1.0, x, 0.0)
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add(&self, y : &mut Y, x : &X){
+    fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){
         self.gemv(y, 1.0, x, 1.0)
     }
 }
 
 
 /// Bounded linear operators
-pub trait BoundedLinear<X> : Linear<X> {
-    type FloatType : Float;
+pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X>
+where
+    F : Num,
+    X : Space,
+    XExp : NormExponent,
+    CodExp : NormExponent
+{
     /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
     /// This is not expected to be the norm, just any bound on it that can be
-    /// reasonably implemented.
-    fn opnorm_bound(&self) -> Self::FloatType;
+    /// reasonably implemented. The [`NormExponent`] `xexp`  indicates the norm
+    /// in `X`, and `codexp` in the codomain.
+    fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F;
 }
 
 // Linear operator application into mutable target. The [`AsRef`] bound
@@ -70,16 +80,16 @@
 }*/
 
 /// Trait for forming the adjoint operator of `Self`.
-pub trait Adjointable<X,Yʹ> : Linear<X> {
-    type AdjointCodomain;
+pub trait Adjointable<X, Yʹ> : Linear<X>
+where
+    X : Space,
+    Yʹ : Space,
+{
+    type AdjointCodomain : Space;
     type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
 
     /// Form the adjoint operator of `self`.
     fn adjoint(&self) -> Self::Adjoint<'_>;
-
-    /*fn adjoint_apply(&self, y : &Yʹ) -> Self::AdjointCodomain {
-        self.adjoint().apply(y)
-    }*/
 }
 
 /// Trait for forming a preadjoint of an operator.
@@ -89,189 +99,206 @@
 /// operator. The space `Ypre` is the predual of its codomain, and should be the
 /// domain of the adjointed operator. `Self::Preadjoint` should be
 /// [`Adjointable`]`<'a,Ypre,X>`.
-pub trait Preadjointable<X,Ypre> : Linear<X> {
-    type PreadjointCodomain;
-    type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a;
 
-    /// Form the preadjoint operator of `self`.
+pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> {
+    type PreadjointCodomain : Space;
+    type Preadjoint<'a> : Adjointable<
+        Ypre, X, AdjointCodomain=Self::Codomain, Codomain=Self::PreadjointCodomain
+    > where Self : 'a;
+
+    /// Form the adjoint operator of `self`.
     fn preadjoint(&self) -> Self::Preadjoint<'_>;
 }
 
 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
-pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
-impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
+pub trait SimplyAdjointable<X : Space> : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
+impl<'a,X : Space, T> SimplyAdjointable<X> for T
+where T : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
 
 /// The identity operator
 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
 pub struct IdOp<X> (PhantomData<X>);
 
 impl<X> IdOp<X> {
-    fn new() -> IdOp<X> { IdOp(PhantomData) }
+    pub fn new() -> IdOp<X> { IdOp(PhantomData) }
 }
 
-impl<X> Apply<X> for IdOp<X> {
-    type Output = X;
+impl<X : Clone + Space> Mapping<X> for IdOp<X> {
+    type Codomain = X;
 
-    fn apply(&self, x : X) -> X {
-        x
+    fn apply<I : Instance<X>>(&self, x : I) -> X {
+        x.own()
     }
 }
 
-impl<'a, X> Apply<&'a X> for IdOp<X> where X : Clone {
-    type Output = X;
-    
-    fn apply(&self, x : &'a X) -> X {
-        x.clone()
-    }
-}
-
-impl<X> Linear<X> for IdOp<X> where X : Clone {
-    type Codomain = X;
-}
-
+impl<X : Clone + Space> Linear<X> for IdOp<X>
+{ }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone {
+impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X>
+where
+    Y : AXPY<F, X>,
+    X : Clone + Space
+{
     // Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv(&self, y : &mut Y, α : F, x : &X, β : F) {
+    fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
         y.axpy(α, x, β)
     }
 
-    fn apply_mut(&self, y : &mut Y, x : &X){
+    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
         y.copy_from(x);
     }
 }
 
-impl<X> BoundedLinear<X> for IdOp<X> where X : Clone {
-    type FloatType = float;
-    fn opnorm_bound(&self) -> float { 1.0 }
+impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
+where
+    X : Space + Clone,
+    F : Num,
+    E : NormExponent
+{
+    fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE }
 }
 
-impl<X> Adjointable<X,X> for IdOp<X> where X : Clone {
+impl<X : Clone + Space> Adjointable<X,X> for IdOp<X> {
     type AdjointCodomain=X;
     type Adjoint<'a> = IdOp<X> where X : 'a;
+
     fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
 }
 
+impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> {
+    type PreadjointCodomain=X;
+    type Preadjoint<'a> = IdOp<X> where X : 'a;
+
+    fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() }
+}
+
 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
 pub struct RowOp<S, T>(pub S, pub T);
 
 use std::ops::Add;
 
-impl<A, B, S, T> Apply<Pair<A, B>> for RowOp<S, T>
+impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
 where
-    S : Apply<A>,
-    T : Apply<B>,
-    S::Output : Add<T::Output>
+    A : Space,
+    B : Space,
+    S : Mapping<A>,
+    T : Mapping<B>,
+    S::Codomain : Add<T::Codomain>,
+    <S::Codomain as Add<T::Codomain>>::Output : Space,
+
 {
-    type Output = <S::Output as Add<T::Output>>::Output;
+    type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
 
-    fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output {
+    fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
+        let Pair(a, b) = x.decompose();
         self.0.apply(a) + self.1.apply(b)
     }
 }
 
-impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for RowOp<S, T>
+impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
 where
-    S : Apply<&'a A>,
-    T : Apply<&'a B>,
-    S::Output : Add<T::Output>
-{
-    type Output = <S::Output as Add<T::Output>>::Output;
+    A : Space,
+    B : Space,
+    S : Linear<A>,
+    T : Linear<B>,
+    S::Codomain : Add<T::Codomain>,
+    <S::Codomain as Add<T::Codomain>>::Output : Space,
+{ }
 
-    fn apply(&self, Pair(ref a, ref b) : &'a Pair<A, B>) -> Self::Output {
-        self.0.apply(a) + self.1.apply(b)
-    }
-}
 
-impl<A, B, S, T, D> Linear<Pair<A, B>> for RowOp<S, T>
+impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
 where
-    RowOp<S, T> : Apply<Pair<A, B>, Output=D> + for<'a>  Apply<&'a Pair<A, B>, Output=D>,
-{
-    type Codomain = D;
-}
-
-impl<F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
-where
+    U : Space,
+    V : Space,
     S : GEMV<F, U, Y>,
     T : GEMV<F, V, Y>,
     F : Num,
     Self : Linear<Pair<U, V>, Codomain=Y>
 {
-    fn gemv(&self, y : &mut Y, α : F, x : &Pair<U, V>, β : F) {
-        self.0.gemv(y, α, &x.0, β);
-        self.1.gemv(y, α, &x.1, β);
+    fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) {
+        let Pair(u, v) = x.decompose();
+        self.0.gemv(y, α, u, β);
+        self.1.gemv(y, α, v, F::ONE);
     }
 
-    fn apply_mut(&self, y : &mut Y, x : &Pair<U, V>){
-        self.0.apply_mut(y, &x.0);
-        self.1.apply_mut(y, &x.1);
+    fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
+        let Pair(u, v) = x.decompose();
+        self.0.apply_mut(y, u);
+        self.1.apply_mut(y, v);
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add(&self, y : &mut Y, x : &Pair<U, V>){
-        self.0.apply_add(y, &x.0);
-        self.1.apply_add(y, &x.1);
+    fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
+        let Pair(u, v) = x.decompose();
+        self.0.apply_add(y, u);
+        self.1.apply_add(y, v);
     }
 }
 
-
 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
 pub struct ColOp<S, T>(pub S, pub T);
 
-impl<A, S, T, O> Apply<A> for ColOp<S, T>
+impl<A, S, T> Mapping<A> for ColOp<S, T>
 where
-    S : for<'a> Apply<&'a A, Output=O>,
-    T : Apply<A>,
-    A : std::borrow::Borrow<A>,
+    A : Space,
+    S : Mapping<A>,
+    T : Mapping<A>,
 {
-    type Output = Pair<O, T::Output>;
+    type Codomain = Pair<S::Codomain, T::Codomain>;
 
-    fn apply(&self, a : A) -> Self::Output {
-        Pair(self.0.apply(a.borrow()), self.1.apply(a))
+    fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain {
+        Pair(self.0.apply(a.ref_instance()), self.1.apply(a))
     }
 }
 
-impl<A, S, T, D> Linear<A> for ColOp<S, T>
+impl<A, S, T> Linear<A> for ColOp<S, T>
 where
-    ColOp<S, T> : Apply<A, Output=D> + for<'a>  Apply<&'a A, Output=D>,
-{
-    type Codomain = D;
-}
+    A : Space,
+    S : Mapping<A>,
+    T : Mapping<A>,
+{ }
 
 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
 where
+    X : Space,
     S : GEMV<F, X, A>,
     T : GEMV<F, X, B>,
     F : Num,
     Self : Linear<X, Codomain=Pair<A, B>>
 {
-    fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &X, β : F) {
-        self.0.gemv(&mut y.0, α, x, β);
+    fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
+        self.0.gemv(&mut y.0, α, x.ref_instance(), β);
         self.1.gemv(&mut y.1, α, x, β);
     }
 
-    fn apply_mut(&self, y : &mut Pair<A, B>, x : &X){
-        self.0.apply_mut(&mut y.0, x);
+    fn apply_mut<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
+        self.0.apply_mut(&mut y.0, x.ref_instance());
         self.1.apply_mut(&mut y.1, x);
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add(&self, y : &mut Pair<A, B>, x : &X){
-        self.0.apply_add(&mut y.0, x);
+    fn apply_add<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
+        self.0.apply_add(&mut y.0, x.ref_instance());
         self.1.apply_add(&mut y.1, x);
     }
 }
 
 
-impl<A, B, Yʹ, R, S, T> Adjointable<Pair<A,B>,Yʹ> for RowOp<S, T>
+impl<A, B, Yʹ, S, T> Adjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
 where
+    A : Space,
+    B : Space,
+    Yʹ : Space,
     S : Adjointable<A, Yʹ>,
     T : Adjointable<B, Yʹ>,
     Self : Linear<Pair<A, B>>,
-    for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Yʹ, Codomain=R>,
+    // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
+    //     Yʹ,
+    //     Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
+    // >,
 {
-    type AdjointCodomain = R;
+    type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
     type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
 
     fn adjoint(&self) -> Self::Adjoint<'_> {
@@ -279,14 +306,21 @@
     }
 }
 
-impl<A, B, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Yʹ> for RowOp<S, T>
+impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
 where
+    A : Space,
+    B : Space,
+    Yʹ : Space,
     S : Preadjointable<A, Yʹ>,
     T : Preadjointable<B, Yʹ>,
     Self : Linear<Pair<A, B>>,
-    for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, Pair<A,B>, Codomain=R>,
+    for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
+        Yʹ, Pair<A,B>,
+        Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>,
+        AdjointCodomain = Self::Codomain,
+    >,
 {
-    type PreadjointCodomain = R;
+    type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
     type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
 
     fn preadjoint(&self) -> Self::Preadjoint<'_> {
@@ -297,10 +331,17 @@
 
 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
 where
-    S : Adjointable<A, Xʹ>,
-    T : Adjointable<A, Yʹ>,
+    A : Space,
+    Xʹ : Space,
+    Yʹ : Space,
+    R : Space + ClosedAdd,
+    S : Adjointable<A, Xʹ, AdjointCodomain = R>,
+    T : Adjointable<A, Yʹ, AdjointCodomain = R>,
     Self : Linear<A>,
-    for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>,
+    // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
+    //     Pair<Xʹ,Yʹ>,
+    //     Codomain=R,
+    // >,
 {
     type AdjointCodomain = R;
     type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
@@ -312,10 +353,18 @@
 
 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
 where
-    S : Preadjointable<A, Xʹ>,
-    T : Preadjointable<A, Yʹ>,
+    A : Space,
+    Xʹ : Space,
+    Yʹ : Space,
+    R : Space + ClosedAdd,
+    S : Preadjointable<A, Xʹ, PreadjointCodomain = R>,
+    T : Preadjointable<A, Yʹ, PreadjointCodomain = R>,
     Self : Linear<A>,
-    for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, A, Codomain=R>,
+    for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
+        Pair<Xʹ,Yʹ>, A,
+        Codomain = R,
+        AdjointCodomain = Self::Codomain,
+    >,
 {
     type PreadjointCodomain = R;
     type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
@@ -328,67 +377,73 @@
 /// Diagonal operator
 pub struct DiagOp<S, T>(pub S, pub T);
 
-impl<A, B, S, T> Apply<Pair<A, B>> for DiagOp<S, T>
+impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
 where
-    S : Apply<A>,
-    T : Apply<B>,
+    A : Space,
+    B : Space,
+    S : Mapping<A>,
+    T : Mapping<B>,
 {
-    type Output = Pair<S::Output, T::Output>;
+    type Codomain = Pair<S::Codomain, T::Codomain>;
 
-    fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output {
-        Pair(self.0.apply(a), self.1.apply(b))
-    }
-}
-
-impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for DiagOp<S, T>
-where
-    S : Apply<&'a A>,
-    T : Apply<&'a B>,
-{
-    type Output = Pair<S::Output, T::Output>;
-
-    fn apply(&self, Pair(ref a, ref b) : &'a Pair<A, B>) -> Self::Output {
+    fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
+        let Pair(a, b) = x.decompose();
         Pair(self.0.apply(a), self.1.apply(b))
     }
 }
 
-impl<A, B, S, T, D> Linear<Pair<A, B>> for DiagOp<S, T>
+impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
 where
-    DiagOp<S, T> : Apply<Pair<A, B>, Output=D> + for<'a>  Apply<&'a Pair<A, B>, Output=D>,
-{
-    type Codomain = D;
-}
+    A : Space,
+    B : Space,
+    S : Linear<A>,
+    T : Linear<B>,
+{ }
 
 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
 where
+    A : Space,
+    B : Space,
+    U : Space,
+    V : Space,
     S : GEMV<F, U, A>,
     T : GEMV<F, V, B>,
     F : Num,
-    Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>
+    Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>,
 {
-    fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &Pair<U, V>, β : F) {
-        self.0.gemv(&mut y.0, α, &x.0, β);
-        self.1.gemv(&mut y.1, α, &x.1, β);
+    fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
+        let Pair(u, v) = x.decompose();
+        self.0.gemv(&mut y.0, α, u, β);
+        self.1.gemv(&mut y.1, α, v, β);
     }
 
-    fn apply_mut(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){
-        self.0.apply_mut(&mut y.0, &x.0);
-        self.1.apply_mut(&mut y.1, &x.1);
+    fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
+        let Pair(u, v) = x.decompose();
+        self.0.apply_mut(&mut y.0, u);
+        self.1.apply_mut(&mut y.1, v);
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){
-        self.0.apply_add(&mut y.0, &x.0);
-        self.1.apply_add(&mut y.1, &x.1);
+    fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
+        let Pair(u, v) = x.decompose();
+        self.0.apply_add(&mut y.0, u);
+        self.1.apply_add(&mut y.1, v);
     }
 }
 
-impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T>
+impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
 where
+    A : Space,
+    B : Space,
+    Xʹ: Space,
+    Yʹ : Space,
+    R : Space,
     S : Adjointable<A, Xʹ>,
     T : Adjointable<B, Yʹ>,
     Self : Linear<Pair<A, B>>,
-    for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>,
+    for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
+        Pair<Xʹ,Yʹ>, Codomain=R,
+    >,
 {
     type AdjointCodomain = R;
     type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
@@ -398,12 +453,21 @@
     }
 }
 
-impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T>
+impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
 where
+    A : Space,
+    B : Space,
+    Xʹ: Space,
+    Yʹ : Space,
+    R : Space,
     S : Preadjointable<A, Xʹ>,
     T : Preadjointable<B, Yʹ>,
     Self : Linear<Pair<A, B>>,
-    for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, Pair<A, B>, Codomain=R>,
+    for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
+        Pair<Xʹ,Yʹ>, Pair<A, B>,
+        Codomain=R,
+        AdjointCodomain = Self::Codomain,
+    >,
 {
     type PreadjointCodomain = R;
     type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
@@ -416,3 +480,65 @@
 /// Block operator
 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;
 
+
+macro_rules! pairnorm {
+    ($expj:ty) => {
+        impl<F, A, B, S, T, ExpA, ExpB, ExpR>
+        BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F>
+        for RowOp<S, T>
+        where
+            F : Float,
+            A : Space,
+            B : Space,
+            S : BoundedLinear<A, ExpA, ExpR, F>,
+            T : BoundedLinear<B, ExpB, ExpR, F>,
+            S::Codomain : Add<T::Codomain>,
+            <S::Codomain as Add<T::Codomain>>::Output : Space,
+            ExpA : NormExponent,
+            ExpB : NormExponent,
+            ExpR : NormExponent,
+        {
+            fn opnorm_bound(
+                &self,
+                PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>,
+                expr : ExpR
+            ) -> F {
+                // An application of the triangle inequality bounds the norm by the maximum
+                // of the individual norms. A simple observation shows this to be exact.
+                let na = self.0.opnorm_bound(expa, expr);
+                let nb = self.1.opnorm_bound(expb, expr);
+                na.max(nb)
+            }
+        }
+        
+        impl<F, A, S, T, ExpA, ExpS, ExpT>
+        BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
+        for ColOp<S, T>
+        where
+            F : Float,
+            A : Space,
+            S : BoundedLinear<A, ExpA, ExpS, F>,
+            T : BoundedLinear<A, ExpA, ExpT, F>,
+            ExpA : NormExponent,
+            ExpS : NormExponent,
+            ExpT : NormExponent,
+        {
+            fn opnorm_bound(
+                &self,
+                expa : ExpA,
+                PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj>
+            ) -> F {
+                // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence,
+                // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖}
+                let ns = self.0.opnorm_bound(expa, exps);
+                let nt = self.1.opnorm_bound(expa, expt);
+                ns.max(nt)
+            }
+        }
+    }
+}
+
+pairnorm!(L1);
+pairnorm!(L2);
+pairnorm!(Linfinity);
+
diff -r 1a38447a89fa -r 9226980e45a7 src/loc.rs
--- a/src/loc.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/loc.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -10,9 +10,12 @@
 use crate::maputil::{FixedLength,FixedLengthMut,map1,map2,map1_mut,map2_mut};
 use crate::euclidean::*;
 use crate::norms::*;
-use crate::linops::AXPY;
+use crate::linops::{AXPY, Mapping, Linear};
+use crate::instance::{Instance, BasicDecomposition};
+use crate::mapping::Space;
 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
@@ -657,19 +660,35 @@
     }
 }
 
+
+impl<F : Num, const N : usize> Space for Loc<F, N> {
+    type Decomp = BasicDecomposition;
+}
+
+impl<F : Num, const N : usize> Mapping<Loc<F, N>> for Loc<F, N> {
+    type Codomain = F;
+
+    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+        x.eval(|x̃| self.dot(x̃))
+    }
+}
+
+impl<F : Num, const N : usize> Linear<Loc<F, N>> for Loc<F, N> { }
+
 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) })
-        }
+    fn axpy<I : Instance<Loc<F, N>>>(&mut self, α : F, x : I, β : F) {
+        x.eval(|x̃| {
+            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 )
+    fn copy_from<I : Instance<Loc<F, N>>>(&mut self, x : I) {
+        x.eval(|x̃| map2_mut(self, x̃, |yi, xi| *yi = *xi ))
     }
 }
diff -r 1a38447a89fa -r 9226980e45a7 src/mapping.rs
--- a/src/mapping.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/mapping.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -3,43 +3,20 @@
 */
 
 use std::marker::PhantomData;
+use std::borrow::Cow;
 use crate::types::Float;
 use serde::Serialize;
 use crate::loc::Loc;
-
-/// Trait for application of `Self` as a mathematical function or operator on `X`.
-pub trait Apply<X> {
-    type Output;
-
-    /// Compute the value of `self` at `x`.
-    fn apply(&self, x : X) -> Self::Output;
-}
-
-/// This blanket implementation is a workaround helper to Rust trait system limitations.
-///
-/// It is introduced because the trait system does not allow blanket implementations of both
-/// [`Apply<X>`] and [`Apply<&'a X>`]. With this, the latter is implemented automatically for
-/// the reference, which can be sufficient to apply the operation in another blanket implementation.
-impl<'a, T, X> Apply<X> for &'a T where T : Apply<X> {
-    type Output = <T as Apply<X>>::Output;
-
-    #[inline]
-    fn apply(&self, x : X) -> Self::Output {
-        (*self).apply(x)
-    }
-}
+pub use crate::instance::{Instance, Decomposition, BasicDecomposition, Space};
 
 /// A mapping from `Domain` to `Codomain`.
 ///
 /// This is automatically implemented when the relevant [`Apply`] are implemented.
-pub trait Mapping<Domain> : Apply<Domain, Output=Self::Codomain>
-                            + for<'a> Apply<&'a Domain, Output=Self::Codomain> {
-    type Codomain;
-}
+pub trait Mapping<Domain : Space> {
+    type Codomain : Space;
 
-impl<Domain, Codomain, T> Mapping<Domain> for T
-where T : Apply<Domain, Output=Codomain> + for<'a> Apply<&'a Domain, Output=Codomain> {
-    type Codomain = Codomain;
+    /// Compute the value of `self` at `x`.
+    fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain;
 }
 
 /// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc<F, N>`] to `F`.
@@ -49,15 +26,6 @@
 impl<F : Float, T, const N : usize> RealMapping<F, N> for T
 where T : Mapping<Loc<F, N>, Codomain = F> {}
 
-/// Automatically implemented shorthand for referring to differentiable [`Mapping`]s from
-/// [`Loc<F, N>`] to `F`.
-pub trait DifferentiableRealMapping<F : Float, const N : usize>
-: DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain=Loc<F, N>> {}
-
-impl<F : Float, T, const N : usize> DifferentiableRealMapping<F, N> for T
-where T : DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain=Loc<F, N>> {}
-
-
 /// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to [`Loc<F, M>`].
 pub trait RealVectorField<F : Float, const N : usize, const M : usize>
 : Mapping<Loc<F, N>, Codomain = Loc<F, M>> {}
@@ -65,69 +33,70 @@
 impl<F : Float, T, const N : usize, const M : usize> RealVectorField<F, N, M> for T
 where T : Mapping<Loc<F, N>, Codomain = Loc<F, M>> {}
 
-
-/// Trait for calculation the differential of `Self` as a mathematical function on `X`.
-pub trait Differentiable<X> : Sized {
-    type Derivative;
-
-    /// Compute the differential of `self` at `x`.
-    fn differential(&self, x : X) -> Self::Derivative;
-}
-
-impl<'g, X, G : Differentiable<X>> Differentiable<X> for &'g G {
-    type Derivative = G::Derivative;
-    #[inline]
-    fn differential(&self, x : X) -> Self::Derivative {
-        (*self).differential(x)
-    }
-}
-
 /// A differentiable mapping from `Domain` to [`Mapping::Codomain`], with differentials
 /// `Differential`.
 ///
-/// This is automatically implemented when the relevant [`Differentiable`] are implemented.
-pub trait DifferentiableMapping<Domain>
-: Mapping<Domain>
-  + Differentiable<Domain, Derivative=Self::DerivativeDomain>
-  + for<'a> Differentiable<&'a Domain, Derivative=Self::DerivativeDomain> {
-    type DerivativeDomain;
-    type Differential : Mapping<Domain, Codomain=Self::DerivativeDomain>;
-    type DifferentialRef<'b> : Mapping<Domain, Codomain=Self::DerivativeDomain> where Self : 'b;
+/// This is automatically implemented when [`DifferentiableImpl`] is.
+pub trait DifferentiableMapping<Domain : Space> : Mapping<Domain> {
+    type DerivativeDomain : Space;
+    type Differential<'b> : Mapping<Domain, Codomain=Self::DerivativeDomain> where Self : 'b;
+
+    /// Calculate differential at `x`
+    fn differential<I : Instance<Domain>>(&self, x : I) -> Self::DerivativeDomain;
 
     /// Form the differential mapping of `self`.
-    fn diff(self) -> Self::Differential;
+    fn diff(self) -> Self::Differential<'static>;
+
     /// Form the differential mapping of `self`.
-    fn diff_ref(&self) -> Self::DifferentialRef<'_>;
+    fn diff_ref(&self) -> Self::Differential<'_>;
 }
 
+/// Automatically implemented shorthand for referring to differentiable [`Mapping`]s from
+/// [`Loc<F, N>`] to `F`.
+pub trait DifferentiableRealMapping<F : Float, const N : usize>
+: DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain = Loc<F, N>> {}
 
-impl<Domain, Derivative, T> DifferentiableMapping<Domain> for T
-where T : Mapping<Domain>
-          + Differentiable<Domain, Derivative=Derivative>
-          + for<'a> Differentiable<&'a Domain,Derivative=Derivative> {
-    type DerivativeDomain = Derivative;
-    type Differential = Differential<Domain, Self>;
-    type DifferentialRef<'b> = Differential<Domain, &'b Self> where Self : 'b;
+impl<F : Float, T, const N : usize> DifferentiableRealMapping<F, N> for T
+where T : DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain = Loc<F, N>> {}
+
+/// Helper trait for implementing [`DifferentiableMapping`]
+pub trait DifferentiableImpl<X : Space> : Sized {
+    type Derivative : Space;
+
+    /// Compute the differential of `self` at `x`, consuming the input.
+    fn differential_impl<I : Instance<X>>(&self, x : I) -> Self::Derivative;
+}
+
+impl<T, Domain> DifferentiableMapping<Domain> for T
+where
+    Domain : Space,
+    T : Clone + Mapping<Domain> + DifferentiableImpl<Domain>
+{
+    type DerivativeDomain = T::Derivative;
+    type Differential<'b> = Differential<'b, Domain, Self> where Self : 'b;
     
-    /// Form the differential mapping of `self`.
-    fn diff(self) -> Self::Differential {
-        Differential{ g : self, _space : PhantomData }
+    #[inline]
+    fn differential<I : Instance<Domain>>(&self, x : I) -> Self::DerivativeDomain {
+        self.differential_impl(x)
     }
 
-    /// Form the differential mapping of `self`.
-    fn diff_ref(&self) -> Self::DifferentialRef<'_> {
-        Differential{ g : self, _space : PhantomData }
+    fn diff(self) -> Differential<'static, Domain, Self> {
+        Differential{ g : Cow::Owned(self), _space : PhantomData }
+    }
+
+    fn diff_ref(&self) -> Differential<'_, Domain, Self> {
+        Differential{ g : Cow::Borrowed(self), _space : PhantomData }
     }
 }
 
 /// A sum of [`Mapping`]s.
 #[derive(Serialize, Debug, Clone)]
-pub struct Sum<Domain, M : Mapping<Domain>> {
+pub struct Sum<Domain, M> {
     components : Vec<M>,
     _domain : PhantomData<Domain>,
 }
 
-impl<Domain, M : Mapping<Domain>> Sum<Domain, M> {
+impl<Domain, M> Sum<Domain, M> {
     /// Construct from an iterator.
     pub fn new<I : Iterator<Item = M>>(iter : I) -> Self {
         Sum { components : iter.collect(), _domain : PhantomData }
@@ -140,123 +109,107 @@
 }
 
 
-impl<Domain, M> Apply<Domain> for Sum<Domain, M>
-where M : Mapping<Domain>,
-      M::Codomain : std::iter::Sum {
-    type Output = M::Codomain;
+impl<Domain, M> Mapping<Domain> for Sum<Domain, M>
+where
+    Domain : Space + Clone,
+    M : Mapping<Domain>,
+    M::Codomain : std::iter::Sum + Clone
+{
+    type Codomain = M::Codomain;
 
-    fn apply(&self, x : Domain) -> Self::Output {
-        self.components.iter().map(|c| c.apply(&x)).sum()
+    fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain {
+        let xr = x.ref_instance();
+        self.components.iter().map(|c| c.apply(xr)).sum()
     }
 }
 
-impl<'a, Domain, M> Apply<&'a Domain> for Sum<Domain, M>
-where M : Mapping<Domain>,
-      M::Codomain : std::iter::Sum {
-    type Output = M::Codomain;
-
-    fn apply(&self, x : &'a Domain) -> Self::Output {
-        self.components.iter().map(|c| c.apply(x)).sum()
-    }
-}
-
-impl<Domain, M> Differentiable<Domain> for Sum<Domain, M>
-where M : DifferentiableMapping<Domain>,
-      M :: Codomain : std::iter::Sum,
-      M :: DerivativeDomain : std::iter::Sum,
-      Domain : Copy {
-
+impl<Domain, M> DifferentiableImpl<Domain> for Sum<Domain, M>
+where
+    Domain : Space + Clone,
+    M : DifferentiableMapping<Domain>,
+    M :: DerivativeDomain : std::iter::Sum
+{
     type Derivative = M::DerivativeDomain;
 
-    fn differential(&self, x : Domain) -> Self::Derivative {
-        self.components.iter().map(|c| c.differential(x)).sum()
+    fn differential_impl<I : Instance<Domain>>(&self, x : I) -> Self::Derivative {
+        let xr = x.ref_instance();
+        self.components.iter().map(|c| c.differential(xr)).sum()
     }
 }
 
 /// Container for the differential [`Mapping`] of a [`Differentiable`] mapping.
-pub struct Differential<X, G : DifferentiableMapping<X>> {
-    g : G,
+pub struct Differential<'a, X, G : Clone> {
+    g : Cow<'a, G>,
     _space : PhantomData<X>
 }
 
-impl<X, G : DifferentiableMapping<X>> Differential<X, G> {
+impl<'a, X, G : Clone> Differential<'a, X, G> {
     pub fn base_fn(&self) -> &G {
         &self.g
     }
 }
 
-impl<X, G : DifferentiableMapping<X>> Apply<X> for Differential<X, G> {
-    type Output = G::DerivativeDomain;
+impl<'a, X, G> Mapping<X> for Differential<'a, X, G>
+where
+    X : Space,
+    G : Clone + DifferentiableMapping<X>
+{
+    type Codomain = G::DerivativeDomain;
 
     #[inline]
-    fn apply(&self, x : X) -> Self::Output {
-        self.g.differential(x)
+    fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain {
+        (*self.g).differential(x)
     }
 }
 
-impl<'a, X, G : DifferentiableMapping<X>> Apply<&'a X> for Differential<X, G> {
-    type Output = G::DerivativeDomain;
-
-    #[inline]
-    fn apply(&self, x : &'a X) -> Self::Output {
-        self.g.differential(x)
-    }
-}
-
-
 /// Container for flattening [`Loc`]`<F, 1>` codomain of a [`Mapping`] to `F`.
-pub struct FlattenedCodomain<X, F, G : Mapping<X, Codomain=Loc<F, 1>>> {
+pub struct FlattenedCodomain<X, F, G> {
     g : G,
     _phantoms : PhantomData<(X, F)>
 }
 
-impl<X, F, G : Mapping<X, Codomain=Loc<F, 1>>> Apply<X> for FlattenedCodomain<X, F, G> {
-    type Output = F;
+impl<F : Space, X, G> Mapping<X> for FlattenedCodomain<X, F, G>
+where
+    X : Space,
+    G: Mapping<X, Codomain=Loc<F, 1>>
+{
+    type Codomain = F;
 
     #[inline]
-    fn apply(&self, x : X) -> Self::Output {
+    fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain {
         self.g.apply(x).flatten1d()
     }
 }
 
 /// An auto-trait for constructing a [`FlattenCodomain`] structure for
 /// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`.
-pub trait FlattenCodomain<X, F> : Mapping<X, Codomain=Loc<F, 1>> + Sized {
+pub trait FlattenCodomain<X : Space, F> : Mapping<X, Codomain=Loc<F, 1>> + Sized {
     /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`.
     fn flatten_codomain(self) -> FlattenedCodomain<X, F, Self> {
         FlattenedCodomain{ g : self, _phantoms : PhantomData }
     }
 }
 
-impl<X, F, G : Sized + Mapping<X, Codomain=Loc<F, 1>>> FlattenCodomain<X, F> for G {}
-
+impl<X : Space, F, G : Sized + Mapping<X, Codomain=Loc<F, 1>>> FlattenCodomain<X, F> for G {}
 
 /// Container for dimensional slicing [`Loc`]`<F, N>` codomain of a [`Mapping`] to `F`.
-pub struct SlicedCodomain<X, F, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> {
-    g : G,
+pub struct SlicedCodomain<'a, X, F, G : Clone, const N : usize> {
+    g : Cow<'a, G>,
     slice : usize,
     _phantoms : PhantomData<(X, F)>
 }
 
-impl<X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> Apply<X>
-for SlicedCodomain<X, F, G, N> {
-    type Output = F;
+impl<'a, X, F, G, const N : usize> Mapping<X> for SlicedCodomain<'a, X, F, G, N>
+where
+    X : Space,
+    F : Copy + Space,
+    G : Mapping<X, Codomain=Loc<F, N>> + Clone,
+{
+    type Codomain = F;
 
     #[inline]
-    fn apply(&self, x : X) -> Self::Output {
-        let tmp : [F; N] = self.g.apply(x).into();
-        // Safety: `slice_codomain` below checks the range.
-        unsafe { *tmp.get_unchecked(self.slice) }
-    }
-}
-
-impl<'a, X, F : Copy, G : Mapping<X, Codomain=Loc<F, N>>, const N : usize> Apply<&'a X>
-for SlicedCodomain<X, F, G, N> {
-    type Output = F;
-
-    #[inline]
-    fn apply(&self, x : &'a X) -> Self::Output {
-        let tmp : [F; N] = self.g.apply(x).into();
+    fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain {
+        let tmp : [F; N] = (*self.g).apply(x).into();
         // Safety: `slice_codomain` below checks the range.
         unsafe { *tmp.get_unchecked(self.slice) }
     }
@@ -264,21 +217,22 @@
 
 /// An auto-trait for constructing a [`FlattenCodomain`] structure for
 /// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`.
-pub trait SliceCodomain<X, F : Copy, const N : usize> : Mapping<X, Codomain=Loc<F, N>> + Sized {
+pub trait SliceCodomain<X : Space, F : Copy, const N : usize>
+    : Mapping<X, Codomain=Loc<F, N>> + Clone + Sized
+{
     /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`.
-    fn slice_codomain(self, slice : usize) -> SlicedCodomain<X, F, Self, N> {
+    fn slice_codomain(self, slice : usize) -> SlicedCodomain<'static, X, F, Self, N> {
         assert!(slice < N);
-        SlicedCodomain{ g : self, slice, _phantoms : PhantomData }
+        SlicedCodomain{ g : Cow::Owned(self), slice, _phantoms : PhantomData }
     }
 
     /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`.
-    fn slice_codomain_ref(&self, slice : usize) -> SlicedCodomain<X, F, &'_ Self, N> {
+    fn slice_codomain_ref(&self, slice : usize) -> SlicedCodomain<'_, X, F, Self, N> {
         assert!(slice < N);
-        SlicedCodomain{ g : self, slice, _phantoms : PhantomData }
+        SlicedCodomain{ g : Cow::Borrowed(self), slice, _phantoms : PhantomData }
     }
 }
 
-impl<X, F : Copy, G : Sized + Mapping<X, Codomain=Loc<F, N>>, const N : usize>
+impl<X : Space, F : Copy, G : Sized + Mapping<X, Codomain=Loc<F, N>> + Clone, const N : usize>
 SliceCodomain<X, F, N>
 for G {}
-
diff -r 1a38447a89fa -r 9226980e45a7 src/metaprogramming.rs
--- a/src/metaprogramming.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/metaprogramming.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -41,63 +41,3 @@
     (noref, &$x:ty) => { $x };
 }
 
-/// Use as
-/// ```ignore
-/// impl_vectorspace_ops!(impl_binop, impl_assignop, impl_scalarop, impl_scalar_assignop, impl_unaryop);
-/// ```
-/// with `impl_binop`, `impl_assignop`, `impl_scalarop`, and `impl_unaryop` macros provided.
-/// For example usage see the [power][crate::vectorspace::powerspace] and
-/// [product spaces][crate::vectorspace::productspace] implementations.
-#[macro_export]
-macro_rules! impl_vectorspace_ops {
-    ($impl_binop:ident, $impl_assignop:ident, $impl_scalarop:ident, $impl_scalarlhs_op:ident,
-        $impl_scalar_assignop:ident, $impl_unaryop:ident) => {
-        impl_vectorspace_ops!($impl_binop, $impl_assignop, $impl_scalarop, $impl_scalarlhs_op,
-                                $impl_scalar_assignop, $impl_unaryop,
-                                (f32, f64,
-                                num::complex::Complex<f32>,
-                                num::complex::Complex<f64>));
-    };
-    ($impl_binop:ident, $impl_assignop:ident, $impl_scalarop:ident, $impl_scalarlhs_op:ident,
-        $impl_scalar_assignop:ident, $impl_unaryop:ident, ($($field:ty),+)) => {
-        impl_vectorspace_ops!(@binary, $impl_binop, Add, add);
-        impl_vectorspace_ops!(@binary, $impl_binop, Sub, sub);
-        impl_vectorspace_ops!(@assign, $impl_assignop, AddAssign, add_assign);
-        impl_vectorspace_ops!(@assign, $impl_assignop, SubAssign, sub_assign);
-        impl_vectorspace_ops!(@scalar, $impl_scalarop, Mul, mul);
-        impl_vectorspace_ops!(@scalar, $impl_scalarop, Div, div);
-        // Compiler overflow
-        // $(
-        //     impl_vectorspace_ops!(@scalar_lhs, $impl_scalarlhs_op, Mul, mul, $field);
-        // )*
-        impl_vectorspace_ops!(@scalar_assign, $impl_scalar_assignop, MulAssign, mul_assign);
-        impl_vectorspace_ops!(@scalar_assign, $impl_scalar_assignop, DivAssign, div_assign);
-        impl_vectorspace_ops!(@unary, $impl_unaryop, Neg, neg);
-    };
-    (@binary, $impl:ident, $trait : ident, $fn : ident) => {
-        $impl!($trait, $fn, ref, ref);
-        $impl!($trait, $fn, ref, noref);
-        $impl!($trait, $fn, noref, ref);
-        $impl!($trait, $fn, noref, noref);
-    };
-    (@assign, $impl:ident, $trait : ident, $fn :ident) => {
-        $impl!($trait, $fn, ref);
-        $impl!($trait, $fn, noref);
-    };
-    (@scalar, $impl:ident, $trait : ident, $fn :ident) => {
-        $impl!($trait, $fn, ref);
-        $impl!($trait, $fn, noref);
-    };
-    (@scalar_lhs, $impl:ident, $trait : ident, $fn : ident, $field: ty) => {
-        // These operators need workarounds
-        $impl!($trait, $fn, ref, $field);
-        $impl!($trait, $fn, noref, $field);
-    };
-    (@scalar_assign, $impl:ident, $trait : ident, $fn :ident) => {
-        $impl!($trait, $fn);
-    };
-    (@unary, $impl:ident, $trait : ident, $fn :ident) => {
-        $impl!($trait, $fn, ref);
-        $impl!($trait, $fn, noref);
-    };
-}
diff -r 1a38447a89fa -r 9226980e45a7 src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/nalgebra_support.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -23,51 +23,49 @@
 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 + Zero + One + ClosedAddAssign + ClosedMulAssign,
-        DefaultAllocator : Allocator<N,K>,
-        DefaultAllocator : Allocator<M,K>,
-        DefaultAllocator : Allocator<N,M>,
-        DefaultAllocator : Allocator<M,N> {
-    type Output = OMatrix<E,N,K>;
-
-    #[inline]
-    fn apply(&self, x : Matrix<E,M,K,SV>) -> Self::Output {
-        self.mul(x)
-    }
+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<'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 + Zero + One + ClosedAddAssign + ClosedMulAssign,
-        DefaultAllocator : Allocator<N,K>,
-        DefaultAllocator : Allocator<M,K>,
-        DefaultAllocator : Allocator<N,M>,
-        DefaultAllocator : Allocator<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>,
+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<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> Linear<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> {
 }
 
 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>,
+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>,
@@ -75,34 +73,36 @@
       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>,
+where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone,
       M : Dim, E : Scalar + Zero + One + Float,
       DefaultAllocator : Allocator<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(&mut self, y : &Vector<E,M,SV1>) {
-        Matrix::copy_from(self, y)
+    fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) {
+        y.eval(|ỹ| Matrix::copy_from(self, ỹ))
     }
 }
 
 impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM>
-where SM: StorageMut<E,M>,
+where SM: StorageMut<E,M> + Clone,
       M : Dim, E : Scalar + Zero + One + Float + RealField,
       DefaultAllocator : Allocator<M> {
     #[inline]
@@ -111,9 +111,9 @@
     }
 }
 
-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>,
+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>,
@@ -170,7 +170,7 @@
 impl<E,M,S> Euclidean<E>
 for Vector<E,M,S>
 where M : Dim,
-      S : StorageMut<E,M>,
+      S : StorageMut<E,M> + Clone,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
 
@@ -195,7 +195,7 @@
 impl<E,M,S> StaticEuclidean<E>
 for Vector<E,M,S>
 where M : DimName,
-      S : StorageMut<E,M>,
+      S : StorageMut<E,M> + Clone,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
 
diff -r 1a38447a89fa -r 9226980e45a7 src/norms.rs
--- a/src/norms.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/norms.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -3,18 +3,31 @@
 */
 
 use serde::Serialize;
+use std::marker::PhantomData;
 use crate::types::*;
 use crate::euclidean::*;
+use crate::mapping::{Mapping, Space, Instance};
 
 //
 // Abstract norms
 //
 
+#[derive(Copy,Clone,Debug)]
+/// Helper structure to convert a [`NormExponent`] into a [`Mapping`]
+pub struct NormMapping<F : Float, E : NormExponent>{
+    pub(crate) exponent : E,
+    _phantoms : PhantomData<F>
+}
+
 /// An exponent for norms.
 ///
 // Just a collection of desirable attributes for a marker type
-pub trait NormExponent : Copy + Send + Sync + 'static {}
-
+pub trait NormExponent : Copy + Send + Sync + 'static {
+    /// Return the norm as a mappin
+    fn as_mapping<F : Float>(self) -> NormMapping<F, Self> {
+        NormMapping{ exponent : self, _phantoms : PhantomData }
+    }
+}
 
 /// Exponent type for the 1-[`Norm`].
 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
@@ -37,6 +50,15 @@
 pub struct L21;
 impl NormExponent for L21 {}
 
+/// Norms for pairs (a, b). ‖(a,b)‖ = ‖(‖a‖_A, ‖b‖_B)‖_J
+/// For use with [`crate::direct_product::Pair`]
+#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
+pub struct PairNorm<A, B, J>(pub A, pub B, pub J);
+
+impl<A, B, J> NormExponent for PairNorm<A, B, J>
+where A : NormExponent, B : NormExponent, J : NormExponent {}
+
+
 /// A Huber/Moreau–Yosida smoothed [`L1`] norm. (Not a norm itself.)
 ///
 /// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher
@@ -154,3 +176,29 @@
     }
 }
 
+impl<F : Float, E : Norm<F, L2>> Norm<F, L21> for Vec<E> {
+    fn norm(&self, _l21 : L21) -> F {
+        self.iter().map(|e| e.norm(L2)).sum()
+    }
+}
+
+impl<F : Float, E : Dist<F, L2>> Dist<F, L21> for Vec<E> {
+    fn dist(&self, other : &Self, _l21 : L21) -> F {
+        self.iter().zip(other.iter()).map(|(e, g)| e.dist(g, L2)).sum()
+    }
+}
+
+impl<E, F, Domain> Mapping<Domain> for NormMapping<F, E>
+where
+    F : Float,
+    E : NormExponent,
+    Domain : Space + Norm<F, E>,
+{
+    type Codomain = F;
+
+    #[inline]
+    fn apply<I : Instance<Domain>>(&self, x : I) -> F {
+        x.eval(|r| r.norm(self.exponent))
+    }
+}
+
diff -r 1a38447a89fa -r 9226980e45a7 src/types.rs
--- a/src/types.rs	Sat Dec 14 09:31:27 2024 -0500
+++ b/src/types.rs	Tue Dec 31 08:30:02 2024 -0500
@@ -12,6 +12,13 @@
 //use trait_set::trait_set;
 pub use num_traits::Float as NumTraitsFloat; // needed to re-export functions.
 pub use num_traits::cast::AsPrimitive;
+pub use simba::scalar::{
+    ClosedAdd, ClosedAddAssign,
+    ClosedSub, ClosedSubAssign,
+    ClosedMul, ClosedMulAssign,
+    ClosedDiv, ClosedDivAssign,
+    ClosedNeg
+};
 
 /// Typical integer type
 #[allow(non_camel_case_types)]
@@ -57,7 +64,8 @@
                  + CastFrom<u128> + CastFrom<usize>
                  + CastFrom<i8>   + CastFrom<i16> + CastFrom<i32> + CastFrom<i64>
                  + CastFrom<i128> + CastFrom<isize>
-                 + CastFrom<f32>  + CastFrom<f64> {
+                 + CastFrom<f32>  + CastFrom<f64>
+                 + crate::instance::Space {
 
     const ZERO : Self;
     const ONE : Self;