diff -r 39c5e6c7759d -r 0693cc9ba9f0 src/regularisation.rs
--- a/src/regularisation.rs	Mon Feb 17 13:45:11 2025 -0500
+++ b/src/regularisation.rs	Mon Feb 17 13:51:50 2025 -0500
@@ -2,53 +2,41 @@
 Regularisation terms
 */
 
-use numeric_literals::replace_float_literals;
-use serde::{Serialize, Deserialize};
-use alg_tools::norms::Norm;
-use alg_tools::linops::Mapping;
-use alg_tools::instance::Instance;
-use alg_tools::loc::Loc;
-use crate::types::*;
-use crate::measures::{
-    RNDM,
-    DeltaMeasure,
-    Radon,
-};
-use crate::fb::FBGenericConfig;
 #[allow(unused_imports)] // Used by documentation.
 use crate::fb::pointsource_fb_reg;
+use crate::fb::FBGenericConfig;
+use crate::measures::{DeltaMeasure, Radon, RNDM};
 #[allow(unused_imports)] // Used by documentation.
 use crate::sliding_fb::pointsource_sliding_fb_reg;
+use crate::types::*;
+use alg_tools::instance::Instance;
+use alg_tools::linops::Mapping;
+use alg_tools::loc::Loc;
+use alg_tools::norms::Norm;
+use numeric_literals::replace_float_literals;
+use serde::{Deserialize, Serialize};
 
-use nalgebra::{DVector, DMatrix};
-use alg_tools::nalgebra_support::ToNalgebraRealField;
+use crate::subproblem::{
+    l1squared_nonneg::l1squared_nonneg, l1squared_unconstrained::l1squared_unconstrained,
+    nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained,
+};
 use alg_tools::bisection_tree::{
-    BTFN,
-    Bounds,
-    BTSearch,
-    P2Minimise,
-    SupportGenerator,
-    LocalAnalysis,
-    Bounded,
-};
-use crate::subproblem::{
-    nonneg::quadratic_nonneg,
-    unconstrained::quadratic_unconstrained,
-    l1squared_unconstrained::l1squared_unconstrained,
-    l1squared_nonneg::l1squared_nonneg
+    BTSearch, Bounded, Bounds, LocalAnalysis, P2Minimise, SupportGenerator, BTFN,
 };
 use alg_tools::iterate::AlgIteratorFactory;
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use nalgebra::{DMatrix, DVector};
 
-use std::cmp::Ordering::{Greater, Less, Equal};
+use std::cmp::Ordering::{Equal, Greater, Less};
 
 /// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other
 /// algorithms.
 ///
 /// The only member of the struct is the regularisation parameter α.
 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
-pub struct NonnegRadonRegTerm<F : Float>(pub F /* α */);
+pub struct NonnegRadonRegTerm<F: Float>(pub F /* α */);
 
-impl<'a, F : Float> NonnegRadonRegTerm<F> {
+impl<'a, F: Float> NonnegRadonRegTerm<F> {
     /// Returns the regularisation parameter
     pub fn α(&self) -> F {
         let &NonnegRadonRegTerm(α) = self;
@@ -56,25 +44,24 @@
     }
 }
 
