diff -r aec67cdd6b14 -r efa60bc4f743 src/regularisation.rs
--- a/src/regularisation.rs	Tue Aug 01 10:32:12 2023 +0300
+++ b/src/regularisation.rs	Thu Aug 29 00:00:00 2024 -0500
@@ -34,9 +34,13 @@
 use crate::subproblem::{
     nonneg::quadratic_nonneg,
     unconstrained::quadratic_unconstrained,
+    l1squared_unconstrained::l1squared_unconstrained,
+    l1squared_nonneg::l1squared_nonneg
 };
 use alg_tools::iterate::AlgIteratorFactory;
 
+use std::cmp::Ordering::{Greater, Less, Equal};
+
 /// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other
 /// algorithms.
 ///
@@ -131,6 +135,23 @@
         config : &FBGenericConfig<F>
     ) -> usize;
 
+    /// Approximately solve the problem
+    /// <div>$$
+    ///     \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x)
+    /// $$</div>
+    /// for $G$ depending on the trait implementation.
+    ///
+    /// 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>
+    ) -> usize;
+
     /// Find a point where `d` may violate the tolerance `ε`.
     ///
     /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
@@ -151,8 +172,37 @@
     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`].
+    ///
+    /// 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
+    /// regulariser.
+    ///
+    /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check
+    /// terminating early. Otherwise returns a possibly violating point, the value of `d` there,
+    /// 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,
+    ) -> 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>;
 
+
     /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
     ///
     /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
@@ -169,6 +219,28 @@
           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`].
+    /// The [`DiscreteMeasure`]s `μ` 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>,
+        μ : &DiscreteMeasure<Loc<F, N>, F>,
+        τ : F,
+        ε : F,
+        config : &FBGenericConfig<F>,
+        radon_μ :&DiscreteMeasure<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>;
+
+
     /// TODO: document this
     fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>>;
 
@@ -197,6 +269,9 @@
           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;
 }
 
 #[replace_float_literals(F::cast_from(literal))]
@@ -215,30 +290,47 @@
     ) -> usize {
         let inner_tolerance = ε * config.inner.tolerance_mult;
         let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-        let inner_τ = config.inner.τ0 / mA_normest;
-        quadratic_nonneg(config.inner.method, mA, g, τ * self.α(), x,
-                         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>
+    ) -> 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)
+    }
+
+
     #[inline]
-    fn find_tolerance_violation<G, BT>(
+    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
     ) -> 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> {
         let τα = τ * self.α();
-        let keep_below = τα + ε;
-        let maximise_above = τα + ε * config.insertion_cutoff_factor;
+        let keep_below = τα + slack + ε;
+        let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
         let refinement_tolerance = ε * config.refinement.tolerance_mult;
 
-        // If preliminary check indicates that we are in bonds, and if it otherwise matches
+        // If preliminary check indicates that we are in bounds, and if it otherwise matches
         // the insertion strategy, skip insertion.
         if skip_by_rough_check && d.bounds().upper() <= keep_below {
             None
@@ -271,9 +363,9 @@
         return (
             bnd.lower() >= keep_supp_above
             ||
-            μ.iter_spikes().map(|&DeltaMeasure{ α : β, ref x }| {
-                (β == 0.0) || d.apply(x) >= keep_supp_above
-            }).all(std::convert::identity)
+            μ.iter_spikes().all(|&DeltaMeasure{ α, ref x }| {
+                (α == 0.0) || d.apply(x) >= keep_supp_above
+            })
          ) && (
             bnd.upper() <= keep_below
             ||
@@ -281,6 +373,50 @@
         )
     }
 
+    fn verify_merge_candidate_radonsq<G, BT>(
+        &self,
+        d : &mut BTFN<F, G, BT, N>,
+        μ : &DiscreteMeasure<Loc<F, N>, F>,
+        τ : F,
+        ε : F,
+        config : &FBGenericConfig<F>,
+        radon_μ :&DiscreteMeasure<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> {
+        let τα = τ * self.α();
+        let refinement_tolerance = ε * config.refinement.tolerance_mult;
+        let merge_tolerance = config.merge_tolerance_mult * ε;
+        let slack = radon_μ.norm(Radon);
+        let bnd = d.bounds();
+
+        return {
+            μ.both_matching(radon_μ)
+             .all(|(α, rα, x)| {
+                let v = d.apply(x);
+                let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
+                    Greater => (τα, τα),
+                    _ => (F::NEG_INFINITY, τα),
+                    // Less should not happen; treated as Equal
+                };
+                let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) {
+                    Greater => (slack, slack),
+                    Equal => (-slack, slack),
+                    Less => (-slack, -slack),
+                };
+                // TODO: both fail.
+                (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)
+         }
+    }
+
     fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>> {
         let τα = τ * self.α();
         Some(Bounds(τα - ε,  τα + ε))
@@ -310,9 +446,12 @@
           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.α()));
-        let τw = |x| τ.min((ε + d.apply(x))/self.α());
-        τw(z) - τw(y)
+        let w = |x| 1.0.min((ε + d.apply(x))/(τ * self.α()));
+        w(z) - w(y)
+    }
+
+    fn radon_norm_bound(&self, b : F) -> F {
+        b / self.α()
     }
 }
 
