# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1745846764 18000
# Node ID 9e5b9fc81c52c9641d8226c6a1ca4ff1b335f12d
# Parent  2e4517b55442207ff49474d456df3e61b8fc0633
Quadratic Mappings; Lipschitz trait (moved from pointsource_algs).

diff -r 2e4517b55442 -r 9e5b9fc81c52 src/linops.rs
--- a/src/linops.rs	Sun Apr 27 20:41:36 2025 -0500
+++ b/src/linops.rs	Mon Apr 28 08:26:04 2025 -0500
@@ -2,38 +2,46 @@
 Abstract linear operators.
 */
 
-use numeric_literals::replace_float_literals;
-use std::marker::PhantomData;
-use serde::Serialize;
-use crate::types::*;
-pub use crate::mapping::{Mapping, Space, Composition};
 use crate::direct_product::Pair;
 use crate::instance::Instance;
-use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm};
+pub use crate::mapping::{Composition, DifferentiableImpl, Mapping, Space};
+use crate::norms::{Linfinity, Norm, NormExponent, PairNorm, L1, L2};
+use crate::types::*;
+use numeric_literals::replace_float_literals;
+use serde::Serialize;
+use std::marker::PhantomData;
 
 /// Trait for linear operators on `X`.
-pub trait Linear<X : Space> : Mapping<X>
-{ }
+pub trait Linear<X: Space>: Mapping<X> {}
+
+// impl<X: Space, A: Linear<X>> DifferentiableImpl<X> for A {
+//     type Derivative = <Self as Mapping<X>>::Codomain;
+
+//     /// Compute the differential of `self` at `x`, consuming the input.
+//     fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
+//         self.apply(x)
+//     }
+// }
 
 /// Efficient in-place summation.
 #[replace_float_literals(F::cast_from(literal))]
-pub trait AXPY<F, X = Self> : Space + std::ops::MulAssign<F>
+pub trait AXPY<F, X = Self>: Space + std::ops::MulAssign<F>
 where
-    F : Num,
-    X : Space,
+    F: Num,
+    X: Space,
 {
-    type Owned : AXPY<F, X>;
+    type Owned: AXPY<F, X>;
 
     /// Computes  `y = βy + αx`, where `y` is `Self`.
-    fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F);
+    fn axpy<I: Instance<X>>(&mut self, α: F, x: I, β: F);
 
     /// Copies `x` to `self`.
-    fn copy_from<I : Instance<X>>(&mut self, x : I) {
+    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<I : Instance<X>>(&mut self, α : F, x : I) {
+    fn scale_from<I: Instance<X>>(&mut self, α: F, x: I) {
         self.axpy(α, x, 0.0)
     }
 
@@ -46,37 +54,36 @@
 
 /// Efficient in-place application for [`Linear`] operators.
 #[replace_float_literals(F::cast_from(literal))]
-pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<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<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F);
+    fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F);
 
     #[inline]
     /// Computes `y = Ax`, where `A` is `Self`
-    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
+    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
         self.gemv(y, 1.0, x, 0.0)
     }
 
     #[inline]
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){
+    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, XExp, CodExp, F = f64> : Linear<X>
+pub trait BoundedLinear<X, XExp, CodExp, F = f64>: Linear<X>
 where
-    F : Num,
-    X : Space + Norm<F, XExp>,
-    XExp : NormExponent,
-    CodExp : NormExponent
+    F: Num,
+    X: Space + Norm<F, XExp>,
+    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. The [`NormExponent`] `xexp`  indicates the norm
     /// in `X`, and `codexp` in the codomain.
-    fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F;
+    fn opnorm_bound(&self, xexp: XExp, codexp: CodExp) -> F;
 }
 
 // Linear operator application into mutable target. The [`AsRef`] bound
@@ -90,13 +97,15 @@
 }*/
 
 /// Trait for forming the adjoint operator of `Self`.
-pub trait Adjointable<X, Yʹ> : Linear<X>
+pub trait Adjointable<X, Yʹ>: Linear<X>
 where
