diff -r efa60bc4f743 -r b087e3eab191 src/kernels/base.rs
--- a/src/kernels/base.rs	Thu Aug 29 00:00:00 2024 -0500
+++ b/src/kernels/base.rs	Tue Dec 31 09:25:45 2024 -0500
@@ -14,16 +14,22 @@
     GlobalAnalysis,
     Bounded,
 };
-use alg_tools::mapping::{Apply, Differentiable};
+use alg_tools::mapping::{
+    Mapping,
+    DifferentiableImpl,
+    DifferentiableMapping,
+    Differential,
+};
+use alg_tools::instance::{Instance, Space};
 use alg_tools::maputil::{array_init, map2, map1_indexed};
 use alg_tools::sets::SetOrd;
 
 use crate::fourier::Fourier;
-use crate::types::Lipschitz;
+use crate::types::*;
 
 /// Representation of the product of two kernels.
 ///
-/// The kernels typically implement [`Support`] and [`Mapping`][alg_tools::mapping::Mapping].
+/// The kernels typically implement [`Support`] and [`Mapping`].
 ///
 /// The implementation [`Support`] only uses the [`Support::support_hint`] of the first parameter!
 #[derive(Copy,Clone,Serialize,Debug)]
@@ -34,59 +40,94 @@
     pub B
 );
 
-impl<A, B, F : Float, const N : usize> Apply<Loc<F, N>>
+impl<A, B, F : Float, const N : usize> Mapping<Loc<F, N>>
 for SupportProductFirst<A, B>
-where A : for<'a> Apply<&'a Loc<F, N>, Output=F>,
-      B : for<'a> Apply<&'a Loc<F, N>, Output=F> {
-    type Output = F;
+where
+    A : Mapping<Loc<F, N>, Codomain = F>,
+    B : Mapping<Loc<F, N>, Codomain = F>,
+{
+    type Codomain = F;
+
     #[inline]
-    fn apply(&self, x : Loc<F, N>) -> Self::Output {
-        self.0.apply(&x) * self.1.apply(&x)
+    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+        self.0.apply(x.ref_instance()) * self.1.apply(x)
     }
 }
 
-impl<'a, A, B, F : Float, const N : usize> Apply<&'a Loc<F, N>>
+impl<A, B, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
 for SupportProductFirst<A, B>
-where A : Apply<&'a Loc<F, N>, Output=F>,
-      B : Apply<&'a Loc<F, N>, Output=F> {
-    type Output = F;
+where
+    A : DifferentiableMapping<
+        Loc<F, N>,
+        DerivativeDomain=Loc<F, N>,
+        Codomain = F
+    >,
+    B : DifferentiableMapping<
+        Loc<F, N>,
+        DerivativeDomain=Loc<F, N>,
+        Codomain = F,
+    >
+{
+    type Derivative = Loc<F, N>;
+
     #[inline]
-    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
-        self.0.apply(x) * self.1.apply(x)
+    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
+        let xr = x.ref_instance();
+        self.0.differential(xr) * self.1.apply(xr) + self.1.differential(xr) * self.0.apply(x)
     }
 }
 
-impl<A, B, F : Float, const N : usize> Differentiable<Loc<F, N>>
+impl<A, B, M : Copy, F : Float> Lipschitz<M>
 for SupportProductFirst<A, B>
-where A : for<'a> Apply<&'a Loc<F, N>, Output=F>
-          + for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
-      B : for<'a> Apply<&'a Loc<F, N>, Output=F>
-          + for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
-    type Output = Loc<F, N>;
+where A : Lipschitz<M, FloatType = F> + Bounded<F>,
+      B : Lipschitz<M, FloatType = F> + Bounded<F> {
+    type FloatType = F;
     #[inline]
-    fn differential(&self, x : Loc<F, N>) -> Self::Output {
-        self.0.differential(&x) * self.1.apply(&x) + self.1.differential(&x) * self.0.apply(&x)
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        // f(x)g(x) - f(y)g(y) = f(x)[g(x)-g(y)] - [f(y)-f(x)]g(y)
+        let &SupportProductFirst(ref f, ref g) = self;
+        f.lipschitz_factor(m).map(|l| l * g.bounds().uniform())
+         .zip(g.lipschitz_factor(m).map(|l| l * f.bounds().uniform()))
+         .map(|(a, b)| a + b)
     }
 }
 
