diff -r 9738b51d90d7 -r 4f468d35fa29 src/kernels/base.rs
--- a/src/kernels/base.rs	Sun Apr 27 15:03:51 2025 -0500
+++ b/src/kernels/base.rs	Thu Feb 26 11:38:43 2026 -0500
@@ -1,28 +1,20 @@
-
 //! Things for constructing new kernels from component kernels and traits for analysing them
-use serde::Serialize;
 use numeric_literals::replace_float_literals;
+use serde::Serialize;
 
-use alg_tools::types::*;
+use alg_tools::bisection_tree::Support;
+use alg_tools::bounds::{Bounded, Bounds, GlobalAnalysis, LocalAnalysis};
+use alg_tools::instance::{Instance, Space};
+use alg_tools::loc::Loc;
+use alg_tools::mapping::{
+    DifferentiableImpl, DifferentiableMapping, LipschitzDifferentiableImpl, Mapping,
+};
+use alg_tools::maputil::{array_init, map1_indexed, map2};
 use alg_tools::norms::*;
-use alg_tools::loc::Loc;
 use alg_tools::sets::Cube;
-use alg_tools::bisection_tree::{
-    Support,
-    Bounds,
-    LocalAnalysis,
-    GlobalAnalysis,
-    Bounded,
-};
-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 alg_tools::types::*;
+use anyhow::anyhow;
 
 use crate::fourier::Fourier;
 use crate::types::*;
@@ -32,134 +24,129 @@
 /// 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)]
+#[derive(Copy, Clone, Serialize, Debug)]
 pub struct SupportProductFirst<A, B>(
     /// First kernel
     pub A,
     /// Second kernel
-    pub B
+    pub B,
 );
 
-impl<A, B, F : Float, const N : usize> Mapping<Loc<F, N>>
-for SupportProductFirst<A, B>
+impl<A, B, F: Float, const N: usize> Mapping<Loc<N, F>> for SupportProductFirst<A, B>
 where
-    A : Mapping<Loc<F, N>, Codomain = F>,
-    B : Mapping<Loc<F, N>, Codomain = F>,
+    A: Mapping<Loc<N, F>, Codomain = F>,
+    B: Mapping<Loc<N, F>, Codomain = F>,
 {
     type Codomain = F;
 
     #[inline]
-    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
-        self.0.apply(x.ref_instance()) * self.1.apply(x)
+    fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain {
+        x.eval_ref(|r| self.0.apply(r)) * self.1.apply(x)
     }
 }
 
-impl<A, B, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
-for SupportProductFirst<A, B>
+impl<A, B, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for SupportProductFirst<A, B>
 where
-    A : DifferentiableMapping<
-        Loc<F, N>,
-        DerivativeDomain=Loc<F, N>,
-        Codomain = F
-    >,
-    B : DifferentiableMapping<
-        Loc<F, N>,
-        DerivativeDomain=Loc<F, N>,
-        Codomain = F,
-    >
+    A: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>, Codomain = F>,
+    B: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>, Codomain = F>,
 {
-    type Derivative = Loc<F, N>;
+    type Derivative = Loc<N, F>;
 
     #[inline]
-    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)
+    fn differential_impl<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Derivative {
+        x.eval_ref(|xr| {
+            self.0.differential(xr) * self.1.apply(xr) + self.1.differential(xr) * self.0.apply(xr)
+        })
     }
 }
 
-impl<A, B, M : Copy, F : Float> Lipschitz<M>
-for SupportProductFirst<A, B>
-where A : Lipschitz<M, FloatType = F> + Bounded<F>,
-      B : Lipschitz<M, FloatType = F> + Bounded<F> {
+impl<A, B, M: Copy, F: Float> Lipschitz<M> for SupportProductFirst<A, B>
+where
+    A: Lipschitz<M, FloatType = F> + Bounded<F>,
+    B: Lipschitz<M, FloatType = F> + Bounded<F>,
+{
     type FloatType = F;
     #[inline]
-    fn lipschitz_factor(&self, m : M) -> Option<F> {
+    fn lipschitz_factor(&self, m: M) -> DynResult<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)
+        Ok(f.lipschitz_factor(m)? * g.bounds().uniform()
+            + g.lipschitz_factor(m)? * f.bounds().uniform())
     }
 }
 
