pointsource_algs: comparison src/fb.rs

-:aec67cdd6b14
+:efa60bc4f743
 calculate_residual,
 L2Squared,
 DataTerm,
 };
-/// Method for constructing $μ$ on each iteration
-#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
-#[allow(dead_code)]
-pub enum InsertionStyle {
-/// Resuse previous $μ$ from previous iteration, optimising weights
-/// before inserting new spikes.
-Reuse,
-/// Start each iteration with $μ=0$.
-Zero,
-}
 /// Settings for [`pointsource_fb_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
 pub struct FBConfig<F : Float> {
 /// Step length scaling
 pub τ0 : F,
 /// Generic parameters
-pub insertion : FBGenericConfig<F>,
+pub generic : FBGenericConfig<F>,
 }
 /// Settings for the solution of the stepwise optimality condition in algorithms based on
 /// [`generic_pointsource_fb_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
 pub struct FBGenericConfig<F : Float> {
-/// Method for constructing $μ$ on each iteration; see [`InsertionStyle`].
-pub insertion_style : InsertionStyle,
 /// Tolerance for point insertion.
 pub tolerance : Tolerance<F>,
 /// Stop looking for predual maximum (where to isert a new point) below
 /// `tolerance` multiplied by this factor.
+///
+/// Not used by [`super::radon_fb`].
 pub insertion_cutoff_factor : F,
 /// Settings for branch and bound refinement when looking for predual maxima
 pub refinement : RefinementSettings<F>,
 /// Maximum insertions within each outer iteration
+///
+/// Not used by [`super::radon_fb`].
 pub max_insertions : usize,
 /// Pair `(n, m)` for maximum insertions `m` on first `n` iterations.
+///
+/// Not used by [`super::radon_fb`].
 pub bootstrap_insertions : Option<(usize, usize)>,
 /// Inner method settings
 pub inner : InnerSettings<F>,
 /// Spike merging method
 pub merging : SpikeMergingMethod<F>,
 /// Tolerance multiplier for merges
 pub merge_tolerance_mult : F,
 /// Spike merging method after the last step
 pub final_merging : SpikeMergingMethod<F>,
 /// Iterations between merging heuristic tries
 pub merge_every : usize,
 /// Save $μ$ for postprocessing optimisation
 pub postprocessing : bool
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float> Default for FBConfig<F> {
 fn default() -> Self {
 FBConfig {
 τ0 : 0.99,
-insertion : Default::default()
+generic : Default::default(),
 }
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float> Default for FBGenericConfig<F> {
 fn default() -> Self {
 FBGenericConfig {
-insertion_style : InsertionStyle::Reuse,
 tolerance : Default::default(),
 insertion_cutoff_factor : 1.0,
 refinement : Default::default(),
 max_insertions : 100,
 //bootstrap_insertions : None,
 bootstrap_insertions : Some((10, 1)),
 inner : InnerSettings {
-method : InnerMethod::SSN,
+method : InnerMethod::Default,
 .. Default::default()
 },
 merging : SpikeMergingMethod::None,
 //merging : Default::default(),
 final_merging : Default::default(),
 postprocessing : false,
 }
 }
 }
-#[replace_float_literals(F::cast_from(literal))]
+/// TODO: document.
-pub(crate) fn μ_diff<F : Float, const N : usize>(
+/// `μ_base + ν_delta` is the base point, where `μ` and `μ_base` are assumed to have the same spike
-μ_new : &DiscreteMeasure<Loc<F, N>, F>,
+/// locations, while `ν_delta` may have different locations.
-μ_base : &DiscreteMeasure<Loc<F, N>, F>,
-ν_delta : Option<&DiscreteMeasure<Loc<F, N>, F>>,
-config : &FBGenericConfig<F>
-) -> DiscreteMeasure<Loc<F, N>, F> {
-let mut ν : DiscreteMeasure<Loc<F, N>, F> = match config.insertion_style {
-InsertionStyle::Reuse => {
-μ_new.iter_spikes()
-.zip(μ_base.iter_masses().chain(std::iter::repeat(0.0)))
-.map(|(δ, α_base)| (δ.x, α_base - δ.α))
-.collect()
-},
-InsertionStyle::Zero => {
-μ_new.iter_spikes()
-.map(|δ| -δ)
-.chain(μ_base.iter_spikes().copied())
-.collect()
-}
-};
-ν.prune(); // Potential small performance improvement
-// Add ν_delta if given
-match ν_delta {
-None => ν,
-Some(ν_d) => ν + ν_d,
-}
-}
 #[replace_float_literals(F::cast_from(literal))]
 pub(crate) fn insert_and_reweigh<
 'a, F, GA, 𝒟, BTA, G𝒟, S, K, Reg, State, const N : usize
 >(
 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
 Reg : RegTerm<F, N>,
 State : AlgIteratorState {
 // Maximum insertion count and measure difference calculation depend on insertion style.
-let (m, warn_insertions) = match (state.iteration(), config.bootstrap_insertions) {
+let (max_insertions, warn_insertions) = match (state.iteration(), config.bootstrap_insertions) {
 (i, Some((l, k))) if i <= l => (k, false),
 _ => (config.max_insertions, !state.is_quiet()),
 };
-let max_insertions = match config.insertion_style {
-InsertionStyle::Zero => {
+// TODO: should avoid a copy of μ_base here.
-todo!("InsertionStyle::Zero does not currently work with FISTA, so diabled.");
-// let n = μ.len();
-// μ = DiscreteMeasure::new();
-// n + m
-},
-InsertionStyle::Reuse => m,
-};
-// TODO: should avoid a second copy of μ here; μ_base already stores a copy.
 let ω0 = op𝒟.apply(match ν_delta {
-None => μ.clone(),
+None => μ_base.clone(),
-Some(ν_d) => &*μ + ν_d,
+Some(ν_d) => &*μ_base + ν_d,
 });
 // Add points to support until within error tolerance or maximum insertion count reached.
 let mut count = 0;
 let (within_tolerances, d) = 'insertion: loop {
 μ.set_masses_dvector(&x);
 }
 // Form d = ω0 - τv - 𝒟μ = -𝒟(μ - μ^k) - τv for checking the proximate optimality
 // conditions in the predual space, and finding new points for insertion, if necessary.
-let mut d = minus_τv + op𝒟.preapply(μ_diff(μ, μ_base, ν_delta, config));
+let mut d = minus_τv + op𝒟.preapply(match ν_delta {
+None => μ_base.sub_matching(μ),
+Some(ν) =>  μ_base.sub_matching(μ) + ν
+});
 // If no merging heuristic is used, let's be more conservative about spike insertion,
 // and skip it after first round. If merging is done, being more greedy about spike
 // insertion also seems to improve performance.
 let skip_by_rough_check = if let SpikeMergingMethod::None = config.merging {
 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
 Reg : RegTerm<F, N>,
 State : AlgIteratorState {
 if state.iteration() % config.merge_every == 0 {
-let n_before_merge = μ.len();
+stats.merged += μ.merge_spikes(config.merging, |μ_candidate| {
-μ.merge_spikes(config.merging, |μ_candidate| {
+let mut d = minus_τv + op𝒟.preapply(μ_base.sub_matching(&μ_candidate));
-let μd = μ_diff(&μ_candidate, &μ_base, None, config);
-let mut d = minus_τv + op𝒟.preapply(μd);
 reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, &config)
-.then_some(())
 });
-debug_assert!(μ.len() >= n_before_merge);
-stats.merged += μ.len() - n_before_merge;
 }
 let n_before_prune = μ.len();
 μ.prune();
 debug_assert!(μ.len() <= n_before_prune);
 PlotLookup : Plotting<N>,
 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
 Reg : RegTerm<F, N> {
 // Set up parameters
-let config = &fbconfig.insertion;
+let config = &fbconfig.generic;
 let op𝒟norm = op𝒟.opnorm_bound();
 let τ = fbconfig.τ0/opA.lipschitz_factor(&op𝒟).unwrap();
 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
 // by τ compared to the conditional gradient approach.
 let tolerance = config.tolerance * τ * reg.tolerance_scaling();
 PlotLookup : Plotting<N>,
 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
 Reg : RegTerm<F, N> {
 // Set up parameters
-let config = &fbconfig.insertion;
+let config = &fbconfig.generic;
 let op𝒟norm = op𝒟.opnorm_bound();
 let τ = fbconfig.τ0/opA.lipschitz_factor(&op𝒟).unwrap();
 let mut λ = 1.0;
 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
 // by τ compared to the conditional gradient approach.

Mercurial > repos > pointsource_algs / file comparison

comparison: src/fb.rs

src/fb.rs