diff -r efa60bc4f743 -r b087e3eab191 src/kernels/hat_convolution.rs
--- a/src/kernels/hat_convolution.rs	Thu Aug 29 00:00:00 2024 -0500
+++ b/src/kernels/hat_convolution.rs	Tue Dec 31 09:25:45 2024 -0500
@@ -14,7 +14,12 @@
     GlobalAnalysis,
     Bounded,
 };
-use alg_tools::mapping::{Apply, Differentiable};
+use alg_tools::mapping::{
+    Mapping,
+    Instance,
+    DifferentiableImpl,
+    Differential,
+};
 use alg_tools::maputil::array_init;
 
 use crate::types::Lipschitz;
@@ -39,6 +44,31 @@
 ///         -\frac{2}{3} (y-1)^3 & \frac{1}{2}\leq y<1. \\\\
 ///     \end{cases}
 /// $$
+// Hence
+// $$
+//     (h\*h)'(y) =
+//     \begin{cases}
+//         2 (y+1)^2 & -1<y\leq -\frac{1}{2}, \\\\
+//         -6 y^2-4 y & -\frac{1}{2}<y\leq 0, \\\\
+//         6 y^2-4 y & 0<y<\frac{1}{2}, \\\\
+//         -2 (y-1)^2 & \frac{1}{2}\leq y<1. \\\\
+//     \end{cases}
+// $$
+// as well as
+// $$
+//     (h\*h)''(y) =
+//     \begin{cases}
+//         4 (y+1) & -1<y\leq -\frac{1}{2}, \\\\
+//         -12 y-4 & -\frac{1}{2}<y\leq 0, \\\\
+//         12 y-4 & 0<y<\frac{1}{2}, \\\\
+//         -4 (y-1) & \frac{1}{2}\leq y<1. \\\\
+//     \end{cases}
+// $$
+// This is maximised at y=±1/2 with value 2, and minimised at y=0 with value -4.
+// Now observe that
+// $$
+//     [∇f(x\_1, …, x\_n)]_j = \frac{4}{σ} (h\*h)'(x\_j/σ) \prod\_{j ≠ i} \frac{4}{σ} (h\*h)(x\_i/σ)
+// $$
 #[derive(Copy,Clone,Debug,Serialize,Eq)]
 pub struct HatConv<S : Constant, const N : usize> {
     /// The parameter $σ$ of the kernel.
@@ -61,27 +91,19 @@
     }
 }
 