-    X : Space,
-    Yʹ : Space,
+    X: Space,
+    Yʹ: Space,
 {
-    type AdjointCodomain : Space;
-    type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
+    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<'_>;
@@ -112,144 +121,168 @@
 /// We do not make additional restrictions on `Self::Preadjoint` (in particular, it
 /// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint
 /// have the full space as the codomain, etc.
-pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> {
-    type PreadjointCodomain : Space;
-    type Preadjoint<'a> : Linear<
-        Ypre, Codomain=Self::PreadjointCodomain
-    > where Self : 'a;
+pub trait Preadjointable<X: Space, Ypre: Space = <Self as Mapping<X>>::Codomain>:
+    Linear<X>
+{
+    type PreadjointCodomain: Space;
+    type Preadjoint<'a>: Linear<Ypre, 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 : 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> {}
+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>);
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
+pub struct IdOp<X>(PhantomData<X>);
 
 impl<X> IdOp<X> {
-    pub fn new() -> IdOp<X> { IdOp(PhantomData) }
+    pub fn new() -> IdOp<X> {
+        IdOp(PhantomData)
+    }
 }
 
-impl<X : Clone + Space> Mapping<X> for IdOp<X> {
+impl<X: Clone + Space> Mapping<X> for IdOp<X> {
     type Codomain = X;
 
-    fn apply<I : Instance<X>>(&self, x : I) -> X {
+    fn apply<I: Instance<X>>(&self, x: I) -> X {
         x.own()
     }
 }
 
-impl<X : Clone + Space> Linear<X> for IdOp<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>
+impl<F: Num, X, Y> GEMV<F, X, Y> for IdOp<X>
 where
-    Y : AXPY<F, X>,
-    X : Clone + Space
+    Y: AXPY<F, X>,
+    X: Clone + Space,
 {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
+    fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
         y.axpy(α, x, β)
     }
 
-    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
+    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
         y.copy_from(x);
     }
 }
 
 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
 where
-    X : Space + Clone + Norm<F, E>,
-    F : Num,
-    E : NormExponent
+    X: Space + Clone + Norm<F, E>,
+    F: Num,
+    E: NormExponent,
 {
-    fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE }
-}
-
-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() }
+    fn opnorm_bound(&self, _xexp: E, _codexp: E) -> F {
+        F::ONE
+    }
 }
 
-impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> {
-    type PreadjointCodomain=X;
-    type Preadjoint<'a> = IdOp<X> where X : 'a;
+impl<X: Clone + Space> Adjointable<X, X> for IdOp<X> {
+    type AdjointCodomain = X;
+    type Adjoint<'a>
+        = IdOp<X>
+    where
+        X: 'a;
 
-    fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() }
+    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()
+    }
+}
 
 /// The zero operator
-#[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
 pub struct ZeroOp<'a, X, XD, Y, F> {
-    zero : &'a Y, // TODO: don't pass this in `new`; maybe not even store.
-    dual_or_predual_zero : XD,
-    _phantoms : PhantomData<(X, Y, F)>,
+    zero: &'a Y, // TODO: don't pass this in `new`; maybe not even store.
+    dual_or_predual_zero: XD,
+    _phantoms: PhantomData<(X, Y, F)>,
 }
 
 // TODO: Need to make Zero in Instance.
 
