Mercurial > repos > pointsource_algs / file comparison

comparison: src/prox_penalty/radon_squared.rs

src/prox_penalty/radon_squared.rs

changeset 70
ed16d0f10d08
parent 63
7a8a55fd41c0
equal deleted inserted replaced
58:6099ba025aac 70:ed16d0f10d08
2 Solver for the point source localisation problem using a simplified forward-backward splitting method. 2 Solver for the point source localisation problem using a simplified forward-backward splitting method.
3 3
4 Instead of the $𝒟$-norm of `fb.rs`, this uses a standard Radon norm for the proximal map. 4 Instead of the $𝒟$-norm of `fb.rs`, this uses a standard Radon norm for the proximal map.
5 */ 5 */
6 6
7 use super::{InsertionConfig, ProxPenalty, ProxTerm, StepLengthBound, StepLengthBoundPD};
8 use crate::dataterm::QuadraticDataTerm;
9 use crate::forward_model::ForwardModel;
10 use crate::measures::merging::SpikeMerging;
11 use crate::measures::{DiscreteMeasure, Radon};
12 use crate::regularisation::RadonSquaredRegTerm;
13 use crate::types::*;
14 use alg_tools::bounds::MinMaxMapping;
15 use alg_tools::error::DynResult;
16 use alg_tools::instance::{Instance, Space};
17 use alg_tools::iterate::{AlgIterator, AlgIteratorIteration};
18 use alg_tools::linops::BoundedLinear;
19 use alg_tools::nalgebra_support::ToNalgebraRealField;
20 use alg_tools::norms::{Norm, L2};
7 use numeric_literals::replace_float_literals; 21 use numeric_literals::replace_float_literals;
8 use serde::{Serialize, Deserialize}; 22 use serde::{Deserialize, Serialize};
9 use nalgebra::DVector;
10
11 use alg_tools::iterate::{
12 AlgIteratorIteration,
13 AlgIterator
14 };
15 use alg_tools::norms::{L2, Norm};
16 use alg_tools::linops::Mapping;
17 use alg_tools::bisection_tree::{
18 BTFN,
19 Bounds,
20 BTSearch,
21 SupportGenerator,
22 LocalAnalysis,
23 };
24 use alg_tools::mapping::RealMapping;
25 use alg_tools::nalgebra_support::ToNalgebraRealField;
26
27 use crate::types::*;
28 use crate::measures::{
29 RNDM,
30 DeltaMeasure,
31 Radon,
32 };
33 use crate::measures::merging::SpikeMerging;
34 use crate::regularisation::RegTerm;
35 use crate::forward_model::{
36 ForwardModel,
37 AdjointProductBoundedBy
38 };
39 use super::{
40 FBGenericConfig,
41 ProxPenalty,
42 };
43 23
44 /// Radon-norm squared proximal penalty 24 /// Radon-norm squared proximal penalty
45 25
46 #[derive(Copy,Clone,Serialize,Deserialize)] 26 #[derive(Copy, Clone, Serialize, Deserialize)]
47 pub struct RadonSquared; 27 pub struct RadonSquared;
48 28
49 #[replace_float_literals(F::cast_from(literal))] 29 #[replace_float_literals(F::cast_from(literal))]
50 impl<F, GA, BTA, S, Reg, const N : usize> 30 impl<Domain, F, M, Reg> ProxPenalty<Domain, M, Reg, F> for RadonSquared
51 ProxPenalty<F, BTFN<F, GA, BTA, N>, Reg, N> for RadonSquared
52 where 31 where
53 F : Float + ToNalgebraRealField, 32 Domain: Space + Clone + PartialEq + 'static,
54 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone, 33 for<'a> &'a Domain: Instance<Domain>,
55 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>, 34 F: Float + ToNalgebraRealField,
56 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 35 M: MinMaxMapping<Domain, F>,
57 Reg : RegTerm<F, N>, 36 Reg: RadonSquaredRegTerm<Domain, F>,
58 RNDM<F, N> : SpikeMerging<F>, 37 DiscreteMeasure<Domain, F>: SpikeMerging<F>,
59 { 38 {
60 type ReturnMapping = BTFN<F, GA, BTA, N>; 39 type ReturnMapping = M;
40
41 fn prox_type() -> ProxTerm {
42 ProxTerm::RadonSquared
43 }
61 44
62 fn insert_and_reweigh<I>( 45 fn insert_and_reweigh<I>(
63 &self, 46 &self,
64 μ : &mut RNDM<F, N>, 47 μ: &mut DiscreteMeasure<Domain, F>,
65 τv : &mut BTFN<F, GA, BTA, N>, 48 τv: &mut M,
66 μ_base : &RNDM<F, N>, 49 τ: F,
67 ν_delta: Option<&RNDM<F, N>>, 50 ε: F,
68 τ : F, 51 config: &InsertionConfig<F>,
69 ε : F, 52 reg: &Reg,
70 config : &FBGenericConfig<F>, 53 _state: &AlgIteratorIteration<I>,
71 reg : &Reg, 54 stats: &mut IterInfo<F>,
72 _state : &AlgIteratorIteration<I>, 55 ) -> DynResult<(Option<Self::ReturnMapping>, bool)>
73 stats : &mut IterInfo<F, N>,
74 ) -> (Option<Self::ReturnMapping>, bool)
75 where 56 where
76 I : AlgIterator 57 I: AlgIterator,
77 { 58 {
78 let mut y = μ_base.masses_dvector(); 59 let violation = reg.find_tolerance_violation(τv, τ, ε, true, config);
60 reg.solve_oc_radonsq(μ, τv, τ, ε, violation, config, stats);
79 61
80 assert!(μ_base.len() <= μ.len()); 62 Ok((None, true))
81
82 'i_and_w: for i in 0..=1 {
83 // Optimise weights
84 if μ.len() > 0 {
85 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
86 // from the beginning of the iteration are all contained in the immutable c and g.
87 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
88 // problems have not yet been updated to sign change.
89 let g̃ = DVector::from_iterator(μ.len(),
90 μ.iter_locations()
91 .map(|ζ| - F::to_nalgebra_mixed(τv.apply(ζ))));
92 let mut x = μ.masses_dvector();
93 y.extend(std::iter::repeat(0.0.to_nalgebra_mixed()).take(0.max(x.len()-y.len())));
94 assert_eq!(y.len(), x.len());
95 // Solve finite-dimensional subproblem.
96 // TODO: This assumes that ν_delta has no common locations with μ-μ_base, to
97 // ignore it.
98 stats.inner_iters += reg.solve_findim_l1squared(&y, &g̃, τ, &mut x, ε, config);
99
100 // Update masses of μ based on solution of finite-dimensional subproblem.
101 μ.set_masses_dvector(&x);
102 }
103
104 if i>0 {
105 // Simple debugging test to see if more inserts would be needed. Doesn't seem so.
106 //let n = μ.dist_matching(μ_base);
107 //println!("{:?}", reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n));
108 break 'i_and_w
109 }
110
111 // Calculate ‖μ - μ_base‖_ℳ
112 // TODO: This assumes that ν_delta has no common locations with μ-μ_base.
113 let n = μ.dist_matching(μ_base) + ν_delta.map_or(0.0, |ν| ν.norm(Radon));
114
115 // Find a spike to insert, if needed.
116 // This only check the overall tolerances, not tolerances on support of μ-μ_base or μ,
117 // which are supposed to have been guaranteed by the finite-dimensional weight optimisation.
118 match reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n) {
119 None => { break 'i_and_w },
120 Some((ξ, _v_ξ, _in_bounds)) => {
121 // Weight is found out by running the finite-dimensional optimisation algorithm
122 // above
123 *μ += DeltaMeasure { x : ξ, α : 0.0 };
124 stats.inserted += 1;
125 }
126 };
127 }
128
129 (None, true)
130 } 63 }
131 64
132 fn merge_spikes( 65 fn merge_spikes(
133 &self, 66 &self,
134 μ : &mut RNDM<F, N>, 67 μ: &mut DiscreteMeasure<Domain, F>,
135 τv : &mut BTFN<F, GA, BTA, N>, 68 τv: &mut M,
136 μ_base : &RNDM<F, N>, 69 μ_base: &DiscreteMeasure<Domain, F>,
137 ν_delta: Option<&RNDM<F, N>>, 70 τ: F,
138 τ : F, 71 ε: F,
139 ε : F, 72 config: &InsertionConfig<F>,
140 config : &FBGenericConfig<F>, 73 reg: &Reg,
141 reg : &Reg, 74 fitness: Option<impl Fn(&DiscreteMeasure<Domain, F>) -> F>,
142 fitness : Option<impl Fn(&RNDM<F, N>) -> F>, 75 ) -> usize {
143 ) -> usize
144 {
145 if config.fitness_merging { 76 if config.fitness_merging {
146 if let Some(f) = fitness { 77 if let Some(f) = fitness {
147 return μ.merge_spikes_fitness(config.merging, f, |&v| v) 78 return μ.merge_spikes_fitness(config.merging, f, |&v| v).1;
148 .1
149 } 79 }
150 } 80 }
151 μ.merge_spikes(config.merging, |μ_candidate| { 81 μ.merge_spikes(config.merging, |μ_candidate| {
152 // Important: μ_candidate's new points are afterwards, 82 // Important: μ_candidate's new points are afterwards,
153 // and do not conflict with μ_base. 83 // and do not conflict with μ_base.
154 // TODO: could simplify to requiring μ_base instead of μ_radon. 84 // TODO: could simplify to requiring μ_base instead of μ_radon.
155 // but may complicate with sliding base's exgtra points that need to be 85 // but may complicate with sliding base's exgtra points that need to be
156 // after μ_candidate's extra points. 86 // after μ_candidate's extra points.
157 // TODO: doesn't seem to work, maybe need to merge μ_base as well? 87 // TODO: doesn't seem to work, maybe need to merge μ_base as well?
158 // Although that doesn't seem to make sense. 88 // Although that doesn't seem to make sense.
159 let μ_radon = match ν_delta { 89 let μ_radon = μ_candidate.sub_matching(μ_base);
160 None => μ_candidate.sub_matching(μ_base),
161 Some(ν) => μ_candidate.sub_matching(μ_base) - ν,
162 };
163 reg.verify_merge_candidate_radonsq(τv, μ_candidate, τ, ε, &config, &μ_radon) 90 reg.verify_merge_candidate_radonsq(τv, μ_candidate, τ, ε, &config, &μ_radon)
164 //let n = μ_candidate.dist_matching(μ_base); 91 //let n = μ_candidate.dist_matching(μ_base);
165 //reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n).is_none() 92 //reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n).is_none()
166 }) 93 })
167 } 94 }
168 } 95 }
169 96
170 97 #[replace_float_literals(F::cast_from(literal))]
171 impl<F, A, const N : usize> AdjointProductBoundedBy<RNDM<F, N>, RadonSquared> 98 impl<'a, F, A, Domain> StepLengthBound<F, QuadraticDataTerm<F, Domain, A>> for RadonSquared
172 for A
173 where 99 where
174 F : Float, 100 F: Float + ToNalgebraRealField,
175 A : ForwardModel<RNDM<F, N>, F> 101 Domain: Space + Norm<Radon, F>,
102 A: ForwardModel<Domain, F> + BoundedLinear<Domain, Radon, L2, F>,
176 { 103 {
177 type FloatType = F; 104 fn step_length_bound(&self, f: &QuadraticDataTerm<F, Domain, A>) -> DynResult<F> {
178 105 // TODO: direct squared calculation
179 fn adjoint_product_bound(&self, _ : &RadonSquared) -> Option<Self::FloatType> { 106 Ok(f.operator().opnorm_bound(Radon, L2)?.powi(2))
180 self.opnorm_bound(Radon, L2).powi(2).into()
181 } 107 }
182 } 108 }
109
110 #[replace_float_literals(F::cast_from(literal))]
111 impl<'a, F, A, Domain> StepLengthBoundPD<F, A, DiscreteMeasure<Domain, F>> for RadonSquared
112 where
113 Domain: Space + Clone + PartialEq + 'static,
114 F: Float + ToNalgebraRealField,
115 A: BoundedLinear<DiscreteMeasure<Domain, F>, Radon, L2, F>,
116 {
117 fn step_length_bound_pd(&self, opA: &A) -> DynResult<F> {
118 opA.opnorm_bound(Radon, L2)
119 }
120 }

Mercurial Repository: pointsource_algs

mercurial