-impl<'a, A, B, F : Float, const N : usize> Differentiable<&'a Loc<F, N>>
-for SupportProductFirst<A, B>
-where A : Apply<&'a Loc<F, N>, Output=F>
-          + Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
-      B : Apply<&'a Loc<F, N>, Output=F>
-          + Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
-    type Output = Loc<F, N>;
+impl<'a, A, B, M : Copy, Domain, F : Float> Lipschitz<M>
+for Differential<'a, Domain, SupportProductFirst<A, B>>
+where
+    Domain : Space,
+    A : Clone + DifferentiableMapping<Domain> + Lipschitz<M, FloatType = F> + Bounded<F>,
+    B : Clone + DifferentiableMapping<Domain> + Lipschitz<M, FloatType = F> + Bounded<F>,
+    SupportProductFirst<A, B> :  DifferentiableMapping<Domain>,
+    for<'b> A::Differential<'b> : Lipschitz<M, FloatType = F> + NormBounded<L2, FloatType=F>,
+    for<'b> B::Differential<'b> : Lipschitz<M, FloatType = F> + NormBounded<L2, FloatType=F>
+{
+    type FloatType = F;
     #[inline]
-    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
-        self.0.differential(&x) * self.1.apply(&x) + self.1.differential(&x) * self.0.apply(&x)
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        // ∇[gf] = f∇g + g∇f
+        // ⟹ ∇[gf](x) - ∇[gf](y) = f(x)∇g(x) + g(x)∇f(x) - f(y)∇g(y) + g(y)∇f(y)
+        //                        = f(x)[∇g(x)-∇g(y)] + g(x)∇f(x) - [f(y)-f(x)]∇g(y) + g(y)∇f(y)
+        //                        = f(x)[∇g(x)-∇g(y)] + g(x)[∇f(x)-∇f(y)]
+        //                          - [f(y)-f(x)]∇g(y) + [g(y)-g(x)]∇f(y)
+        let &SupportProductFirst(ref f, ref g) = self.base_fn();
+        let (df, dg) = (f.diff_ref(), g.diff_ref());
+        [
+            df.lipschitz_factor(m).map(|l| l * g.bounds().uniform()),
+            dg.lipschitz_factor(m).map(|l| l * f.bounds().uniform()),
+            f.lipschitz_factor(m).map(|l| l * dg.norm_bound(L2)),
+            g.lipschitz_factor(m).map(|l| l * df.norm_bound(L2))
+        ].into_iter().sum()
     }
 }
 
 
 impl<'a, A, B, F : Float, const N : usize> Support<F, N>
 for SupportProductFirst<A, B>