-impl<'a, A, B, M : Copy, Domain, F : Float> Lipschitz<M>
-for Differential<'a, Domain, SupportProductFirst<A, B>>
+impl<'a, A, B, M: Copy, Domain, F: Float> LipschitzDifferentiableImpl<Domain, M>
+    for 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>
+    Domain: Space,
+    Self: DifferentiableImpl<Domain>,
+    A: DifferentiableMapping<Domain> + Lipschitz<M, FloatType = F> + Bounded<F>,
+    B: 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 lipschitz_factor(&self, m : M) -> Option<F> {
+    fn diff_lipschitz_factor(&self, m: M) -> DynResult<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 &SupportProductFirst(ref f, ref g) = self;
         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()
+        Ok([
+            df.lipschitz_factor(m)? * g.bounds().uniform(),
+            dg.lipschitz_factor(m)? * f.bounds().uniform(),
+            f.lipschitz_factor(m)? * dg.norm_bound(L2),
+            g.lipschitz_factor(m)? * df.norm_bound(L2),
+        ]
+        .into_iter()
+        .sum())
     }
 }
 
-
-impl<'a, A, B, F : Float, const N : usize> Support<F, N>
-for SupportProductFirst<A, B>
+impl<'a, A, B, F: Float, const N: usize> Support<N, F> for SupportProductFirst<A, B>
 where
-    A : Support<F, N>,
-    B : Support<F, N>
+    A: Support<N, F>,
+    B: Support<N, F>,
 {
     #[inline]
-    fn support_hint(&self) -> Cube<F, N> {
+    fn support_hint(&self) -> Cube<N, F> {
         self.0.support_hint()
     }
 
     #[inline]
-    fn in_support(&self, x : &Loc<F, N>) -> bool {
+    fn in_support(&self, x: &Loc<N, F>) -> bool {
         self.0.in_support(x)
     }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
         self.0.bisection_hint(cube)
     }
 }
 
-impl<'a, A, B, F : Float> GlobalAnalysis<F, Bounds<F>>
-for SupportProductFirst<A, B>
-where A : GlobalAnalysis<F, Bounds<F>>,
-      B : GlobalAnalysis<F, Bounds<F>> {
+impl<'a, A, B, F: Float> GlobalAnalysis<F, Bounds<F>> for SupportProductFirst<A, B>
+where
+    A: GlobalAnalysis<F, Bounds<F>>,
+    B: GlobalAnalysis<F, Bounds<F>>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         self.0.global_analysis() * self.1.global_analysis()
     }
 }
 
-impl<'a, A, B, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N>
-for SupportProductFirst<A, B>
-where A : LocalAnalysis<F, Bounds<F>, N>,
-      B : LocalAnalysis<F, Bounds<F>, N> {
+impl<'a, A, B, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N>
+    for SupportProductFirst<A, B>
+where
+    A: LocalAnalysis<F, Bounds<F>, N>,
+    B: LocalAnalysis<F, Bounds<F>, N>,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
         self.0.local_analysis(cube) * self.1.local_analysis(cube)
     }
 }
@@ -169,125 +156,116 @@
 /// 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)]
+#[derive(Copy, Clone, Serialize, Debug)]
 pub struct SupportSum<A, B>(
     /// First kernel
     pub A,
     /// Second kernel
-    pub B
+    pub B,
 );
 
-impl<'a, A, B, F : Float, const N : usize> Mapping<Loc<F, N>>
-for SupportSum<A, B>
+impl<'a, A, B, F: Float, const N: usize> Mapping<Loc<N, F>> for SupportSum<A, B>
 where
