Mercurial > repos > pointsource_algs / changeset
changeset
Implement non-negativity constraints for the conditional gradient methods
Tue, 21 Mar 2023 20:31:01 +0200
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Tue, 21 Mar 2023 20:31:01 +0200
- changeset 25
- 79943be70720
- parent 24
- d29d1fcf5423
- child 26
- acf57c458740
Implement non-negativity constraints for the conditional gradient methods
| src/experiments.rs | file | annotate | diff | comparison | revisions | |
| src/frank_wolfe.rs | file | annotate | diff | comparison | revisions | |
| src/run.rs | file | annotate | diff | comparison | revisions |
--- 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(),
--- 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) +}
--- 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();