-where A : Support<F, N>,
-      B : Support<F, N> {
+where
+    A : Support<F, N>,
+    B : Support<F, N>
+{
     #[inline]
     fn support_hint(&self) -> Cube<F, N> {
         self.0.support_hint()
@@ -125,7 +166,7 @@
 
 /// Representation of the sum of two kernels
 ///
-/// The kernels typically implement [`Support`] and [`Mapping`][alg_tools::mapping::Mapping].
+/// The kernels typically implement [`Support`] and [`Mapping`].
 ///
 /// The implementation [`Support`] only uses the [`Support::support_hint`] of the first parameter!
 #[derive(Copy,Clone,Serialize,Debug)]
@@ -136,55 +177,48 @@
     pub B
 );
 
-impl<'a, A, B, F : Float, const N : usize> Apply<&'a Loc<F, N>>
+impl<'a, A, B, F : Float, const N : usize> Mapping<Loc<F, N>>
 for SupportSum<A, B>
-where A : Apply<&'a Loc<F, N>, Output=F>,
-      B : Apply<&'a Loc<F, N>, Output=F> {
-    type Output = F;
+where
+    A : Mapping<Loc<F, N>, Codomain = F>,
+    B : Mapping<Loc<F, N>, Codomain = F>,
+{
+    type Codomain = F;
+
     #[inline]
-    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
-        self.0.apply(x) + self.1.apply(x)
-    }
-}
-
-impl<A, B, F : Float, const N : usize> Apply<Loc<F, N>>
-for SupportSum<A, B>
-where A : for<'a> Apply<&'a Loc<F, N>, Output=F>,
-      B : for<'a> Apply<&'a Loc<F, N>, Output=F> {
-    type Output = F;
-    #[inline]
-    fn apply(&self, x : Loc<F, N>) -> Self::Output {
-        self.0.apply(&x) + self.1.apply(&x)
+    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+        self.0.apply(x.ref_instance()) + self.1.apply(x)
     }
 }
 
-impl<'a, A, B, F : Float, const N : usize> Differentiable<&'a Loc<F, N>>
+impl<'a, A, B, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
 for SupportSum<A, B>
-where A : Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
-      B : Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
-    type Output = Loc<F, N>;
+where
+    A : DifferentiableMapping<
+        Loc<F, N>,
+        DerivativeDomain = Loc<F, N>
+    >,
+    B : DifferentiableMapping<
+        Loc<F, N>,
+        DerivativeDomain = Loc<F, N>,
+    >
+{
+
+    type Derivative = Loc<F, N>;
+
     #[inline]
-    fn differential(&self, x : &'a Loc<F, N>) -> Self::Output {
-        self.0.differential(x) + self.1.differential(x)
+    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
+        self.0.differential(x.ref_instance()) + self.1.differential(x)
     }
 }
 
-impl<A, B, F : Float, const N : usize> Differentiable<Loc<F, N>>
-for SupportSum<A, B>
-where A : for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>>,
-      B : for<'a> Differentiable<&'a Loc<F, N>, Output=Loc<F, N>> {
-    type Output = Loc<F, N>;
-    #[inline]
-    fn differential(&self, x : Loc<F, N>) -> Self::Output {
-        self.0.differential(&x) + self.1.differential(&x)
-    }
-}
 
 impl<'a, A, B, F : Float, const N : usize> Support<F, N>
 for SupportSum<A, B>
 where A : Support<F, N>,
       B : Support<F, N>,
       Cube<F, N> : SetOrd {
+
     #[inline]
     fn support_hint(&self) -> Cube<F, N> {
         self.0.support_hint().common(&self.1.support_hint())
@@ -237,10 +271,29 @@
     }
 }
 
+impl<'b, F : Float, M : Copy, A, B, Domain> Lipschitz<M>
+for Differential<'b, Domain, SupportSum<A, B>>
+where
+    Domain : Space,
+    A : Clone + DifferentiableMapping<Domain, Codomain=F>,
+    B : Clone + DifferentiableMapping<Domain, Codomain=F>,
+    SupportSum<A, B> : DifferentiableMapping<Domain, Codomain=F>,
+    for<'a> A :: Differential<'a> : Lipschitz<M, FloatType = F>,
+    for<'a> B :: Differential<'a> : Lipschitz<M, FloatType = F>
+{
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        let base = self.base_fn();
+        base.0.diff_ref().lipschitz_factor(m)
+            .zip(base.1.diff_ref().lipschitz_factor(m))
+            .map(|(a, b)| a + b)
+    }
+}
 
 /// Representation of the convolution of two kernels.
 ///
-/// The kernels typically implement [`Support`]s and [`Mapping`][alg_tools::mapping::Mapping].
+/// The kernels typically implement [`Support`]s and [`Mapping`].
 //
 /// Trait implementations have to be on a case-by-case basis.
 #[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
@@ -252,18 +305,45 @@
 );
 
 impl<F : Float, M, A, B> Lipschitz<M> for Convolution<A, B>