@@ -331,28 +470,43 @@
     ) -> usize {
         let inner_tolerance = ε * config.inner.tolerance_mult;
         let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-        let inner_τ = config.inner.τ0 / mA_normest;
-        quadratic_unconstrained(config.inner.method, mA, g, τ * self.α(), x,
-                                inner_τ, inner_it)
+        quadratic_unconstrained(mA, g, τ * self.α(), x,
+                                mA_normest, &config.inner, inner_it)
     }
 
-   fn find_tolerance_violation<G, BT>(
+    fn solve_findim_l1squared(
+        &self,
+        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)
+    }
+
+    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,
     ) -> 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> {
         let τα = τ * self.α();
-        let keep_below = τα + ε;
-        let keep_above = -τα - ε;
-        let maximise_above = τα + ε * config.insertion_cutoff_factor;
-        let minimise_below = -τα - ε * config.insertion_cutoff_factor;
+        let keep_below = τα + slack + ε;
+        let keep_above = -(τα + slack) - ε;
+        let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
+        let minimise_below = -(τα + slack) - ε * config.insertion_cutoff_factor;
         let refinement_tolerance = ε * config.refinement.tolerance_mult;
 
         // If preliminary check indicates that we are in bonds, and if it otherwise matches
@@ -405,14 +559,13 @@
         return (
             (bnd.lower() >= keep_supp_pos_above && bnd.upper() <= keep_supp_neg_below)
             ||
-            μ.iter_spikes().map(|&DeltaMeasure{ α : β, ref x }| {
-                use std::cmp::Ordering::*;
+            μ.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,
                 }
-            }).all(std::convert::identity)
+            })
         ) && (
             bnd.upper() <= keep_below
             ||
@@ -426,6 +579,56 @@
         )
     }
 
+    fn verify_merge_candidate_radonsq<G, BT>(
+        &self,
+        d : &mut BTFN<F, G, BT, N>,
+        μ : &DiscreteMeasure<Loc<F, N>, F>,
+        τ : F,
+        ε : F,
+        config : &FBGenericConfig<F>,
+        radon_μ : &DiscreteMeasure<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> {
+        let τα = τ * self.α();
+        let refinement_tolerance = ε * config.refinement.tolerance_mult;
+        let merge_tolerance = config.merge_tolerance_mult * ε;
+        let slack = radon_μ.norm(Radon);
+        let bnd = d.bounds();
+
+        return {
+            μ.both_matching(radon_μ)
+             .all(|(α, rα, x)| {
+                let v = d.apply(x);
+                let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
+                    Greater => (τα, τα),
+                    Equal => (-τα, τα),
+                    Less => (-τα, -τα),
+                };
+                let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) {
+                    Greater => (slack, slack),
+                    Equal => (-slack, slack),
+                    Less => (-slack, -slack),
+                };
+                (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)
+        } && {
+            let keep_above = -τα - slack - merge_tolerance;
+            bnd.lower() >= keep_above
+            ||
+            d.has_lower_bound(keep_above, refinement_tolerance,
+                              config.refinement.max_steps)
+        }
+    }
+
     fn target_bounds(&self, τ : F, ε : F) -> Option<Bounds<F>> {
         let τα = τ * self.α();
         Some(Bounds(-τα - ε,  τα + ε))
@@ -457,14 +660,14 @@
                            + LocalAnalysis<F, Bounds<F>, N> {
 
         let α = self.α();
-        // let w = |x| {
-        //     let dx = d.apply(x);
-        //     ((-ε + dx)/(τ * α)).max(1.0.min(ε + dx)/(τ * α))
-        // };
-        let τw = |x| {
+        let w = |x| {
             let dx = d.apply(x);
-            ((-ε + dx)/α).max(τ.min(ε + dx)/α)
+            ((-ε + dx)/(τ * α)).max(1.0.min(ε + dx)/(τ * α))
         };
-        τw(z) - τw(y)
+        w(z) - w(y)
+    }
+
+    fn radon_norm_bound(&self, b : F) -> F {
+        b / self.α()
     }
 }
\ No newline at end of file