-impl<'a, F : Num, X : Space, XD, Y : Space + Clone> ZeroOp<'a, X, XD, Y, F> {
-    pub fn new(zero : &'a Y, dual_or_predual_zero : XD) -> Self {
-        ZeroOp{ zero, dual_or_predual_zero, _phantoms : PhantomData }
+impl<'a, F: Num, X: Space, XD, Y: Space + Clone> ZeroOp<'a, X, XD, Y, F> {
+    pub fn new(zero: &'a Y, dual_or_predual_zero: XD) -> Self {
+        ZeroOp {
+            zero,
+            dual_or_predual_zero,
+            _phantoms: PhantomData,
+        }
     }
 }
 
-impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> {
+impl<'a, F: Num, X: Space, XD, Y: AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> {
     type Codomain = Y;
 
-    fn apply<I : Instance<X>>(&self, _x : I) -> Y {
+    fn apply<I: Instance<X>>(&self, _x: I) -> Y {
         self.zero.clone()
     }
 }
 
-impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F>
-{ }
+impl<'a, F: Num, X: Space, XD, Y: AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F> {}
 
 #[replace_float_literals(F::cast_from(literal))]
 impl<'a, F, X, XD, Y> GEMV<F, X, Y> for ZeroOp<'a, X, XD, Y, F>
 where
-    F : Num,
-    Y : AXPY<F, Y> + Clone,
-    X : Space
+    F: Num,
+    Y: AXPY<F, Y> + Clone,
+    X: Space,
 {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv<I : Instance<X>>(&self, y : &mut Y, _α : F, _x : I, β : F) {
+    fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
         *y *= β;
     }
 
-    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, _x : I){
+    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
         y.set_zero();
     }
 }
 
 impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F>
 where
-    X : Space + Norm<F, E1>,
-    Y : AXPY<F> + Clone + Norm<F, E2>,
-    F : Num,
-    E1 : NormExponent,
-    E2 : NormExponent,
+    X: Space + Norm<F, E1>,
+    Y: AXPY<F> + Clone + Norm<F, E2>,
+    F: Num,
+    E1: NormExponent,
+    E2: NormExponent,
 {
-    fn opnorm_bound(&self, _xexp : E1, _codexp : E2) -> F { F::ZERO }
+    fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> F {
+        F::ZERO
+    }
 }
 
-impl<'a, F : Num, X, XD, Y, Yprime : Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F>
+impl<'a, F: Num, X, XD, Y, Yprime: Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F>
 where
-     X : Space,
-     Y :  AXPY<F> + Clone + 'static,
-     XD : AXPY<F> + Clone + 'static,
+    X: Space,
+    Y: AXPY<F> + Clone + 'static,
+    XD: AXPY<F> + Clone + 'static,
 {
     type AdjointCodomain = XD;
-    type Adjoint<'b> = ZeroOp<'b, Yprime, (), XD, F> where Self : 'b;
-        // () means not (pre)adjointable.
+    type Adjoint<'b>
+        = ZeroOp<'b, Yprime, (), XD, F>
+    where
+        Self: 'b;
+    // () means not (pre)adjointable.
 
     fn adjoint(&self) -> Self::Adjoint<'_> {
         ZeroOp::new(&self.dual_or_predual_zero, ())
@@ -258,15 +291,18 @@
 
 impl<'a, F, X, XD, Y, Ypre> Preadjointable<X, Ypre> for ZeroOp<'a, X, XD, Y, F>
 where
-    F : Num,
-    X : Space,
-    Y : AXPY<F> + Clone,
-    Ypre : Space,
-    XD : AXPY<F> + Clone + 'static,
+    F: Num,
+    X: Space,
+    Y: AXPY<F> + Clone,
+    Ypre: Space,
+    XD: AXPY<F> + Clone + 'static,
 {
     type PreadjointCodomain = XD;
-    type Preadjoint<'b> = ZeroOp<'b, Ypre, (), XD, F> where Self : 'b;
-        // () means not (pre)adjointable.
+    type Preadjoint<'b>
+        = ZeroOp<'b, Ypre, (), XD, F>
+    where
+        Self: 'b;
+    // () means not (pre)adjointable.
 
     fn preadjoint(&self) -> Self::Preadjoint<'_> {
         ZeroOp::new(&self.dual_or_predual_zero, ())
@@ -275,45 +311,46 @@
 
 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
 where
-    X : Space,
-    T : Linear<X>,
-    S : Linear<T::Codomain>
-{ }
+    X: Space,
+    T: Linear<X>,
+    S: Linear<T::Codomain>,
+{
+}
 
 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
 where
-    F : Num,
-    X : Space,
-    T : Linear<X>,
-    S : GEMV<F, T::Codomain, Y>,
+    F: Num,
+    X: Space,
+    T: Linear<X>,
+    S: GEMV<F, T::Codomain, Y>,
 {
-    fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
+    fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
         self.outer.gemv(y, α, self.inner.apply(x), β)
     }
 
     /// Computes `y = Ax`, where `A` is `Self`
-    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
+    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
         self.outer.apply_mut(y, self.inner.apply(x))
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){
+    fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) {
         self.outer.apply_add(y, self.inner.apply(x))
     }
 }
 
 impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp>
 where
-    F : Num,
-    X : Space + Norm<F, Xexp>,
-    Z : Space + Norm<F, Zexp>,
-    Xexp : NormExponent,
-    Yexp : NormExponent,
-    Zexp : NormExponent,
-    T : BoundedLinear<X, Xexp, Zexp, F, Codomain=Z>,
-    S : BoundedLinear<Z, Zexp, Yexp, F>,
+    F: Num,
+    X: Space + Norm<F, Xexp>,
+    Z: Space + Norm<F, Zexp>,
+    Xexp: NormExponent,
+    Yexp: NormExponent,
+    Zexp: NormExponent,
+    T: BoundedLinear<X, Xexp, Zexp, F, Codomain = Z>,
+    S: BoundedLinear<Z, Zexp, Yexp, F>,
 {
-    fn opnorm_bound(&self, xexp : Xexp, yexp : Yexp) -> F {
+    fn opnorm_bound(&self, xexp: Xexp, yexp: Yexp) -> F {
         let zexp = self.intermediate_norm_exponent;
         self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp)
     }
@@ -326,17 +363,16 @@
 
 impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
 where
-    A : Space,
-    B : Space,
-    S : Mapping<A>,
-    T : Mapping<B>,
-    S::Codomain : Add<T::Codomain>,
-    <S::Codomain as Add<T::Codomain>>::Output : Space,
-
+    A: Space,
+    B: Space,
+    S: Mapping<A>,
+    T: Mapping<B>,
+    S::Codomain: Add<T::Codomain>,
+    <S::Codomain as Add<T::Codomain>>::Output: Space,
 {
     type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
 
-    fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
+    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)
     }
@@ -344,38 +380,38 @@
 
 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
 where
-    A : Space,
-    B : Space,
-    S : Linear<A>,
-    T : Linear<B>,
-    S::Codomain : Add<T::Codomain>,
-    <S::Codomain as Add<T::Codomain>>::Output : Space,
-{ }
-
+    A: Space,
+    B: Space,
+    S: Linear<A>,
+    T: Linear<B>,
+    S::Codomain: Add<T::Codomain>,
+    <S::Codomain as Add<T::Codomain>>::Output: Space,
+{
+}
 
 impl<'b, 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>
+    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<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) {
+    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<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
+    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_add(y, v);
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
+    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);
@@ -387,129 +423,137 @@
 
 impl<A, S, T> Mapping<A> for ColOp<S, T>
 where
-    A : Space,
-    S : Mapping<A>,
-    T : Mapping<A>,
+    A: Space,
+    S: Mapping<A>,
+    T: Mapping<A>,
 {
     type Codomain = Pair<S::Codomain, T::Codomain>;
 
-    fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain {
+    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> Linear<A> for ColOp<S, T>
 where
-    A : Space,
-    S : Mapping<A>,
-    T : Mapping<A>,
-{ }
+    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>>
+    X: Space,
+    S: GEMV<F, X, A>,
+    T: GEMV<F, X, B>,
+    F: Num,
+    Self: Linear<X, Codomain = Pair<A, B>>,
 {
-    fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
+    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<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
+    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<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
+    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ʹ, 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>>,
+    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=Pair<S::AdjointCodomain, T::AdjointCodomain>
     // >,
 {
     type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
-    type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+    type Adjoint<'a>
+        = ColOp<S::Adjoint<'a>, T::Adjoint<'a>>
+    where
+        Self: 'a;
 
     fn adjoint(&self) -> Self::Adjoint<'_> {
         ColOp(self.0.adjoint(), self.1.adjoint())
     }
 }
 
-impl<A, B, Yʹ, 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>> : Linear<
-        Yʹ, Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>,
-    >,
+    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>>:
+        Linear<Yʹ, Codomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>>,
 {
     type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
-    type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+    type Preadjoint<'a>
+        = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
+    where
+        Self: 'a;
 
     fn preadjoint(&self) -> Self::Preadjoint<'_> {
         ColOp(self.0.preadjoint(), self.1.preadjoint())
     }
 }
 
-
-impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
+impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
 where
-    A : Space,
-    Xʹ : Space,
-    Yʹ : Space,
-    R : Space + ClosedAdd,
-    S : Adjointable<A, Xʹ, AdjointCodomain = R>,
-    T : Adjointable<A, Yʹ, AdjointCodomain = R>,
-    Self : Linear<A>,
+    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,
     // >,
 {
     type AdjointCodomain = R;
-    type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+    type Adjoint<'a>
+        = RowOp<S::Adjoint<'a>, T::Adjoint<'a>>
+    where
+        Self: 'a;
 
     fn adjoint(&self) -> Self::Adjoint<'_> {
         RowOp(self.0.adjoint(), self.1.adjoint())
     }
 }
 
-impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
+impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
 where
-    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>> : Linear<
-        Pair<Xʹ,Yʹ>, Codomain = R,
-    >,
+    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>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
 {
     type PreadjointCodomain = R;
-    type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+    type Preadjoint<'a>
+        = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
+    where
+        Self: 'a;
 
     fn preadjoint(&self) -> Self::Preadjoint<'_> {
         RowOp(self.0.preadjoint(), self.1.preadjoint())
@@ -521,14 +565,14 @@
 
 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
 where
-    A : Space,
-    B : Space,
-    S : Mapping<A>,
-    T : Mapping<B>,
+    A: Space,
+    B: Space,
+    S: Mapping<A>,
+    T: Mapping<B>,
 {
     type Codomain = Pair<S::Codomain, T::Codomain>;
 
-    fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
+    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))
     }