-where A : Bounded<F> ,
+where A : Norm<F, L1> ,
       B : Lipschitz<M, FloatType = F> {
     type FloatType = F;
 
     fn lipschitz_factor(&self, m : M) -> Option<F> {
-        self.1.lipschitz_factor(m).map(|l| l * self.0.bounds().uniform())
+        // For [f * g](x) = ∫ f(x-y)g(y) dy we have
+        // [f * g](x) - [f * g](z) = ∫ [f(x-y)-f(z-y)]g(y) dy.
+        // Hence |[f * g](x) - [f * g](z)| ≤ ∫ |f(x-y)-f(z-y)|g(y)| dy.
+        //                                 ≤ L|x-z| ∫ |g(y)| dy,
+        // where L is the Lipschitz factor of f.
+        self.1.lipschitz_factor(m).map(|l| l * self.0.norm(L1))
+    }
+}
+
+impl<'b, F : Float, M, A, B, Domain> Lipschitz<M>
+for Differential<'b, Domain, Convolution<A, B>>
+where
+    Domain : Space,
+    A : Clone + Norm<F, L1> ,
+    Convolution<A, B> : DifferentiableMapping<Domain, Codomain=F>,
+    B : Clone + DifferentiableMapping<Domain, Codomain=F>,
+    for<'a> B :: Differential<'a> : Lipschitz<M, FloatType = F>
+{
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        // For [f * g](x) = ∫ f(x-y)g(y) dy we have
+        // ∇[f * g](x) - ∇[f * g](z) = ∫ [∇f(x-y)-∇f(z-y)]g(y) dy.
+        // Hence |∇[f * g](x) - ∇[f * g](z)| ≤ ∫ |∇f(x-y)-∇f(z-y)|g(y)| dy.
+        //                                 ≤ L|x-z| ∫ |g(y)| dy,
+        // where L is the Lipschitz factor of ∇f.
+        let base = self.base_fn();
+        base.1.diff_ref().lipschitz_factor(m).map(|l| l * base.0.norm(L1))
     }
 }
 
 /// Representation of the autoconvolution of a kernel.
 ///
-/// The kernel typically implements [`Support`] and [`Mapping`][alg_tools::mapping::Mapping].
+/// The kernel typically implements [`Support`] and [`Mapping`].
 ///
 /// Trait implementations have to be on a case-by-case basis.
 #[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
@@ -273,11 +353,27 @@
 );
 
 impl<F : Float, M, C> Lipschitz<M> for AutoConvolution<C>
