# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1679423461 -7200
# Node ID 79943be707208eba8e9cc4f301ae7e9755fb69fa
# Parent  d29d1fcf542310e46b5d479b4dabab8b133f2c5a
Implement non-negativity constraints for the conditional gradient methods

diff -r d29d1fcf5423 -r 79943be70720 src/experiments.rs
--- a/src/experiments.rs	Sun Dec 11 23:25:53 2022 +0200
+++ b/src/experiments.rs	Tue Mar 21 20:31:01 2023 +0200
@@ -140,7 +140,7 @@
                 Box::new(Named { name, data : ExperimentV2 {
                     domain : [[0.0, 1.0]].into(),
                     sensor_count : [N_SENSORS_1D],
-                    regularisation : Regularisation::Radon(cli.alpha.unwrap_or(0.09)),
+                    regularisation : Regularisation::NonnegRadon(cli.alpha.unwrap_or(0.09)),
                     noise_distr : SerializableNormal::new(0.0, cli.variance.unwrap_or(0.2))?,
                     dataterm : DataTerm::L2Squared,
                     μ_hat : MU_TRUE_1D_BASIC.into(),
@@ -157,7 +157,7 @@
                 Box::new(Named { name, data : ExperimentV2 {
                     domain : [[0.0, 1.0]].into(),
                     sensor_count : [N_SENSORS_1D],
-                    regularisation : Regularisation::Radon(cli.alpha.unwrap_or(0.06)),
+                    regularisation : Regularisation::NonnegRadon(cli.alpha.unwrap_or(0.06)),
                     noise_distr : SerializableNormal::new(0.0, cli.variance.unwrap_or(0.2))?,
                     dataterm : DataTerm::L2Squared,
                     μ_hat : MU_TRUE_1D_BASIC.into(),
@@ -175,7 +175,7 @@
                 Box::new(Named { name, data : ExperimentV2 {
                     domain : [[0.0, 1.0]; 2].into(),
                     sensor_count : [N_SENSORS_2D; 2],
-                    regularisation : Regularisation::Radon(cli.alpha.unwrap_or(0.19)),
+                    regularisation : Regularisation::NonnegRadon(cli.alpha.unwrap_or(0.19)),
                     noise_distr : SerializableNormal::new(0.0, cli.variance.unwrap_or(0.25))?,
                     dataterm : DataTerm::L2Squared,
                     μ_hat : MU_TRUE_2D_BASIC.into(),
@@ -194,7 +194,7 @@
                 Box::new(Named { name, data : ExperimentV2 {
                     domain : [[0.0, 1.0]; 2].into(),
                     sensor_count : [N_SENSORS_2D; 2],
-                    regularisation : Regularisation::Radon(cli.alpha.unwrap_or(0.12)),
+                    regularisation : Regularisation::NonnegRadon(cli.alpha.unwrap_or(0.12)),
                     noise_distr : SerializableNormal::new(0.0, cli.variance.unwrap_or(0.15))?, //0.25
                     dataterm : DataTerm::L2Squared,
                     μ_hat : MU_TRUE_2D_BASIC.into(),
diff -r d29d1fcf5423 -r 79943be70720 src/frank_wolfe.rs
--- a/src/frank_wolfe.rs	Sun Dec 11 23:25:53 2022 +0200
+++ b/src/frank_wolfe.rs	Tue Mar 21 20:31:01 2023 +0200
@@ -21,6 +21,7 @@
     AlgIteratorFactory,
     AlgIteratorState,
     AlgIteratorOptions,
+    ValueIteratorFactory,
 };
 use alg_tools::euclidean::Euclidean;
 use alg_tools::norms::Norm;
@@ -54,6 +55,7 @@
 #[allow(unused_imports)] // Used in documentation
 use crate::subproblem::{
     unconstrained::quadratic_unconstrained,
+    nonneg::quadratic_nonneg,
     InnerSettings,
     InnerMethod,
 };
@@ -63,6 +65,11 @@
     Plotting,
     PlotLookup
 };
+use crate::regularisation::{
+    NonnegRadonRegTerm,
+    RadonRegTerm,
+};
+use crate::fb::RegTerm;
 
 /// Settings for [`pointsource_fw`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
@@ -109,85 +116,295 @@
 ///
 /// The pre-initialisation is done by [`prepare_optimise_weights`].
 pub struct FindimData<F : Float> {
-    opAnorm_squared : F
+    /// ‖A‖^2
+    opAnorm_squared : F,
+    /// Bound $M_0$ from the Bredies–Pikkarainen article.
+    m0 : F
+}
+
+/// Trait for finite dimensional weight optimisation.
+pub trait WeightOptim<
+    F : Float + ToNalgebraRealField,
+    A : ForwardModel<Loc<F, N>, 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>;
+
+    /// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
+    /// point source localisation problem.
+    ///
+    /// That is, we minimise
+    /// <div>$$
+    ///     μ ↦ \frac{1}{2}\|Aμ-b\|_w^2 + G(μ)
+    /// $$</div>
+    /// only with respect to the weights of $μ$. Here $G$ is a regulariser modelled by `Self`.
+    ///
+    /// The parameter `μ` is the discrete measure whose weights are to be optimised.
+    /// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
+    /// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
+    /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
+    /// save intermediate iteration states as images. The parameter `findim_data` should be
+    /// prepared using [`Self::prepare_optimise_weights`]:
+    ///
+    /// Returns the number of iterations taken by the method configured in `inner`.
+    fn optimise_weights<'a>(
+        &self,
+        μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+        opA : &'a A,
+        b : &A::Observable,
+        findim_data : &FindimData<F>,
+        inner : &InnerSettings<F>,
+        iterator : I
+    ) -> usize;
 }
 
-/// Return a pre-initialisation struct for [`prepare_optimise_weights`].
-///
-/// The parameter `opA` is the forward operator $A$.
-pub fn prepare_optimise_weights<F, A, const N : usize>(opA : &A) -> FindimData<F>
-where F : Float + ToNalgebraRealField,
+/// Trait for regularisation terms supported by [`pointsource_fw_reg`].
+pub trait RegTermFW<
+    F : Float + ToNalgebraRealField,
+    A : ForwardModel<Loc<F, N>, F>,
+    I : AlgIteratorFactory<F>,
+    const N : usize
+> : RegTerm<F, N>
+    + WeightOptim<F, A, I, N>
+    + for<'a> Apply<&'a DiscreteMeasure<Loc<F, N>, F>, Output = 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.
+    fn find_insertion(
+        &self,
+        g : &mut A::PreadjointCodomain,
+        refinement_tolerance : F,
+        max_steps : usize
+    ) -> (Loc<F, N>, F);
+
+    /// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen.
+    fn relaxed_insert<'a>(
+        &self,
+        μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+        g : &A::PreadjointCodomain,
+        opA : &'a A,
+        ξ : Loc<F, N>,
+        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 : ForwardModel<Loc<F, N>, F> {
-    FindimData{
-        opAnorm_squared : opA.opnorm_bound().powi(2)
+
+    fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> {
+        FindimData{
+            opAnorm_squared : opA.opnorm_bound().powi(2),
+            m0 : b.norm2_squared() / (2.0 * self.α()),
+        }
+    }
+
+    fn optimise_weights<'a>(
+        &self,
+        μ : &mut DiscreteMeasure<Loc<F, 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();
+
+        // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
+        // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
+        // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
+        // 2-norm, we estimate
+        //      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
+        //                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
+        // 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 inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
+        let iters = quadratic_unconstrained(inner.method, &Ã, &g̃, self.α(),
+                                            &mut x, inner_τ, iterator);
+        // Update masses of μ based on solution of finite-dimensional subproblem.
+        μ.set_masses_dvector(&x);
+
+        iters
     }
 }
 
-/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
-/// point source localisation problem.
-///
-/// That is, we minimise
-/// <div>$$
-///     μ ↦ \frac{1}{2}\|Aμ-b\|_w^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
-/// $$</div>
-/// only with respect to the weights of $μ$.
-///
-/// The parameter `μ` is the discrete measure whose weights are to be optimised.
-/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
-/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
-/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
-/// save intermediate iteration states as images. The parameter `findim_data` should be
-/// prepared using [`prepare_optimise_weights`]:
-///
-/// Returns the number of iterations taken by the method configured in `inner`.
-pub fn optimise_weights<'a, F, A, I, const N : usize>(
-    μ : &mut DiscreteMeasure<Loc<F, N>, F>,
-    opA : &'a A,
-    b : &A::Observable,
-    α : F,
-    findim_data : &FindimData<F>,
-    inner : &InnerSettings<F>,
-    iterator : I
-) -> usize
-where F : Float + ToNalgebraRealField,
+#[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>
+where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
       I : AlgIteratorFactory<F>,
-      A : ForwardModel<Loc<F, N>, F>
-{
-    // Form and solve finite-dimensional subproblem.
-    let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
-    let mut x = μ.masses_dvector();
+      S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
+      GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
+      A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>,
+      BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>> {
+
+    fn find_insertion(
+        &self,
+        g : &mut A::PreadjointCodomain,
+        refinement_tolerance : F,
+        max_steps : usize
+    ) -> (Loc<F, N>, 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 {
+            (ξmin, v_ξmin)
+        } else {
+            (ξmax, v_ξmax)
+        }
+    }
+
+    fn relaxed_insert<'a>(
+        &self,
+        μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+        g : &A::PreadjointCodomain,
+        opA : &'a A,
+        ξ : Loc<F, N>,
+        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 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 - s;
+        *μ += δ * s;
+    }
+}
+
+#[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 : ForwardModel<Loc<F, N>, F> {
+
+    fn prepare_optimise_weights(&self, opA : &A, b : &A::Observable) -> FindimData<F> {
+        FindimData{
+            opAnorm_squared : opA.opnorm_bound().powi(2),
+            m0 : b.norm2_squared() / (2.0 * self.α()),
+        }
+    }
+
+    fn optimise_weights<'a>(
+        &self,
+        μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+        opA : &'a A,
+        b : &A::Observable,
+        findim_data : &FindimData<F>,
+        inner : &InnerSettings<F>,
+        iterator : I
+    ) -> usize {
 
-    // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
-    // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
-    // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
-    // 2-norm, we estimate
-    //      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
-    //                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
-    // 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 inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
-    let iters = quadratic_unconstrained(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator);
-    // Update masses of μ based on solution of finite-dimensional subproblem.
-    μ.set_masses_dvector(&x);
+        // Form and solve finite-dimensional subproblem.
+        let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
+        let mut x = μ.masses_dvector();
+
+        // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
+        // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
+        // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
+        // 2-norm, we estimate
+        //      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
+        //                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
+        // 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 inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
+        let iters = quadratic_nonneg(inner.method, &Ã, &g̃, self.α(),
+                                     &mut x, inner_τ, iterator);
+        // Update masses of μ based on solution of finite-dimensional subproblem.
+        μ.set_masses_dvector(&x);
+
+        iters
+    }
+}
+
+#[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>
+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 : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>,
+      BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>> {
 
-    iters
+    fn find_insertion(
+        &self,
+        g : &mut A::PreadjointCodomain,
+        refinement_tolerance : F,
+        max_steps : usize
+    ) -> (Loc<F, N>, F) {
+        g.maximise(refinement_tolerance, max_steps)
+    }
+
+
+    fn relaxed_insert<'a>(
+        &self,
+        μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+        g : &A::PreadjointCodomain,
+        opA : &'a A,
+        ξ : Loc<F, N>,
+        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 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 - s;
+        *μ += δ * s;
+    }
 }
 
+
 /// Solve point source localisation problem using a conditional gradient method
 /// for the 2-norm-squared data fidelity, i.e., the problem
 /// <div>$$
-///     \min_μ \frac{1}{2}\|Aμ-b\|_w^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ).
-/// $$</div>
+///     \min_μ \frac{1}{2}\|Aμ-b\|_w^2 + G(μ),
+/// $$
+/// where $G$ is the regularisation term modelled by `reg`.
+/// </div>
 ///
-/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
+/// The `opA` parameter is the forward operator $A$, while `b`$ is as in the
 /// objective above. The method parameter are set in `config` (see [`FWConfig`]), while
 /// `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<'a, F, I, A, GA, BTA, S, const N : usize>(
+pub fn pointsource_fw_reg<'a, F, I, A, GA, BTA, S, Reg, const N : usize>(
     opA : &'a A,
     b : &A::Observable,
-    α : F,
+    reg : Reg,
     //domain : Cube<F, N>,
     config : &FWConfig<F>,
     iterator : I,
@@ -205,15 +422,14 @@
       BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
       Cube<F, N>: P2Minimise<Loc<F, N>, F>,
       PlotLookup : Plotting<N>,
-      DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
+      DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
+      Reg : RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N> {
 
     // Set up parameters
     // We multiply tolerance by α for all algoritms.
-    let tolerance = config.tolerance * α;
+    let tolerance = config.tolerance * reg.tolerance_scaling();
     let mut ε = tolerance.initial();
-    let findim_data = prepare_optimise_weights(opA);
-    let m0 = b.norm2_squared() / (2.0 * α);
-    let φ = |t| if t <= m0 { α * t } else { α / (2.0 * m0) * (t*t + m0 * m0) };
+    let findim_data = reg.prepare_optimise_weights(opA, b);
 
     // Initialise operators
     let preadjA = opA.preadjoint();
@@ -244,15 +460,8 @@
         let mut g = -preadjA.apply(r);
 
         // Find absolute value maximising point
-        let (ξmax, v_ξmax) = g.maximise(refinement_tolerance,
-                                        config.refinement.max_steps);
-        let (ξmin, v_ξmin) = g.minimise(refinement_tolerance,
-                                        config.refinement.max_steps);
-        let (ξ, v_ξ) = if v_ξmin < 0.0 && -v_ξmin > v_ξmax {
-            (ξmin, v_ξmin)
-        } else {
-            (ξmax, v_ξmax)
-        };
+        let (ξ, v_ξ) = reg.find_insertion(&mut g, refinement_tolerance,
+                                          config.refinement.max_steps);
 
         let inner_it = match config.variant {
             FWVariant::FullyCorrective => {
@@ -262,24 +471,13 @@
             },
             FWVariant::Relaxed => {
                 // Perform a relaxed initialisation of μ
-                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 - s;
-                μ += δ * s;
+                reg.relaxed_insert(&mut μ, &g, opA, ξ, v_ξ, &findim_data);
                 // The stop_target is only needed for the type system.
                 AlgIteratorOptions{ max_iter : 1, .. config.inner.iterator_options}.stop_target(0.0)
             }
         };
 
-        inner_iters += optimise_weights(&mut μ, opA, b, α, &findim_data, &config.inner, inner_it);
+        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 n_before_merge = μ.len();
@@ -306,7 +504,7 @@
                 None, &μ
             );
             let res = IterInfo {
-                value : residual.norm2_squared_div2() + α * μ.norm(Radon),
+                value : residual.norm2_squared_div2() + reg.apply(&μ),
                 n_spikes : μ.len(),
                 inner_iters,
                 this_iters,
@@ -327,6 +525,50 @@
     μ
 }
 
-
+//
+// Deprecated interface
+//
 
+#[deprecated(note = "Use `pointsource_fw_reg`")]
+pub fn pointsource_fw<'a, F, I, A, GA, BTA, S, const N : usize>(
+    opA : &'a A,
+    b : &A::Observable,
+    α : F,
+    //domain : Cube<F, N>,
+    config : &FWConfig<F>,
+    iterator : I,
+    plotter : SeqPlotter<F, N>,
+) -> DiscreteMeasure<Loc<F, N>, F>
+where F : Float + ToNalgebraRealField,
+      I : AlgIteratorFactory<IterInfo<F, N>>,
+      for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
+                                  //+ std::ops::Mul<F, Output=A::Observable>,  <-- FIXME: compiler overflow
+      A::Observable : std::ops::MulAssign<F>,
+      GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
+      A : ForwardModel<Loc<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>,
+      DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F> {
 
+      pointsource_fw_reg(opA, b, NonnegRadonRegTerm(α), config, iterator, plotter)
+}
+
+#[deprecated(note = "Use `WeightOptim::optimise_weights`")]
+pub fn optimise_weights<'a, F, A, I, const N : usize>(
+    μ : &mut DiscreteMeasure<Loc<F, N>, F>,
+    opA : &'a A,
+    b : &A::Observable,
+    α : F,
+    findim_data : &FindimData<F>,
+    inner : &InnerSettings<F>,
+    iterator : I
+) -> usize
+where F : Float + ToNalgebraRealField,
+      I : AlgIteratorFactory<F>,
+      A : ForwardModel<Loc<F, N>, F>
+{
+     NonnegRadonRegTerm(α).optimise_weights(μ, opA, b, findim_data, inner, iterator)
+}
diff -r d29d1fcf5423 -r 79943be70720 src/run.rs
--- a/src/run.rs	Sun Dec 11 23:25:53 2022 +0200
+++ b/src/run.rs	Tue Mar 21 20:31:01 2023 +0200
@@ -57,9 +57,8 @@
 use crate::frank_wolfe::{
     FWConfig,
     FWVariant,
-    pointsource_fw,
-    prepare_optimise_weights,
-    optimise_weights,
+    pointsource_fw_reg,
+    WeightOptim,
 };
 use crate::subproblem::InnerSettings;
 use crate::seminorms::*;
@@ -436,7 +435,11 @@
             };
             // Create Logger and IteratorFactory
             let mut logger = Logger::new();
-            let findim_data = prepare_optimise_weights(&opA);
+            let reg : Box<dyn WeightOptim<_, _, _, N>> = match regularisation {
+                Regularisation::Radon(α) => Box::new(RadonRegTerm(α)),
+                Regularisation::NonnegRadon(α) => Box::new(NonnegRadonRegTerm(α)),
+            };
+            let findim_data = reg.prepare_optimise_weights(&opA, &b);
             let inner_config : InnerSettings<F> = Default::default();
             let inner_it = inner_config.iterator_options;
             let logmap = |iter, Timed { cpu_time, data }| {
@@ -450,12 +453,12 @@
                     this_iters,
                     ..
                 } = data;
-                let post_value = match (postprocessing, dataterm, regularisation) {
-                    (Some(mut μ), DataTerm::L2Squared, Regularisation::Radon(α)) => {
+                let post_value = match (postprocessing, dataterm) {
+                    (Some(mut μ), DataTerm::L2Squared) => {
                         // Comparison postprocessing is only implemented for the case handled
                         // by the FW variants.
-                        optimise_weights(
-                            &mut μ, &opA, &b, α, &findim_data, &inner_config,
+                        reg.optimise_weights(
+                            &mut μ, &opA, &b, &findim_data, &inner_config,
                             inner_it
                         );
                         dataterm.value_at_residual(opA.apply(&μ) - &b)
@@ -544,7 +547,13 @@
                     match (regularisation, dataterm) {
                         (Regularisation::Radon(α), DataTerm::L2Squared) => {
                             running();
-                            pointsource_fw(&opA, &b, α, algconfig, iterator, plotter)
+                            pointsource_fw_reg(&opA, &b, RadonRegTerm(α),
+                                               algconfig, iterator, plotter)
+                        },
+                        (Regularisation::NonnegRadon(α), DataTerm::L2Squared) => {
+                            running();
+                            pointsource_fw_reg(&opA, &b, NonnegRadonRegTerm(α),
+                                               algconfig, iterator, plotter)
                         },
                         _ => {
                             not_implemented();