Mercurial > repos > pointsource_algs / file comparison

comparison: src/regularisation.rs

src/regularisation.rs

changeset 70
ed16d0f10d08
parent 63
7a8a55fd41c0
equal deleted inserted replaced
58:6099ba025aac 70:ed16d0f10d08
2 Regularisation terms 2 Regularisation terms
3 */ 3 */
4 4
5 #[allow(unused_imports)] // Used by documentation. 5 #[allow(unused_imports)] // Used by documentation.
6 use crate::fb::pointsource_fb_reg; 6 use crate::fb::pointsource_fb_reg;
7 use crate::fb::FBGenericConfig; 7 use crate::fb::InsertionConfig;
8 use crate::measures::{DeltaMeasure, Radon, RNDM}; 8 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
9 #[allow(unused_imports)] // Used by documentation. 9 #[allow(unused_imports)] // Used by documentation.
10 use crate::sliding_fb::pointsource_sliding_fb_reg; 10 use crate::sliding_fb::pointsource_sliding_fb_reg;
11 use crate::types::*; 11 use crate::types::*;
12 use alg_tools::instance::Instance; 12 use alg_tools::instance::{Instance, Space};
13 use alg_tools::linops::Mapping; 13 use alg_tools::linops::Mapping;
14 use alg_tools::loc::Loc; 14 use alg_tools::loc::Loc;
15 use alg_tools::norms::Norm; 15 use alg_tools::norms::Norm;
16 use numeric_literals::replace_float_literals; 16 use numeric_literals::replace_float_literals;
17 use serde::{Deserialize, Serialize}; 17 use serde::{Deserialize, Serialize};
18 18
19 use crate::subproblem::{ 19 use crate::subproblem::{
20 l1squared_nonneg::l1squared_nonneg, l1squared_unconstrained::l1squared_unconstrained, 20 l1squared_nonneg::l1squared_nonneg, l1squared_unconstrained::l1squared_unconstrained,
21 nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained, 21 nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained,
22 }; 22 };
23 use alg_tools::bisection_tree::{ 23 use alg_tools::bounds::{Bounds, MinMaxMapping};
24 BTSearch, Bounded, Bounds, LocalAnalysis, P2Minimise, SupportGenerator, BTFN,
25 };
26 use alg_tools::iterate::AlgIteratorFactory; 24 use alg_tools::iterate::AlgIteratorFactory;
27 use alg_tools::nalgebra_support::ToNalgebraRealField; 25 use alg_tools::nalgebra_support::ToNalgebraRealField;
28 use nalgebra::{DMatrix, DVector}; 26 use nalgebra::{DMatrix, DVector};
29 27
30 use std::cmp::Ordering::{Equal, Greater, Less}; 28 use std::cmp::Ordering::{Equal, Greater, Less};
32 /// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other 30 /// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other
33 /// algorithms. 31 /// algorithms.
34 /// 32 ///
35 /// The only member of the struct is the regularisation parameter α. 33 /// The only member of the struct is the regularisation parameter α.
36 #[derive(Copy, Clone, Debug, Serialize, Deserialize)] 34 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
37 pub struct NonnegRadonRegTerm<F: Float>(pub F /* α */); 35 pub struct NonnegRadonRegTerm<F: Float = f64>(pub F /* α */);
38 36
39 impl<'a, F: Float> NonnegRadonRegTerm<F> { 37 impl<'a, F: Float> NonnegRadonRegTerm<F> {
40 /// Returns the regularisation parameter 38 /// Returns the regularisation parameter
41 pub fn α(&self) -> F { 39 pub fn α(&self) -> F {
42 let &NonnegRadonRegTerm(α) = self; 40 let &NonnegRadonRegTerm(α) = self;
43 α 41 α
44 } 42 }
45 } 43 }
46 44
47 impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for NonnegRadonRegTerm<F> { 45 impl<'a, F: Float, const N: usize> Mapping<RNDM<N, F>> for NonnegRadonRegTerm<F> {
48 type Codomain = F; 46 type Codomain = F;
49 47
50 fn apply<I>(&self, μ: I) -> F 48 fn apply<I>(&self, μ: I) -> F
51 where 49 where
52 I: Instance<RNDM<F, N>>, 50 I: Instance<RNDM<N, F>>,
53 { 51 {
54 self.α() * μ.eval(|x| x.norm(Radon)) 52 self.α() * μ.eval(|x| x.norm(Radon))
55 } 53 }
56 } 54 }
57 55
58 /// The regularisation term $α\|μ\|_{ℳ(Ω)}$ for [`pointsource_fb_reg`]. 56 /// The regularisation term $α\|μ\|_{ℳ(Ω)}$ for [`pointsource_fb_reg`].
59 /// 57 ///
60 /// The only member of the struct is the regularisation parameter α. 58 /// The only member of the struct is the regularisation parameter α.
61 #[derive(Copy, Clone, Debug, Serialize, Deserialize)] 59 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
62 pub struct RadonRegTerm<F: Float>(pub F /* α */); 60 pub struct RadonRegTerm<F: Float = f64>(pub F /* α */);
63 61
64 impl<'a, F: Float> RadonRegTerm<F> { 62 impl<'a, F: Float> RadonRegTerm<F> {
65 /// Returns the regularisation parameter 63 /// Returns the regularisation parameter
66 pub fn α(&self) -> F { 64 pub fn α(&self) -> F {
67 let &RadonRegTerm(α) = self; 65 let &RadonRegTerm(α) = self;
68 α 66 α
69 } 67 }
70 } 68 }
71 69
72 impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for RadonRegTerm<F> { 70 impl<'a, F: Float, const N: usize> Mapping<RNDM<N, F>> for RadonRegTerm<F> {
73 type Codomain = F; 71 type Codomain = F;
74 72
75 fn apply<I>(&self, μ: I) -> F 73 fn apply<I>(&self, μ: I) -> F
76 where 74 where
77 I: Instance<RNDM<F, N>>, 75 I: Instance<RNDM<N, F>>,
78 { 76 {
79 self.α() * μ.eval(|x| x.norm(Radon)) 77 self.α() * μ.eval(|x| x.norm(Radon))
80 } 78 }
81 } 79 }
82 80
83 /// Regularisation term configuration 81 /// Regularisation term configuration
84 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] 82 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
85 pub enum Regularisation<F: Float> { 83 pub enum Regularisation<F: Float = f64> {
86 /// $α \\|μ\\|\_{ℳ(Ω)}$ 84 /// $α \\|μ\\|\_{ℳ(Ω)}$
87 Radon(F), 85 Radon(F),
88 /// $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ 86 /// $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$
89 NonnegRadon(F), 87 NonnegRadon(F),
90 } 88 }
91 89
92 impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for Regularisation<F> { 90 impl<'a, F: Float, const N: usize> Mapping<RNDM<N, F>> for Regularisation<F> {
93 type Codomain = F; 91 type Codomain = F;
94 92
95 fn apply<I>(&self, μ: I) -> F 93 fn apply<I>(&self, μ: I) -> F
96 where 94 where
97 I: Instance<RNDM<F, N>>, 95 I: Instance<RNDM<N, F>>,
98 { 96 {
99 match *self { 97 match *self {
100 Self::Radon(α) => RadonRegTerm(α).apply(μ), 98 Self::Radon(α) => RadonRegTerm(α).apply(μ),
101 Self::NonnegRadon(α) => NonnegRadonRegTerm(α).apply(μ), 99 Self::NonnegRadon(α) => NonnegRadonRegTerm(α).apply(μ),
102 } 100 }
103 } 101 }
104 } 102 }
105 103
106 /// Abstraction of regularisation terms. 104 /// Abstraction of regularisation terms.
107 pub trait RegTerm<F: Float + ToNalgebraRealField, const N: usize>: 105 pub trait RegTerm<Domain, F = f64>: Mapping<DiscreteMeasure<Domain, F>, Codomain = F>
108 Mapping<RNDM<F, N>, Codomain = F> 106 where
107 Domain: Space + Clone,
108 F: Float + ToNalgebraRealField,
109 { 109 {
110 /// Approximately solve the problem 110 /// Approximately solve the problem
111 /// <div>$$ 111 /// <div>$$
112 /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + τ G(x) 112 /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + τ G(x)
113 /// $$</div> 113 /// $$</div>
123 g: &DVector<F::MixedType>, 123 g: &DVector<F::MixedType>,
124 τ: F, 124 τ: F,
125 x: &mut DVector<F::MixedType>, 125 x: &mut DVector<F::MixedType>,
126 mA_normest: F, 126 mA_normest: F,
127 ε: F, 127 ε: F,
128 config: &FBGenericConfig<F>, 128 config: &InsertionConfig<F>,
129 ) -> usize;
130
131 /// Approximately solve the problem
132 /// <div>$$
133 /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x)
134 /// $$</div>
135 /// for $G$ depending on the trait implementation.
136 ///
137 /// Returns the number of iterations taken.
138 fn solve_findim_l1squared(
139 &self,
140 y: &DVector<F::MixedType>,
141 g: &DVector<F::MixedType>,
142 τ: F,
143 x: &mut DVector<F::MixedType>,
144 ε: F,
145 config: &FBGenericConfig<F>,
146 ) -> usize; 129 ) -> usize;
147 130
148 /// Find a point where `d` may violate the tolerance `ε`. 131 /// Find a point where `d` may violate the tolerance `ε`.
149 /// 132 ///
150 /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we 133 /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
152 /// regulariser. 135 /// regulariser.
153 /// 136 ///
154 /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check 137 /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check
155 /// terminating early. Otherwise returns a possibly violating point, the value of `d` there, 138 /// terminating early. Otherwise returns a possibly violating point, the value of `d` there,
156 /// and a boolean indicating whether the found point is in bounds. 139 /// and a boolean indicating whether the found point is in bounds.
157 fn find_tolerance_violation<G, BT>( 140 fn find_tolerance_violation<M>(
158 &self, 141 &self,
159 d: &mut BTFN<F, G, BT, N>, 142 d: &mut M,
160 τ: F, 143 τ: F,
161 ε: F, 144 ε: F,
162 skip_by_rough_check: bool, 145 skip_by_rough_check: bool,
163 config: &FBGenericConfig<F>, 146 config: &InsertionConfig<F>,
164 ) -> Option<(Loc<F, N>, F, bool)> 147 ) -> Option<(Domain, F, bool)>
165 where 148 where
166 BT: BTSearch<F, N, Agg = Bounds<F>>, 149 M: MinMaxMapping<Domain, F>;
167 G: SupportGenerator<F, N, Id = BT::Data>,
168 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
169 {
170 self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO)
171 }
172
173 /// Find a point where `d` may violate the tolerance `ε`.
174 ///
175 /// This version includes a `slack` parameter to expand the tolerances.
176 /// It is used for Radon-norm squared proximal term in [`crate::prox_penalty::radon_squared`].
177 ///
178 /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
179 /// are in bounds. `ε` is the current main tolerance and `τ` a scaling factor for the
180 /// regulariser.
181 ///
182 /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check
183 /// terminating early. Otherwise returns a possibly violating point, the value of `d` there,
184 /// and a boolean indicating whether the found point is in bounds.
185 fn find_tolerance_violation_slack<G, BT>(
186 &self,
187 d: &mut BTFN<F, G, BT, N>,
188 τ: F,
189 ε: F,
190 skip_by_rough_check: bool,
191 config: &FBGenericConfig<F>,
192 slack: F,
193 ) -> Option<(Loc<F, N>, F, bool)>
194 where
195 BT: BTSearch<F, N, Agg = Bounds<F>>,
196 G: SupportGenerator<F, N, Id = BT::Data>,
197 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
198 150
199 /// Verify that `d` is in bounds `ε` for a merge candidate `μ` 151 /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
200 /// 152 ///
201 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser. 153 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
202 fn verify_merge_candidate<G, BT>( 154 fn verify_merge_candidate<M>(
203 &self, 155 &self,
204 d: &mut BTFN<F, G, BT, N>, 156 d: &mut M,
205 μ: &RNDM<F, N>, 157 μ: &DiscreteMeasure<Domain, F>,
206 τ: F, 158 τ: F,
207 ε: F, 159 ε: F,
208 config: &FBGenericConfig<F>, 160 config: &InsertionConfig<F>,
209 ) -> bool 161 ) -> bool
210 where 162 where
211 BT: BTSearch<F, N, Agg = Bounds<F>>, 163 M: MinMaxMapping<Domain, F>;
212 G: SupportGenerator<F, N, Id = BT::Data>, 164
213 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>; 165 /// TODO: document this
166 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
167
168 /// Returns a scaling factor for the tolerance sequence.
169 ///
170 /// Typically this is the regularisation parameter.
171 fn tolerance_scaling(&self) -> F;
172 }
173
174 /// Regularisation term that can be used with [`crate::prox_penalty::radon_squared::RadonSquared`]
175 /// proximal penalty.
176 pub trait RadonSquaredRegTerm<Domain, F = f64>: RegTerm<Domain, F>
177 where
178 Domain: Space + Clone,
179 F: Float + ToNalgebraRealField,
180 {
181 /// Adapt weights of μ, possibly insertion a new point at tolerance_violation (which should
182 /// be returned by [`RegTerm::find_tolerance_violation`].
183 fn solve_oc_radonsq<M>(
184 &self,
185 μ: &mut DiscreteMeasure<Domain, F>,
186 τv: &mut M,
187 τ: F,
188 ε: F,
189 tolerance_violation: Option<(Domain, F, bool)>,
190 config: &InsertionConfig<F>,
191 stats: &mut IterInfo<F>,
192 ) where
193 M: Mapping<Domain, Codomain = F>;
214 194
215 /// Verify that `d` is in bounds `ε` for a merge candidate `μ` 195 /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
216 /// 196 ///
217 /// This version is s used for Radon-norm squared proximal term in 197 /// This version is s used for Radon-norm squared proximal term in
218 /// [`crate::prox_penalty::radon_squared`]. 198 /// [`crate::prox_penalty::radon_squared`].
219 /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have 199 /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have
220 /// same coordinates at same agreeing indices. 200 /// same coordinates at same agreeing indices.
221 /// 201 ///
222 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser. 202 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
223 fn verify_merge_candidate_radonsq<G, BT>( 203 fn verify_merge_candidate_radonsq<M>(
224 &self, 204 &self,
225 d: &mut BTFN<F, G, BT, N>, 205 d: &mut M,
226 μ: &RNDM<F, N>, 206 μ: &DiscreteMeasure<Domain, F>,
227 τ: F, 207 τ: F,
228 ε: F, 208 ε: F,
229 config: &FBGenericConfig<F>, 209 config: &InsertionConfig<F>,
230 radon_μ: &RNDM<F, N>, 210 radon_μ: &DiscreteMeasure<Domain, F>,
231 ) -> bool 211 ) -> bool
232 where 212 where
233 BT: BTSearch<F, N, Agg = Bounds<F>>, 213 M: MinMaxMapping<Domain, F>;
234 G: SupportGenerator<F, N, Id = BT::Data>,
235 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
236
237 /// TODO: document this
238 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
239
240 /// Returns a scaling factor for the tolerance sequence.
241 ///
242 /// Typically this is the regularisation parameter.
243 fn tolerance_scaling(&self) -> F;
244 } 214 }
245 215
246 /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`]. 216 /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`].
247 pub trait SlidingRegTerm<F: Float + ToNalgebraRealField, const N: usize>: RegTerm<F, N> { 217 pub trait SlidingRegTerm<Domain, F = f64>: RegTerm<Domain, F>
218 where
219 Domain: Space + Clone,
220 F: Float + ToNalgebraRealField,
221 {
248 /// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation 222 /// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation
249 /// term, such that $-ε ≤ τw - d ≤ ε$. 223 /// term, such that $-ε ≤ τw - d ≤ ε$.
250 fn goodness<G, BT>( 224 fn goodness<M>(
251 &self, 225 &self,
252 d: &mut BTFN<F, G, BT, N>, 226 d: &mut M,
253 μ: &RNDM<F, N>, 227 μ: &DiscreteMeasure<Domain, F>,
254 y: &Loc<F, N>, 228 y: &Domain,
255 z: &Loc<F, N>, 229 z: &Domain,
256 τ: F, 230 τ: F,
257 ε: F, 231 ε: F,
258 config: &FBGenericConfig<F>, 232 config: &InsertionConfig<F>,
259 ) -> F 233 ) -> F
260 where 234 where
261 BT: BTSearch<F, N, Agg = Bounds<F>>, 235 M: MinMaxMapping<Domain, F>;
262 G: SupportGenerator<F, N, Id = BT::Data>,
263 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
264 236
265 /// Convert bound on the regulariser to a bond on the Radon norm 237 /// Convert bound on the regulariser to a bond on the Radon norm
266 fn radon_norm_bound(&self, b: F) -> F; 238 fn radon_norm_bound(&self, b: F) -> F;
267 } 239 }
268 240
269 #[replace_float_literals(F::cast_from(literal))] 241 #[replace_float_literals(F::cast_from(literal))]
270 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for NonnegRadonRegTerm<F> 242 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<Loc<N, F>, F>
271 where 243 for NonnegRadonRegTerm<F>
272 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
273 { 244 {
274 fn solve_findim( 245 fn solve_findim(
275 &self, 246 &self,
276 mA: &DMatrix<F::MixedType>, 247 mA: &DMatrix<F::MixedType>,
277 g: &DVector<F::MixedType>, 248 g: &DVector<F::MixedType>,
278 τ: F, 249 τ: F,
279 x: &mut DVector<F::MixedType>, 250 x: &mut DVector<F::MixedType>,
280 mA_normest: F, 251 mA_normest: F,
281 ε: F, 252 ε: F,
282 config: &FBGenericConfig<F>, 253 config: &InsertionConfig<F>,
283 ) -> usize { 254 ) -> usize {
284 let inner_tolerance = ε * config.inner.tolerance_mult; 255 let inner_tolerance = ε * config.inner.tolerance_mult;
285 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); 256 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
286 quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) 257 quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
287 } 258 }
288 259
289 fn solve_findim_l1squared(
290 &self,
291 y: &DVector<F::MixedType>,
292 g: &DVector<F::MixedType>,
293 τ: F,
294 x: &mut DVector<F::MixedType>,
295 ε: F,
296 config: &FBGenericConfig<F>,
297 ) -> usize {
298 let inner_tolerance = ε * config.inner.tolerance_mult;
299 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
300 l1squared_nonneg(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
301 }
302
303 #[inline] 260 #[inline]
304 fn find_tolerance_violation_slack<G, BT>( 261 fn find_tolerance_violation<M>(
305 &self, 262 &self,
306 d: &mut BTFN<F, G, BT, N>, 263 d: &mut M,
307 τ: F, 264 τ: F,
308 ε: F, 265 ε: F,
309 skip_by_rough_check: bool, 266 skip_by_rough_check: bool,
310 config: &FBGenericConfig<F>, 267 config: &InsertionConfig<F>,
311 slack: F, 268 ) -> Option<(Loc<N, F>, F, bool)>
312 ) -> Option<(Loc<F, N>, F, bool)> 269 where
313 where 270 M: MinMaxMapping<Loc<N, F>, F>,
314 BT: BTSearch<F, N, Agg = Bounds<F>>, 271 {
315 G: SupportGenerator<F, N, Id = BT::Data>, 272 let τα = τ * self.α();
316 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, 273 let keep_above = -τα - ε;
317 { 274 let minimise_below = -τα - ε * config.insertion_cutoff_factor;
318 let τα = τ * self.α();
319 let keep_above = -τα - slack - ε;
320 let minimise_below = -τα - slack - ε * config.insertion_cutoff_factor;
321 let refinement_tolerance = ε * config.refinement.tolerance_mult; 275 let refinement_tolerance = ε * config.refinement.tolerance_mult;
322 276
323 // If preliminary check indicates that we are in bounds, and if it otherwise matches 277 // If preliminary check indicates that we are in bounds, and if it otherwise matches
324 // the insertion strategy, skip insertion. 278 // the insertion strategy, skip insertion.
325 if skip_by_rough_check && d.bounds().lower() >= keep_above { 279 if skip_by_rough_check && d.bounds().lower() >= keep_above {
326 None 280 None
327 } else { 281 } else {
328 // If the rough check didn't indicate no insertion needed, find minimising point. 282 // If the rough check didn't indicate no insertion needed, find minimising point.
329 d.minimise_below( 283 let res = d.minimise_below(
330 minimise_below, 284 minimise_below,
331 refinement_tolerance, 285 refinement_tolerance,
332 config.refinement.max_steps, 286 config.refinement.max_steps,
333 ) 287 );
334 .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above)) 288
289 res.map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
335 } 290 }
336 } 291 }
337 292
338 fn verify_merge_candidate<G, BT>( 293 fn verify_merge_candidate<M>(
339 &self, 294 &self,
340 d: &mut BTFN<F, G, BT, N>, 295 d: &mut M,
341 μ: &RNDM<F, N>, 296 μ: &RNDM<N, F>,
342 τ: F, 297 τ: F,
343 ε: F, 298 ε: F,
344 config: &FBGenericConfig<F>, 299 config: &InsertionConfig<F>,
345 ) -> bool 300 ) -> bool
346 where 301 where
347 BT: BTSearch<F, N, Agg = Bounds<F>>, 302 M: MinMaxMapping<Loc<N, F>, F>,
348 G: SupportGenerator<F, N, Id = BT::Data>,
349 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
350 { 303 {
351 let τα = τ * self.α(); 304 let τα = τ * self.α();
352 let refinement_tolerance = ε * config.refinement.tolerance_mult; 305 let refinement_tolerance = ε * config.refinement.tolerance_mult;
353 let merge_tolerance = config.merge_tolerance_mult * ε; 306 let merge_tolerance = config.merge_tolerance_mult * ε;
354 let keep_above = -τα - merge_tolerance; 307 let keep_above = -τα - merge_tolerance;
365 refinement_tolerance, 318 refinement_tolerance,
366 config.refinement.max_steps, 319 config.refinement.max_steps,
367 )); 320 ));
368 } 321 }
369 322
370 fn verify_merge_candidate_radonsq<G, BT>( 323 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
371 &self, 324 let τα = τ * self.α();
372 d: &mut BTFN<F, G, BT, N>, 325 Some(Bounds(τα - ε, τα + ε))
373 μ: &RNDM<F, N>, 326 }
374 τ: F, 327
375 ε: F, 328 fn tolerance_scaling(&self) -> F {
376 config: &FBGenericConfig<F>, 329 self.α()
377 radon_μ: &RNDM<F, N>, 330 }
331 }
332
333 #[replace_float_literals(F::cast_from(literal))]
334 impl<F: Float + ToNalgebraRealField, const N: usize> RadonSquaredRegTerm<Loc<N, F>, F>
335 for NonnegRadonRegTerm<F>
336 {
337 fn solve_oc_radonsq<M>(
338 &self,
339 μ: &mut DiscreteMeasure<Loc<N, F>, F>,
340 τv: &mut M,
341 τ: F,
342 ε: F,
343 tolerance_violation: Option<(Loc<N, F>, F, bool)>,
344 config: &InsertionConfig<F>,
345 stats: &mut IterInfo<F>,
346 ) where
347 M: Mapping<Loc<N, F>, Codomain = F>,
348 {
349 let τα = τ * self.α();
350 let mut g: Vec<_> = μ
351 .iter_locations()
352 .map(|ζ| F::to_nalgebra_mixed(-τv.apply(ζ)))
353 .collect();
354
355 let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation {
356 // Don't insert if existing spikes are almost as good
357 if g.iter().all(|minus_τv| {
358 -F::from_nalgebra_mixed(*minus_τv) > v_ξ + ε * config.refinement.tolerance_mult
359 }) {
360 // Weight is found out by running the finite-dimensional optimisation algorithm
361 // above
362 // NOTE: cannot set α here before y is extracted
363 *μ += DeltaMeasure { x: ξ, α: 0.0 /*-v_ξ - τα*/ };
364 g.push(F::to_nalgebra_mixed(-v_ξ));
365 Some(-v_ξ - τα)
366 } else {
367 None
368 }
369 } else {
370 None
371 };
372
373 // Optimise weights
374 if μ.len() > 0 {
375 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
376 // from the beginning of the iteration are all contained in the immutable c and g.
377 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
378 // problems have not yet been updated to sign change.
379 let y = μ.masses_dvector();
380 let mut x = y.clone();
381 let g_na = DVector::from_vec(g);
382 if let (Some(β), Some(dest)) = (new_spike_initial_weight, x.as_mut_slice().last_mut())
383 {
384 *dest = F::to_nalgebra_mixed(β);
385 }
386 // Solve finite-dimensional subproblem.
387 let inner_tolerance = ε * config.inner.tolerance_mult;
388 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
389 stats.inner_iters +=
390 l1squared_nonneg(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it);
391
392 // Update masses of μ based on solution of finite-dimensional subproblem.
393 μ.set_masses_dvector(&x);
394 }
395 }
396
397 fn verify_merge_candidate_radonsq<M>(
398 &self,
399 d: &mut M,
400 μ: &RNDM<N, F>,
401 τ: F,
402 ε: F,
403 config: &InsertionConfig<F>,
404 radon_μ: &RNDM<N, F>,
378 ) -> bool 405 ) -> bool
379 where 406 where
380 BT: BTSearch<F, N, Agg = Bounds<F>>, 407 M: MinMaxMapping<Loc<N, F>, F>,
381 G: SupportGenerator<F, N, Id = BT::Data>,
382 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
383 { 408 {
384 let τα = τ * self.α(); 409 let τα = τ * self.α();
385 let refinement_tolerance = ε * config.refinement.tolerance_mult; 410 let refinement_tolerance = ε * config.refinement.tolerance_mult;
386 let merge_tolerance = config.merge_tolerance_mult * ε; 411 let merge_tolerance = config.merge_tolerance_mult * ε;
387 let slack = radon_μ.norm(Radon); 412 let slack = radon_μ.norm(Radon);
412 refinement_tolerance, 437 refinement_tolerance,
413 config.refinement.max_steps, 438 config.refinement.max_steps,
414 ) 439 )
415 }; 440 };
416 } 441 }
417
418 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
419 let τα = τ * self.α();
420 Some(Bounds(τα - ε, τα + ε))
421 }
422
423 fn tolerance_scaling(&self) -> F {
424 self.α()
425 }
426 } 442 }
427 443
428 #[replace_float_literals(F::cast_from(literal))] 444 #[replace_float_literals(F::cast_from(literal))]
429 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for NonnegRadonRegTerm<F> 445 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<Loc<N, F>, F>
430 where 446 for NonnegRadonRegTerm<F>
431 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
432 { 447 {
433 fn goodness<G, BT>( 448 fn goodness<M>(
434 &self, 449 &self,
435 d: &mut BTFN<F, G, BT, N>, 450 d: &mut M,
436 _μ: &RNDM<F, N>, 451 _μ: &RNDM<N, F>,
437 y: &Loc<F, N>, 452 y: &Loc<N, F>,
438 z: &Loc<F, N>, 453 z: &Loc<N, F>,
439 τ: F, 454 τ: F,
440 ε: F, 455 ε: F,
441 _config: &FBGenericConfig<F>, 456 _config: &InsertionConfig<F>,
442 ) -> F 457 ) -> F
443 where 458 where
444 BT: BTSearch<F, N, Agg = Bounds<F>>, 459 M: MinMaxMapping<Loc<N, F>, F>,
445 G: SupportGenerator<F, N, Id = BT::Data>,
446 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
447 { 460 {
448 let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α())); 461 let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α()));
449 w(z) - w(y) 462 w(z) - w(y)
450 } 463 }
451 464
453 b / self.α() 466 b / self.α()
454 } 467 }
455 } 468 }
456 469
457 #[replace_float_literals(F::cast_from(literal))] 470 #[replace_float_literals(F::cast_from(literal))]
458 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for RadonRegTerm<F> 471 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<Loc<N, F>, F> for RadonRegTerm<F> {
459 where
460 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
461 {
462 fn solve_findim( 472 fn solve_findim(
463 &self, 473 &self,
464 mA: &DMatrix<F::MixedType>, 474 mA: &DMatrix<F::MixedType>,
465 g: &DVector<F::MixedType>, 475 g: &DVector<F::MixedType>,
466 τ: F, 476 τ: F,
467 x: &mut DVector<F::MixedType>, 477 x: &mut DVector<F::MixedType>,
468 mA_normest: F, 478 mA_normest: F,
469 ε: F, 479 ε: F,
470 config: &FBGenericConfig<F>, 480 config: &InsertionConfig<F>,
471 ) -> usize { 481 ) -> usize {
472 let inner_tolerance = ε * config.inner.tolerance_mult; 482 let inner_tolerance = ε * config.inner.tolerance_mult;
473 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); 483 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
474 quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) 484 quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
475 } 485 }
476 486
477 fn solve_findim_l1squared( 487 fn find_tolerance_violation<M>(
478 &self, 488 &self,
479 y: &DVector<F::MixedType>, 489 d: &mut M,
480 g: &DVector<F::MixedType>,
481 τ: F,
482 x: &mut DVector<F::MixedType>,
483 ε: F,
484 config: &FBGenericConfig<F>,
485 ) -> usize {
486 let inner_tolerance = ε * config.inner.tolerance_mult;
487 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
488 l1squared_unconstrained(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
489 }
490
491 fn find_tolerance_violation_slack<G, BT>(
492 &self,
493 d: &mut BTFN<F, G, BT, N>,
494 τ: F, 490 τ: F,
495 ε: F, 491 ε: F,
496 skip_by_rough_check: bool, 492 skip_by_rough_check: bool,
497 config: &FBGenericConfig<F>, 493 config: &InsertionConfig<F>,
498 slack: F, 494 ) -> Option<(Loc<N, F>, F, bool)>
499 ) -> Option<(Loc<F, N>, F, bool)> 495 where
500 where 496 M: MinMaxMapping<Loc<N, F>, F>,
501 BT: BTSearch<F, N, Agg = Bounds<F>>, 497 {
502 G: SupportGenerator<F, N, Id = BT::Data>, 498 let τα = τ * self.α();
503 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, 499 let keep_below = τα + ε;
504 { 500 let keep_above = -τα - ε;
505 let τα = τ * self.α(); 501 let maximise_above = τα + ε * config.insertion_cutoff_factor;
506 let keep_below = τα + slack + ε; 502 let minimise_below = -τα - ε * config.insertion_cutoff_factor;
507 let keep_above = -(τα + slack) - ε;
508 let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
509 let minimise_below = -(τα + slack) - ε * config.insertion_cutoff_factor;
510 let refinement_tolerance = ε * config.refinement.tolerance_mult; 503 let refinement_tolerance = ε * config.refinement.tolerance_mult;
511 504
512 // If preliminary check indicates that we are in bonds, and if it otherwise matches 505 // If preliminary check indicates that we are in bonds, and if it otherwise matches
513 // the insertion strategy, skip insertion. 506 // the insertion strategy, skip insertion.
514 if skip_by_rough_check && Bounds(keep_above, keep_below).superset(&d.bounds()) { 507 if skip_by_rough_check && Bounds(keep_above, keep_below).superset(&d.bounds()) {
539 } 532 }
540 } 533 }
541 } 534 }
542 } 535 }
543 536
544 fn verify_merge_candidate<G, BT>( 537 fn verify_merge_candidate<M>(
545 &self, 538 &self,
546 d: &mut BTFN<F, G, BT, N>, 539 d: &mut M,
547 μ: &RNDM<F, N>, 540 μ: &RNDM<N, F>,
548 τ: F, 541 τ: F,
549 ε: F, 542 ε: F,
550 config: &FBGenericConfig<F>, 543 config: &InsertionConfig<F>,
551 ) -> bool 544 ) -> bool
552 where 545 where
553 BT: BTSearch<F, N, Agg = Bounds<F>>, 546 M: MinMaxMapping<Loc<N, F>, F>,
554 G: SupportGenerator<F, N, Id = BT::Data>,
555 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
556 { 547 {
557 let τα = τ * self.α(); 548 let τα = τ * self.α();
558 let refinement_tolerance = ε * config.refinement.tolerance_mult; 549 let refinement_tolerance = ε * config.refinement.tolerance_mult;
559 let merge_tolerance = config.merge_tolerance_mult * ε; 550 let merge_tolerance = config.merge_tolerance_mult * ε;
560 let keep_below = τα + merge_tolerance; 551 let keep_below = τα + merge_tolerance;
583 refinement_tolerance, 574 refinement_tolerance,
584 config.refinement.max_steps, 575 config.refinement.max_steps,
585 )); 576 ));
586 } 577 }
587 578
588 fn verify_merge_candidate_radonsq<G, BT>( 579 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
589 &self, 580 let τα = τ * self.α();
590 d: &mut BTFN<F, G, BT, N>, 581 Some(Bounds(-τα - ε, τα + ε))
591 μ: &RNDM<F, N>, 582 }
592 τ: F, 583
593 ε: F, 584 fn tolerance_scaling(&self) -> F {
594 config: &FBGenericConfig<F>, 585 self.α()
595 radon_μ: &RNDM<F, N>, 586 }
587 }
588
589 #[replace_float_literals(F::cast_from(literal))]
590 impl<F: Float + ToNalgebraRealField, const N: usize> RadonSquaredRegTerm<Loc<N, F>, F>
591 for RadonRegTerm<F>
592 {
593 fn solve_oc_radonsq<M>(
594 &self,
595 μ: &mut DiscreteMeasure<Loc<N, F>, F>,
596 τv: &mut M,
597 τ: F,
598 ε: F,
599 tolerance_violation: Option<(Loc<N, F>, F, bool)>,
600 config: &InsertionConfig<F>,
601 stats: &mut IterInfo<F>,
602 ) where
603 M: Mapping<Loc<N, F>, Codomain = F>,
604 {
605 let τα = τ * self.α();
606 let mut g: Vec<_> = μ
607 .iter_locations()
608 .map(|ζ| F::to_nalgebra_mixed(τv.apply(-ζ)))
609 .collect();
610
611 let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation {
612 // Don't insert if existing spikes are almost as good
613 let n = v_ξ.abs();
614 if g.iter().all(|minus_τv| {
615 F::from_nalgebra_mixed(*minus_τv).abs() < n - ε * config.refinement.tolerance_mult
616 }) {
617 // Weight is found out by running the finite-dimensional optimisation algorithm
618 // above
619 // NOTE: cannot initialise α before y is extracted.
620 *μ += DeltaMeasure { x: ξ, α: 0.0 /*-(n + τα) * v_ξ.signum()*/ };
621 g.push(F::to_nalgebra_mixed(-v_ξ));
622 Some(-(n + τα) * v_ξ.signum())
623 } else {
624 None
625 }
626 } else {
627 None
628 };
629
630 // Optimise weights
631 if μ.len() > 0 {
632 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
633 // from the beginning of the iteration are all contained in the immutable c and g.
634 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
635 // problems have not yet been updated to sign change.
636 let y = μ.masses_dvector();
637 let mut x = y.clone();
638 if let (Some(β), Some(dest)) = (new_spike_initial_weight, x.as_mut_slice().last_mut())
639 {
640 *dest = F::to_nalgebra_mixed(β);
641 }
642 let g_na = DVector::from_vec(g);
643 // Solve finite-dimensional subproblem.
644 let inner_tolerance = ε * config.inner.tolerance_mult;
645 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
646 stats.inner_iters +=
647 l1squared_unconstrained(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it);
648
649 // Update masses of μ based on solution of finite-dimensional subproblem.
650 μ.set_masses_dvector(&x);
651 }
652 }
653
654 fn verify_merge_candidate_radonsq<M>(
655 &self,
656 d: &mut M,
657 μ: &RNDM<N, F>,
658 τ: F,
659 ε: F,
660 config: &InsertionConfig<F>,
661 radon_μ: &RNDM<N, F>,
596 ) -> bool 662 ) -> bool
597 where 663 where
598 BT: BTSearch<F, N, Agg = Bounds<F>>, 664 M: MinMaxMapping<Loc<N, F>, F>,
599 G: SupportGenerator<F, N, Id = BT::Data>,
600 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
601 { 665 {
602 let τα = τ * self.α(); 666 let τα = τ * self.α();
603 let refinement_tolerance = ε * config.refinement.tolerance_mult; 667 let refinement_tolerance = ε * config.refinement.tolerance_mult;
604 let merge_tolerance = config.merge_tolerance_mult * ε; 668 let merge_tolerance = config.merge_tolerance_mult * ε;
605 let slack = radon_μ.norm(Radon); 669 let slack = radon_μ.norm(Radon);
636 refinement_tolerance, 700 refinement_tolerance,
637 config.refinement.max_steps, 701 config.refinement.max_steps,
638 ) 702 )
639 }; 703 };
640 } 704 }
641
642 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
643 let τα = τ * self.α();
644 Some(Bounds(-τα - ε, τα + ε))
645 }
646
647 fn tolerance_scaling(&self) -> F {
648 self.α()
649 }
650 } 705 }
651 706
652 #[replace_float_literals(F::cast_from(literal))] 707 #[replace_float_literals(F::cast_from(literal))]
653 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for RadonRegTerm<F> 708 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<Loc<N, F>, F>
654 where 709 for RadonRegTerm<F>
655 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
656 { 710 {
657 fn goodness<G, BT>( 711 fn goodness<M>(
658 &self, 712 &self,
659 d: &mut BTFN<F, G, BT, N>, 713 d: &mut M,
660 _μ: &RNDM<F, N>, 714 _μ: &RNDM<N, F>,
661 y: &Loc<F, N>, 715 y: &Loc<N, F>,
662 z: &Loc<F, N>, 716 z: &Loc<N, F>,
663 τ: F, 717 τ: F,
664 ε: F, 718 ε: F,
665 _config: &FBGenericConfig<F>, 719 _config: &InsertionConfig<F>,
666 ) -> F 720 ) -> F
667 where 721 where
668 BT: BTSearch<F, N, Agg = Bounds<F>>, 722 M: MinMaxMapping<Loc<N, F>, F>,
669 G: SupportGenerator<F, N, Id = BT::Data>,
670 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
671 { 723 {
672 let α = self.α(); 724 let α = self.α();
673 let w = |x| { 725 let w = |x| {
674 let dx = d.apply(x); 726 let dx = d.apply(x);
675 ((-ε + dx) / (τ * α)).max(1.0.min(ε + dx) / (τ * α)) 727 ((-ε + dx) / (τ * α)).max(1.0.min(ε + dx) / (τ * α))

Mercurial Repository: pointsource_algs

mercurial