diff -r 9738b51d90d7 -r 4f468d35fa29 src/kernels/hat_convolution.rs
--- a/src/kernels/hat_convolution.rs	Sun Apr 27 15:03:51 2025 -0500
+++ b/src/kernels/hat_convolution.rs	Thu Feb 26 11:38:43 2026 -0500
@@ -1,30 +1,22 @@
 //! Implementation of the convolution of two hat functions,
 //! and its convolution with a [`CubeIndicator`].
+
+use alg_tools::bisection_tree::{Constant, Support};
+use alg_tools::bounds::{Bounded, Bounds, GlobalAnalysis, LocalAnalysis};
+use alg_tools::error::DynResult;
+use alg_tools::loc::Loc;
+use alg_tools::mapping::{DifferentiableImpl, Instance, LipschitzDifferentiableImpl, Mapping};
+use alg_tools::maputil::array_init;
+use alg_tools::norms::*;
+use alg_tools::sets::Cube;
+use alg_tools::types::*;
+use anyhow::anyhow;
 use numeric_literals::replace_float_literals;
 use serde::Serialize;
-use alg_tools::types::*;
-use alg_tools::norms::*;
-use alg_tools::loc::Loc;
-use alg_tools::sets::Cube;
-use alg_tools::bisection_tree::{
-    Support,
-    Constant,
-    Bounds,
-    LocalAnalysis,
-    GlobalAnalysis,
-    Bounded,
-};
-use alg_tools::mapping::{
-    Mapping,
-    Instance,
-    DifferentiableImpl,
-    Differential,
-};
-use alg_tools::maputil::array_init;
 
-use crate::types::Lipschitz;
+use super::ball_indicator::CubeIndicator;
 use super::base::*;
-use super::ball_indicator::CubeIndicator;
+use crate::types::Lipschitz;
 
 /// Hat convolution kernel.
 ///
@@ -69,21 +61,26 @@
 // $$
 //     [∇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> {
+#[derive(Copy, Clone, Debug, Serialize, Eq)]
+pub struct HatConv<S: Constant, const N: usize> {
     /// The parameter $σ$ of the kernel.
-    pub radius : S,
+    pub radius: S,
 }
 
