Mercurial > repos > pointsource_algs / file comparison

comparison: src/regularisation.rs

src/regularisation.rs

changeset 52
f0e8704d3f0e
parent 51
0693cc9ba9f0
equal deleted inserted replaced
31:6105b5cd8d89 52:f0e8704d3f0e
1 /*! 1 /*!
2 Regularisation terms 2 Regularisation terms
3 */ 3 */
4 4
5 use serde::{Serialize, Deserialize}; 5 #[allow(unused_imports)] // Used by documentation.
6 use crate::fb::pointsource_fb_reg;
7 use crate::fb::FBGenericConfig;
8 use crate::measures::{DeltaMeasure, Radon, RNDM};
9 #[allow(unused_imports)] // Used by documentation.
10 use crate::sliding_fb::pointsource_sliding_fb_reg;
11 use crate::types::*;
12 use alg_tools::instance::Instance;
13 use alg_tools::linops::Mapping;
14 use alg_tools::loc::Loc;
6 use alg_tools::norms::Norm; 15 use alg_tools::norms::Norm;
7 use alg_tools::linops::Apply; 16 use numeric_literals::replace_float_literals;
8 use alg_tools::loc::Loc; 17 use serde::{Deserialize, Serialize};
9 use crate::types::*; 18
10 use crate::measures::{ 19 use crate::subproblem::{
11 DiscreteMeasure, 20 l1squared_nonneg::l1squared_nonneg, l1squared_unconstrained::l1squared_unconstrained,
12 Radon 21 nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained,
13 }; 22 };
14 #[allow(unused_imports)] // Used by documentation. 23 use alg_tools::bisection_tree::{
15 use crate::fb::generic_pointsource_fb_reg; 24 BTSearch, Bounded, Bounds, LocalAnalysis, P2Minimise, SupportGenerator, BTFN,
16 25 };
17 /// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`generic_pointsource_fb_reg`]. 26 use alg_tools::iterate::AlgIteratorFactory;
27 use alg_tools::nalgebra_support::ToNalgebraRealField;
28 use nalgebra::{DMatrix, DVector};
29
30 use std::cmp::Ordering::{Equal, Greater, Less};
31
32 /// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other
33 /// algorithms.
18 /// 34 ///
19 /// The only member of the struct is the regularisation parameter α. 35 /// The only member of the struct is the regularisation parameter α.
20 #[derive(Copy, Clone, Debug, Serialize, Deserialize)] 36 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
21 pub struct NonnegRadonRegTerm<F : Float>(pub F /* α */); 37 pub struct NonnegRadonRegTerm<F: Float>(pub F /* α */);
22 38
23 impl<'a, F : Float> NonnegRadonRegTerm<F> { 39 impl<'a, F: Float> NonnegRadonRegTerm<F> {
24 /// Returns the regularisation parameter 40 /// Returns the regularisation parameter
25 pub fn α(&self) -> F { 41 pub fn α(&self) -> F {
26 let &NonnegRadonRegTerm(α) = self; 42 let &NonnegRadonRegTerm(α) = self;
27 α 43 α
28 } 44 }
29 } 45 }
30 46
31 impl<'a, F : Float, const N : usize> Apply<&'a DiscreteMeasure<Loc<F, N>, F>> 47 impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for NonnegRadonRegTerm<F> {
32 for NonnegRadonRegTerm<F> { 48 type Codomain = F;
33 type Output = F; 49
34 50 fn apply<I>(&self, μ: I) -> F
35 fn apply(&self, μ : &'a DiscreteMeasure<Loc<F, N>, F>) -> F { 51 where
36 self.α() * μ.norm(Radon) 52 I: Instance<RNDM<F, N>>,
37 } 53 {
38 } 54 self.α() * μ.eval(|x| x.norm(Radon))
39 55 }
40 56 }
41 /// The regularisation term $α\|μ\|_{ℳ(Ω)}$ for [`generic_pointsource_fb_reg`]. 57
58 /// The regularisation term $α\|μ\|_{ℳ(Ω)}$ for [`pointsource_fb_reg`].
42 /// 59 ///
43 /// The only member of the struct is the regularisation parameter α. 60 /// The only member of the struct is the regularisation parameter α.
44 #[derive(Copy, Clone, Debug, Serialize, Deserialize)] 61 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
45 pub struct RadonRegTerm<F : Float>(pub F /* α */); 62 pub struct RadonRegTerm<F: Float>(pub F /* α */);
46 63
47 64 impl<'a, F: Float> RadonRegTerm<F> {
48 impl<'a, F : Float> RadonRegTerm<F> {
49 /// Returns the regularisation parameter 65 /// Returns the regularisation parameter
50 pub fn α(&self) -> F { 66 pub fn α(&self) -> F {
51 let &RadonRegTerm(α) = self; 67 let &RadonRegTerm(α) = self;
52 α 68 α
53 } 69 }
54 } 70 }
55 71
56 impl<'a, F : Float, const N : usize> Apply<&'a DiscreteMeasure<Loc<F, N>, F>> 72 impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for RadonRegTerm<F> {
57 for RadonRegTerm<F> { 73 type Codomain = F;
58 type Output = F; 74
59 75 fn apply<I>(&self, μ: I) -> F
60 fn apply(&self, μ : &'a DiscreteMeasure<Loc<F, N>, F>) -> F { 76 where
61 self.α() * μ.norm(Radon) 77 I: Instance<RNDM<F, N>>,
78 {
79 self.α() * μ.eval(|x| x.norm(Radon))
62 } 80 }
63 } 81 }
64 82
65 /// Regularisation term configuration 83 /// Regularisation term configuration
66 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] 84 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
67 pub enum Regularisation<F : Float> { 85 pub enum Regularisation<F: Float> {
68 /// $α \\|μ\\|\_{ℳ(Ω)}$ 86 /// $α \\|μ\\|\_{ℳ(Ω)}$
69 Radon(F), 87 Radon(F),
70 /// $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ 88 /// $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$
71 NonnegRadon(F), 89 NonnegRadon(F),
72 } 90 }
73 91
74 impl<'a, F : Float, const N : usize> Apply<&'a DiscreteMeasure<Loc<F, N>, F>> 92 impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for Regularisation<F> {
75 for Regularisation<F> { 93 type Codomain = F;
76 type Output = F; 94
77 95 fn apply<I>(&self, μ: I) -> F
78 fn apply(&self, μ : &'a DiscreteMeasure<Loc<F, N>, F>) -> F { 96 where
97 I: Instance<RNDM<F, N>>,
98 {
79 match *self { 99 match *self {
80 Self::Radon(α) => RadonRegTerm(α).apply(μ), 100 Self::Radon(α) => RadonRegTerm(α).apply(μ),
81 Self::NonnegRadon(α) => NonnegRadonRegTerm(α).apply(μ), 101 Self::NonnegRadon(α) => NonnegRadonRegTerm(α).apply(μ),
82 } 102 }
83 } 103 }
84 } 104 }
105
106 /// Abstraction of regularisation terms.
107 pub trait RegTerm<F: Float + ToNalgebraRealField, const N: usize>:
108 Mapping<RNDM<F, N>, Codomain = F>
109 {
110 /// Approximately solve the problem
111 /// <div>$$
112 /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + τ G(x)
113 /// $$</div>
114 /// for $G$ depending on the trait implementation.
115 ///
116 /// The parameter `mA` is $A$. An estimate for its opeator norm should be provided in
117 /// `mA_normest`. The initial iterate and output is `x`. The current main tolerance is `ε`.
118 ///
119 /// Returns the number of iterations taken.
120 fn solve_findim(
121 &self,
122 mA: &DMatrix<F::MixedType>,
123 g: &DVector<F::MixedType>,
124 τ: F,
125 x: &mut DVector<F::MixedType>,
126 mA_normest: F,
127 ε: F,
128 config: &FBGenericConfig<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;
147
148 /// Find a point where `d` may violate the tolerance `ε`.
149 ///
150 /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
151 /// are in bounds. `ε` is the current main tolerance and `τ` a scaling factor for the
152 /// regulariser.
153 ///
154 /// 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,
156 /// and a boolean indicating whether the found point is in bounds.
157 fn find_tolerance_violation<G, BT>(
158 &self,
159 d: &mut BTFN<F, G, BT, N>,
160 τ: F,
161 ε: F,
162 skip_by_rough_check: bool,
163 config: &FBGenericConfig<F>,
164 ) -> Option<(Loc<F, N>, F, bool)>
165 where
166 BT: BTSearch<F, N, Agg = Bounds<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
199 /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
200 ///
201 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
202 fn verify_merge_candidate<G, BT>(
203 &self,
204 d: &mut BTFN<F, G, BT, N>,
205 μ: &RNDM<F, N>,
206 τ: F,
207 ε: F,
208 config: &FBGenericConfig<F>,
209 ) -> bool
210 where
211 BT: BTSearch<F, N, Agg = Bounds<F>>,
212 G: SupportGenerator<F, N, Id = BT::Data>,
213 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
214
215 /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
216 ///
217 /// This version is s used for Radon-norm squared proximal term in
218 /// [`crate::prox_penalty::radon_squared`].
219 /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have
220 /// same coordinates at same agreeing indices.
221 ///
222 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
223 fn verify_merge_candidate_radonsq<G, BT>(
224 &self,
225 d: &mut BTFN<F, G, BT, N>,
226 μ: &RNDM<F, N>,
227 τ: F,
228 ε: F,
229 config: &FBGenericConfig<F>,
230 radon_μ: &RNDM<F, N>,
231 ) -> bool
232 where
233 BT: BTSearch<F, N, Agg = Bounds<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 }
245
246 /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`].
247 pub trait SlidingRegTerm<F: Float + ToNalgebraRealField, const N: usize>: RegTerm<F, N> {
248 /// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation
249 /// term, such that $-ε ≤ τw - d ≤ ε$.
250 fn goodness<G, BT>(
251 &self,
252 d: &mut BTFN<F, G, BT, N>,
253 μ: &RNDM<F, N>,
254 y: &Loc<F, N>,
255 z: &Loc<F, N>,
256 τ: F,
257 ε: F,
258 config: &FBGenericConfig<F>,
259 ) -> F
260 where
261 BT: BTSearch<F, N, Agg = Bounds<F>>,
262 G: SupportGenerator<F, N, Id = BT::Data>,
263 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
264
265 /// Convert bound on the regulariser to a bond on the Radon norm
266 fn radon_norm_bound(&self, b: F) -> F;
267 }
268
269 #[replace_float_literals(F::cast_from(literal))]
270 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for NonnegRadonRegTerm<F>
271 where
272 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
273 {
274 fn solve_findim(
275 &self,
276 mA: &DMatrix<F::MixedType>,
277 g: &DVector<F::MixedType>,
278 τ: F,
279 x: &mut DVector<F::MixedType>,
280 mA_normest: F,
281 ε: F,
282 config: &FBGenericConfig<F>,
283 ) -> usize {
284 let inner_tolerance = ε * config.inner.tolerance_mult;
285 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
286 quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
287 }
288
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]
304 fn find_tolerance_violation_slack<G, BT>(
305 &self,
306 d: &mut BTFN<F, G, BT, N>,
307 τ: F,
308 ε: F,
309 skip_by_rough_check: bool,
310 config: &FBGenericConfig<F>,
311 slack: F,
312 ) -> Option<(Loc<F, N>, F, bool)>
313 where
314 BT: BTSearch<F, N, Agg = Bounds<F>>,
315 G: SupportGenerator<F, N, Id = BT::Data>,
316 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
317 {
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;
322
323 // If preliminary check indicates that we are in bounds, and if it otherwise matches
324 // the insertion strategy, skip insertion.
325 if skip_by_rough_check && d.bounds().lower() >= keep_above {
326 None
327 } else {
328 // If the rough check didn't indicate no insertion needed, find minimising point.
329 d.minimise_below(
330 minimise_below,
331 refinement_tolerance,
332 config.refinement.max_steps,
333 )
334 .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
335 }
336 }
337
338 fn verify_merge_candidate<G, BT>(
339 &self,
340 d: &mut BTFN<F, G, BT, N>,
341 μ: &RNDM<F, N>,
342 τ: F,
343 ε: F,
344 config: &FBGenericConfig<F>,
345 ) -> bool
346 where
347 BT: BTSearch<F, N, Agg = Bounds<F>>,
348 G: SupportGenerator<F, N, Id = BT::Data>,
349 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
350 {
351 let τα = τ * self.α();
352 let refinement_tolerance = ε * config.refinement.tolerance_mult;
353 let merge_tolerance = config.merge_tolerance_mult * ε;
354 let keep_above = -τα - merge_tolerance;
355 let keep_supp_below = -τα + merge_tolerance;
356 let bnd = d.bounds();
357
358 return (bnd.upper() <= keep_supp_below
359 || μ
360 .iter_spikes()
361 .all(|&DeltaMeasure { α, ref x }| (α == 0.0) || d.apply(x) <= keep_supp_below))
362 && (bnd.lower() >= keep_above
363 || d.has_lower_bound(
364 keep_above,
365 refinement_tolerance,
366 config.refinement.max_steps,
367 ));
368 }
369
370 fn verify_merge_candidate_radonsq<G, BT>(
371 &self,
372 d: &mut BTFN<F, G, BT, N>,
373 μ: &RNDM<F, N>,
374 τ: F,
375 ε: F,
376 config: &FBGenericConfig<F>,
377 radon_μ: &RNDM<F, N>,
378 ) -> bool
379 where
380 BT: BTSearch<F, N, Agg = Bounds<F>>,
381 G: SupportGenerator<F, N, Id = BT::Data>,
382 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
383 {
384 let τα = τ * self.α();
385 let refinement_tolerance = ε * config.refinement.tolerance_mult;
386 let merge_tolerance = config.merge_tolerance_mult * ε;
387 let slack = radon_μ.norm(Radon);
388 let bnd = d.bounds();
389
390 return {
391 μ.both_matching(radon_μ).all(|(α, rα, x)| {
392 let v = -d.apply(x); // TODO: observe ad hoc negation here, after minus_τv
393 // switch to τv.
394 let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
395 Greater => (τα, τα),
396 _ => (F::NEG_INFINITY, τα),
397 // Less should not happen; treated as Equal
398 };
399 let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) {
400 Greater => (slack, slack),
401 Equal => (-slack, slack),
402 Less => (-slack, -slack),
403 };
404 // TODO: both fail.
405 (l1 + l2 - merge_tolerance <= v) && (v <= u1 + u2 + merge_tolerance)
406 })
407 } && {
408 let keep_above = -τα - slack - merge_tolerance;
409 bnd.lower() <= keep_above
410 || d.has_lower_bound(
411 keep_above,
412 refinement_tolerance,
413 config.refinement.max_steps,
414 )
415 };
416 }
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 }
427
428 #[replace_float_literals(F::cast_from(literal))]
429 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for NonnegRadonRegTerm<F>
430 where
431 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
432 {
433 fn goodness<G, BT>(
434 &self,
435 d: &mut BTFN<F, G, BT, N>,
436 _μ: &RNDM<F, N>,
437 y: &Loc<F, N>,
438 z: &Loc<F, N>,
439 τ: F,
440 ε: F,
441 _config: &FBGenericConfig<F>,
442 ) -> F
443 where
444 BT: BTSearch<F, N, Agg = Bounds<F>>,
445 G: SupportGenerator<F, N, Id = BT::Data>,
446 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
447 {
448 let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α()));
449 w(z) - w(y)
450 }
451
452 fn radon_norm_bound(&self, b: F) -> F {
453 b / self.α()
454 }
455 }
456
457 #[replace_float_literals(F::cast_from(literal))]
458 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for RadonRegTerm<F>
459 where
460 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
461 {
462 fn solve_findim(
463 &self,
464 mA: &DMatrix<F::MixedType>,
465 g: &DVector<F::MixedType>,
466 τ: F,
467 x: &mut DVector<F::MixedType>,
468 mA_normest: F,
469 ε: F,
470 config: &FBGenericConfig<F>,
471 ) -> usize {
472 let inner_tolerance = ε * config.inner.tolerance_mult;
473 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
474 quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
475 }
476
477 fn solve_findim_l1squared(
478 &self,
479 y: &DVector<F::MixedType>,
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,
495 ε: F,
496 skip_by_rough_check: bool,
497 config: &FBGenericConfig<F>,
498 slack: F,
499 ) -> Option<(Loc<F, N>, F, bool)>
500 where
501 BT: BTSearch<F, N, Agg = Bounds<F>>,
502 G: SupportGenerator<F, N, Id = BT::Data>,
503 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
504 {
505 let τα = τ * self.α();
506 let keep_below = τα + slack + ε;
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;
511
512 // If preliminary check indicates that we are in bonds, and if it otherwise matches
513 // the insertion strategy, skip insertion.
514 if skip_by_rough_check && Bounds(keep_above, keep_below).superset(&d.bounds()) {
515 None
516 } else {
517 // If the rough check didn't indicate no insertion needed, find maximising point.
518 let mx = d.maximise_above(
519 maximise_above,
520 refinement_tolerance,
521 config.refinement.max_steps,
522 );
523 let mi = d.minimise_below(
524 minimise_below,
525 refinement_tolerance,
526 config.refinement.max_steps,
527 );
528
529 match (mx, mi) {
530 (None, None) => None,
531 (Some((ξ, v_ξ)), None) => Some((ξ, v_ξ, keep_below >= v_ξ)),
532 (None, Some((ζ, v_ζ))) => Some((ζ, v_ζ, keep_above <= v_ζ)),
533 (Some((ξ, v_ξ)), Some((ζ, v_ζ))) => {
534 if v_ξ - τα > τα - v_ζ {
535 Some((ξ, v_ξ, keep_below >= v_ξ))
536 } else {
537 Some((ζ, v_ζ, keep_above <= v_ζ))
538 }
539 }
540 }
541 }
542 }
543
544 fn verify_merge_candidate<G, BT>(
545 &self,
546 d: &mut BTFN<F, G, BT, N>,
547 μ: &RNDM<F, N>,
548 τ: F,
549 ε: F,
550 config: &FBGenericConfig<F>,
551 ) -> bool
552 where
553 BT: BTSearch<F, N, Agg = Bounds<F>>,
554 G: SupportGenerator<F, N, Id = BT::Data>,
555 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
556 {
557 let τα = τ * self.α();
558 let refinement_tolerance = ε * config.refinement.tolerance_mult;
559 let merge_tolerance = config.merge_tolerance_mult * ε;
560 let keep_below = τα + merge_tolerance;
561 let keep_above = -τα - merge_tolerance;
562 let keep_supp_pos_above = τα - merge_tolerance;
563 let keep_supp_neg_below = -τα + merge_tolerance;
564 let bnd = d.bounds();
565
566 return ((bnd.lower() >= keep_supp_pos_above && bnd.upper() <= keep_supp_neg_below)
567 || μ
568 .iter_spikes()
569 .all(|&DeltaMeasure { α: β, ref x }| match β.partial_cmp(&0.0) {
570 Some(Greater) => d.apply(x) >= keep_supp_pos_above,
571 Some(Less) => d.apply(x) <= keep_supp_neg_below,
572 _ => true,
573 }))
574 && (bnd.upper() <= keep_below
575 || d.has_upper_bound(
576 keep_below,
577 refinement_tolerance,
578 config.refinement.max_steps,
579 ))
580 && (bnd.lower() >= keep_above
581 || d.has_lower_bound(
582 keep_above,
583 refinement_tolerance,
584 config.refinement.max_steps,
585 ));
586 }
587
588 fn verify_merge_candidate_radonsq<G, BT>(
589 &self,
590 d: &mut BTFN<F, G, BT, N>,
591 μ: &RNDM<F, N>,
592 τ: F,
593 ε: F,
594 config: &FBGenericConfig<F>,
595 radon_μ: &RNDM<F, N>,
596 ) -> bool
597 where
598 BT: BTSearch<F, N, Agg = Bounds<F>>,
599 G: SupportGenerator<F, N, Id = BT::Data>,
600 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
601 {
602 let τα = τ * self.α();
603 let refinement_tolerance = ε * config.refinement.tolerance_mult;
604 let merge_tolerance = config.merge_tolerance_mult * ε;
605 let slack = radon_μ.norm(Radon);
606 let bnd = d.bounds();
607
608 return {
609 μ.both_matching(radon_μ).all(|(α, rα, x)| {
610 let v = d.apply(x);
611 let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
612 Greater => (τα, τα),
613 Equal => (-τα, τα),
614 Less => (-τα, -τα),
615 };
616 let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) {
617 Greater => (slack, slack),
618 Equal => (-slack, slack),
619 Less => (-slack, -slack),
620 };
621 (l1 + l2 - merge_tolerance <= v) && (v <= u1 + u2 + merge_tolerance)
622 })
623 } && {
624 let keep_below = τα + slack + merge_tolerance;
625 bnd.upper() <= keep_below
626 || d.has_upper_bound(
627 keep_below,
628 refinement_tolerance,
629 config.refinement.max_steps,
630 )
631 } && {
632 let keep_above = -τα - slack - merge_tolerance;
633 bnd.lower() >= keep_above
634 || d.has_lower_bound(
635 keep_above,
636 refinement_tolerance,
637 config.refinement.max_steps,
638 )
639 };
640 }
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 }
651
652 #[replace_float_literals(F::cast_from(literal))]
653 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for RadonRegTerm<F>
654 where
655 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
656 {
657 fn goodness<G, BT>(
658 &self,
659 d: &mut BTFN<F, G, BT, N>,
660 _μ: &RNDM<F, N>,
661 y: &Loc<F, N>,
662 z: &Loc<F, N>,
663 τ: F,
664 ε: F,
665 _config: &FBGenericConfig<F>,
666 ) -> F
667 where
668 BT: BTSearch<F, N, Agg = Bounds<F>>,
669 G: SupportGenerator<F, N, Id = BT::Data>,
670 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
671 {
672 let α = self.α();
673 let w = |x| {
674 let dx = d.apply(x);
675 ((-ε + dx) / (τ * α)).max(1.0.min(ε + dx) / (τ * α))
676 };
677 w(z) - w(y)
678 }
679
680 fn radon_norm_bound(&self, b: F) -> F {
681 b / self.α()
682 }
683 }

Mercurial Repository: pointsource_algs

mercurial