@@ -536,81 +580,84 @@
 
 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
 where
-    A : Space,
-    B : Space,
-    S : Linear<A>,
-    T : Linear<B>,
-{ }
+    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>>,
+    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>>,
 {
-    fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
+    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<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
+    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<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
+    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,
+    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,
-    >,
+    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>,
 {
     type AdjointCodomain = R;
-    type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+    type Adjoint<'a>
+        = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>
+    where
+        Self: 'a;
 
     fn adjoint(&self) -> Self::Adjoint<'_> {
         DiagOp(self.0.adjoint(), self.1.adjoint())
     }
 }
 
-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,
+    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>> : Linear<
-        Pair<Xʹ,Yʹ>, Codomain=R,
-    >,
+    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>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
 {
     type PreadjointCodomain = R;
-    type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+    type Preadjoint<'a>
+        = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
+    where
+        Self: 'a;
 
     fn preadjoint(&self) -> Self::Preadjoint<'_> {
         DiagOp(self.0.preadjoint(), self.1.preadjoint())
@@ -620,28 +667,26 @@
 /// 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>
+            BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> for RowOp<S, T>
         where
-            F : Float,
-            A : Space + Norm<F, ExpA>,
-            B : Space + Norm<F, ExpB>,
-            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,
+            F: Float,
+            A: Space + Norm<F, ExpA>,
+            B: Space + Norm<F, ExpB>,
+            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
+                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.
@@ -650,23 +695,22 @@
                 na.max(nb)
             }
         }
