diff -r 9738b51d90d7 -r 4f468d35fa29 src/frank_wolfe.rs
--- a/src/frank_wolfe.rs	Sun Apr 27 15:03:51 2025 -0500
+++ b/src/frank_wolfe.rs	Thu Feb 26 11:38:43 2026 -0500
@@ -13,82 +13,51 @@
     DOI: [10.1051/cocv/2011205](https://doi.org/0.1051/cocv/2011205).
 */
 
+use nalgebra::{DMatrix, DVector};
 use numeric_literals::replace_float_literals;
-use nalgebra::{DMatrix, DVector};
-use serde::{Serialize, Deserialize};
+use serde::{Deserialize, Serialize};
 //use colored::Colorize;
-
-use alg_tools::iterate::{
-    AlgIteratorFactory,
-    AlgIteratorOptions,
-    ValueIteratorFactory,
-};
-use alg_tools::euclidean::Euclidean;
-use alg_tools::norms::Norm;
-use alg_tools::linops::Mapping;
-use alg_tools::sets::Cube;
-use alg_tools::loc::Loc;
-use alg_tools::bisection_tree::{
-    BTFN,
-    Bounds,
-    BTNodeLookup,
-    BTNode,
-    BTSearch,
-    P2Minimise,
-    SupportGenerator,
-    LocalAnalysis,
-};
-use alg_tools::mapping::RealMapping;
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::L2;
-
-use crate::types::*;
-use crate::measures::{
-    RNDM,
-    DiscreteMeasure,
-    DeltaMeasure,
-    Radon,
-};
-use crate::measures::merging::{
-    SpikeMergingMethod,
-    SpikeMerging,
-};
+use crate::dataterm::QuadraticDataTerm;
 use crate::forward_model::ForwardModel;
+use crate::measures::merging::{SpikeMerging, SpikeMergingMethod};
+use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
+use crate::plot::Plotter;
+use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, RegTerm};
 #[allow(unused_imports)] // Used in documentation
 use crate::subproblem::{
-    unconstrained::quadratic_unconstrained,
-    nonneg::quadratic_nonneg,
-    InnerSettings,
-    InnerMethod,
+    nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained, InnerMethod, InnerSettings,
 };
 use crate::tolerance::Tolerance;
-use crate::plot::{
-    SeqPlotter,
-    Plotting,
-    PlotLookup
-};
-use crate::regularisation::{
-    NonnegRadonRegTerm,
-    RadonRegTerm,
-    RegTerm
-};
+use crate::types::*;
+use alg_tools::bisection_tree::P2Minimise;
+use alg_tools::bounds::MinMaxMapping;
+use alg_tools::error::DynResult;
+use alg_tools::euclidean::Euclidean;
+use alg_tools::instance::Instance;
+use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorOptions, ValueIteratorFactory};
+use alg_tools::linops::Mapping;
+use alg_tools::loc::Loc;
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::norms::Norm;
+use alg_tools::norms::L2;
+use alg_tools::sets::Cube;
 
 /// Settings for [`pointsource_fw_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
-pub struct FWConfig<F : Float> {
+pub struct FWConfig<F: Float> {
     /// Tolerance for branch-and-bound new spike location discovery
-    pub tolerance : Tolerance<F>,
+    pub tolerance: Tolerance<F>,
     /// Inner problem solution configuration. Has to have `method` set to [`InnerMethod::FB`]
     /// as the conditional gradient subproblems' optimality conditions do not in general have an
     /// invertible Newton derivative for SSN.
-    pub inner : InnerSettings<F>,
+    pub inner: InnerSettings<F>,
     /// Variant of the conditional gradient method
-    pub variant : FWVariant,
+    pub variant: FWVariant,
     /// Settings for branch and bound refinement when looking for predual maxima
-    pub refinement : RefinementSettings<F>,
+    pub refinement: RefinementSettings<F>,
     /// Spike merging heuristic
-    pub merging : SpikeMergingMethod<F>,
+    pub merging: SpikeMergingMethod<F>,
 }
 
 /// Conditional gradient method variant; see also [`FWConfig`].
@@ -101,51 +70,51 @@
     Relaxed,
 }
 