-    A : Mapping<Loc<F, N>, Codomain = F>,
-    B : Mapping<Loc<F, N>, Codomain = F>,
+    A: Mapping<Loc<N, F>, Codomain = F>,
+    B: Mapping<Loc<N, F>, Codomain = F>,
 {
     type Codomain = F;
 
     #[inline]
-    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
-        self.0.apply(x.ref_instance()) + self.1.apply(x)
+    fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain {
+        x.eval_ref(|r| self.0.apply(r)) + self.1.apply(x)
     }
 }
 
-impl<'a, A, B, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
-for SupportSum<A, B>
+impl<'a, A, B, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for SupportSum<A, B>
 where
-    A : DifferentiableMapping<
-        Loc<F, N>,
-        DerivativeDomain = Loc<F, N>
-    >,
-    B : DifferentiableMapping<
-        Loc<F, N>,
-        DerivativeDomain = Loc<F, N>,
-    >
+    A: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>>,
+    B: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>>,
 {
-
-    type Derivative = Loc<F, N>;
+    type Derivative = Loc<N, F>;
 
     #[inline]
-    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
-        self.0.differential(x.ref_instance()) + self.1.differential(x)
+    fn differential_impl<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Derivative {
+        x.eval_ref(|r| self.0.differential(r)) + 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 {
-
+impl<'a, A, B, F: Float, const N: usize> Support<N, F> for SupportSum<A, B>
+where
+    A: Support<N, F>,
+    B: Support<N, F>,
+    Cube<N, F>: SetOrd,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F, N> {
+    fn support_hint(&self) -> Cube<N, F> {
         self.0.support_hint().common(&self.1.support_hint())
     }
 
     #[inline]
-    fn in_support(&self, x : &Loc<F, N>) -> bool {
+    fn in_support(&self, x: &Loc<N, F>) -> bool {
         self.0.in_support(x) || self.1.in_support(x)
     }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
-        map2(self.0.bisection_hint(cube),
-             self.1.bisection_hint(cube),
-             |a, b| a.or(b))
+    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
+        map2(
+            self.0.bisection_hint(cube),
+            self.1.bisection_hint(cube),
+            |a, b| a.or(b),
+        )
     }
 }
 
-impl<'a, A, B, F : Float> GlobalAnalysis<F, Bounds<F>>
-for SupportSum<A, B>
-where A : GlobalAnalysis<F, Bounds<F>>,
-      B : GlobalAnalysis<F, Bounds<F>> {
+impl<'a, A, B, F: Float> GlobalAnalysis<F, Bounds<F>> for SupportSum<A, B>
+where
+    A: GlobalAnalysis<F, Bounds<F>>,
+    B: GlobalAnalysis<F, Bounds<F>>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         self.0.global_analysis() + self.1.global_analysis()
     }
 }
 
-impl<'a, A, B, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N>
-for SupportSum<A, B>
-where A : LocalAnalysis<F, Bounds<F>, N>,
-      B : LocalAnalysis<F, Bounds<F>, N>,
-      Cube<F, N> : SetOrd {
+impl<'a, A, B, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for SupportSum<A, B>
+where
+    A: LocalAnalysis<F, Bounds<F>, N>,
+    B: LocalAnalysis<F, Bounds<F>, N>,
+    Cube<N, F>: SetOrd,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
         self.0.local_analysis(cube) + self.1.local_analysis(cube)
     }
 }
 
-impl<F : Float, M : Copy, A, B> Lipschitz<M> for SupportSum<A, B>
-where A : Lipschitz<M, FloatType = F>,
-      B : Lipschitz<M, FloatType = F> {
+impl<F: Float, M: Copy, A, B> Lipschitz<M> for SupportSum<A, B>
+where
+    A: Lipschitz<M, FloatType = F>,
+    B: Lipschitz<M, FloatType = F>,
+{
     type FloatType = F;
 
-    fn lipschitz_factor(&self, m : M) -> Option<F> {
-        match (self.0.lipschitz_factor(m), self.1.lipschitz_factor(m)) {
-            (Some(l0), Some(l1)) => Some(l0 + l1),
-            _ => None
-        }
+    fn lipschitz_factor(&self, m: M) -> DynResult<F> {
+        Ok(self.0.lipschitz_factor(m)? * self.1.lipschitz_factor(m)?)
     }
 }
 
-impl<'b, F : Float, M : Copy, A, B, Domain> Lipschitz<M>
-for Differential<'b, Domain, SupportSum<A, B>>
+impl<'b, F: Float, M: Copy, A, B, Domain> LipschitzDifferentiableImpl<Domain, M>
+    for 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>
+    Domain: Space,
+    Self: DifferentiableImpl<Domain>,
+    A: DifferentiableMapping<Domain, Codomain = F>,
+    B: 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)
+    fn diff_lipschitz_factor(&self, m: M) -> DynResult<F> {
+        Ok(self.0.diff_ref().lipschitz_factor(m)? + self.1.diff_ref().lipschitz_factor(m)?)
     }
 }
 
@@ -296,48 +274,52 @@
 /// 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)]
+#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
 pub struct Convolution<A, B>(
     /// First kernel
     pub A,
     /// Second kernel
-    pub B
+    pub B,
 );
 
-impl<F : Float, M, A, B> Lipschitz<M> for Convolution<A, B>
-where A : Norm<F, L1> ,
-      B : Lipschitz<M, FloatType = F> {
+impl<F: Float, M, A, B> Lipschitz<M> for Convolution<A, B>
+where
+    A: Norm<L1, F>,
+    B: Lipschitz<M, FloatType = F>,
+{
     type FloatType = F;
 
-    fn lipschitz_factor(&self, m : M) -> Option<F> {
+    fn lipschitz_factor(&self, m: M) -> DynResult<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.
-        self.1.lipschitz_factor(m).map(|l| l * self.0.norm(L1))
+        Ok(self.1.lipschitz_factor(m)? * self.0.norm(L1))
     }
 }
 
-impl<'b, F : Float, M, A, B, Domain> Lipschitz<M>
-for Differential<'b, Domain, Convolution<A, B>>
+impl<'b, F: Float, M, A, B, Domain> LipschitzDifferentiableImpl<Domain, M> for 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>
+    Domain: Space,
+    Self: DifferentiableImpl<Domain>,
+    A: Norm<L1, F>,
+    Convolution<A, B>: DifferentiableMapping<Domain, Codomain = F>,
+    B: DifferentiableMapping<Domain, Codomain = F>,
+    for<'a> B::Differential<'a>: Lipschitz<M, FloatType = F>,
 {
     type FloatType = F;
 
-    fn lipschitz_factor(&self, m : M) -> Option<F> {
+    fn diff_lipschitz_factor(&self, m: M) -> DynResult<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))
+        self.1
+            .diff_ref()
+            .lipschitz_factor(m)
+            .map(|l| l * self.0.norm(L1))
     }
 }
 