-where C : Lipschitz<M, FloatType = F> + Bounded<F> {
+where C : Lipschitz<M, FloatType = F> + Norm<F, L1> {
     type FloatType = F;
 
     fn lipschitz_factor(&self, m : M) -> Option<F> {
-        self.0.lipschitz_factor(m).map(|l| l * self.0.bounds().uniform())
+        self.0.lipschitz_factor(m).map(|l| l * self.0.norm(L1))
+    }
+}
+
+impl<'b, F : Float, M, C, Domain> Lipschitz<M>
+for Differential<'b, Domain, AutoConvolution<C>>
+where
+    Domain : Space,
+    C : Clone + Norm<F, L1> + DifferentiableMapping<Domain, Codomain=F>,
+    AutoConvolution<C> : DifferentiableMapping<Domain, Codomain=F>,
+    for<'a> C :: Differential<'a> : Lipschitz<M, FloatType = F>
+{
+    type FloatType = F;
+
+    fn lipschitz_factor(&self, m : M) -> Option<F> {
+        let base = self.base_fn();
+        base.0.diff_ref().lipschitz_factor(m).map(|l| l * base.0.norm(L1))
     }
 }
 
@@ -285,42 +381,47 @@
 /// Representation a multi-dimensional product of a one-dimensional kernel.
 ///
 /// For $G: ℝ → ℝ$, this is the function $F(x\_1, …, x\_n) := \prod_{i=1}^n G(x\_i)$.
-/// The kernel $G$ typically implements [`Support`] and [`Mapping`][alg_tools::mapping::Mapping]
+/// The kernel $G$ typically implements [`Support`] and [`Mapping`]
 /// on [`Loc<F, 1>`]. Then the product implements them on [`Loc<F, N>`].
 #[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
+#[allow(dead_code)]
 struct UniformProduct<G, const N : usize>(
     /// The one-dimensional kernel
     G
 );
 
-impl<'a, G, F : Float, const N : usize> Apply<&'a Loc<F, N>>
+impl<'a, G, F : Float, const N : usize> Mapping<Loc<F, N>>
 for UniformProduct<G, N>
-where G : Apply<Loc<F, 1>, Output=F> {
-    type Output = F;
+where
+    G : Mapping<Loc<F, 1>, Codomain = F>
+{
+    type Codomain = F;
+
     #[inline]
-    fn apply(&self, x : &'a Loc<F, N>) -> F {
-        x.iter().map(|&y| self.0.apply(Loc([y]))).product()
+    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> F {
+        x.cow().iter().map(|&y| self.0.apply(Loc([y]))).product()
     }
 }
 
-impl<G, F : Float, const N : usize> Apply<Loc<F, N>>
+
+
+impl<'a, G, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
 for UniformProduct<G, N>
-where G : Apply<Loc<F, 1>, Output=F> {
-    type Output = F;
+where
+    G : DifferentiableMapping<
+        Loc<F, 1>,
+        DerivativeDomain = F,
+        Codomain = F,
+    >
+{
+    type Derivative = Loc<F, N>;
+
     #[inline]
-    fn apply(&self, x : Loc<F, N>) -> F {
-        x.into_iter().map(|y| self.0.apply(Loc([y]))).product()
-    }
-}
-
-impl<'a, G, F : Float, const N : usize> Differentiable<&'a Loc<F, N>>
-for UniformProduct<G, N>
-where G : Apply<Loc<F, 1>, Output=F> + Differentiable<Loc<F, 1>, Output=F> {
-    type Output = Loc<F, N>;
-    #[inline]
-    fn differential(&self, x : &'a Loc<F, N>) -> Loc<F, N> {
-        let vs = x.map(|y| self.0.apply(Loc([y])));
-        product_differential(x, &vs, |y| self.0.differential(Loc([y])))
+    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x0 : I) -> Loc<F, N> {
+        x0.eval(|x| {
+            let vs = x.map(|y| self.0.apply(Loc([y])));
+            product_differential(x, &vs, |y| self.0.differential(Loc([y])))
+        })
     }
 }
 
@@ -342,13 +443,39 @@
     }).into()
 }
 
-impl<G, F : Float, const N : usize> Differentiable<Loc<F, N>>
-for UniformProduct<G, N>
-where G : Apply<Loc<F, 1>, Output=F> + Differentiable<Loc<F, 1>, Output=F> {
-    type Output = Loc<F, N>;
-    #[inline]
-    fn differential(&self, x : Loc<F, N>) -> Loc<F, N> {
-        self.differential(&x)
+/// Helper function to calulate the Lipschitz factor of $∇f$ for $f(x)=∏_{i=1}^N g(x_i)$.
+///
+/// The parameter `bound` is a bound on $|g|_∞$, `lip` is a Lipschitz factor for $g$,
+/// `dbound` is a bound on $|∇g|_∞$, and `dlip` a Lipschitz factor for $∇g$.
+#[inline]
+pub(crate) fn product_differential_lipschitz_factor<F : Float, const N : usize>(
+    bound : F,
+    lip : F,
+    dbound : F,
+    dlip : F
+) -> F {
+    // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have
+    // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j ≠ i} ψ_j(x_j)
+    // by a simple recursive argument. In particular, if ψ_i=g for all i, j, we have
+    // |ψ(x) - ψ(y)| ≤ ∑_i L_g M_g^{n-1}|x-y|, where L_g is the Lipschitz factor of g, and
+    // M_g a bound on it.
+    //
+    // We also have in the general case ∇ψ(x) = ∑_i ∇ψ_i(x_i) ∏_{j ≠ i} ψ_j(x_j), whence
+    // using the previous formula for each i with f_i=∇ψ_i and f_j=ψ_j for j ≠ i, we get
+    //  ∇ψ(x) - ∇ψ(y) = ∑_i[ ∇ψ_i(x_i)∏_{j ≠ i} ψ_j(x_j) - ∇ψ_i(y_i)∏_{j ≠ i} ψ_j(y_j)]
+    //                = ∑_i[ [∇ψ_i(x_i) - ∇ψ_j(x_j)] ∏_{j ≠ i}ψ_j(x_j)
+    //                       + [∑_{k ≠ i} [ψ_k(x_k) - ∇ψ_k(x_k)] ∏_{j ≠ i, k}ψ_j(x_j)]∇ψ_i(x_i)].
+    // With $ψ_i=g for all i, j, it follows that
+    // |∇ψ(x) - ∇ψ(y)| ≤ ∑_i L_{∇g} M_g^{n-1} + ∑_{k ≠ i} L_g M_g^{n-2} M_{∇g}
+    //                 = n [L_{∇g} M_g^{n-1} + (n-1) L_g M_g^{n-2} M_{∇g}].
+    //                 = n M_g^{n-2}[L_{∇g} M_g + (n-1) L_g M_{∇g}].
+    if N >= 2 {
+        F::cast_from(N) * bound.powi((N-2) as i32)
+                        * (dlip * bound  + F::cast_from(N-1) * lip * dbound)
+    } else if N==1 {
+        dlip
+    } else {
+        panic!("Invalid dimension")
     }
 }