-impl<F : Float> Default for FWConfig<F> {
+impl<F: Float> Default for FWConfig<F> {
     fn default() -> Self {
         FWConfig {
-            tolerance : Default::default(),
-            refinement : Default::default(),
-            inner : Default::default(),
-            variant : FWVariant::FullyCorrective,
-            merging : SpikeMergingMethod { enabled : true, ..Default::default() },
+            tolerance: Default::default(),
+            refinement: Default::default(),
+            inner: Default::default(),
+            variant: FWVariant::FullyCorrective,
+            merging: SpikeMergingMethod { enabled: true, ..Default::default() },
         }
     }
 }
 
-pub trait FindimQuadraticModel<Domain, F> : ForwardModel<DiscreteMeasure<Domain, F>, F>
+pub trait FindimQuadraticModel<Domain, F>: ForwardModel<DiscreteMeasure<Domain, F>, F>
 where
-    F : Float + ToNalgebraRealField,
-    Domain : Clone + PartialEq,
+    F: Float + ToNalgebraRealField,
+    Domain: Clone + PartialEq,
 {
     /// Return A_*A and A_* b
     fn findim_quadratic_model(
         &self,
-        μ : &DiscreteMeasure<Domain, F>,
-        b : &Self::Observable
+        μ: &DiscreteMeasure<Domain, F>,
+        b: &Self::Observable,
     ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>);
 }
 
 /// Helper struct for pre-initialising the finite-dimensional subproblem solver.
-pub struct FindimData<F : Float> {
+pub struct FindimData<F: Float> {
     /// ‖A‖^2
-    opAnorm_squared : F,
+    opAnorm_squared: F,
     /// Bound $M_0$ from the Bredies–Pikkarainen article.
-    m0 : F
+    m0: F,
 }
 
 /// Trait for finite dimensional weight optimisation.
 pub trait WeightOptim<
-    F : Float + ToNalgebraRealField,
-    A : ForwardModel<RNDM<F, N>, F>,
-    I : AlgIteratorFactory<F>,
-    const N : usize
-> {
-
+    F: Float + ToNalgebraRealField,
+    A: ForwardModel<RNDM<N, F>, F>,
+    I: AlgIteratorFactory<F>,
+    const N: usize,
+>
+{
     /// Return a pre-initialisation struct for [`Self::optimise_weights`].
     ///
     /// The parameter `opA` is the forward operator $A$.
-    fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F>;
+    fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>>;
 
     /// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
     /// point source localisation problem.
@@ -166,72 +135,70 @@
     /// Returns the number of iterations taken by the method configured in `inner`.
     fn optimise_weights<'a>(
         &self,
-        μ : &mut RNDM<F, N>,
-        opA : &'a A,
-        b : &A::Observable,
-        findim_data : &FindimData<F>,
-        inner : &InnerSettings<F>,
-        iterator : I
+        μ: &mut RNDM<N, F>,
+        opA: &'a A,
+        b: &A::Observable,
+        findim_data: &FindimData<F>,
+        inner: &InnerSettings<F>,
+        iterator: I,
     ) -> usize;
 }
 
 /// Trait for regularisation terms supported by [`pointsource_fw_reg`].
 pub trait RegTermFW<
-    F : Float + ToNalgebraRealField,
-    A : ForwardModel<RNDM<F, N>, F>,
-    I : AlgIteratorFactory<F>,
-    const N : usize
-> : RegTerm<F, N>
-    + WeightOptim<F, A, I, N>
-    + Mapping<RNDM<F, N>, Codomain = F> {
-
+    F: Float + ToNalgebraRealField,
+    A: ForwardModel<RNDM<N, F>, F>,
+    I: AlgIteratorFactory<F>,
+    const N: usize,
+>: RegTerm<Loc<N, F>, F> + WeightOptim<F, A, I, N> + Mapping<RNDM<N, F>, Codomain = F>
+{
     /// With $g = A\_\*(Aμ-b)$, returns $(x, g(x))$ for $x$ a new point to be inserted
     /// into $μ$, as determined by the regulariser.
     ///
     /// The parameters `refinement_tolerance` and `max_steps` are passed to relevant
-    /// [`BTFN`] minimisation and maximisation routines.
+    /// [`MinMaxMapping`] minimisation and maximisation routines.
     fn find_insertion(
         &self,
-        g : &mut A::PreadjointCodomain,
-        refinement_tolerance : F,
-        max_steps : usize
-    ) -> (Loc<F, N>, F);
+        g: &mut A::PreadjointCodomain,
+        refinement_tolerance: F,
+        max_steps: usize,
+    ) -> (Loc<N, F>, F);
 
     /// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen.
     fn relaxed_insert<'a>(
         &self,
-        μ : &mut RNDM<F, N>,
-        g : &A::PreadjointCodomain,
-        opA : &'a A,
-        ξ : Loc<F, N>,
-        v_ξ : F,
-        findim_data : &FindimData<F>
+        μ: &mut RNDM<N, F>,
+        g: &A::PreadjointCodomain,
+        opA: &'a A,
+        ξ: Loc<N, F>,
+        v_ξ: F,
+        findim_data: &FindimData<F>,
     );
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, A, I, const N : usize> WeightOptim<F, A, I, N>
-for RadonRegTerm<F>
-where I : AlgIteratorFactory<F>,
-      A : FindimQuadraticModel<Loc<F, N>, F>  {
-
-    fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> {
-        FindimData{
-            opAnorm_squared : opA.opnorm_bound(Radon, L2).powi(2),
-            m0 : b.norm2_squared() / (2.0 * self.α()),
-        }
+impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N>
+    for RadonRegTerm<F>
+where
+    I: AlgIteratorFactory<F>,
+    A: FindimQuadraticModel<Loc<N, F>, F>,
+{
+    fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
+        Ok(FindimData {
+            opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2),
+            m0: b.norm2_squared() / (2.0 * self.α()),
+        })
     }
 
     fn optimise_weights<'a>(
         &self,
-        μ : &mut RNDM<F, N>,
-        opA : &'a A,
-        b : &A::Observable,
-        findim_data : &FindimData<F>,
-        inner : &InnerSettings<F>,
-        iterator : I
+        μ: &mut RNDM<N, F>,
+        opA: &'a A,
+        b: &A::Observable,
+        findim_data: &FindimData<F>,
+        inner: &InnerSettings<F>,
+        iterator: I,
     ) -> usize {
-
         // Form and solve finite-dimensional subproblem.
         let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
         let mut x = μ.masses_dvector();
@@ -245,8 +212,7 @@
         // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
         // square root is needed when we scale:
         let normest = findim_data.opAnorm_squared * F::cast_from(μ.len());
-        let iters = quadratic_unconstrained(&Ã, &g̃, self.α(), &mut x,
-                                            normest, inner, iterator);
+        let iters = quadratic_unconstrained(&Ã, &g̃, self.α(), &mut x, normest, inner, iterator);
         // Update masses of μ based on solution of finite-dimensional subproblem.
         μ.set_masses_dvector(&x);
 
@@ -255,28 +221,23 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, A, I, S, GA, BTA, const N : usize> RegTermFW<F, A, I, N>
-for RadonRegTerm<F>
+impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N> for RadonRegTerm<F>
 where
-    Cube<F, N> : P2Minimise<Loc<F, N>, F>,
-    I : AlgIteratorFactory<F>,
-    S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
-    GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
-    A : FindimQuadraticModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>,
-    BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
+    Cube<N, F>: P2Minimise<Loc<N, F>, F>,
+    I: AlgIteratorFactory<F>,
+    A: FindimQuadraticModel<Loc<N, F>, F>,
+    A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
+    for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
     // FIXME: the following *should not* be needed, they are already implied
-    RNDM<F, N> : Mapping<A::PreadjointCodomain, Codomain = F>,
-    DeltaMeasure<Loc<F, N>, F> : Mapping<A::PreadjointCodomain, Codomain = F>,
-    //A : Mapping<RNDM<F, N>, Codomain = A::Observable>,
-    //A : Mapping<DeltaMeasure<Loc<F, N>, F>, Codomain = A::Observable>,
+    RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
+    DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
 {
-
     fn find_insertion(
         &self,
-        g : &mut A::PreadjointCodomain,
-        refinement_tolerance : F,
-        max_steps : usize
-    ) -> (Loc<F, N>, F) {
+        g: &mut A::PreadjointCodomain,
+        refinement_tolerance: F,
+        max_steps: usize,
+    ) -> (Loc<N, F>, F) {
         let (ξmax, v_ξmax) = g.maximise(refinement_tolerance, max_steps);
         let (ξmin, v_ξmin) = g.minimise(refinement_tolerance, max_steps);
         if v_ξmin < 0.0 && -v_ξmin > v_ξmax {
@@ -288,25 +249,35 @@
 
     fn relaxed_insert<'a>(
         &self,
-        μ : &mut RNDM<F, N>,
-        g : &A::PreadjointCodomain,
-        opA : &'a A,
-        ξ : Loc<F, N>,
-        v_ξ : F,
-        findim_data : &FindimData<F>
+        μ: &mut RNDM<N, F>,
+        g: &A::PreadjointCodomain,
+        opA: &'a A,
+        ξ: Loc<N, F>,
+        v_ξ: F,
+        findim_data: &FindimData<F>,
     ) {
         let α = self.0;
         let m0 = findim_data.m0;
-        let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) };
-        let v = if v_ξ.abs() <= α { 0.0 } else { m0 / α * v_ξ };
-        let δ = DeltaMeasure { x : ξ, α : v };
+        let φ = |t| {
+            if t <= m0 {
+                α * t
+            } else {
+                α / (2.0 * m0) * (t * t + m0 * m0)
+            }
+        };
+        let v = if v_ξ.abs() <= α {
+            0.0
+        } else {
+            m0 / α * v_ξ
+        };
+        let δ = DeltaMeasure { x: ξ, α: v };
         let dp = μ.apply(g) - δ.apply(g);
         let d = opA.apply(&*μ) - opA.apply(δ);
         let r = d.norm2_squared();
         let s = if r == 0.0 {
             1.0
         } else {
-            1.0.min( (α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
+            1.0.min((α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
         };
         *μ *= 1.0 - s;
         *μ += δ * s;
@@ -314,28 +285,28 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, A, I, const N : usize> WeightOptim<F, A, I, N>
-for NonnegRadonRegTerm<F>
-where I : AlgIteratorFactory<F>,
-      A : FindimQuadraticModel<Loc<F, N>, F> {
-
-    fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> {
-        FindimData{
-            opAnorm_squared : opA.opnorm_bound(Radon, L2).powi(2),
-            m0 : b.norm2_squared() / (2.0 * self.α()),
-        }
+impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N>
+    for NonnegRadonRegTerm<F>
+where
+    I: AlgIteratorFactory<F>,
+    A: FindimQuadraticModel<Loc<N, F>, F>,
+{
+    fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
+        Ok(FindimData {
+            opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2),
+            m0: b.norm2_squared() / (2.0 * self.α()),
+        })
     }
 
     fn optimise_weights<'a>(
         &self,
-        μ : &mut RNDM<F, N>,
-        opA : &'a A,
-        b : &A::Observable,
-        findim_data : &FindimData<F>,
-        inner : &InnerSettings<F>,
-        iterator : I
+        μ: &mut RNDM<N, F>,
+        opA: &'a A,
+        b: &A::Observable,
+        findim_data: &FindimData<F>,
+        inner: &InnerSettings<F>,
+        iterator: I,
     ) -> usize {
-
         // Form and solve finite-dimensional subproblem.
         let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
         let mut x = μ.masses_dvector();
@@ -349,8 +320,7 @@
         // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
         // square root is needed when we scale:
         let normest = findim_data.opAnorm_squared * F::cast_from(μ.len());
-        let iters = quadratic_nonneg(&Ã, &g̃, self.α(), &mut x,
-                                     normest, inner, iterator);
+        let iters = quadratic_nonneg(&Ã, &g̃, self.α(), &mut x, normest, inner, iterator);
         // Update masses of μ based on solution of finite-dimensional subproblem.
         μ.set_masses_dvector(&x);
 
@@ -359,59 +329,65 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float + ToNalgebraRealField, A, I, S, GA, BTA, const N : usize> RegTermFW<F, A, I, N>
-for NonnegRadonRegTerm<F>
+impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N>
+    for NonnegRadonRegTerm<F>
 where
-    Cube<F, N> : P2Minimise<Loc<F, N>, F>,
-    I : AlgIteratorFactory<F>,
-    S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
-    GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
-    A : FindimQuadraticModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>,
-    BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
+    Cube<N, F>: P2Minimise<Loc<N, F>, F>,
+    I: AlgIteratorFactory<F>,
+    A: FindimQuadraticModel<Loc<N, F>, F>,
+    A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
+    for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
     // FIXME: the following *should not* be needed, they are already implied
-    RNDM<F, N> : Mapping<A::PreadjointCodomain, Codomain = F>,
-    DeltaMeasure<Loc<F, N>, F> : Mapping<A::PreadjointCodomain, Codomain = F>,
+    RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
+    DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
 {
-
     fn find_insertion(
         &self,
-        g : &mut A::PreadjointCodomain,
-        refinement_tolerance : F,
-        max_steps : usize
-    ) -> (Loc<F, N>, F) {
+        g: &mut A::PreadjointCodomain,
+        refinement_tolerance: F,
+        max_steps: usize,
+    ) -> (Loc<N, F>, F) {
         g.maximise(refinement_tolerance, max_steps)
     }
 
-
     fn relaxed_insert<'a>(
         &self,
-        μ : &mut RNDM<F, N>,
-        g : &A::PreadjointCodomain,
-        opA : &'a A,
-        ξ : Loc<F, N>,
-        v_ξ : F,
-        findim_data : &FindimData<F>
+        μ: &mut RNDM<N, F>,
+        g: &A::PreadjointCodomain,
+        opA: &'a A,
+        ξ: Loc<N, F>,
+        v_ξ: F,
+        findim_data: &FindimData<F>,
     ) {
         // This is just a verbatim copy of RadonRegTerm::relaxed_insert.
         let α = self.0;
         let m0 = findim_data.m0;
-        let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) };
-        let v = if v_ξ.abs() <= α { 0.0 } else { m0 / α * v_ξ };
-        let δ = DeltaMeasure { x : ξ, α : v };
+        let φ = |t| {
+            if t <= m0 {
+                α * t
+            } else {
+                α / (2.0 * m0) * (t * t + m0 * m0)
+            }
+        };
+        let v = if v_ξ.abs() <= α {
+            0.0
+        } else {
+            m0 / α * v_ξ
+        };
+        let δ = DeltaMeasure { x: ξ, α: v };
         let dp = μ.apply(g) - δ.apply(g);
         let d = opA.apply(&*μ) - opA.apply(&δ);
         let r = d.norm2_squared();
         let s = if r == 0.0 {
             1.0
         } else {
-            1.0.min( (α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
+            1.0.min((α * μ.norm(Radon) - φ(v.abs()) - dp) / r)
         };
         *μ *= 1.0 - s;
         *μ += δ * s;
     }
 }
 
-
 /// Solve point source localisation problem using a conditional gradient method
 /// for the 2-norm-squared data fidelity, i.e., the problem
 /// <div>$$
@@ -425,49 +401,48 @@
 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
 /// save intermediate iteration states as images.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fw_reg<F, I, A, GA, BTA, S, Reg, const N : usize>(
-    opA : &A,
-    b : &A::Observable,
-    reg : Reg,
-    //domain : Cube<F, N>,
-    config : &FWConfig<F>,
-    iterator : I,
-    mut plotter : SeqPlotter<F, N>,
-) -> RNDM<F, N>
-where F : Float + ToNalgebraRealField,
-      I : AlgIteratorFactory<IterInfo<F, N>>,
-      for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
-      GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
-      A : ForwardModel<RNDM<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>,
-      BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
-      S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
-      BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
-      Cube<F, N>: P2Minimise<Loc<F, N>, F>,
-      PlotLookup : Plotting<N>,
-      RNDM<F, N> : SpikeMerging<F>,
-      Reg : RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N> {
+pub fn pointsource_fw_reg<'a, F, I, A, Reg, Plot, const N: usize>(
+    f: &'a QuadraticDataTerm<F, RNDM<N, F>, A>,
+    reg: &Reg,
+    //domain : Cube<N, F>,
+    config: &FWConfig<F>,
+    iterator: I,
+    mut plotter: Plot,
+    μ0 : Option<RNDM<N, F>>,
+) -> DynResult<RNDM<N, F>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    A: ForwardModel<RNDM<N, F>, F>,
+    A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
+    &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
+    Cube<N, F>: P2Minimise<Loc<N, F>, F>,
+    RNDM<N, F>: SpikeMerging<F>,
+    Reg: RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N>,
+    Plot: Plotter<A::PreadjointCodomain, A::PreadjointCodomain, RNDM<N, F>>,
+{
+    let opA = f.operator();
+    let b = f.data();
 
     // Set up parameters
     // We multiply tolerance by α for all algoritms.
     let tolerance = config.tolerance * reg.tolerance_scaling();
     let mut ε = tolerance.initial();
-    let findim_data = reg.prepare_optimise_weights(opA, b);
+    let findim_data = reg.prepare_optimise_weights(opA, b)?;
 
     // Initialise operators
     let preadjA = opA.preadjoint();
 
     // Initialise iterates
-    let mut μ = DiscreteMeasure::new();
-    let mut residual = -b;
+    let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
+    let mut residual = f.residual(&μ);
 
     // Statistics
-    let full_stats = |residual : &A::Observable,
-                      ν : &RNDM<F, N>,
-                      ε, stats| IterInfo {
-        value : residual.norm2_squared_div2() + reg.apply(ν),
-        n_spikes : ν.len(),
+    let full_stats = |residual: &A::Observable, ν: &RNDM<N, F>, ε, stats| IterInfo {
+        value: residual.norm2_squared_div2() + reg.apply(ν),
+        n_spikes: ν.len(),
         ε,
-        .. stats
+        ..stats
     };
     let mut stats = IterInfo::new();
 
@@ -480,32 +455,34 @@
         let mut g = preadjA.apply(residual * (-1.0));
 
         // Find absolute value maximising point
-        let (ξ, v_ξ) = reg.find_insertion(&mut g, refinement_tolerance,
-                                          config.refinement.max_steps);
+        let (ξ, v_ξ) =
+            reg.find_insertion(&mut g, refinement_tolerance, config.refinement.max_steps);
 
         let inner_it = match config.variant {
             FWVariant::FullyCorrective => {
                 // No point in optimising the weight here: the finite-dimensional algorithm is fast.
-                μ += DeltaMeasure { x : ξ, α : 0.0 };
+                μ += DeltaMeasure { x: ξ, α: 0.0 };
                 stats.inserted += 1;
                 config.inner.iterator_options.stop_target(inner_tolerance)
-            },
+            }
             FWVariant::Relaxed => {
                 // Perform a relaxed initialisation of μ
                 reg.relaxed_insert(&mut μ, &g, opA, ξ, v_ξ, &findim_data);
                 stats.inserted += 1;
                 // The stop_target is only needed for the type system.
-                AlgIteratorOptions{ max_iter : 1, .. config.inner.iterator_options}.stop_target(0.0)
+                AlgIteratorOptions { max_iter: 1, ..config.inner.iterator_options }.stop_target(0.0)
             }
         };
 
-        stats.inner_iters += reg.optimise_weights(&mut μ, opA, b, &findim_data,
-                                                  &config.inner, inner_it);
-   
+        stats.inner_iters +=
+            reg.optimise_weights(&mut μ, opA, b, &findim_data, &config.inner, inner_it);
+
         // Merge spikes and update residual for next step and `if_verbose` below.
-        let (r, count) = μ.merge_spikes_fitness(config.merging,
-                                                |μ̃| opA.apply(μ̃) - b,
-                                                A::Observable::norm2_squared);
+        let (r, count) = μ.merge_spikes_fitness(
+            config.merging,
+            |μ̃| f.residual(μ̃),
+            A::Observable::norm2_squared,
+        );
         residual = r;
         stats.merged += count;
 
@@ -520,8 +497,13 @@
 
         // Give statistics if needed
         state.if_verbose(|| {
-            plotter.plot_spikes(iter, Some(&g), Option::<&S>::None, &μ);
-            full_stats(&residual, &μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
+            plotter.plot_spikes(iter, Some(&g), Option::<&A::PreadjointCodomain>::None, &μ);
+            full_stats(
+                &residual,
+                &μ,
+                ε,
+                std::mem::replace(&mut stats, IterInfo::new()),
+            )
         });
 
         // Update tolerance
@@ -529,5 +511,5 @@
     }
 
     // Return final iterate
-    μ
+    Ok(μ)
 }