-impl<'a, S, const N : usize> Apply<&'a Loc<S::Type, N>> for HatConv<S, N>
+impl<'a, S, const N : usize> Mapping<Loc<S::Type, N>> for HatConv<S, N>
 where S : Constant {
-    type Output = S::Type;
+    type Codomain = S::Type;
+
     #[inline]
-    fn apply(&self, y : &'a Loc<S::Type, N>) -> Self::Output {
+    fn apply<I : Instance<Loc<S::Type, N>>>(&self, y : I) -> Self::Codomain {
         let σ = self.radius();
-        y.product_map(|x| {
+        y.cow().product_map(|x| {
             self.value_1d_σ1(x  / σ) / σ
         })
     }
 }
 
-impl<'a, S, const N : usize> Apply<Loc<S::Type, N>> for HatConv<S, N>
-where S : Constant {
-    type Output = S::Type;
-    #[inline]
-    fn apply(&self, y : Loc<S::Type, N>) -> Self::Output {
-        self.apply(&y)
-    }
-}
-
 #[replace_float_literals(S::Type::cast_from(literal))]
 impl<S, const N : usize> Lipschitz<L1> for HatConv<S, N>
 where S : Constant {
@@ -95,7 +117,7 @@
         // = ∑_{j=1}^N [ψ_j(x_j)-ψ_j(y_j)]∏_{i > j} ψ_i(x_i) ∏_{i < j} ψ_i(y_i)
         // Thus
         // |∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)|
-        // ≤ ∑_{j=1}^N |ψ_j(x_j)-ψ_j(y_j)| ∏_{j ≠ i} \max_i |ψ_i|
+        // ≤ ∑_{j=1}^N |ψ_j(x_j)-ψ_j(y_j)| ∏_{j ≠ i} \max_j |ψ_j|
         let σ = self.radius();
         let l1d = self.lipschitz_1d_σ1() / (σ*σ);
         let m1d = self.value_1d_σ1(0.0) / σ;
@@ -113,31 +135,45 @@
 }
 
 
-impl<'a, S, const N : usize> Differentiable<&'a Loc<S::Type, N>> for HatConv<S, N>
+impl<'a, S, const N : usize> DifferentiableImpl<Loc<S::Type, N>> for HatConv<S, N>
 where S : Constant {
-    type Output = Loc<S::Type, N>;
+    type Derivative = Loc<S::Type, N>;
+
     #[inline]
-    fn differential(&self, y : &'a Loc<S::Type, N>) -> Self::Output {
+    fn differential_impl<I : Instance<Loc<S::Type, N>>>(&self, y0 : I) -> Self::Derivative {
+        let y = y0.cow();
         let σ = self.radius();
         let σ2 = σ * σ;
         let vs = y.map(|x| {
             self.value_1d_σ1(x  / σ) / σ
         });
-        product_differential(y, &vs, |x| {
+        product_differential(&*y, &vs, |x| {
             self.diff_1d_σ1(x  / σ) / σ2
         })
     }
 }
 
-impl<'a, S, const N : usize> Differentiable<Loc<S::Type, N>> for HatConv<S, N>
-where S : Constant {
-    type Output = Loc<S::Type, N>;
+
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<'a, F : Float, S, const N : usize> Lipschitz<L2>
+for Differential<'a, Loc<F, N>, HatConv<S, N>>
+where S : Constant<Type=F> {
+    type FloatType = F;
+
     #[inline]
-    fn differential(&self, y : Loc<S::Type, N>) -> Self::Output {
-        self.differential(&y)
+    fn lipschitz_factor(&self, _l2 : L2) -> Option<F> {
+        let h = self.base_fn();
+        let σ = h.radius();
+        Some(product_differential_lipschitz_factor::<F, N>(
+            h.value_1d_σ1(0.0) / σ,
+            h.lipschitz_1d_σ1() / (σ*σ),
+            h.maxabsdiff_1d_σ1() / (σ*σ),
+            h.lipschitz_diff_1d_σ1() / (σ*σ),
+        ))
     }
 }
 
+
 #[replace_float_literals(S::Type::cast_from(literal))]
 impl<'a, F : Float, S, const N : usize> HatConv<S, N>
 where S : Constant<Type=F> {
@@ -173,6 +209,34 @@
         // Maximal absolute differential achieved at ±0.5 by diff_1d_σ1 analysis
         2.0
     }
+
+    /// Computes the maximum absolute differential of the kernel for $n=1$ with $σ=1$.
+    #[inline]
+    fn maxabsdiff_1d_σ1(&self) -> F {
+        // Maximal absolute differential achieved at ±0.5 by diff_1d_σ1 analysis
+        2.0
+    }
+
+    /// Computes the second differential of the kernel for $n=1$ with $σ=1$.
+    #[inline]
+    #[allow(dead_code)]
+    fn diff2_1d_σ1(&self, x : F) -> F {
+        let y = x.abs();
+        if y >= 1.0 {
+            0.0
+        } else if y > 0.5 {
+            - 16.0 * (y - 1.0)
+        } else /* 0 ≤ y ≤ 0.5 */ {
+            48.0 * y - 16.0
+        }
+    }
+
+    /// Computes the differential of the kernel for $n=1$ with $σ=1$.
+    #[inline]
+    fn lipschitz_diff_1d_σ1(&self) -> F {
+        // Maximal absolute second differential achieved at 0 by diff2_1d_σ1 analysis
+        16.0
+    }
 }
 
 impl<'a, S, const N : usize> Support<S::Type, N> for HatConv<S, N>
@@ -235,21 +299,21 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, R, C, const N : usize> Apply<&'a Loc<F, N>>
+impl<'a, F : Float, R, C, const N : usize> Mapping<Loc<F, N>>
 for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
 where R : Constant<Type=F>,
       C : Constant<Type=F> {
 
-    type Output = F;
+    type Codomain = F;
 
     #[inline]
-    fn apply(&self, y : &'a Loc<F, N>) -> F {
+    fn apply<I : Instance<Loc<F, N>>>(&self, y : I) -> F {
         let Convolution(ref ind, ref hatconv) = self;
         let β = ind.r.value();
         let σ = hatconv.radius();
 
         // This is just a product of one-dimensional versions
-        y.product_map(|x| {
+        y.cow().product_map(|x| {
             // With $u_σ(x) = u_1(x/σ)/σ$ the normalised hat convolution
             // we have
             // $$
@@ -264,29 +328,17 @@
     }
 }
 
-impl<'a, F : Float, R, C, const N : usize> Apply<Loc<F, N>>
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, const N : usize> DifferentiableImpl<Loc<F, N>>
 for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
 where R : Constant<Type=F>,
       C : Constant<Type=F> {
 
-    type Output = F;
+    type Derivative = Loc<F, N>;
 
     #[inline]
-    fn apply(&self, y : Loc<F, N>) -> F {
-        self.apply(&y)
-    }
-}
-
-#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, R, C, const N : usize> Differentiable<&'a Loc<F, N>>
-for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
-where R : Constant<Type=F>,
-      C : Constant<Type=F> {
-
-    type Output = Loc<F, N>;
-
-    #[inline]
-    fn differential(&self, y : &'a Loc<F, N>) -> Loc<F, N> {
+    fn differential_impl<I : Instance<Loc<F, N>>>(&self, y0 : I) -> Loc<F, N> {
+        let y = y0.cow();
         let Convolution(ref ind, ref hatconv) = self;
         let β = ind.r.value();
         let σ = hatconv.radius();
@@ -295,24 +347,12 @@
         let vs = y.map(|x| {
             self.value_1d_σ1(x / σ, β / σ)
         });
-        product_differential(y, &vs, |x| {
+        product_differential(&*y, &vs, |x| {
             self.diff_1d_σ1(x  / σ, β / σ) / σ2
         })
     }
 }
 
-impl<'a, F : Float, R, C, const N : usize> Differentiable<Loc<F, N>>
-for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
-where R : Constant<Type=F>,
-      C : Constant<Type=F> {
-
-    type Output = Loc<F, N>;
-
-    #[inline]
-    fn differential(&self, y : Loc<F, N>) -> Loc<F, N> {
-        self.differential(&y)
-    }
-}
 
 /// Integrate $f$, whose support is $[c, d]$, on $[a, b]$.
 /// If $b > d$, add $g()$ to the result.
@@ -415,6 +455,22 @@
     }
 }
 
+/*
+impl<'a, F : Float, R, C, const N : usize> Lipschitz<L2>
+for Differential<Loc<F, N>, Convolution<CubeIndicator<R, N>, HatConv<C, N>>>
+where R : Constant<Type=F>,
+      C : Constant<Type=F> {
+
+    type FloatType = F;
+
+    #[inline]
+    fn lipschitz_factor(&self, _l2 : L2) -> Option<F> {
+        dbg!("unimplemented");
+        None
+    }
+}
+*/
+
 impl<F : Float, R, C, const N : usize>
 Convolution<CubeIndicator<R, N>, HatConv<C, N>>
 where R : Constant<Type=F>,
@@ -556,7 +612,7 @@
 #[cfg(test)]
 mod tests {
     use alg_tools::lingrid::linspace;
-    use alg_tools::mapping::Apply;
+    use alg_tools::mapping::Mapping;
     use alg_tools::norms::Linfinity;
     use alg_tools::loc::Loc;
     use crate::kernels::{BallIndicator, CubeIndicator, Convolution};