diff -r 9738b51d90d7 -r 4f468d35fa29 src/kernels/ball_indicator.rs
--- a/src/kernels/ball_indicator.rs	Sun Apr 27 15:03:51 2025 -0500
+++ b/src/kernels/ball_indicator.rs	Thu Feb 26 11:38:43 2026 -0500
@@ -1,56 +1,44 @@
-
 //! Implementation of the indicator function of a ball with respect to various norms.
+use super::base::*;
+use crate::types::*;
+use alg_tools::bisection_tree::{Bounds, Constant, GlobalAnalysis, LocalAnalysis, Support};
+use alg_tools::coefficients::factorial;
+use alg_tools::euclidean::StaticEuclidean;
+use alg_tools::instance::Instance;
+use alg_tools::loc::Loc;
+use alg_tools::mapping::{DifferentiableImpl, Differential, LipschitzDifferentiableImpl, Mapping};
+use alg_tools::maputil::array_init;
+use alg_tools::norms::*;
+use alg_tools::sets::Cube;
+use anyhow::anyhow;
 use float_extras::f64::tgamma as gamma;
 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,
-};
-use alg_tools::mapping::{
-    Mapping,
-    Differential,
-    DifferentiableImpl,
-};
-use alg_tools::instance::Instance;
-use alg_tools::euclidean::StaticEuclidean;
-use alg_tools::maputil::array_init;
-use alg_tools::coefficients::factorial;
-use crate::types::*;
-use super::base::*;
 
 /// Representation of the indicator of the ball $𝔹_q = \\{ x ∈ ℝ^N \mid \\|x\\|\_q ≤ r \\}$,
 /// where $q$ is the `Exponent`, and $r$ is the radius [`Constant`] `C`.