-        
-        impl<F, A, S, T, ExpA, ExpS, ExpT>
-        BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
-        for ColOp<S, T>
+
+        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 + Norm<F, ExpA>,
-            S : BoundedLinear<A, ExpA, ExpS, F>,
-            T : BoundedLinear<A, ExpA, ExpT, F>,
-            ExpA : NormExponent,
-            ExpS : NormExponent,
-            ExpT : NormExponent,
+            F: Float,
+            A: Space + Norm<F, ExpA>,
+            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>
+                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‖}
@@ -675,10 +719,9 @@
                 ns.max(nt)
             }
         }
-    }
+    };
 }
 
 pairnorm!(L1);
 pairnorm!(L2);
 pairnorm!(Linfinity);
-
diff -r 2e4517b55442 -r 9e5b9fc81c52 src/mapping.rs
--- a/src/mapping.rs	Sun Apr 27 20:41:36 2025 -0500
+++ b/src/mapping.rs	Mon Apr 28 08:26:04 2025 -0500
@@ -2,84 +2,95 @@
 Traits for mathematical functions.
 */
 
-use std::marker::PhantomData;
-use std::borrow::Cow;
-use crate::types::{Num, Float, ClosedMul};
+pub use crate::instance::{BasicDecomposition, Decomposition, Instance, Space};
 use crate::loc::Loc;
-pub use crate::instance::{Instance, Decomposition, BasicDecomposition, Space};
 use crate::norms::{Norm, NormExponent};
-use crate::operator_arithmetic::{Weighted, Constant};
+use crate::operator_arithmetic::{Constant, Weighted};
+use crate::types::{ClosedMul, Float, Num};
+use std::borrow::Cow;
+use std::marker::PhantomData;
 
 /// A mapping from `Domain` to `Self::Codomain`.
