pointsource_algs: comparison src/prox_penalty/radon

-:9738b51d90d7
+:4f468d35fa29
 Solver for the point source localisation problem using a simplified forward-backward splitting method.
 Instead of the $𝒟$-norm of `fb.rs`, this uses a standard Radon norm for the proximal map.
 */
+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::types::*;
+use alg_tools::bounds::MinMaxMapping;
+use alg_tools::error::DynResult;
+use alg_tools::instance::{Instance, Space};
+use alg_tools::iterate::{AlgIterator, AlgIteratorIteration};
+use alg_tools::linops::BoundedLinear;
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::norms::{Norm, L2};
+use anyhow::ensure;
+use nalgebra::DVector;
 use numeric_literals::replace_float_literals;
-use serde::{Serialize, Deserialize};
+use serde::{Deserialize, Serialize};
-use nalgebra::DVector;
-use alg_tools::iterate::{
-AlgIteratorIteration,
-AlgIterator
-};
-use alg_tools::norms::{L2, Norm};
-use alg_tools::linops::Mapping;
-use alg_tools::bisection_tree::{
-BTFN,
-Bounds,
-BTSearch,
-SupportGenerator,
-LocalAnalysis,
-};
-use alg_tools::mapping::RealMapping;
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-use crate::types::*;
-use crate::measures::{
-RNDM,
-DeltaMeasure,
-Radon,
-};
-use crate::measures::merging::SpikeMerging;
-use crate::regularisation::RegTerm;
-use crate::forward_model::{
-ForwardModel,
-AdjointProductBoundedBy
-};
-use super::{
-FBGenericConfig,
-ProxPenalty,
-};
 /// Radon-norm squared proximal penalty
-#[derive(Copy,Clone,Serialize,Deserialize)]
+#[derive(Copy, Clone, Serialize, Deserialize)]
 pub struct RadonSquared;
 #[replace_float_literals(F::cast_from(literal))]
-impl<F, GA, BTA, S, Reg, const N : usize>
+impl<Domain, F, M, Reg> ProxPenalty<Domain, M, Reg, F> for RadonSquared
-ProxPenalty<F, BTFN<F, GA, BTA, N>, Reg, N> for RadonSquared
 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>,
+M: 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, GA, BTA, N>;
+type ReturnMapping = M;
+fn prox_type() -> ProxTerm {
+ProxTerm::RadonSquared
+}
 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>,
+μ_base: &DiscreteMeasure<Domain, F>,
-ν_delta: Option<&RNDM<F, N>>,
+ν_delta: Option<&DiscreteMeasure<Domain, F>>,
-τ : F,
+τ: F,
-ε : F,
+ε: F,
-config : &FBGenericConfig<F>,
+config: &InsertionConfig<F>,
-reg : &Reg,
+reg: &Reg,
-_state : &AlgIteratorIteration<I>,
+_state: &AlgIteratorIteration<I>,
-stats : &mut IterInfo<F, N>,
+stats: &mut IterInfo<F>,
-) -> (Option<Self::ReturnMapping>, bool)
+) -> DynResult<(Option<Self::ReturnMapping>, bool)>
 where
-I : AlgIterator
+I: AlgIterator,
 {
 let mut y = μ_base.masses_dvector();
-assert!(μ_base.len() <= μ.len());
+ensure!(μ_base.len() <= μ.len());
 'i_and_w: for i in 0..=1 {
 // Optimise weights
 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 g̃ = DVector::from_iterator(μ.len(),
+let g̃ = DVector::from_iterator(
-μ.iter_locations()
+μ.len(),
-.map(|ζ| - F::to_nalgebra_mixed(τv.apply(ζ))));
+μ.iter_locations()
+.map(|ζ| -F::to_nalgebra_mixed(τv.apply(ζ))),
+);
 let mut x = μ.masses_dvector();
-y.extend(std::iter::repeat(0.0.to_nalgebra_mixed()).take(0.max(x.len()-y.len())));
+y.extend(std::iter::repeat(0.0.to_nalgebra_mixed()).take(0.max(x.len() - y.len())));
 assert_eq!(y.len(), x.len());
 // Solve finite-dimensional subproblem.
 // TODO: This assumes that ν_delta has no common locations with μ-μ_base, to
 // ignore it.
 stats.inner_iters += reg.solve_findim_l1squared(&y, &g̃, τ, &mut x, ε, config);
 // Update masses of μ based on solution of finite-dimensional subproblem.
 μ.set_masses_dvector(&x);
 }