-#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
-pub struct BallIndicator<C : Constant, Exponent : NormExponent, const N : usize> {
+#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
+pub struct BallIndicator<C: Constant, Exponent: NormExponent, const N: usize> {
     /// The radius of the ball.
-    pub r : C,
+    pub r: C,
     /// The exponent $q$ of the norm creating the ball
-    pub exponent : Exponent,
+    pub exponent: Exponent,
 }
 
 /// Alias for the representation of the indicator of the $∞$-norm-ball
 /// $𝔹_∞ = \\{ x ∈ ℝ^N \mid \\|x\\|\_∞ ≤ c \\}$.
-pub type CubeIndicator<C, const N : usize> = BallIndicator<C, Linfinity, N>;
+pub type CubeIndicator<C, const N: usize> = BallIndicator<C, Linfinity, N>;
 
 #[replace_float_literals(C::Type::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Mapping<Loc<C::Type, N>>
-for BallIndicator<C, Exponent, N>
+impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
+    Mapping<Loc<N, C::Type>> for BallIndicator<C, Exponent, N>
 where
-    Loc<F, N> : Norm<F, Exponent>
+    Loc<N, F>: Norm<Exponent, F>,
 {
     type Codomain = C::Type;
 
     #[inline]
-    fn apply<I : Instance<Loc<C::Type, N>>>(&self, x : I) -> Self::Codomain {
+    fn apply<I: Instance<Loc<N, C::Type>>>(&self, x: I) -> Self::Codomain {
         let r = self.r.value();
         let n = x.eval(|x| x.norm(self.exponent));
         if n <= r {
@@ -61,114 +49,105 @@
     }
 }
 
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-DifferentiableImpl<Loc<C::Type, N>>
-for BallIndicator<C, Exponent, N>
+impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
+    DifferentiableImpl<Loc<N, C::Type>> for BallIndicator<C, Exponent, N>
 where
-    C : Constant,
-     Loc<F, N> : Norm<F, Exponent>
+    C: Constant,
+    Loc<N, F>: Norm<Exponent, F>,
 {
-    type Derivative = Loc<C::Type, N>;
+    type Derivative = Loc<N, C::Type>;
 
     #[inline]
-    fn differential_impl<I : Instance<Loc<C::Type, N>>>(&self, _x : I) -> Self::Derivative {
+    fn differential_impl<I: Instance<Loc<N, C::Type>>>(&self, _x: I) -> Self::Derivative {
         Self::Derivative::origin()
     }
 }
 
-impl<F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Lipschitz<L2>
-for BallIndicator<C, Exponent, N>
-where C : Constant,
-      Loc<F, N> : Norm<F, Exponent> {
+impl<F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize> Lipschitz<L2>
+    for BallIndicator<C, Exponent, N>
+where
+    C: Constant,
+    Loc<N, F>: Norm<Exponent, F>,
+{
     type FloatType = C::Type;
 
-    fn lipschitz_factor(&self, _l2 : L2) -> Option<C::Type> {
-        None
+    fn lipschitz_factor(&self, _l2: L2) -> DynResult<C::Type> {
+        Err(anyhow!("Not a Lipschitz function"))
     }
 }
 
-impl<'b, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Lipschitz<L2>
-for Differential<'b, Loc<F, N>, BallIndicator<C, Exponent, N>>
-where C : Constant,
-      Loc<F, N> : Norm<F, Exponent> {
+impl<'b, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
+    LipschitzDifferentiableImpl<Loc<N, F>, L2> for BallIndicator<C, Exponent, N>
+where
+    C: Constant,
+    Loc<N, F>: Norm<Exponent, F>,
+{
     type FloatType = C::Type;
 
-    fn lipschitz_factor(&self, _l2 : L2) -> Option<C::Type> {
-        None
+    fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<C::Type> {
+        Err(anyhow!("Not a Lipschitz-differentiable function"))
     }
 }
 
-impl<'a, 'b, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Lipschitz<L2>
-for Differential<'b, Loc<F, N>, &'a BallIndicator<C, Exponent, N>>
-where C : Constant,
-      Loc<F, N> : Norm<F, Exponent> {
+impl<'b, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize> NormBounded<L2>
+    for Differential<'b, Loc<N, F>, BallIndicator<C, Exponent, N>>
+where
+    C: Constant,
+    Loc<N, F>: Norm<Exponent, F>,
+{
     type FloatType = C::Type;
 
-    fn lipschitz_factor(&self, _l2 : L2) -> Option<C::Type> {
-        None
-    }
-}
-
-
-impl<'b, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-NormBounded<L2>
-for Differential<'b, Loc<F, N>, BallIndicator<C, Exponent, N>>
-where C : Constant,
-      Loc<F, N> : Norm<F, Exponent> {
-    type FloatType = C::Type;
-
-    fn norm_bound(&self, _l2 : L2) -> C::Type {
+    fn norm_bound(&self, _l2: L2) -> C::Type {
         F::INFINITY
     }
 }
 
-impl<'a, 'b, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-NormBounded<L2>
-for Differential<'b, Loc<F, N>, &'a BallIndicator<C, Exponent, N>>
-where C : Constant,
-      Loc<F, N> : Norm<F, Exponent> {
+impl<'a, 'b, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
+    NormBounded<L2> for Differential<'b, Loc<N, F>, &'a BallIndicator<C, Exponent, N>>
+where
+    C: Constant,
+    Loc<N, F>: Norm<Exponent, F>,
+{
     type FloatType = C::Type;
 
-    fn norm_bound(&self, _l2 : L2) -> C::Type {
+    fn norm_bound(&self, _l2: L2) -> C::Type {
         F::INFINITY
     }
 }
 
-
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Support<C::Type, N>
-for BallIndicator<C, Exponent, N>
-where Loc<F, N> : Norm<F, Exponent>,
-      Linfinity : Dominated<F, Exponent, Loc<F, N>> {
-
+impl<'a, F: Float, C: Constant<Type = F>, Exponent, const N: usize> Support<N, C::Type>
+    for BallIndicator<C, Exponent, N>
+where
+    Exponent: NormExponent + Sync + Send + 'static,
+    Loc<N, F>: Norm<Exponent, F>,
+    Linfinity: Dominated<F, Exponent, Loc<N, F>>,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F,N> {
+    fn support_hint(&self) -> Cube<N, F> {
         let r = Linfinity.from_norm(self.r.value(), self.exponent);
         array_init(|| [-r, r]).into()
     }
 
     #[inline]
-    fn in_support(&self, x : &Loc<F,N>) -> bool {
+    fn in_support(&self, x: &Loc<N, F>) -> bool {
         let r = Linfinity.from_norm(self.r.value(), self.exponent);
         x.norm(self.exponent) <= r
     }
 
     /// This can only really work in a reasonable fashion for N=1.
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
         let r = Linfinity.from_norm(self.r.value(), self.exponent);
         cube.map(|a, b| symmetric_interval_hint(r, a, b))
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-GlobalAnalysis<F, Bounds<F>>
-for BallIndicator<C, Exponent, N>
-where Loc<F, N> : Norm<F, Exponent> {
+impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
+    GlobalAnalysis<F, Bounds<F>> for BallIndicator<C, Exponent, N>
+where
+    Loc<N, F>: Norm<Exponent, F>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         Bounds(0.0, 1.0)
@@ -176,29 +155,28 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Norm<F, Linfinity>
-for BallIndicator<C, Exponent, N>
-where Loc<F, N> : Norm<F, Exponent> {
+impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize> Norm<Linfinity, F>
+    for BallIndicator<C, Exponent, N>
+where
+    Loc<N, F>: Norm<Exponent, F>,
+{
     #[inline]
-    fn norm(&self, _ : Linfinity) -> F {
+    fn norm(&self, _: Linfinity) -> F {
         1.0
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
-Norm<F, L1>
-for BallIndicator<C, L1, N> {
+impl<'a, F: Float, C: Constant<Type = F>, const N: usize> Norm<L1, F> for BallIndicator<C, L1, N> {
     #[inline]
-    fn norm(&self, _ : L1) -> F {
+    fn norm(&self, _: L1) -> F {
         // Using https://en.wikipedia.org/wiki/Volume_of_an_n-ball#Balls_in_Lp_norms,
         // we have V_N^1(r) = (2r)^N / N!
         let r = self.r.value();
-        if N==1 {
+        if N == 1 {
             2.0 * r
-        } else if N==2 {
-            r*r
+        } else if N == 2 {
+            r * r
         } else {
             (2.0 * r).powi(N as i32) * F::cast_from(factorial(N))
         }
@@ -206,17 +184,15 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
-Norm<F, L1>
-for BallIndicator<C, L2, N> {
+impl<'a, F: Float, C: Constant<Type = F>, const N: usize> Norm<L1, F> for BallIndicator<C, L2, N> {
     #[inline]
-    fn norm(&self, _ : L1) -> F {
+    fn norm(&self, _: L1) -> F {
         // See https://en.wikipedia.org/wiki/Volume_of_an_n-ball#The_volume.
         let r = self.r.value();
         let π = F::PI;
-        if N==1 {
+        if N == 1 {
             2.0 * r
-        } else if N==2 {
+        } else if N == 2 {
             π * (r * r)
         } else {
             let ndiv2 = F::cast_from(N) / 2.0;
@@ -227,94 +203,100 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
-Norm<F, L1>
-for BallIndicator<C, Linfinity, N> {
+impl<'a, F: Float, C: Constant<Type = F>, const N: usize> Norm<L1, F>
+    for BallIndicator<C, Linfinity, N>
+{
     #[inline]
-    fn norm(&self, _ : L1) -> F {
+    fn norm(&self, _: L1) -> F {
         let two_r = 2.0 * self.r.value();
         two_r.powi(N as i32)
     }
 }
 
-
 macro_rules! indicator_local_analysis {
     ($exponent:ident) => {
-        impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
-        LocalAnalysis<F, Bounds<F>, N>
-        for BallIndicator<C, $exponent, N>
-        where Loc<F, N> : Norm<F, $exponent>,
-            Linfinity : Dominated<F, $exponent, Loc<F, N>> {
+        impl<'a, F: Float, C: Constant<Type = F>, const N: usize> LocalAnalysis<F, Bounds<F>, N>
+            for BallIndicator<C, $exponent, N>
+        where
+            Loc<N, F>: Norm<$exponent, F>,
+            Linfinity: Dominated<F, $exponent, Loc<N, 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 2-norm is minimised/maximised.
                 let lower = self.apply(cube.maxnorm_point());
                 let upper = self.apply(cube.minnorm_point());
                 Bounds(lower, upper)
             }
         }
-    }
+    };
 }
 
 indicator_local_analysis!(L1);
 indicator_local_analysis!(L2);
 indicator_local_analysis!(Linfinity);
 
-
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, R, const N : usize> Mapping<Loc<F, N>>
-for AutoConvolution<CubeIndicator<R, N>>
-where R : Constant<Type=F> {
+impl<'a, F: Float, R, const N: usize> Mapping<Loc<N, F>> for AutoConvolution<CubeIndicator<R, N>>
+where
+    R: 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 two_r = 2.0 * self.0.r.value();
         // This is just a product of one-dimensional versions
-        y.cow().iter().map(|&x| {
-            0.0.max(two_r - x.abs())
-        }).product()
+        y.decompose()
+            .iter()
+            .map(|&x| 0.0.max(two_r - x.abs()))
+            .product()
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float, R, const N : usize> Support<F, N>
-for AutoConvolution<CubeIndicator<R, N>>
-where R : Constant<Type=F> {
+impl<F: Float, R, const N: usize> Support<N, F> for AutoConvolution<CubeIndicator<R, N>>
+where
+    R: Constant<Type = F>,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F, N> {
+    fn support_hint(&self) -> Cube<N, F> {
         let two_r = 2.0 * self.0.r.value();
         array_init(|| [-two_r, two_r]).into()
     }
 
     #[inline]
-    fn in_support(&self, y : &Loc<F, N>) -> bool {
+    fn in_support(&self, y: &Loc<N, F>) -> bool {
         let two_r = 2.0 * self.0.r.value();
         y.iter().all(|x| x.abs() <= two_r)
     }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
         let two_r = 2.0 * self.0.r.value();
         cube.map(|c, d| symmetric_interval_hint(two_r, c, d))
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float, R, const N : usize> GlobalAnalysis<F, Bounds<F>>
-for AutoConvolution<CubeIndicator<R, N>>
-where R : Constant<Type=F> {
+impl<F: Float, R, const N: usize> GlobalAnalysis<F, Bounds<F>>
+    for AutoConvolution<CubeIndicator<R, N>>
+where
+    R: Constant<Type = F>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         Bounds(0.0, self.apply(Loc::ORIGIN))
     }
 }
 
-impl<F : Float, R, const N : usize> LocalAnalysis<F, Bounds<F>, N>
-for AutoConvolution<CubeIndicator<R, N>>
-where R : Constant<Type=F> {
+impl<F: Float, R, const N: usize> LocalAnalysis<F, Bounds<F>, N>
+    for AutoConvolution<CubeIndicator<R, N>>
+where
+    R: 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());