-pub trait Mapping<Domain : Space> {
-    type Codomain : Space;
+pub trait Mapping<Domain: Space> {
+    type Codomain: Space;
 
     /// Compute the value of `self` at `x`.
-    fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain;
+    fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain;
 
     #[inline]
     /// Form the composition `self ∘ other`
-    fn compose<X : Space, T : Mapping<X, Codomain=Domain>>(self, other : T)
-        -> Composition<Self, T>
+    fn compose<X: Space, T: Mapping<X, Codomain = Domain>>(self, other: T) -> Composition<Self, T>
     where
-        Self : Sized
+        Self: Sized,
     {
-        Composition{ outer : self, inner : other, intermediate_norm_exponent : () }
+        Composition {
+            outer: self,
+            inner: other,
+            intermediate_norm_exponent: (),
+        }
     }
 
-
     #[inline]
     /// Form the composition `self ∘ other`, assigning a norm to the inermediate space
-    fn compose_with_norm<F, X, T, E>(
-        self, other : T, norm : E
-    )  -> Composition<Self, T, E>
+    fn compose_with_norm<F, X, T, E>(self, other: T, norm: E) -> Composition<Self, T, E>
     where
-        Self : Sized,
-        X : Space,
-        T : Mapping<X, Codomain=Domain>,
-        E : NormExponent,
-        Domain : Norm<F, E>,
-        F : Num
+        Self: Sized,
+        X: Space,
+        T: Mapping<X, Codomain = Domain>,
+        E: NormExponent,
+        Domain: Norm<F, E>,
+        F: Num,
     {
-        Composition{ outer : self, inner : other, intermediate_norm_exponent : norm }
+        Composition {
+            outer: self,
+            inner: other,
+            intermediate_norm_exponent: norm,
+        }
     }
 
     /// Multiply `self` by the scalar `a`.
     #[inline]
-    fn weigh<C>(self, a : C) -> Weighted<Self, C>
+    fn weigh<C>(self, a: C) -> Weighted<Self, C>
     where
-        Self : Sized,
-        C : Constant,
-        Self::Codomain : ClosedMul<C::Type>,
+        Self: Sized,
+        C: Constant,
+        Self::Codomain: ClosedMul<C::Type>,
     {
-        Weighted { weight : a, base_fn : self }
+        Weighted {
+            weight: a,
+            base_fn: self,
+        }
     }
 }
 
 /// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc<F, N>`] to `F`.
-pub trait RealMapping<F : Float, const N : usize>
-: Mapping<Loc<F, N>, Codomain = F> {}
+pub trait RealMapping<F: Float, const N: usize>: Mapping<Loc<F, N>, Codomain = F> {}
 
-impl<F : Float, T, const N : usize> RealMapping<F, N> for T
-where T : Mapping<Loc<F, N>, Codomain = F> {}
+impl<F: Float, T, const N: usize> RealMapping<F, N> for T where T: Mapping<Loc<F, N>, Codomain = F> {}
 
 /// 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>> {}
+pub trait RealVectorField<F: Float, const N: usize, const M: usize>:
+    Mapping<Loc<F, N>, Codomain = Loc<F, M>>
+{
+}
 
-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>> {}
+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>>
+{
+}
 
 /// A differentiable mapping from `Domain` to [`Mapping::Codomain`], with differentials
 /// `Differential`.
 ///
 /// 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;
+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;
+    fn differential<I: Instance<Domain>>(&self, x: I) -> Self::DerivativeDomain;
 
     /// Form the differential mapping of `self`.
     fn diff(self) -> Self::Differential<'static>;
@@ -90,50 +101,62 @@
 
 /// 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>> {}
+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>> {}
+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;
+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;
+    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>
+    Domain: Space,
+    T: Clone + Mapping<Domain> + DifferentiableImpl<Domain>,
 {
     type DerivativeDomain = T::Derivative;
-    type Differential<'b> = Differential<'b, Domain, Self> where Self : 'b;
-    
+    type Differential<'b>
+        = Differential<'b, Domain, Self>
+    where
+        Self: 'b;
+
     #[inline]
-    fn differential<I : Instance<Domain>>(&self, x : I) -> Self::DerivativeDomain {
+    fn differential<I: Instance<Domain>>(&self, x: I) -> Self::DerivativeDomain {
         self.differential_impl(x)
     }
 
     fn diff(self) -> Differential<'static, Domain, Self> {
-        Differential{ g : Cow::Owned(self), _space : PhantomData }
+        Differential {
+            g: Cow::Owned(self),
+            _space: PhantomData,
+        }
     }
 
     fn diff_ref(&self) -> Differential<'_, Domain, Self> {
-        Differential{ g : Cow::Borrowed(self), _space : PhantomData }
+        Differential {
+            g: Cow::Borrowed(self),
+            _space: PhantomData,
+        }
     }
 }
 
-
 /// Container for the differential [`Mapping`] of a [`DifferentiableMapping`].
-pub struct Differential<'a, X, G : Clone> {
-    g : Cow<'a, G>,
-    _space : PhantomData<X>
+pub struct Differential<'a, X, G: Clone> {
+    g: Cow<'a, G>,
+    _space: PhantomData<X>,
 }
 
