pointsource_algs: comparison src/prox

-:6099ba025aac
+:ed16d0f10d08
 /*!
 Basic proximal penalty based on convolution operators $𝒟$.
 */
+use super::{InsertionConfig, ProxPenalty, ProxTerm, StepLengthBound, StepLengthBoundPD};
+use crate::dataterm::QuadraticDataTerm;
+use crate::forward_model::ForwardModel;
+use crate::measures::merging::SpikeMerging;
+use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon};
+use crate::regularisation::RegTerm;
+use crate::seminorms::DiscreteMeasureOp;
+use crate::types::IterInfo;
+use alg_tools::bounds::MinMaxMapping;
+use alg_tools::error::DynResult;
+use alg_tools::iterate::{AlgIterator, AlgIteratorIteration};
+use alg_tools::linops::BoundedLinear;
+use alg_tools::mapping::{Instance, Mapping, Space};
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::norms::{Linfinity, Norm, NormExponent, L2};
+use alg_tools::types::*;
+use colored::Colorize;
+use nalgebra::DVector;
 use numeric_literals::replace_float_literals;
-use nalgebra::DVector;
-use colored::Colorize;
-use alg_tools::types::*;
-use alg_tools::loc::Loc;
-use alg_tools::mapping::{Mapping, RealMapping};
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::Linfinity;
-use alg_tools::iterate::{
-AlgIteratorIteration,
-AlgIterator,
-};
-use alg_tools::bisection_tree::{
-BTFN,
-PreBTFN,
-Bounds,
-BTSearch,
-SupportGenerator,
-LocalAnalysis,
-BothGenerators,
-};
-use crate::measures::{
-RNDM,
-DeltaMeasure,
-Radon,
-};
-use crate::measures::merging::{
-SpikeMerging,
-};
-use crate::seminorms::DiscreteMeasureOp;
-use crate::types::{
-IterInfo,
-};
-use crate::regularisation::RegTerm;
-use super::{ProxPenalty, FBGenericConfig};
 #[replace_float_literals(F::cast_from(literal))]
-impl<F, GA, BTA, S, Reg, 𝒟, G𝒟, K, const N : usize>
+impl<F, M, Reg, 𝒟, O, Domain> ProxPenalty<Domain, M, Reg, F> for 𝒟
-ProxPenalty<F, BTFN<F, GA, BTA, N>, Reg, N> for 𝒟
 where
-F : Float + ToNalgebraRealField,
+Domain: Space + Clone + PartialEq + 'static,
-GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
+for<'a> &'a Domain: Instance<Domain>,
-BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
+F: Float + ToNalgebraRealField,
-S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
+𝒟: DiscreteMeasureOp<Domain, F>,
-G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone,
+𝒟::Codomain: Mapping<Domain, Codomain = F>,
-𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>,
+M: Mapping<Domain, Codomain = F>,
-𝒟::Codomain : RealMapping<F, N>,
+for<'a> &'a M: std::ops::Add<𝒟::PreCodomain, Output = O>,
-K : RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
+O: MinMaxMapping<Domain, F>,
-Reg : RegTerm<F, N>,
+Reg: RegTerm<Domain, F>,
-RNDM<F, N> : SpikeMerging<F>,
+DiscreteMeasure<Domain, F>: SpikeMerging<F>,
 {
-type ReturnMapping = BTFN<F, BothGenerators<GA, G𝒟>, BTA, N>;
+type ReturnMapping = O;
+fn prox_type() -> ProxTerm {
+ProxTerm::Wave
+}
 fn insert_and_reweigh<I>(
 &self,
-μ : &mut RNDM<F, N>,
+μ: &mut DiscreteMeasure<Domain, F>,
-τv : &mut BTFN<F, GA, BTA, N>,
+τv: &mut M,
-μ_base : &RNDM<F, N>,
+τ: F,
-ν_delta: Option<&RNDM<F, N>>,
+ε: F,
-τ : F,
+config: &InsertionConfig<F>,
-ε : F,
+reg: &Reg,
-config : &FBGenericConfig<F>,
+state: &AlgIteratorIteration<I>,
-reg : &Reg,
+stats: &mut IterInfo<F>,
-state : &AlgIteratorIteration<I>,
+) -> DynResult<(Option<Self::ReturnMapping>, bool)>
-stats : &mut IterInfo<F, N>,
-) -> (Option<BTFN<F, BothGenerators<GA, G𝒟>, BTA, N>>, bool)
 where
