pointsource_algs: comparison src/fb.rs

-:51dbb2cc25d9
+:8d1c959432c8
 We solve this with either SSN or FB as determined by
 [`crate::subproblem::InnerSettings`] in [`InsertionConfig::inner`].
 */
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::measures::merging::SpikeMerging;
 use crate::measures::{DiscreteMeasure, RNDM};
 use crate::plot::Plotter;
 pub use crate::prox_penalty::{InsertionConfig, ProxPenalty, StepLengthBound};
 use crate::regularisation::RegTerm;
 use alg_tools::instance::Instance;
 use alg_tools::iterate::AlgIteratorFactory;
 use alg_tools::mapping::DifferentiableMapping;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use colored::Colorize;
+use measures::DeltaMeasure;
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
 /// Settings for [`pointsource_fb_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 ///
 /// For details on the mathematical formulation, see the [module level](self) documentation.
 ///
 /// Returns the final iterate.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fb_reg<F, I, Dat, Reg, P, Plot, const N: usize>(
+pub fn pointsource_fb_dom<F, I, Dat, Reg, P, Plot, const N: usize>(
 f: &Dat,
 reg: &Reg,
 prox_penalty: &P,
 fbconfig: &FBConfig<F>,
 iterator: I,
 mut plotter: Plot,
 μ0: Option<RNDM<N, F>>,
+domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<RNDM<N, F>>
 where
 F: Float + ToNalgebraRealField,
 I: AlgIteratorFactory<IterInfo<F>>,
 RNDM<N, F>: SpikeMerging<F>,
 };
 let mut stats = IterInfo::new();
 // Run the algorithm
 for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
-let maybe_μ_base = config.if_merge_now(&state, |_| μ.clone());
+let maybe_μ_base = config.if_merge_now_check(&state, || μ.clone());
 let μ_base_len = μ.len();
 // Calculate smooth part of surrogate model.
 let mut τv = f.differential(&μ) * τ;
 prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
 stats.inserted += μ.len() - μ_base_len;
 // Prune and possibly merge spikes
-if let Some(μ_base) = maybe_μ_base {
+let did_merge = if let Some(μ_base) = maybe_μ_base {
+// Merge spikes.
+// This crucially expects the merge routine to be stable with respect to spike locations,
+// and not to performing any pruning. That is be to done below simultaneously for γ.
 config.if_merge_now(&state, |ins| {
 stats.merged += prox_penalty.merge_spikes(
 &mut μ,
 &mut τv,
 &μ_base,
 τ,
 ε,
 ins,
 &reg,
 Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃)),
-)
+);
-});
+stats.pruned += prune_with_stats(&mut μ);
-}
+if config.random_insertions {
-stats.pruned += prune_with_stats(&mut μ);
+if let Some(x) = domain.random_point() {
+μ += DeltaMeasure { x, α: 0.0 };
-// Do extra weight optimisation step heuristic
+stats.inserted += 1;
-for _ in 0..config.extra_weight_optimisation_steps {
+}
-τv = f.differential(&μ) * τ;
+}
-prox_penalty.reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
-}
+// Do extra weight optimisation step heuristic
+for _ in 1..config.extra_weight_optimisation_steps {
-if config.extra_weight_optimisation_steps > 0 {
+τv = f.differential(&μ) * τ;
+prox_penalty
+.reweigh(&mut μ, &mut τv, τ, ε, &config, &reg, &state, &mut stats)?;
+}
+Ok(())
+})?
+} else {
+None
+};
+if did_merge.is_none() {
 stats.pruned += prune_with_stats(&mut μ);
 }
 let iter = state.iteration();
 stats.this_iters += 1;
 ε = tolerance.update(ε, iter);
 }
 //postprocess(μ_prev, config, f)
 postprocess(μ, config, |μ̃| f.apply(μ̃) + reg.apply(μ̃))
+}
+#[deprecated]
+pub fn pointsource_fb_reg<F, I, Dat, Reg, P, Plot, const N: usize>(
+f: &Dat,
+reg: &Reg,
+prox_penalty: &P,
+fbconfig: &FBConfig<F>,
+iterator: I,
+plotter: Plot,
+μ0: Option<RNDM<N, F>>,
+) -> DynResult<RNDM<N, F>>
+where
+F: Float + ToNalgebraRealField,
+I: AlgIteratorFactory<IterInfo<F>>,
+RNDM<N, F>: SpikeMerging<F>,
+Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F>,
+Dat::DerivativeDomain: ClosedMul<F>,
+Reg: RegTerm<Loc<N, F>, F>,
+P: ProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
+{
+pointsource_fb_dom(
+f,
+reg,
+prox_penalty,
+fbconfig,
+iterator,
+plotter,
+μ0,
+&mut FullSpace,
+)
 }
 /// Iteratively solve the pointsource localisation problem using inertial forward-backward splitting.
 ///
 /// The settings in `config` have their [respective documentation](FBConfig). `opA` is the
 stats.inserted += μ.len() - μ_base_len;
 // (Do not) merge spikes.
 if !warned_merging {
-config.if_merge_now(&state, |_| {
+config.if_merge_now_check(&state, || {
 let err = format!("Merging not supported for μFISTA");
 println!("{}", err.red());
 warned_merging = true;
 });
 }

Mercurial > repos > pointsource_algs / file comparison

comparison: src/fb.rs

src/fb.rs