-impl<'a, X, G : Clone> Differential<'a, X, G> {
+impl<'a, X, G: Clone> Differential<'a, X, G> {
     pub fn base_fn(&self) -> &G {
         &self.g
     }
@@ -141,65 +164,68 @@
 
 impl<'a, X, G> Mapping<X> for Differential<'a, X, G>
 where
-    X : Space,
-    G : Clone + DifferentiableMapping<X>
+    X: Space,
+    G: Clone + DifferentiableMapping<X>,
 {
     type Codomain = G::DerivativeDomain;
 
     #[inline]
-    fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain {
+    fn apply<I: Instance<X>>(&self, x: I) -> Self::Codomain {
         (*self.g).differential(x)
     }
 }
 
 /// Container for flattening [`Loc`]`<F, 1>` codomain of a [`Mapping`] to `F`.
 pub struct FlattenedCodomain<X, F, G> {
-    g : G,
-    _phantoms : PhantomData<(X, F)>
+    g: G,
+    _phantoms: PhantomData<(X, F)>,
 }
 
-impl<F : Space, X, G> Mapping<X> for FlattenedCodomain<X, F, G>
+impl<F: Space, X, G> Mapping<X> for FlattenedCodomain<X, F, G>
 where
-    X : Space,
-    G: Mapping<X, Codomain=Loc<F, 1>>
+    X: Space,
+    G: Mapping<X, Codomain = Loc<F, 1>>,
 {
     type Codomain = F;
 
     #[inline]
-    fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain {
+    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 : Space, 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 }
+        FlattenedCodomain {
+            g: self,
+            _phantoms: PhantomData,
+        }
     }
 }
 
-impl<X : Space, 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<'a, X, F, G : Clone, const N : usize> {
-    g : Cow<'a, G>,
-    slice : usize,
-    _phantoms : PhantomData<(X, F)>
+pub struct SlicedCodomain<'a, X, F, G: Clone, const N: usize> {
+    g: Cow<'a, G>,
+    slice: usize,
+    _phantoms: PhantomData<(X, F)>,
 }
 
-impl<'a, X, F, G, const N : usize> Mapping<X> for SlicedCodomain<'a, X, F, G, N>
+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,
+    X: Space,
+    F: Copy + Space,
+    G: Mapping<X, Codomain = Loc<F, N>> + Clone,
 {
     type Codomain = F;
 
     #[inline]
-    fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain {
-        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) }
     }
@@ -207,44 +233,64 @@
 
 /// 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 : Space, F : Copy, const N : usize>
-    : Mapping<X, Codomain=Loc<F, N>> + Clone + 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<'static, X, F, Self, N> {
+    fn slice_codomain(self, slice: usize) -> SlicedCodomain<'static, X, F, Self, N> {
         assert!(slice < N);
-        SlicedCodomain{ g : Cow::Owned(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 : Cow::Borrowed(self), slice, _phantoms : PhantomData }
+        SlicedCodomain {
+            g: Cow::Borrowed(self),
+            slice,
+            _phantoms: PhantomData,
+        }
     }
 }
 
-impl<X : Space, F : Copy, G : Sized + Mapping<X, Codomain=Loc<F, N>> + Clone, const N : usize>
-SliceCodomain<X, F, N>
-for G {}
-
+impl<X: Space, F: Copy, G: Sized + Mapping<X, Codomain = Loc<F, N>> + Clone, const N: usize>
+    SliceCodomain<X, F, N> for G
+{
+}
 
 /// The composition S ∘ T. `E` is for storing a `NormExponent` for the intermediate space.
 pub struct Composition<S, T, E = ()> {
-    pub outer : S,
-    pub inner : T,
-    pub intermediate_norm_exponent : E
+    pub outer: S,
+    pub inner: T,
+    pub intermediate_norm_exponent: E,
 }
 
 impl<S, T, X, E> Mapping<X> for Composition<S, T, E>
 where