-I : AlgIterator
+I: AlgIterator,
 {
+let op𝒟norm = self.opnorm_bound(Radon, Linfinity)?;
-let op𝒟norm = self.opnorm_bound(Radon, Linfinity);
 // Maximum insertion count and measure difference calculation depend on insertion style.
-let (max_insertions, warn_insertions) = match (state.iteration(), config.bootstrap_insertions) {
+let (max_insertions, warn_insertions) =
-(i, Some((l, k))) if i <= l => (k, false),
+match (state.iteration(), config.bootstrap_insertions) {
-_ => (config.max_insertions, !state.is_quiet()),
+(i, Some((l, k))) if i <= l => (k, false),
-};
+_ => (config.max_insertions, !state.is_quiet()),
+};
-let ω0 = match ν_delta {
-None => self.apply(μ_base),
+let μ_base = μ.clone();
-Some(ν) => self.apply(μ_base + ν),
+let ω0 = self.apply(&μ_base);
-};
 // Add points to support until within error tolerance or maximum insertion count reached.
 let mut count = 0;
 let (within_tolerances, d) = 'insertion: loop {
 if μ.len() > 0 {
 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
 // from the beginning of the iteration are all contained in the immutable c and g.
 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
 // problems have not yet been updated to sign change.
 let Ã = self.findim_matrix(μ.iter_locations());
-let g̃ = DVector::from_iterator(μ.len(),
+let g̃ = DVector::from_iterator(
-μ.iter_locations()
+μ.len(),
-.map(|ζ| ω0.apply(ζ) - τv.apply(ζ))
+μ.iter_locations()
-.map(F::to_nalgebra_mixed));
+.map(|ζ| ω0.apply(ζ) - τv.apply(ζ))
+.map(F::to_nalgebra_mixed),
+);
 let mut x = μ.masses_dvector();
 // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃.
 // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩
 // ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟Cx|_∞ ≤  sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟| |Cx|_ℳ
 μ.set_masses_dvector(&x);
 }
 // Form d = τv + 𝒟μ - ω0 = τv + 𝒟(μ - μ^k) for checking the proximate optimality
 // conditions in the predual space, and finding new points for insertion, if necessary.
-let mut d = &*τv + match ν_delta {
+let mut d = &*τv + self.preapply(μ.sub_matching(&μ_base));
-None => self.preapply(μ.sub_matching(μ_base)),
-Some(ν) => self.preapply(μ.sub_matching(μ_base) - ν)
-};
 // 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 config.merging.enabled {
 } else {
 count > 0
 };
 // Find a spike to insert, if needed
-let (ξ, _v_ξ, in_bounds) =  match reg.find_tolerance_violation(
+let (ξ, _v_ξ, in_bounds) =
-&mut d, τ, ε, skip_by_rough_check, config
+match reg.find_tolerance_violation(&mut d, τ, ε, skip_by_rough_check, config) {
-) {
+None => break 'insertion (true, d),
-None => break 'insertion (true, d),
+Some(res) => res,
-Some(res) => res,
+};
-};
 // Break if maximum insertion count reached
 if count >= max_insertions {
-break 'insertion (in_bounds, d)
+break 'insertion (in_bounds, d);
 }
 // No point in optimising the weight here; the finite-dimensional algorithm is fast.
-*μ += DeltaMeasure { x : ξ, α : 0.0 };
+*μ += DeltaMeasure { x: ξ, α: 0.0 };
 count += 1;