-impl<S1, S2, const N : usize> PartialEq<HatConv<S2, N>> for HatConv<S1, N>
-where S1 : Constant,
-      S2 : Constant<Type=S1::Type> {
-    fn eq(&self, other : &HatConv<S2, N>) -> bool {
+impl<S1, S2, const N: usize> PartialEq<HatConv<S2, N>> for HatConv<S1, N>
+where
+    S1: Constant,
+    S2: Constant<Type = S1::Type>,
+{
+    fn eq(&self, other: &HatConv<S2, N>) -> bool {
         self.radius.value() == other.radius.value()
     }
 }
 
-impl<'a, S, const N : usize> HatConv<S, N> where S : Constant {
+impl<'a, S, const N: usize> HatConv<S, N>
+where
+    S: Constant,
+{
     /// Returns the $σ$ parameter of the kernel.
     #[inline]
     pub fn radius(&self) -> S::Type {
@@ -91,25 +88,27 @@
     }
 }
 
-impl<'a, S, const N : usize> Mapping<Loc<S::Type, N>> for HatConv<S, N>
-where S : Constant {
+impl<'a, S, const N: usize> Mapping<Loc<N, S::Type>> for HatConv<S, N>
+where
+    S: Constant,
+{
     type Codomain = S::Type;
 
     #[inline]
-    fn apply<I : Instance<Loc<S::Type, N>>>(&self, y : I) -> Self::Codomain {
+    fn apply<I: Instance<Loc<N, S::Type>>>(&self, y: I) -> Self::Codomain {
         let σ = self.radius();
-        y.cow().product_map(|x| {
-            self.value_1d_σ1(x  / σ) / σ
-        })
+        y.decompose().product_map(|x| self.value_1d_σ1(x / σ) / σ)
     }
 }
 
 #[replace_float_literals(S::Type::cast_from(literal))]
-impl<S, const N : usize> Lipschitz<L1> for HatConv<S, N>
-where S : Constant {
+impl<S, const N: usize> Lipschitz<L1> for HatConv<S, N>
+where
+    S: Constant,
+{
     type FloatType = S::Type;
     #[inline]
-    fn lipschitz_factor(&self, L1 : L1) -> Option<Self::FloatType> {
+    fn lipschitz_factor(&self, L1: L1) -> DynResult<Self::FloatType> {
         // For any ψ_i, we have
         // ∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)
         // = [ψ_1(x_1)-ψ_1(y_1)] ∏_{i=2}^N ψ_i(x_i)
@@ -119,86 +118,91 @@
         // |∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_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 l1d = self.lipschitz_1d_σ1() / (σ * σ);
         let m1d = self.value_1d_σ1(0.0) / σ;
-        Some(l1d * m1d.powi(N as i32 - 1))
+        Ok(l1d * m1d.powi(N as i32 - 1))
     }
 }
 
-impl<S, const N : usize> Lipschitz<L2> for HatConv<S, N>
-where S : Constant {
+impl<S, const N: usize> Lipschitz<L2> for HatConv<S, N>
+where
+    S: Constant,
+{
     type FloatType = S::Type;
     #[inline]
-    fn lipschitz_factor(&self, L2 : L2) -> Option<Self::FloatType> {
-        self.lipschitz_factor(L1).map(|l1| l1 * <S::Type>::cast_from(N).sqrt())
+    fn lipschitz_factor(&self, L2: L2) -> DynResult<Self::FloatType> {
+        self.lipschitz_factor(L1)
+            .map(|l1| l1 * <S::Type>::cast_from(N).sqrt())
     }
 }
 
-
-impl<'a, S, const N : usize> DifferentiableImpl<Loc<S::Type, N>> for HatConv<S, N>
-where S : Constant {
-    type Derivative = Loc<S::Type, N>;
+impl<'a, S, const N: usize> DifferentiableImpl<Loc<N, S::Type>> for HatConv<S, N>
+where
+    S: Constant,
+{
+    type Derivative = Loc<N, S::Type>;
 
     #[inline]
-    fn differential_impl<I : Instance<Loc<S::Type, N>>>(&self, y0 : I) -> Self::Derivative {
-        let y = y0.cow();
+    fn differential_impl<I: Instance<Loc<N, S::Type>>>(&self, y0: I) -> Self::Derivative {
+        let y = y0.decompose();
         let σ = self.radius();
         let σ2 = σ * σ;
-        let vs = y.map(|x| {
-            self.value_1d_σ1(x  / σ) / σ
-        });
-        product_differential(&*y, &vs, |x| {
-            self.diff_1d_σ1(x  / σ) / σ2
-        })
+        let vs = y.map(|x| self.value_1d_σ1(x / σ) / σ);
+        product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ) / σ2)
     }
 }
 
-
 #[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> {
+impl<'a, F: Float, S, const N: usize> LipschitzDifferentiableImpl<Loc<N, S::Type>, L2>
+    for HatConv<S, N>
+where
+    S: Constant<Type = F>,
+{
     type FloatType = F;
 
     #[inline]
-    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() / (σ*σ),
-        ))
+    fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<F> {
+        let σ = self.radius();
+        product_differential_lipschitz_factor::<F, N>(
+            self.value_1d_σ1(0.0) / σ,
+            self.lipschitz_1d_σ1() / (σ * σ),
+            self.maxabsdiff_1d_σ1() / (σ * σ),
+            self.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> {
+impl<'a, F: Float, S, const N: usize> HatConv<S, N>
+where
+    S: Constant<Type = F>,
+{
     /// Computes the value of the kernel for $n=1$ with $σ=1$.
     #[inline]
-    fn value_1d_σ1(&self, x : F) -> F {
+    fn value_1d_σ1(&self, x: F) -> F {
         let y = x.abs();
         if y >= 1.0 {
             0.0
         } else if y > 0.5 {
-            - (8.0/3.0) * (y - 1.0).powi(3)
-        } else /* 0 ≤ y ≤ 0.5 */ {
-            (4.0/3.0) + 8.0 * y * y * (y - 1.0)
+            -(8.0 / 3.0) * (y - 1.0).powi(3)
+        } else
+        /* 0 ≤ y ≤ 0.5 */
+        {
+            (4.0 / 3.0) + 8.0 * y * y * (y - 1.0)
         }
     }
 
     /// Computes the differential of the kernel for $n=1$ with $σ=1$.
     #[inline]
-    fn diff_1d_σ1(&self, x : F) -> F {
+    fn diff_1d_σ1(&self, x: F) -> F {
         let y = x.abs();
         if y >= 1.0 {
             0.0
         } else if y > 0.5 {
-            - 8.0 * (y - 1.0).powi(2)
-        } else /* 0 ≤ y ≤ 0.5 */ {
+            -8.0 * (y - 1.0).powi(2)
+        } else
+        /* 0 ≤ y ≤ 0.5 */
+        {
             (24.0 * y - 16.0) * y
         }
     }
@@ -220,13 +224,15 @@
     /// Computes the second differential of the kernel for $n=1$ with $σ=1$.
     #[inline]
     #[allow(dead_code)]
-    fn diff2_1d_σ1(&self, x : F) -> F {
+    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 */ {
+            -16.0 * (y - 1.0)
+        } else
+        /* 0 ≤ y ≤ 0.5 */
+        {
             48.0 * y - 16.0
         }
     }
@@ -239,40 +245,46 @@
     }
 }
 
-impl<'a, S, const N : usize> Support<S::Type, N> for HatConv<S, N>
-where S : Constant {
+impl<'a, S, const N: usize> Support<N, S::Type> for HatConv<S, N>
+where
+    S: Constant,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<S::Type,N> {
+    fn support_hint(&self) -> Cube<N, S::Type> {
         let σ = self.radius();
         array_init(|| [-σ, σ]).into()
     }
 
     #[inline]
-    fn in_support(&self, y : &Loc<S::Type,N>) -> bool {
+    fn in_support(&self, y: &Loc<N, S::Type>) -> bool {
         let σ = self.radius();
         y.iter().all(|x| x.abs() <= σ)
     }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<S::Type, N>) -> [Option<S::Type>; N] {
+    fn bisection_hint(&self, cube: &Cube<N, S::Type>) -> [Option<S::Type>; N] {
         let σ = self.radius();
         cube.map(|c, d| symmetric_peak_hint(σ, c, d))
     }
 }
 
 #[replace_float_literals(S::Type::cast_from(literal))]
-impl<S, const N : usize> GlobalAnalysis<S::Type, Bounds<S::Type>>  for HatConv<S, N>
-where S : Constant {
+impl<S, const N: usize> GlobalAnalysis<S::Type, Bounds<S::Type>> for HatConv<S, N>
+where
+    S: Constant,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<S::Type> {
         Bounds(0.0, self.apply(Loc::ORIGIN))
     }
 }
 
-impl<S, const N : usize> LocalAnalysis<S::Type, Bounds<S::Type>, N>  for HatConv<S, N>
-where S : Constant {
+impl<S, const N: usize> LocalAnalysis<S::Type, Bounds<S::Type>, N> for HatConv<S, N>
+where
+    S: Constant,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<S::Type, N>) -> Bounds<S::Type> {
+    fn local_analysis(&self, cube: &Cube<N, S::Type>) -> Bounds<S::Type> {
         // The function is maximised/minimised where the 2-norm is minimised/maximised.
         let lower = self.apply(cube.maxnorm_point());
         let upper = self.apply(cube.minnorm_point());
@@ -281,39 +293,38 @@
 }
 
 #[replace_float_literals(C::Type::cast_from(literal))]
-impl<'a, C : Constant, const N : usize> Norm<C::Type, L1>
-for HatConv<C, N> {
+impl<'a, C: Constant, const N: usize> Norm<L1, C::Type> for HatConv<C, N> {
     #[inline]
-    fn norm(&self, _ : L1) -> C::Type {
+    fn norm(&self, _: L1) -> C::Type {
         1.0
     }
 }
 
 #[replace_float_literals(C::Type::cast_from(literal))]
-impl<'a, C : Constant, const N : usize> Norm<C::Type, Linfinity>
-for HatConv<C, N> {
+impl<'a, C: Constant, const N: usize> Norm<Linfinity, C::Type> for HatConv<C, N> {
     #[inline]
-    fn norm(&self, _ : Linfinity) -> C::Type {
+    fn norm(&self, _: Linfinity) -> C::Type {
         self.bounds().upper()
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-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> {
-
+impl<'a, F: Float, R, C, const N: usize> Mapping<Loc<N, F>>
+    for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
+where
+    R: Constant<Type = F>,
+    C: Constant<Type = F>,
+{
     type Codomain = F;
 
     #[inline]
-    fn apply<I : Instance<Loc<F, N>>>(&self, y : I) -> F {
+    fn apply<I: Instance<Loc<N, F>>>(&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.cow().product_map(|x| {
+        y.decompose().product_map(|x| {
             // With $u_σ(x) = u_1(x/σ)/σ$ the normalised hat convolution
             // we have
             // $$
@@ -329,37 +340,32 @@
 }
 
 #[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 Derivative = Loc<F, N>;
+impl<'a, F: Float, R, C, const N: usize> DifferentiableImpl<Loc<N, F>>
+    for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
+where
+    R: Constant<Type = F>,
+    C: Constant<Type = F>,
+{
+    type Derivative = Loc<N, F>;
 
     #[inline]
-    fn differential_impl<I : Instance<Loc<F, N>>>(&self, y0 : I) -> Loc<F, N> {
-        let y = y0.cow();
+    fn differential_impl<I: Instance<Loc<N, F>>>(&self, y0: I) -> Loc<N, F> {
+        let y = y0.decompose();
         let Convolution(ref ind, ref hatconv) = self;
         let β = ind.r.value();
         let σ = hatconv.radius();
         let σ2 = σ * σ;
 
-        let vs = y.map(|x| {
-            self.value_1d_σ1(x / σ, β / σ)
-        });
-        product_differential(&*y, &vs, |x| {
-            self.diff_1d_σ1(x  / σ, β / σ) / σ2
-        })
+        let vs = y.map(|x| self.value_1d_σ1(x / σ, β / σ));
+        product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ) / σ2)
     }
 }
 
-
 /// Integrate $f$, whose support is $[c, d]$, on $[a, b]$.
 /// If $b > d$, add $g()$ to the result.
 #[inline]
 #[replace_float_literals(F::cast_from(literal))]
-fn i<F: Float>(a : F, b : F, c : F, d : F, f : impl Fn(F) -> F,
-                g : impl Fn() -> F) -> F {
+fn i<F: Float>(a: F, b: F, c: F, d: F, f: impl Fn(F) -> F, g: impl Fn() -> F) -> F {
     if b < c {
         0.0
     } else if b <= d {
@@ -368,7 +374,9 @@
         } else {
             f(b) - f(a)
         }
-    } else /* b > d */ {
+    } else
+    /* b > d */
+    {
         g() + if a <= c {
             f(d) - f(c)
         } else if a < d {
@@ -380,84 +388,130 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float, C, R, const N : usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>>
-where R : Constant<Type=F>,
-      C : Constant<Type=F> {
-      
+impl<F: Float, C, R, const N: usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>>
+where
+    R: Constant<Type = F>,
+    C: Constant<Type = F>,
+{
     /// Calculates the value of the 1D hat convolution further convolved by a interval indicator.
     /// As both functions are piecewise polynomials, this is implemented by explicit integral over
     /// all subintervals of polynomiality of the cube indicator, using easily formed
     /// antiderivatives.
     #[inline]
-    pub fn value_1d_σ1(&self, x : F, β : F) -> F {
+    pub fn value_1d_σ1(&self, x: F, β: F) -> F {
         // The integration interval
         let a = x - β;
         let b = x + β;
 
         #[inline]
-        fn pow4<F : Float>(x : F) -> F {
+        fn pow4<F: Float>(x: F) -> F {
             let y = x * x;
             y * y
         }
-        
+
         // Observe the factor 1/6 at the front from the antiderivatives below.
         // The factor 4 is from normalisation of the original function.
-        (4.0/6.0) * i(a, b, -1.0, -0.5,
+        (4.0 / 6.0)
+            * i(
+                a,
+                b,
+                -1.0,
+                -0.5,
                 // (2/3) (y+1)^3  on  -1 < y ≤ -1/2
                 // The antiderivative is  (2/12)(y+1)^4 = (1/6)(y+1)^4
-                |y| pow4(y+1.0),
-                || i(a, b, -0.5, 0.0,
-                    // -2 y^3 - 2 y^2 + 1/3  on  -1/2 < y ≤ 0
-                    // The antiderivative is -1/2 y^4 - 2/3 y^3 + 1/3 y
-                    |y| y*(-y*y*(y*3.0 + 4.0) + 2.0),
-                    || i(a, b, 0.0, 0.5,
-                            // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
-                            // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y
-                            |y| y*(y*y*(y*3.0 - 4.0) + 2.0),
-                            || i(a, b, 0.5, 1.0,
-                                // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
-                                // The antiderivative is  -(2/12)(y-1)^4 = -(1/6)(y-1)^4
-                                |y| -pow4(y-1.0),
-                                || 0.0
+                |y| pow4(y + 1.0),
+                || {
+                    i(
+                        a,
+                        b,
+                        -0.5,
+                        0.0,
+                        // -2 y^3 - 2 y^2 + 1/3  on  -1/2 < y ≤ 0
+                        // The antiderivative is -1/2 y^4 - 2/3 y^3 + 1/3 y
+                        |y| y * (-y * y * (y * 3.0 + 4.0) + 2.0),
+                        || {
+                            i(
+                                a,
+                                b,
+                                0.0,
+                                0.5,
+                                // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
+                                // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y
+                                |y| y * (y * y * (y * 3.0 - 4.0) + 2.0),
+                                || {
+                                    i(
+                                        a,
+                                        b,
+                                        0.5,
+                                        1.0,
+                                        // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
+                                        // The antiderivative is  -(2/12)(y-1)^4 = -(1/6)(y-1)^4
+                                        |y| -pow4(y - 1.0),
+                                        || 0.0,
+                                    )
+                                },
                             )
+                        },
                     )
-                )
-        )
+                },
+            )
     }
 
     /// Calculates the derivative of the 1D hat convolution further convolved by a interval
     /// indicator. The implementation is similar to [`Self::value_1d_σ1`], using the fact that
     /// $(θ * ψ)' = θ * ψ'$.
     #[inline]
-    pub fn diff_1d_σ1(&self, x : F, β : F) -> F {
+    pub fn diff_1d_σ1(&self, x: F, β: F) -> F {
         // The integration interval
         let a = x - β;
         let b = x + β;
 
         // The factor 4 is from normalisation of the original function.
-        4.0 * i(a, b, -1.0, -0.5,
-                // (2/3) (y+1)^3  on  -1 < y ≤ -1/2
-                |y| (2.0/3.0) * (y + 1.0).powi(3),
-                || i(a, b, -0.5, 0.0,
+        4.0 * i(
+            a,
+            b,
+            -1.0,
+            -0.5,
+            // (2/3) (y+1)^3  on  -1 < y ≤ -1/2
+            |y| (2.0 / 3.0) * (y + 1.0).powi(3),
+            || {
+                i(
+                    a,
+                    b,
+                    -0.5,
+                    0.0,
                     // -2 y^3 - 2 y^2 + 1/3  on  -1/2 < y ≤ 0
-                    |y| -2.0*(y + 1.0) * y * y + (1.0/3.0),
-                    || i(a, b, 0.0, 0.5,
+                    |y| -2.0 * (y + 1.0) * y * y + (1.0 / 3.0),
+                    || {
+                        i(
+                            a,
+                            b,
+                            0.0,
+                            0.5,
                             // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
-                            |y| 2.0*(y - 1.0) * y * y + (1.0/3.0),
-                            || i(a, b, 0.5, 1.0,
-                                // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
-                                |y| -(2.0/3.0) * (y - 1.0).powi(3),
-                                || 0.0
-                            )
-                    )
+                            |y| 2.0 * (y - 1.0) * y * y + (1.0 / 3.0),
+                            || {
+                                i(
+                                    a,
+                                    b,
+                                    0.5,
+                                    1.0,
+                                    // -(2/3) (y-1)^3  on  1/2 < y ≤ 1
+                                    |y| -(2.0 / 3.0) * (y - 1.0).powi(3),
+                                    || 0.0,
+                                )
+                            },
+                        )
+                    },
                 )
+            },
         )
     }
 }
 
 /*
 impl<'a, F : Float, R, C, const N : usize> Lipschitz<L2>
-for Differential<Loc<F, N>, Convolution<CubeIndicator<R, N>, HatConv<C, N>>>
+for Differential<Loc<N, F>, Convolution<CubeIndicator<R, N>, HatConv<C, N>>>
 where R : Constant<Type=F>,
       C : Constant<Type=F> {
 
@@ -471,11 +525,11 @@
 }
 */
 
-impl<F : Float, R, C, const N : usize>
-Convolution<CubeIndicator<R, N>, HatConv<C, N>>
-where R : Constant<Type=F>,
-      C : Constant<Type=F> {
-
+impl<F: Float, R, C, const N: usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>>
+where
+    R: Constant<Type = F>,
+    C: Constant<Type = F>,
+{
     #[inline]
     fn get_r(&self) -> F {
         let Convolution(ref ind, ref hatconv) = self;
@@ -483,25 +537,26 @@
     }
 }
 
-impl<F : Float, R, C, const N : usize> Support<F, N>
-for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
-where R : Constant<Type=F>,
-      C : Constant<Type=F> {
-    
+impl<F: Float, R, C, const N: usize> Support<N, F>
+    for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
+where
+    R: Constant<Type = F>,
+    C: Constant<Type = F>,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F, N> {
+    fn support_hint(&self) -> Cube<N, F> {
         let r = self.get_r();
         array_init(|| [-r, r]).into()
     }
 
     #[inline]
-    fn in_support(&self, y : &Loc<F, N>) -> bool {
+    fn in_support(&self, y: &Loc<N, F>) -> bool {
         let r = self.get_r();
         y.iter().all(|x| x.abs() <= r)
     }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
         // It is not difficult to verify that [`HatConv`] is C^2.
         // Therefore, so is [`Convolution<CubeIndicator<R, N>, HatConv<C, N>>`] so that a finer
         // subdivision for the hint than this is not particularly useful.
@@ -510,22 +565,26 @@
     }
 }
 
-impl<F : Float, R, C, const N : usize> GlobalAnalysis<F, Bounds<F>>
-for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
-where R : Constant<Type=F>,
-      C : Constant<Type=F> {
+impl<F: Float, R, C, const N: usize> GlobalAnalysis<F, Bounds<F>>
+    for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
+where
+    R: Constant<Type = F>,
+    C: Constant<Type = F>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         Bounds(F::ZERO, self.apply(Loc::ORIGIN))
     }
 }
 
-impl<F : Float, R, C, const N : usize> LocalAnalysis<F, Bounds<F>, N>
-for Convolution<CubeIndicator<R, N>, HatConv<C,  N>>
-where R : Constant<Type=F>,
-      C : Constant<Type=F> {
+impl<F: Float, R, C, const N: usize> LocalAnalysis<F, Bounds<F>, N>
+    for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
+where
+    R: Constant<Type = F>,
+    C: Constant<Type = F>,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
         // The function is maximised/minimised where the absolute value is minimised/maximised.
         let lower = self.apply(cube.maxnorm_point());
         let upper = self.apply(cube.minnorm_point());
@@ -542,7 +601,6 @@
     }
 }
 
-
 /// This [`BoundedBy`] implementation bounds $u * u$ by $(ψ * ψ) u$ for $u$ a hat convolution and
 /// $ψ = χ_{[-a,a]^N}$ for some $a>0$.
 ///
@@ -550,17 +608,18 @@
 /// where we take $ψ = χ_{[-a,a]^N}$ and $χ = χ_{[-σ,σ]^N}$ for $σ$ the width of the hat
 /// convolution.
 #[replace_float_literals(F::cast_from(literal))]
-impl<F, C, S, const N : usize>
-BoundedBy<F, SupportProductFirst<AutoConvolution<CubeIndicator<S, N>>, HatConv<C, N>>>
-for AutoConvolution<HatConv<C, N>>
-where F : Float,
-      C : Constant<Type=F>,
-      S : Constant<Type=F> {
-
+impl<F, C, S, const N: usize>
+    BoundedBy<F, SupportProductFirst<AutoConvolution<CubeIndicator<S, N>>, HatConv<C, N>>>
+    for AutoConvolution<HatConv<C, N>>
+where
+    F: Float,
+    C: Constant<Type = F>,
+    S: Constant<Type = F>,
+{
     fn bounding_factor(
         &self,
-        kernel : &SupportProductFirst<AutoConvolution<CubeIndicator<S, N>>, HatConv<C, N>>
-    ) -> Option<F> {
+        kernel: &SupportProductFirst<AutoConvolution<CubeIndicator<S, N>>, HatConv<C, N>>,
+    ) -> DynResult<F> {
         // We use the comparison $ℱ[𝒜(ψ v)] ≤ L_1 ℱ[𝒜(ψ)u] ⟺ I_{v̂} v̂ ≤ L_1 û$ with
         // $ψ = χ_{[-w, w]}$ satisfying $supp v ⊂ [-w, w]$, i.e. $w ≥ σ$. Here $v̂ = ℱ[v]$ and
         // $I_{v̂} = ∫ v̂ d ξ. For this relationship to be valid, we need $v̂ ≥ 0$, which is guaranteed
@@ -574,10 +633,10 @@
         // `SupportProductFirst<AutoConvolution<CubeIndicator<S, N>>, HatConv<C, N>>`
         // is wide enough, and that the hat convolution has the same radius as ours.
         if σ <= a && hatconv2 == &self.0 {
-            Some(bounding_1d.powi(N as i32))
+            Ok(bounding_1d.powi(N as i32))
         } else {
             // We cannot compare
-            None
+            Err(anyhow!("Incomparable factors"))
         }
     }
 }
@@ -586,46 +645,44 @@
 ///
 /// This is based on Example 3.3 in the manuscript.
 #[replace_float_literals(F::cast_from(literal))]
-impl<F, C, const N : usize>
-BoundedBy<F, HatConv<C, N>>
-for AutoConvolution<HatConv<C, N>>
-where F : Float,
-      C : Constant<Type=F> {
-
+impl<F, C, const N: usize> BoundedBy<F, HatConv<C, N>> for AutoConvolution<HatConv<C, N>>
+where
+    F: Float,
+    C: Constant<Type = F>,
+{
     /// Returns an estimate of the factor $L_1$.
     ///
     /// Returns  `None` if `kernel` does not have the same width as hat convolution that `self`
     /// is based on.
-    fn bounding_factor(
-        &self,
-        kernel : &HatConv<C, N>
-    ) -> Option<F> {
+    fn bounding_factor(&self, kernel: &HatConv<C, N>) -> DynResult<F> {
         if kernel == &self.0 {
-            Some(1.0)
+            Ok(1.0)
         } else {
             // We cannot compare
-            None
+            Err(anyhow!("Incomparable kernels"))
         }
     }
 }
 
 #[cfg(test)]
 mod tests {
+    use super::HatConv;
+    use crate::kernels::{BallIndicator, Convolution, CubeIndicator};
     use alg_tools::lingrid::linspace;
+    use alg_tools::loc::Loc;
     use alg_tools::mapping::Mapping;
     use alg_tools::norms::Linfinity;
-    use alg_tools::loc::Loc;
-    use crate::kernels::{BallIndicator, CubeIndicator, Convolution};
-    use super::HatConv;
 
     /// Tests numerically that [`HatConv<f64, 1>`] is monotone.
     #[test]
     fn hatconv_monotonicity() {
         let grid = linspace(0.0, 1.0, 100000);
-        let hatconv : HatConv<f64, 1> = HatConv{ radius : 1.0 };
+        let hatconv: HatConv<f64, 1> = HatConv { radius: 1.0 };
         let mut vals = grid.into_iter().map(|t| hatconv.apply(Loc::from(t)));
         let first = vals.next().unwrap();
-        let monotone = vals.fold((first, true), |(prev, ok), t| (prev, ok && prev >= t)).1;
+        let monotone = vals
+            .fold((first, true), |(prev, ok), t| (prev, ok && prev >= t))
+            .1;
         assert!(monotone);
     }
 
@@ -633,21 +690,27 @@
     #[test]
     fn convolution_cubeind_hatconv_monotonicity() {
         let grid = linspace(-2.0, 0.0, 100000);
-        let hatconv : Convolution<CubeIndicator<f64, 1>, HatConv<f64, 1>>
-            = Convolution(BallIndicator { r : 0.5, exponent : Linfinity },
-                          HatConv{ radius : 1.0 } );
+        let hatconv: Convolution<CubeIndicator<f64, 1>, HatConv<f64, 1>> =
+            Convolution(BallIndicator { r: 0.5, exponent: Linfinity }, HatConv {
+                radius: 1.0,
+            });
         let mut vals = grid.into_iter().map(|t| hatconv.apply(Loc::from(t)));
         let first = vals.next().unwrap();
-        let monotone = vals.fold((first, true), |(prev, ok), t| (prev, ok && prev <= t)).1;
+        let monotone = vals
+            .fold((first, true), |(prev, ok), t| (prev, ok && prev <= t))
+            .1;
         assert!(monotone);
 
         let grid = linspace(0.0, 2.0, 100000);
-        let hatconv : Convolution<CubeIndicator<f64, 1>, HatConv<f64, 1>>
-            = Convolution(BallIndicator { r : 0.5, exponent : Linfinity },
-                          HatConv{ radius : 1.0 } );
+        let hatconv: Convolution<CubeIndicator<f64, 1>, HatConv<f64, 1>> =
+            Convolution(BallIndicator { r: 0.5, exponent: Linfinity }, HatConv {
+                radius: 1.0,
+            });
         let mut vals = grid.into_iter().map(|t| hatconv.apply(Loc::from(t)));
         let first = vals.next().unwrap();
-        let monotone = vals.fold((first, true), |(prev, ok), t| (prev, ok && prev >= t)).1;
+        let monotone = vals
+            .fold((first, true), |(prev, ok), t| (prev, ok && prev >= t))
+            .1;
         assert!(monotone);
     }
 }