-    X : Space,
-    T : Mapping<X>,
-    S : Mapping<T::Codomain>
+    X: Space,
+    T: Mapping<X>,
+    S: Mapping<T::Codomain>,
 {
     type Codomain = S::Codomain;
 
     #[inline]
-    fn apply<I : Instance<X>>(&self, x : I) -> Self::Codomain {
+    fn apply<I: Instance<X>>(&self, x: I) -> Self::Codomain {
         self.outer.apply(self.inner.apply(x))
     }
 }
+
+mod quadratic;
+pub use quadratic::Quadratic;
+
+/// Trait for indicating that `Self` is Lipschitz with respect to the (semi)norm `D`.
+pub trait Lipschitz<M> {
+    /// The type of floats
+    type FloatType: Float;
+
+    /// Returns the Lipschitz factor of `self` with respect to the (semi)norm `D`.
+    fn lipschitz_factor(&self, seminorm: M) -> Option<Self::FloatType>;
+}
diff -r 2e4517b55442 -r 9e5b9fc81c52 src/mapping/quadratic.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/mapping/quadratic.rs	Mon Apr 28 08:26:04 2025 -0500
@@ -0,0 +1,98 @@
+/*!
+Quadratic functions  of the form $\frac{1}{2}\|Ax-b\|_2^2$ for an operator $A$
+to a [`Euclidean`] space.
+*/
+
+#![allow(non_snake_case)]
+
+use super::{DifferentiableImpl, Differential, Lipschitz, Mapping};
+use crate::convex::ConvexMapping;
+use crate::euclidean::Euclidean;
+use crate::instance::{Instance, Space};
+use crate::linops::{BoundedLinear, Linear, Preadjointable};
+use crate::norms::{Norm, NormExponent, L2};
+use crate::types::Float;
+use std::marker::PhantomData;
+
+/// Functions of the form $\frac{1}{2}\|Ax-b\|_2^2$ for an operator $A$
+/// to a [`Euclidean`] space.
+#[derive(Clone, Copy)]
+pub struct Quadratic<'a, F: Float, Domain: Space, A: Mapping<Domain>> {
+    opA: &'a A,
+    b: &'a <A as Mapping<Domain>>::Codomain,
+    _phantoms: PhantomData<F>,
+}
+
+#[allow(non_snake_case)]
+impl<'a, F: Float, Domain: Space, A: Mapping<Domain>> Quadratic<'a, F, Domain, A> {
+    pub fn new(opA: &'a A, b: &'a A::Codomain) -> Self {
+        Quadratic {
+            opA,
+            b,
+            _phantoms: PhantomData,
+        }
+    }
+
+    pub fn operator(&self) -> &'a A {
+        self.opA
+    }
+
+    pub fn data(&self) -> &'a <A as Mapping<Domain>>::Codomain {
+        self.b
+    }
+}
+
+//+ AdjointProductBoundedBy<RNDM<F, N>, P, FloatType = F>,
+
+impl<'a, F: Float, X: Space, A: Mapping<X>> Mapping<X> for Quadratic<'a, F, X, A>
+where
+    A::Codomain: Euclidean<F>,
+{
+    type Codomain = F;
+
+    fn apply<I: Instance<X>>(&self, x: I) -> F {
+        // TODO: possibly (if at all more effcient) use GEMV once generalised
+        // to not require preallocation. However, Rust should be pretty efficient
+        // at not doing preallocations or anything here, as the result of self.opA.apply()
+        // can be consumed, so maybe GEMV is no use.
+        (self.opA.apply(x) - self.b).norm2_squared() / F::TWO
+    }
+}
+
+impl<'a, F: Float, X: Space, A: Linear<X>> ConvexMapping<X, F> for Quadratic<'a, F, X, A> where
+    A::Codomain: Euclidean<F>
+{
+}
+
+impl<'a, F, X, A> DifferentiableImpl<X> for Quadratic<'a, F, X, A>
+where
+    F: Float,
+    X: Space,
+    <A as Mapping<X>>::Codomain: Euclidean<F>,
+    A: Linear<X> + Preadjointable<X>,
+    <<A as Mapping<X>>::Codomain as Euclidean<F>>::Output: Instance<<A as Mapping<X>>::Codomain>,
+{
+    type Derivative = A::PreadjointCodomain;
+
+    fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
+        // TODO: possibly (if at all more effcient) use GEMV once generalised
+        // to not require preallocation. However, Rust should be pretty efficient
+        // at not doing preallocations or anything here, as the result of self.opA.apply()
+        // can be consumed, so maybe GEMV is no use.
+        self.opA.preadjoint().apply(self.opA.apply(x) - self.b)
+    }
+}
+
+impl<'a, 'b, F, X, ExpX, A> Lipschitz<ExpX> for Differential<'b, X, Quadratic<'a, F, X, A>>
+where
+    F: Float,
+    X: Space + Clone + Norm<F, ExpX>,
+    ExpX: NormExponent,
+    A: Clone + BoundedLinear<X, ExpX, L2, F>,
+{
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, seminorm: ExpX) -> Option<F> {
+        Some((*self.g).opA.opnorm_bound(seminorm, L2).powi(2))
+    }
+}