-stats.inserted += 1;
 };
 if !within_tolerances && warn_insertions {
 // Complain (but continue) if we failed to get within tolerances
 // by inserting more points.
-let err = format!("Maximum insertions reached without achieving \
+let err = format!(
-subproblem solution tolerance");
+"Maximum insertions reached without achieving \
+subproblem solution tolerance"
+);
 println!("{}", err.red());
 }
-(Some(d), within_tolerances)
+Ok((Some(d), within_tolerances))
 }
 fn merge_spikes(
 &self,
-μ : &mut RNDM<F, N>,
+μ: &mut DiscreteMeasure<Domain, F>,
-τv : &mut BTFN<F, GA, BTA, N>,
+τv: &mut M,
-μ_base : &RNDM<F, N>,
+μ_base: &DiscreteMeasure<Domain, F>,
-ν_delta: Option<&RNDM<F, N>>,
+τ: F,
-τ : F,
+ε: F,
-ε : F,
+config: &InsertionConfig<F>,
-config : &FBGenericConfig<F>,
+reg: &Reg,
-reg : &Reg,
+fitness: Option<impl Fn(&DiscreteMeasure<Domain, F>) -> F>,
-fitness : Option<impl Fn(&RNDM<F, N>) -> F>,
+) -> usize {
-) -> usize
-{
 if config.fitness_merging {
 if let Some(f) = fitness {
-return μ.merge_spikes_fitness(config.merging, f, |&v| v)
+return μ.merge_spikes_fitness(config.merging, f, |&v| v).1;
-.1
 }
 }
 μ.merge_spikes(config.merging, |μ_candidate| {
-let mut d = &*τv + self.preapply(match ν_delta {
+let mut d = &*τv + self.preapply(μ_candidate.sub_matching(μ_base));
-None => μ_candidate.sub_matching(μ_base),
-Some(ν) => μ_candidate.sub_matching(μ_base) - ν,
-});
 reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, config)
 })
 }
 }
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F, A, 𝒟, Domain> StepLengthBound<F, QuadraticDataTerm<F, DiscreteMeasure<Domain, F>, A>>
+for 𝒟
+where
+Domain: Space + Clone + PartialEq + 'static,
+F: Float + ToNalgebraRealField,
+𝒟: DiscreteMeasureOp<Domain, F>,
+A: ForwardModel<DiscreteMeasure<Domain, F>, F>
++ for<'b> BoundedLinear<DiscreteMeasure<Domain, F>, &'b 𝒟, L2, F>,
+DiscreteMeasure<Domain, F>: for<'b> Norm<&'b 𝒟, F>,
+for<'b> &'b 𝒟: NormExponent,
+{
+fn step_length_bound(
+&self,
+f: &QuadraticDataTerm<F, DiscreteMeasure<Domain, F>, A>,
+) -> DynResult<F> {
+// TODO: direct squared calculation
+Ok(f.operator().opnorm_bound(self, L2)?.powi(2))
+}
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F, A, 𝒟, Domain> StepLengthBoundPD<F, A, DiscreteMeasure<Domain, F>> for 𝒟
+where
+Domain: Space + Clone + PartialEq + 'static,
+F: Float + ToNalgebraRealField,
+𝒟: DiscreteMeasureOp<Domain, F>,
+A: for<'a> BoundedLinear<DiscreteMeasure<Domain, F>, &'a 𝒟, L2, F>,
+DiscreteMeasure<Domain, F>: for<'a> Norm<&'a 𝒟, F>,
+for<'b> &'b 𝒟: NormExponent,
+{
+fn step_length_bound_pd(&self, opA: &A) -> DynResult<F> {
+opA.opnorm_bound(self, L2)
+}
+}

comparison: src/prox_penalty/wave.rs

src/prox_penalty/wave.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/prox_penalty/wave.rs

src/prox_penalty/wave.rs