-impl<'a, F : Float, const N : usize> Mapping<RNDM<F, N>>
-for NonnegRadonRegTerm<F> {
+impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for NonnegRadonRegTerm<F> {
     type Codomain = F;
-    
-    fn apply<I>(&self, μ : I) -> F
-    where I : Instance<RNDM<F, N>> {
+
+    fn apply<I>(&self, μ: I) -> F
+    where
+        I: Instance<RNDM<F, N>>,
+    {
         self.α() * μ.eval(|x| x.norm(Radon))
     }
 }
 
-
 /// The regularisation term $α\|μ\|_{ℳ(Ω)}$ for [`pointsource_fb_reg`].
 ///
 /// The only member of the struct is the regularisation parameter α.
 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
-pub struct RadonRegTerm<F : Float>(pub F /* α */);
+pub struct RadonRegTerm<F: Float>(pub F /* α */);
 
-
-impl<'a, F : Float> RadonRegTerm<F> {
+impl<'a, F: Float> RadonRegTerm<F> {
     /// Returns the regularisation parameter
     pub fn α(&self) -> F {
         let &RadonRegTerm(α) = self;
@@ -82,31 +69,33 @@
     }
 }
 
-impl<'a, F : Float, const N : usize> Mapping<RNDM<F, N>>
-for RadonRegTerm<F> {
+impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for RadonRegTerm<F> {
     type Codomain = F;
-    
-    fn apply<I>(&self, μ : I) -> F
-    where I : Instance<RNDM<F, N>> {
+
+    fn apply<I>(&self, μ: I) -> F
+    where
+        I: Instance<RNDM<F, N>>,
+    {
         self.α() * μ.eval(|x| x.norm(Radon))
     }
 }
 
 /// Regularisation term configuration
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
-pub enum Regularisation<F : Float> {
+pub enum Regularisation<F: Float> {
     /// $α \\|μ\\|\_{ℳ(Ω)}$
     Radon(F),
     /// $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$
     NonnegRadon(F),
 }
 
-impl<'a, F : Float, const N : usize> Mapping<RNDM<F, N>>
-for Regularisation<F> {
+impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for Regularisation<F> {
     type Codomain = F;
 
-    fn apply<I>(&self, μ : I) -> F
-    where I : Instance<RNDM<F, N>> {
+    fn apply<I>(&self, μ: I) -> F
+    where
+        I: Instance<RNDM<F, N>>,
+    {
         match *self {
             Self::Radon(α) => RadonRegTerm(α).apply(μ),
             Self::NonnegRadon(α) => NonnegRadonRegTerm(α).apply(μ),
@@ -115,8 +104,9 @@
 }
 
 /// Abstraction of regularisation terms.
-pub trait RegTerm<F : Float + ToNalgebraRealField, const N : usize>
-: Mapping<RNDM<F, N>, Codomain = F> {
+pub trait RegTerm<F: Float + ToNalgebraRealField, const N: usize>:
+    Mapping<RNDM<F, N>, Codomain = F>
+{
     /// Approximately solve the problem
     /// <div>$$
     ///     \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + τ G(x)
@@ -129,13 +119,13 @@
     /// Returns the number of iterations taken.
     fn solve_findim(
         &self,
-        mA : &DMatrix<F::MixedType>,
-        g : &DVector<F::MixedType>,
-        τ : F,
-        x : &mut DVector<F::MixedType>,
-        mA_normest : F,
-        ε : F,
-        config : &FBGenericConfig<F>
+        mA: &DMatrix<F::MixedType>,
+        g: &DVector<F::MixedType>,
+        τ: F,
+        x: &mut DVector<F::MixedType>,
+        mA_normest: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> usize;
 
     /// Approximately solve the problem
@@ -147,12 +137,12 @@
     /// Returns the number of iterations taken.
     fn solve_findim_l1squared(
         &self,
-        y : &DVector<F::MixedType>,
-        g : &DVector<F::MixedType>,
-        τ : F,
-        x : &mut DVector<F::MixedType>,
-        ε : F,
-        config : &FBGenericConfig<F>
+        y: &DVector<F::MixedType>,
+        g: &DVector<F::MixedType>,
+        τ: F,
+        x: &mut DVector<F::MixedType>,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> usize;
 
     /// Find a point where `d` may violate the tolerance `ε`.
@@ -166,22 +156,24 @@
     /// and a boolean indicating whether the found point is in bounds.
     fn find_tolerance_violation<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        τ : F,
-        ε : F,
-        skip_by_rough_check : bool,
-        config : &FBGenericConfig<F>,
+        d: &mut BTFN<F, G, BT, N>,
+        τ: F,
+        ε: F,
+        skip_by_rough_check: bool,
+        config: &FBGenericConfig<F>,
     ) -> Option<(Loc<F, N>, F, bool)>
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO)
     }
 
     /// Find a point where `d` may violate the tolerance `ε`.
     ///
     /// This version includes a `slack` parameter to expand the tolerances.
-    /// It is used for Radon-norm squared proximal term in [`crate::radon_fb`].
+    /// It is used for Radon-norm squared proximal term in [`crate::prox_penalty::radon_squared`].
     ///
     /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
     /// are in bounds. `ε` is the current main tolerance and `τ` a scaling factor for the
@@ -192,56 +184,58 @@
     /// and a boolean indicating whether the found point is in bounds.
     fn find_tolerance_violation_slack<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        τ : F,
-        ε : F,
-        skip_by_rough_check : bool,
-        config : &FBGenericConfig<F>,
-        slack : F,
+        d: &mut BTFN<F, G, BT, N>,
+        τ: F,
+        ε: F,
+        skip_by_rough_check: bool,
+        config: &FBGenericConfig<F>,
+        slack: F,
     ) -> Option<(Loc<F, N>, F, bool)>
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>;
-
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
 
     /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
     ///
     /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
     fn verify_merge_candidate<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        μ : &RNDM<F, N>,
-        τ : F,
-        ε : F,
-        config : &FBGenericConfig<F>,
+        d: &mut BTFN<F, G, BT, N>,
+        μ: &RNDM<F, N>,
+        τ: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> bool
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>;
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
 
     /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
     ///
-    /// This version is s used for Radon-norm squared proximal term in [`crate::radon_fb`].
+    /// This version is s used for Radon-norm squared proximal term in
+    /// [`crate::prox_penalty::radon_squared`].
     /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have
     /// same coordinates at same agreeing indices.
     ///
     /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
     fn verify_merge_candidate_radonsq<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        μ : &RNDM<F, N>,
-        τ : F,
-        ε : F,
-        config : &FBGenericConfig<F>,
-        radon_μ :&RNDM<F, N>,
+        d: &mut BTFN<F, G, BT, N>,
+        μ: &RNDM<F, N>,
+        τ: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
+        radon_μ: &RNDM<F, N>,
     ) -> bool
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>;
-
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
 
     /// TODO: document this
-    fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>>;
+    fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
 
     /// Returns a scaling factor for the tolerance sequence.
     ///
@@ -250,78 +244,77 @@
 }
 
 /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`].
-pub trait SlidingRegTerm<F : Float + ToNalgebraRealField, const N : usize>
-: RegTerm<F, N> {
+pub trait SlidingRegTerm<F: Float + ToNalgebraRealField, const N: usize>: RegTerm<F, N> {
     /// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation
     /// term, such that $-ε ≤ τw - d ≤ ε$.
     fn goodness<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        μ : &RNDM<F, N>,
-        y : &Loc<F, N>,
-        z : &Loc<F, N>,
-        τ : F,
-        ε : F,
-        config : &FBGenericConfig<F>,
+        d: &mut BTFN<F, G, BT, N>,
+        μ: &RNDM<F, N>,
+        y: &Loc<F, N>,
+        z: &Loc<F, N>,
+        τ: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> F
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>;
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
 
     /// Convert bound on the regulariser to a bond on the Radon norm
-    fn radon_norm_bound(&self, b : F) -> F;
+    fn radon_norm_bound(&self, b: F) -> F;
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, const N : usize> RegTerm<F, N>
-for NonnegRadonRegTerm<F>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F> {
+impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for NonnegRadonRegTerm<F>
+where
+    Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+{
     fn solve_findim(
         &self,
-        mA : &DMatrix<F::MixedType>,
-        g : &DVector<F::MixedType>,
-        τ : F,
-        x : &mut DVector<F::MixedType>,
-        mA_normest : F,
-        ε : F,
-        config : &FBGenericConfig<F>
+        mA: &DMatrix<F::MixedType>,
+        g: &DVector<F::MixedType>,
+        τ: F,
+        x: &mut DVector<F::MixedType>,
+        mA_normest: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> usize {
         let inner_tolerance = ε * config.inner.tolerance_mult;
         let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-        quadratic_nonneg(mA, g, τ * self.α(), x,
-                         mA_normest, &config.inner, inner_it)
-
+        quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
     }
 
     fn solve_findim_l1squared(
         &self,
-        y : &DVector<F::MixedType>,
-        g : &DVector<F::MixedType>,
-        τ : F,
-        x : &mut DVector<F::MixedType>,
-        ε : F,
-        config : &FBGenericConfig<F>
+        y: &DVector<F::MixedType>,
+        g: &DVector<F::MixedType>,
+        τ: F,
+        x: &mut DVector<F::MixedType>,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> usize {
         let inner_tolerance = ε * config.inner.tolerance_mult;
         let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-        l1squared_nonneg(y, g, τ * self.α(), 1.0, x,
-                         &config.inner, inner_it)
+        l1squared_nonneg(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
     }
 
-
     #[inline]
     fn find_tolerance_violation_slack<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        τ : F,
-        ε : F,
-        skip_by_rough_check : bool,
-        config : &FBGenericConfig<F>,
-        slack : F
+        d: &mut BTFN<F, G, BT, N>,
+        τ: F,
+        ε: F,
+        skip_by_rough_check: bool,
+        config: &FBGenericConfig<F>,
+        slack: F,
     ) -> Option<(Loc<F, N>, F, bool)>
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         let τα = τ * self.α();
         let keep_above = -τα - slack - ε;
         let minimise_below = -τα - slack - ε * config.insertion_cutoff_factor;
@@ -333,22 +326,28 @@
             None
         } else {
             // If the rough check didn't indicate no insertion needed, find minimising point.
-            d.minimise_below(minimise_below, refinement_tolerance, config.refinement.max_steps)
-             .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
+            d.minimise_below(
+                minimise_below,
+                refinement_tolerance,
+                config.refinement.max_steps,
+            )
+            .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
         }
     }
 
     fn verify_merge_candidate<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        μ : &RNDM<F, N>,
-        τ : F,
-        ε : F,
-        config : &FBGenericConfig<F>,
+        d: &mut BTFN<F, G, BT, N>,
+        μ: &RNDM<F, N>,
+        τ: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> bool
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         let τα = τ * self.α();
         let refinement_tolerance = ε * config.refinement.tolerance_mult;
         let merge_tolerance = config.merge_tolerance_mult * ε;
@@ -356,31 +355,32 @@
         let keep_supp_below = -τα + merge_tolerance;
         let bnd = d.bounds();
 
-        return (
-            bnd.upper() <= keep_supp_below
-            ||
-            μ.iter_spikes().all(|&DeltaMeasure{ α, ref x }| {
-                (α == 0.0) || d.apply(x) <= keep_supp_below
-            })
-         ) && (
-            bnd.lower() >= keep_above
-            ||
-            d.has_lower_bound(keep_above, refinement_tolerance, config.refinement.max_steps)
-        )
+        return (bnd.upper() <= keep_supp_below
+            || μ
+                .iter_spikes()
+                .all(|&DeltaMeasure { α, ref x }| (α == 0.0) || d.apply(x) <= keep_supp_below))
+            && (bnd.lower() >= keep_above
+                || d.has_lower_bound(
+                    keep_above,
+                    refinement_tolerance,
+                    config.refinement.max_steps,
+                ));
     }
 
     fn verify_merge_candidate_radonsq<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        μ : &RNDM<F, N>,
-        τ : F,
-        ε : F,
-        config : &FBGenericConfig<F>,
-        radon_μ :&RNDM<F, N>,
+        d: &mut BTFN<F, G, BT, N>,
+        μ: &RNDM<F, N>,
+        τ: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
+        radon_μ: &RNDM<F, N>,
     ) -> bool
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         let τα = τ * self.α();
         let refinement_tolerance = ε * config.refinement.tolerance_mult;
         let merge_tolerance = config.merge_tolerance_mult * ε;
@@ -388,9 +388,8 @@
         let bnd = d.bounds();
 
         return {
-            μ.both_matching(radon_μ)
-             .all(|(α, rα, x)| {
-                let v = -d.apply(x); // TODO: observe ad hoc negation here, after minus_τv 
+            μ.both_matching(radon_μ).all(|(α, rα, x)| {
+                let v = -d.apply(x); // TODO: observe ad hoc negation here, after minus_τv
                                      // switch to τv.
                 let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
                     Greater => (τα, τα),
@@ -405,17 +404,20 @@
                 // TODO: both fail.
                 (l1 + l2 - merge_tolerance <= v) && (v <= u1 + u2 + merge_tolerance)
             })
-         } && {
+        } && {
             let keep_above = -τα - slack - merge_tolerance;
             bnd.lower() <= keep_above
-            ||
-            d.has_lower_bound(keep_above, refinement_tolerance, config.refinement.max_steps)
-         }
+                || d.has_lower_bound(
+                    keep_above,
+                    refinement_tolerance,
+                    config.refinement.max_steps,
+                )
+        };
     }
 
-    fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>> {
+    fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
         let τα = τ * self.α();
-        Some(Bounds(τα - ε,  τα + ε))
+        Some(Bounds(τα - ε, τα + ε))
     }
 
     fn tolerance_scaling(&self) -> F {
@@ -424,78 +426,82 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, const N : usize> SlidingRegTerm<F, N>
-for NonnegRadonRegTerm<F>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F> {
-
+impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for NonnegRadonRegTerm<F>
+where
+    Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+{
     fn goodness<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        _μ : &RNDM<F, N>,
-        y : &Loc<F, N>,
-        z : &Loc<F, N>,
-        τ : F,
-        ε : F,
-        _config : &FBGenericConfig<F>,
+        d: &mut BTFN<F, G, BT, N>,
+        _μ: &RNDM<F, N>,
+        y: &Loc<F, N>,
+        z: &Loc<F, N>,
+        τ: F,
+        ε: F,
+        _config: &FBGenericConfig<F>,
     ) -> F
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
-        let w = |x| 1.0.min((ε + d.apply(x))/(τ * self.α()));
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
+        let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α()));
         w(z) - w(y)
     }
 
-    fn radon_norm_bound(&self, b : F) -> F {
+    fn radon_norm_bound(&self, b: F) -> F {
         b / self.α()
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, const N : usize> RegTerm<F, N> for RadonRegTerm<F>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F> {
+impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for RadonRegTerm<F>
+where
+    Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+{
     fn solve_findim(
         &self,
-        mA : &DMatrix<F::MixedType>,
-        g : &DVector<F::MixedType>,
-        τ : F,
-        x : &mut DVector<F::MixedType>,
+        mA: &DMatrix<F::MixedType>,
+        g: &DVector<F::MixedType>,
+        τ: F,
+        x: &mut DVector<F::MixedType>,
         mA_normest: F,
-        ε : F,
-        config : &FBGenericConfig<F>
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> usize {
         let inner_tolerance = ε * config.inner.tolerance_mult;
         let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-        quadratic_unconstrained(mA, g, τ * self.α(), x,
-                                mA_normest, &config.inner, inner_it)
+        quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
     }
 
     fn solve_findim_l1squared(
         &self,
-        y : &DVector<F::MixedType>,
-        g : &DVector<F::MixedType>,
-        τ : F,
-        x : &mut DVector<F::MixedType>,
-        ε : F,
-        config : &FBGenericConfig<F>
+        y: &DVector<F::MixedType>,
+        g: &DVector<F::MixedType>,
+        τ: F,
+        x: &mut DVector<F::MixedType>,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> usize {
         let inner_tolerance = ε * config.inner.tolerance_mult;
         let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-        l1squared_unconstrained(y, g, τ * self.α(), 1.0, x,
-                                &config.inner, inner_it)
+        l1squared_unconstrained(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
     }
 
     fn find_tolerance_violation_slack<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        τ : F,
-        ε : F,
-        skip_by_rough_check : bool,
-        config : &FBGenericConfig<F>,
-        slack : F,
+        d: &mut BTFN<F, G, BT, N>,
+        τ: F,
+        ε: F,
+        skip_by_rough_check: bool,
+        config: &FBGenericConfig<F>,
+        slack: F,
     ) -> Option<(Loc<F, N>, F, bool)>
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         let τα = τ * self.α();
         let keep_below = τα + slack + ε;
         let keep_above = -(τα + slack) - ε;
@@ -509,10 +515,16 @@
             None
         } else {
             // If the rough check didn't indicate no insertion needed, find maximising point.
-            let mx = d.maximise_above(maximise_above, refinement_tolerance,
-                                      config.refinement.max_steps);
-            let mi = d.minimise_below(minimise_below, refinement_tolerance,
-                                      config.refinement.max_steps);
+            let mx = d.maximise_above(
+                maximise_above,
+                refinement_tolerance,
+                config.refinement.max_steps,
+            );
+            let mi = d.minimise_below(
+                minimise_below,
+                refinement_tolerance,
+                config.refinement.max_steps,
+            );
 
             match (mx, mi) {
                 (None, None) => None,
@@ -531,15 +543,17 @@
 
     fn verify_merge_candidate<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        μ : &RNDM<F, N>,
-        τ : F,
-        ε : F,
-        config : &FBGenericConfig<F>,
+        d: &mut BTFN<F, G, BT, N>,
+        μ: &RNDM<F, N>,
+        τ: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
     ) -> bool
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>,Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         let τα = τ * self.α();
         let refinement_tolerance = ε * config.refinement.tolerance_mult;
         let merge_tolerance = config.merge_tolerance_mult * ε;
@@ -549,41 +563,42 @@
         let keep_supp_neg_below = -τα + merge_tolerance;
         let bnd = d.bounds();
 
-        return (
-            (bnd.lower() >= keep_supp_pos_above && bnd.upper() <= keep_supp_neg_below)
-            ||
-            μ.iter_spikes().all(|&DeltaMeasure{ α : β, ref x }| {
-                match β.partial_cmp(&0.0) {
+        return ((bnd.lower() >= keep_supp_pos_above && bnd.upper() <= keep_supp_neg_below)
+            || μ
+                .iter_spikes()
+                .all(|&DeltaMeasure { α: β, ref x }| match β.partial_cmp(&0.0) {
                     Some(Greater) => d.apply(x) >= keep_supp_pos_above,
                     Some(Less) => d.apply(x) <= keep_supp_neg_below,
                     _ => true,
-                }
-            })
-        ) && (
-            bnd.upper() <= keep_below
-            ||
-            d.has_upper_bound(keep_below, refinement_tolerance,
-                              config.refinement.max_steps)
-        ) && (
-            bnd.lower() >= keep_above
-            ||
-            d.has_lower_bound(keep_above, refinement_tolerance,
-                              config.refinement.max_steps)
-        )
+                }))
+            && (bnd.upper() <= keep_below
+                || d.has_upper_bound(
+                    keep_below,
+                    refinement_tolerance,
+                    config.refinement.max_steps,
+                ))
+            && (bnd.lower() >= keep_above
+                || d.has_lower_bound(
+                    keep_above,
+                    refinement_tolerance,
+                    config.refinement.max_steps,
+                ));
     }
 
     fn verify_merge_candidate_radonsq<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        μ : &RNDM<F, N>,
-        τ : F,
-        ε : F,
-        config : &FBGenericConfig<F>,
-        radon_μ : &RNDM<F, N>,
+        d: &mut BTFN<F, G, BT, N>,
+        μ: &RNDM<F, N>,
+        τ: F,
+        ε: F,
+        config: &FBGenericConfig<F>,
+        radon_μ: &RNDM<F, N>,
     ) -> bool
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>,Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         let τα = τ * self.α();
         let refinement_tolerance = ε * config.refinement.tolerance_mult;
         let merge_tolerance = config.merge_tolerance_mult * ε;
@@ -591,8 +606,7 @@
         let bnd = d.bounds();
 
         return {
-            μ.both_matching(radon_μ)
-             .all(|(α, rα, x)| {
+            μ.both_matching(radon_μ).all(|(α, rα, x)| {
                 let v = d.apply(x);
                 let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
                     Greater => (τα, τα),
@@ -606,24 +620,28 @@
                 };
                 (l1 + l2 - merge_tolerance <= v) && (v <= u1 + u2 + merge_tolerance)
             })
-         } && {
+        } && {
             let keep_below = τα + slack + merge_tolerance;
             bnd.upper() <= keep_below
-            ||
-            d.has_upper_bound(keep_below, refinement_tolerance,
-                              config.refinement.max_steps)
+                || d.has_upper_bound(
+                    keep_below,
+                    refinement_tolerance,
+                    config.refinement.max_steps,
+                )
         } && {
             let keep_above = -τα - slack - merge_tolerance;
             bnd.lower() >= keep_above
-            ||
-            d.has_lower_bound(keep_above, refinement_tolerance,
-                              config.refinement.max_steps)
-        }
+                || d.has_lower_bound(
+                    keep_above,
+                    refinement_tolerance,
+                    config.refinement.max_steps,
+                )
+        };
     }
 
-    fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>> {
+    fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
         let τα = τ * self.α();
-        Some(Bounds(-τα - ε,  τα + ε))
+        Some(Bounds(-τα - ε, τα + ε))
     }
 
     fn tolerance_scaling(&self) -> F {
@@ -632,34 +650,34 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, const N : usize> SlidingRegTerm<F, N>
-for RadonRegTerm<F>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F> {
-
+impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for RadonRegTerm<F>
+where
+    Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+{
     fn goodness<G, BT>(
         &self,
-        d : &mut BTFN<F, G, BT, N>,
-        _μ : &RNDM<F, N>,
-        y : &Loc<F, N>,
-        z : &Loc<F, N>,
-        τ : F,
-        ε : F,
-        _config : &FBGenericConfig<F>,
+        d: &mut BTFN<F, G, BT, N>,
+        _μ: &RNDM<F, N>,
+        y: &Loc<F, N>,
+        z: &Loc<F, N>,
+        τ: F,
+        ε: F,
+        _config: &FBGenericConfig<F>,
     ) -> F
-    where BT : BTSearch<F, N, Agg=Bounds<F>>,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : Mapping<Loc<F, N>,Codomain=F>
-                           + LocalAnalysis<F, Bounds<F>, N> {
-
+    where
+        BT: BTSearch<F, N, Agg = Bounds<F>>,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    {
         let α = self.α();
         let w = |x| {
             let dx = d.apply(x);
-            ((-ε + dx)/(τ * α)).max(1.0.min(ε + dx)/(τ * α))
+            ((-ε + dx) / (τ * α)).max(1.0.min(ε + dx) / (τ * α))
         };
         w(z) - w(y)
     }
 
-    fn radon_norm_bound(&self, b : F) -> F {
+    fn radon_norm_bound(&self, b: F) -> F {
         b / self.α()
     }
-}
\ No newline at end of file
+}