-if i>0 {
+if i > 0 {
 // Simple debugging test to see if more inserts would be needed. Doesn't seem so.
 //let n = μ.dist_matching(μ_base);
 //println!("{:?}", reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n));
-break 'i_and_w
+break 'i_and_w;
 }
 // Calculate ‖μ - μ_base‖_ℳ
 // TODO: This assumes that ν_delta has no common locations with μ-μ_base.
 let n = μ.dist_matching(μ_base) + ν_delta.map_or(0.0, |ν| ν.norm(Radon));
 // Find a spike to insert, if needed.
 // This only check the overall tolerances, not tolerances on support of μ-μ_base or μ,
 // which are supposed to have been guaranteed by the finite-dimensional weight optimisation.
 match reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n) {
-None => { break 'i_and_w },
+None => break 'i_and_w,
 Some((ξ, _v_ξ, _in_bounds)) => {
 // Weight is found out by running the finite-dimensional optimisation algorithm
 // above
-*μ += DeltaMeasure { x : ξ, α : 0.0 };
+*μ += DeltaMeasure { x: ξ, α: 0.0 };
 stats.inserted += 1;
 }
 };
 }
-(None, true)
+Ok((None, true))
 }
 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>>,
+ν_delta: Option<&DiscreteMeasure<Domain, F>>,
-τ : F,
+τ: F,
-ε : F,
+ε: F,
-config : &FBGenericConfig<F>,
+config: &InsertionConfig<F>,
-reg : &Reg,
+reg: &Reg,
-fitness : Option<impl Fn(&RNDM<F, N>) -> F>,
+fitness: Option<impl Fn(&DiscreteMeasure<Domain, F>) -> 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| {
 // Important: μ_candidate's new points are afterwards,
 // and do not conflict with μ_base.
 //reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n).is_none()
 })
 }
 }
+#[replace_float_literals(F::cast_from(literal))]
-impl<F, A, const N : usize> AdjointProductBoundedBy<RNDM<F, N>, RadonSquared>
+impl<'a, F, A, Domain> StepLengthBound<F, QuadraticDataTerm<F, Domain, A>> for RadonSquared
-for A
 where
-F : Float,
+F: Float + ToNalgebraRealField,
-A : ForwardModel<RNDM<F, N>, F>
+Domain: Space + Norm<Radon, F>,
+A: ForwardModel<Domain, F> + BoundedLinear<Domain, Radon, L2, F>,
 {
-type FloatType = F;
+fn step_length_bound(&self, f: &QuadraticDataTerm<F, Domain, A>) -> DynResult<F> {
+// TODO: direct squared calculation
-fn adjoint_product_bound(&self, _ : &RadonSquared) -> Option<Self::FloatType> {
+Ok(f.operator().opnorm_bound(Radon, L2)?.powi(2))
-self.opnorm_bound(Radon, L2).powi(2).into()
 }
 }
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F, A, Domain> StepLengthBoundPD<F, A, DiscreteMeasure<Domain, F>> for RadonSquared
+where
+Domain: Space + Clone + PartialEq + 'static,
+F: Float + ToNalgebraRealField,
+A: BoundedLinear<DiscreteMeasure<Domain, F>, Radon, L2, F>,
+{
+fn step_length_bound_pd(&self, opA: &A) -> DynResult<F> {
+opA.opnorm_bound(Radon, L2)
+}
+}

comparison: src/prox_penalty/radon_squared.rs

src/prox_penalty/radon_squared.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/prox_penalty/radon_squared.rs

src/prox_penalty/radon_squared.rs