@@ -346,78 +328,76 @@
 /// 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)]
+#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
 pub struct AutoConvolution<A>(
     /// The kernel to be autoconvolved
-    pub A
+    pub A,
 );
 
-impl<F : Float, M, C> Lipschitz<M> for AutoConvolution<C>
-where C : Lipschitz<M, FloatType = F> + Norm<F, L1> {
+impl<F: Float, M, C> Lipschitz<M> for AutoConvolution<C>
+where
+    C: Lipschitz<M, FloatType = F> + Norm<L1, F>,
+{
     type FloatType = F;
 
-    fn lipschitz_factor(&self, m : M) -> Option<F> {
+    fn lipschitz_factor(&self, m: M) -> DynResult<F> {
         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>>
+impl<'b, F: Float, M, C, Domain> LipschitzDifferentiableImpl<Domain, M> for 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>
+    Domain: Space,
+    Self: DifferentiableImpl<Domain>,
+    C: Norm<L1, F> + 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))
+    fn diff_lipschitz_factor(&self, m: M) -> DynResult<F> {
+        self.0
+            .diff_ref()
+            .lipschitz_factor(m)
+            .map(|l| l * self.0.norm(L1))
     }
 }
 
-
 /// 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`]
-/// on [`Loc<F, 1>`]. Then the product implements them on [`Loc<F, N>`].
-#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
+/// on [`Loc<1, F>`]. Then the product implements them on [`Loc<N, F>`].
+#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
 #[allow(dead_code)]
-struct UniformProduct<G, const N : usize>(
+struct UniformProduct<G, const N: usize>(
     /// The one-dimensional kernel
-    G
+    G,
 );
 
-impl<'a, G, F : Float, const N : usize> Mapping<Loc<F, N>>
-for UniformProduct<G, N>
+impl<'a, G, F: Float, const N: usize> Mapping<Loc<N, F>> for UniformProduct<G, N>
 where
-    G : Mapping<Loc<F, 1>, Codomain = F>
+    G: Mapping<Loc<1, F>, Codomain = F>,
 {
     type Codomain = F;
 
     #[inline]
-    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> F {
-        x.cow().iter().map(|&y| self.0.apply(Loc([y]))).product()
+    fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> F {
+        x.decompose()
+            .iter()
+            .map(|&y| self.0.apply(Loc([y])))
+            .product()
     }
 }
 
-
-
-impl<'a, G, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
-for UniformProduct<G, N>
+impl<'a, G, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for UniformProduct<G, N>
 where
-    G : DifferentiableMapping<
-        Loc<F, 1>,
-        DerivativeDomain = F,
-        Codomain = F,
-    >
+    G: DifferentiableMapping<Loc<1, F>, DerivativeDomain = F, Codomain = F>,
 {
-    type Derivative = Loc<F, N>;
+    type Derivative = Loc<N, F>;
 
     #[inline]
-    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x0 : I) -> Loc<F, N> {
+    fn differential_impl<I: Instance<Loc<N, F>>>(&self, x0: I) -> Loc<N, F> {
         x0.eval(|x| {
             let vs = x.map(|y| self.0.apply(Loc([y])));
             product_differential(x, &vs, |y| self.0.differential(Loc([y])))
@@ -430,17 +410,19 @@
 /// The vector `x` is the location, `vs` consists of the values `g(x_i)`, and
 /// `gd` calculates the derivative `g'`.
 #[inline]
-pub(crate) fn product_differential<F : Float, G : Fn(F) -> F, const N : usize>(
-    x : &Loc<F, N>,
-    vs : &Loc<F, N>,
-    gd : G
-) -> Loc<F, N> {
+pub(crate) fn product_differential<F: Float, G: Fn(F) -> F, const N: usize>(
+    x: &Loc<N, F>,
+    vs: &Loc<N, F>,
+    gd: G,
+) -> Loc<N, F> {
     map1_indexed(x, |i, &y| {
-        gd(y) * vs.iter()
-                  .zip(0..)
-                  .filter_map(|(v, j)| (j != i).then_some(*v))
-                  .product()
-    }).into()
+        gd(y)
+            * vs.iter()
+                .zip(0..)
+                .filter_map(|(v, j)| (j != i).then_some(*v))
+                .product()
+    })
+    .into()
 }
 
 /// Helper function to calulate the Lipschitz factor of $∇f$ for $f(x)=∏_{i=1}^N g(x_i)$.
@@ -448,12 +430,12 @@
 /// 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 {
+pub(crate) fn product_differential_lipschitz_factor<F: Float, const N: usize>(
+    bound: F,
+    lip: F,
+    dbound: F,
+    dlip: F,
+) -> DynResult<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
@@ -470,31 +452,33 @@
     //                 = 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
+        Ok(F::cast_from(N)
+            * bound.powi((N - 2) as i32)
+            * (dlip * bound + F::cast_from(N - 1) * lip * dbound))
+    } else if N == 1 {
+        Ok(dlip)
     } else {
-        panic!("Invalid dimension")
+        Err(anyhow!("Invalid dimension"))
     }
 }
 
-impl<G, F : Float, const N : usize> Support<F, N>
-for UniformProduct<G, N>
-where G : Support<F, 1> {
+impl<G, F: Float, const N: usize> Support<N, F> for UniformProduct<G, N>
+where
+    G: Support<1, F>,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F, N> {
-        let [a] : [[F; 2]; 1] = self.0.support_hint().into();
+    fn support_hint(&self) -> Cube<N, F> {
+        let [a]: [[F; 2]; 1] = self.0.support_hint().into();
         array_init(|| a.clone()).into()
     }
 
     #[inline]
-    fn in_support(&self, x : &Loc<F, N>) -> bool {
+    fn in_support(&self, x: &Loc<N, F>) -> bool {
         x.iter().all(|&y| self.0.in_support(&Loc([y])))
     }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
         cube.map(|a, b| {
             let [h] = self.0.bisection_hint(&[[a, b]].into());
             h
@@ -502,9 +486,10 @@
     }
 }
 
-impl<G, F : Float, const N : usize> GlobalAnalysis<F, Bounds<F>>
-for UniformProduct<G, N>
-where G : GlobalAnalysis<F, Bounds<F>> {
+impl<G, F: Float, const N: usize> GlobalAnalysis<F, Bounds<F>> for UniformProduct<G, N>
+where
+    G: GlobalAnalysis<F, Bounds<F>>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         let g = self.0.global_analysis();
@@ -512,88 +497,91 @@
     }
 }
 
-impl<G, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N>
-for UniformProduct<G, N>
-where G : LocalAnalysis<F, Bounds<F>, 1> {
+impl<G, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for UniformProduct<G, N>
+where
+    G: LocalAnalysis<F, Bounds<F>, 1>,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
-        cube.iter_coords().map(
-            |&[a, b]| self.0.local_analysis(&([[a, b]].into()))
-        ).product()
+    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
+        cube.iter_coords()
+            .map(|&[a, b]| self.0.local_analysis(&([[a, b]].into())))
+            .product()
     }
 }
 
 macro_rules! product_lpnorm {
     ($lp:ident) => {
-        impl<G, F : Float, const N : usize> Norm<F, $lp>
-        for UniformProduct<G, N>
-        where G : Norm<F, $lp> {
+        impl<G, F: Float, const N: usize> Norm<$lp, F> for UniformProduct<G, N>
+        where
+            G: Norm<$lp, F>,
+        {
             #[inline]
-            fn norm(&self, lp : $lp) -> F {
+            fn norm(&self, lp: $lp) -> F {
                 self.0.norm(lp).powi(N as i32)
             }
         }
-    }
+    };
 }
 
 product_lpnorm!(L1);
 product_lpnorm!(L2);
 product_lpnorm!(Linfinity);
 
-
 /// Trait for bounding one kernel with respect to another.
 ///
 /// The type `F` is the scalar field, and `T` another kernel to which `Self` is compared.
-pub trait BoundedBy<F : Num, T> {
+pub trait BoundedBy<F: Num, T> {
     /// Calclate a bounding factor $c$ such that the Fourier transforms $ℱ\[v\] ≤ c ℱ\[u\]$ for
     /// $v$ `self` and $u$ `other`.
     ///
     /// If no such factors exits, `None` is returned.
-    fn bounding_factor(&self, other : &T) -> Option<F>;
+    fn bounding_factor(&self, other: &T) -> DynResult<F>;
 }
 
 /// This [`BoundedBy`] implementation bounds $(uv) * (uv)$ by $(ψ * ψ) u$.
 #[replace_float_literals(F::cast_from(literal))]
-impl<F, C, BaseP>
-BoundedBy<F, SupportProductFirst<AutoConvolution<C>, BaseP>>
-for AutoConvolution<SupportProductFirst<C, BaseP>>
-where F : Float,
-      C : Clone + PartialEq,
-      BaseP :  Fourier<F> + PartialOrd, // TODO: replace by BoundedBy,
-      <BaseP as Fourier<F>>::Transformed : Bounded<F> + Norm<F, L1> {
-
-    fn bounding_factor(&self, kernel : &SupportProductFirst<AutoConvolution<C>, BaseP>) -> Option<F> {
+impl<F, C, BaseP> BoundedBy<F, SupportProductFirst<AutoConvolution<C>, BaseP>>
+    for AutoConvolution<SupportProductFirst<C, BaseP>>
+where
+    F: Float,
+    C: PartialEq,
+    BaseP: Fourier<F> + PartialOrd, // TODO: replace by BoundedBy,
+    <BaseP as Fourier<F>>::Transformed: Bounded<F> + Norm<L1, F>,
+{
+    fn bounding_factor(
+        &self,
+        kernel: &SupportProductFirst<AutoConvolution<C>, BaseP>,
+    ) -> DynResult<F> {
         let SupportProductFirst(AutoConvolution(ref cutoff2), base_spread2) = kernel;
         let AutoConvolution(SupportProductFirst(ref cutoff, ref base_spread)) = self;
         let v̂ = base_spread.fourier();
 
         // Verify that the cut-off and ideal physical model (base spread) are the same.
-        if cutoff == cutoff2
-           && base_spread <= base_spread2
-           && v̂.bounds().lower() >= 0.0 {
+        if cutoff == cutoff2 && base_spread <= base_spread2 && v̂.bounds().lower() >= 0.0 {
             // Calculate the factor between the convolution approximation
             // `AutoConvolution<SupportProductFirst<C, BaseP>>` of $A_*A$ and the
             // kernel of the seminorm. This depends on the physical model P being
             // `SupportProductFirst<C, BaseP>` with the kernel `K` being
             // a `SupportSum` involving `SupportProductFirst<AutoConvolution<C>, BaseP>`.
-            Some(v̂.norm(L1))
+            Ok(v̂.norm(L1))
         } else {
             // We cannot compare
-            None
+            Err(anyhow!("Incomprable kernels"))
         }
     }
 }
 
-impl<F : Float, A, B, C> BoundedBy<F, SupportSum<B, C>> for A
-where A : BoundedBy<F, B>,
-      C : Bounded<F> {
-
+impl<F: Float, A, B, C> BoundedBy<F, SupportSum<B, C>> for A
+where
+    A: BoundedBy<F, B>,
+    C: Bounded<F>,
+{
     #[replace_float_literals(F::cast_from(literal))]
-    fn bounding_factor(&self, SupportSum(ref kernel1, kernel2) : &SupportSum<B, C>) -> Option<F> {
+    fn bounding_factor(&self, SupportSum(ref kernel1, kernel2): &SupportSum<B, C>) -> DynResult<F> {
         if kernel2.bounds().lower() >= 0.0 {
             self.bounding_factor(kernel1)
         } else {
-            None
+            Err(anyhow!("Component kernel not lower-bounded by zero"))
         }
     }
 }
@@ -603,7 +591,7 @@
 /// It will attempt to place the subdivision point at $-r$ or $r$.
 /// If neither of these points lies within $[a, b]$, `None` is returned.
 #[inline]
-pub(super) fn symmetric_interval_hint<F : Float>(r : F, a : F, b : F) -> Option<F> {
+pub(super) fn symmetric_interval_hint<F: Float>(r: F, a: F, b: F) -> Option<F> {
     if a < -r && -r < b {
         Some(-r)
     } else if a < r && r < b {
@@ -622,7 +610,7 @@
 /// returned.
 #[replace_float_literals(F::cast_from(literal))]
 #[inline]
-pub(super) fn symmetric_peak_hint<F : Float>(r : F, a : F, b : F) -> Option<F> {
+pub(super) fn symmetric_peak_hint<F: Float>(r: F, a: F, b: F) -> Option<F> {
     let stage1 = if a < -r {
         if b <= -r {
             None
@@ -648,7 +636,7 @@
     // Ignore stage1 hint if either side of subdivision would be just a small fraction of the
     // interval
     match stage1 {
-        Some(h) if (h - a).min(b-h) >= 0.3 * r => Some(h),
-        _ => None
+        Some(h) if (h - a).min(b - h) >= 0.3 * r => Some(h),
+        _ => None,
     }
 }