Mercurial > repos > pointsource_algs / file comparison

comparison: src/regularisation.rs

src/regularisation.rs

branch
dev
changeset 61
4f468d35fa29
parent 51
0693cc9ba9f0
child 63
7a8a55fd41c0
equal deleted inserted replaced
60:9738b51d90d7 61:4f468d35fa29
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; 129 ) -> usize;
130 130
131 /// Approximately solve the problem 131 /// Approximately solve the problem
132 /// <div>$$ 132 /// <div>$$
133 /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x) 133 /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x)
140 y: &DVector<F::MixedType>, 140 y: &DVector<F::MixedType>,
141 g: &DVector<F::MixedType>, 141 g: &DVector<F::MixedType>,
142 τ: F, 142 τ: F,
143 x: &mut DVector<F::MixedType>, 143 x: &mut DVector<F::MixedType>,
144 ε: F, 144 ε: F,
145 config: &FBGenericConfig<F>, 145 config: &InsertionConfig<F>,
146 ) -> usize; 146 ) -> usize;
147 147
148 /// Find a point where `d` may violate the tolerance `ε`. 148 /// Find a point where `d` may violate the tolerance `ε`.
149 /// 149 ///
150 /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we 150 /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
152 /// regulariser. 152 /// regulariser.
153 /// 153 ///
154 /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check 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, 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. 156 /// and a boolean indicating whether the found point is in bounds.
157 fn find_tolerance_violation<G, BT>( 157 fn find_tolerance_violation<M>(
158 &self, 158 &self,
159 d: &mut BTFN<F, G, BT, N>, 159 d: &mut M,
160 τ: F, 160 τ: F,
161 ε: F, 161 ε: F,
162 skip_by_rough_check: bool, 162 skip_by_rough_check: bool,
163 config: &FBGenericConfig<F>, 163 config: &InsertionConfig<F>,
164 ) -> Option<(Loc<F, N>, F, bool)> 164 ) -> Option<(Domain, F, bool)>
165 where 165 where
166 BT: BTSearch<F, N, Agg = Bounds<F>>, 166 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 { 167 {
170 self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO) 168 self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO)
171 } 169 }
172 170
173 /// Find a point where `d` may violate the tolerance `ε`. 171 /// Find a point where `d` may violate the tolerance `ε`.
180 /// regulariser. 178 /// regulariser.
181 /// 179 ///
182 /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check 180 /// 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, 181 /// 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. 182 /// and a boolean indicating whether the found point is in bounds.
185 fn find_tolerance_violation_slack<G, BT>( 183 fn find_tolerance_violation_slack<M>(
186 &self, 184 &self,
187 d: &mut BTFN<F, G, BT, N>, 185 d: &mut M,
188 τ: F, 186 τ: F,
189 ε: F, 187 ε: F,
190 skip_by_rough_check: bool, 188 skip_by_rough_check: bool,
191 config: &FBGenericConfig<F>, 189 config: &InsertionConfig<F>,
192 slack: F, 190 slack: F,
193 ) -> Option<(Loc<F, N>, F, bool)> 191 ) -> Option<(Domain, F, bool)>
194 where 192 where
195 BT: BTSearch<F, N, Agg = Bounds<F>>, 193 M: MinMaxMapping<Domain, F>;
196 G: SupportGenerator<F, N, Id = BT::Data>,
197 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
198 194
199 /// Verify that `d` is in bounds `ε` for a merge candidate `μ` 195 /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
200 /// 196 ///
201 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser. 197 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
202 fn verify_merge_candidate<G, BT>( 198 fn verify_merge_candidate<M>(
203 &self, 199 &self,
204 d: &mut BTFN<F, G, BT, N>, 200 d: &mut M,
205 μ: &RNDM<F, N>, 201 μ: &DiscreteMeasure<Domain, F>,
206 τ: F, 202 τ: F,
207 ε: F, 203 ε: F,
208 config: &FBGenericConfig<F>, 204 config: &InsertionConfig<F>,
209 ) -> bool 205 ) -> bool
210 where 206 where
211 BT: BTSearch<F, N, Agg = Bounds<F>>, 207 M: MinMaxMapping<Domain, F>;
212 G: SupportGenerator<F, N, Id = BT::Data>,
213 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
214 208
215 /// Verify that `d` is in bounds `ε` for a merge candidate `μ` 209 /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
216 /// 210 ///
217 /// This version is s used for Radon-norm squared proximal term in 211 /// This version is s used for Radon-norm squared proximal term in
218 /// [`crate::prox_penalty::radon_squared`]. 212 /// [`crate::prox_penalty::radon_squared`].
219 /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have 213 /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have
220 /// same coordinates at same agreeing indices. 214 /// same coordinates at same agreeing indices.
221 /// 215 ///
222 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser. 216 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
223 fn verify_merge_candidate_radonsq<G, BT>( 217 fn verify_merge_candidate_radonsq<M>(
224 &self, 218 &self,
225 d: &mut BTFN<F, G, BT, N>, 219 d: &mut M,
226 μ: &RNDM<F, N>, 220 μ: &DiscreteMeasure<Domain, F>,
227 τ: F, 221 τ: F,
228 ε: F, 222 ε: F,
229 config: &FBGenericConfig<F>, 223 config: &InsertionConfig<F>,
230 radon_μ: &RNDM<F, N>, 224 radon_μ: &DiscreteMeasure<Domain, F>,
231 ) -> bool 225 ) -> bool
232 where 226 where
233 BT: BTSearch<F, N, Agg = Bounds<F>>, 227 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 228
237 /// TODO: document this 229 /// TODO: document this
238 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>; 230 fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
239 231
240 /// Returns a scaling factor for the tolerance sequence. 232 /// Returns a scaling factor for the tolerance sequence.
242 /// Typically this is the regularisation parameter. 234 /// Typically this is the regularisation parameter.
243 fn tolerance_scaling(&self) -> F; 235 fn tolerance_scaling(&self) -> F;
244 } 236 }
245 237
246 /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`]. 238 /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`].
247 pub trait SlidingRegTerm<F: Float + ToNalgebraRealField, const N: usize>: RegTerm<F, N> { 239 pub trait SlidingRegTerm<Domain, F = f64>: RegTerm<Domain, F>
240 where
241 Domain: Space + Clone,
242 F: Float + ToNalgebraRealField,
243 {
248 /// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation 244 /// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation
249 /// term, such that $-ε ≤ τw - d ≤ ε$. 245 /// term, such that $-ε ≤ τw - d ≤ ε$.
250 fn goodness<G, BT>( 246 fn goodness<M>(
251 &self, 247 &self,
252 d: &mut BTFN<F, G, BT, N>, 248 d: &mut M,
253 μ: &RNDM<F, N>, 249 μ: &DiscreteMeasure<Domain, F>,
254 y: &Loc<F, N>, 250 y: &Domain,
255 z: &Loc<F, N>, 251 z: &Domain,
256 τ: F, 252 τ: F,
257 ε: F, 253 ε: F,
258 config: &FBGenericConfig<F>, 254 config: &InsertionConfig<F>,
259 ) -> F 255 ) -> F
260 where 256 where
261 BT: BTSearch<F, N, Agg = Bounds<F>>, 257 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 258
265 /// Convert bound on the regulariser to a bond on the Radon norm 259 /// Convert bound on the regulariser to a bond on the Radon norm
266 fn radon_norm_bound(&self, b: F) -> F; 260 fn radon_norm_bound(&self, b: F) -> F;
267 } 261 }
268 262
269 #[replace_float_literals(F::cast_from(literal))] 263 #[replace_float_literals(F::cast_from(literal))]
270 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for NonnegRadonRegTerm<F> 264 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<Loc<N, F>, F>
271 where 265 for NonnegRadonRegTerm<F>
272 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
273 { 266 {
274 fn solve_findim( 267 fn solve_findim(
275 &self, 268 &self,
276 mA: &DMatrix<F::MixedType>, 269 mA: &DMatrix<F::MixedType>,
277 g: &DVector<F::MixedType>, 270 g: &DVector<F::MixedType>,
278 τ: F, 271 τ: F,
279 x: &mut DVector<F::MixedType>, 272 x: &mut DVector<F::MixedType>,
280 mA_normest: F, 273 mA_normest: F,
281 ε: F, 274 ε: F,
282 config: &FBGenericConfig<F>, 275 config: &InsertionConfig<F>,
283 ) -> usize { 276 ) -> usize {
284 let inner_tolerance = ε * config.inner.tolerance_mult; 277 let inner_tolerance = ε * config.inner.tolerance_mult;
285 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); 278 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
286 quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) 279 quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
287 } 280 }
291 y: &DVector<F::MixedType>, 284 y: &DVector<F::MixedType>,
292 g: &DVector<F::MixedType>, 285 g: &DVector<F::MixedType>,
293 τ: F, 286 τ: F,
294 x: &mut DVector<F::MixedType>, 287 x: &mut DVector<F::MixedType>,
295 ε: F, 288 ε: F,
296 config: &FBGenericConfig<F>, 289 config: &InsertionConfig<F>,
297 ) -> usize { 290 ) -> usize {
298 let inner_tolerance = ε * config.inner.tolerance_mult; 291 let inner_tolerance = ε * config.inner.tolerance_mult;
299 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); 292 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
300 l1squared_nonneg(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it) 293 l1squared_nonneg(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
301 } 294 }
302 295
303 #[inline] 296 #[inline]
304 fn find_tolerance_violation_slack<G, BT>( 297 fn find_tolerance_violation_slack<M>(
305 &self, 298 &self,
306 d: &mut BTFN<F, G, BT, N>, 299 d: &mut M,
307 τ: F, 300 τ: F,
308 ε: F, 301 ε: F,
309 skip_by_rough_check: bool, 302 skip_by_rough_check: bool,
310 config: &FBGenericConfig<F>, 303 config: &InsertionConfig<F>,
311 slack: F, 304 slack: F,
312 ) -> Option<(Loc<F, N>, F, bool)> 305 ) -> Option<(Loc<N, F>, F, bool)>
313 where 306 where
314 BT: BTSearch<F, N, Agg = Bounds<F>>, 307 M: MinMaxMapping<Loc<N, F>, F>,
315 G: SupportGenerator<F, N, Id = BT::Data>,
316 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
317 { 308 {
318 let τα = τ * self.α(); 309 let τα = τ * self.α();
319 let keep_above = -τα - slack - ε; 310 let keep_above = -τα - slack - ε;
320 let minimise_below = -τα - slack - ε * config.insertion_cutoff_factor; 311 let minimise_below = -τα - slack - ε * config.insertion_cutoff_factor;
321 let refinement_tolerance = ε * config.refinement.tolerance_mult; 312 let refinement_tolerance = ε * config.refinement.tolerance_mult;
313
314 // println!(
315 // "keep_above: {keep_above}, rough lower bound: {}, tolerance: {ε}, slack: {slack}, τα: {τα}",
316 // d.bounds().lower()
317 // );
322 318
323 // If preliminary check indicates that we are in bounds, and if it otherwise matches 319 // If preliminary check indicates that we are in bounds, and if it otherwise matches
324 // the insertion strategy, skip insertion. 320 // the insertion strategy, skip insertion.
325 if skip_by_rough_check && d.bounds().lower() >= keep_above { 321 if skip_by_rough_check && d.bounds().lower() >= keep_above {
326 None 322 None
327 } else { 323 } else {
328 // If the rough check didn't indicate no insertion needed, find minimising point. 324 // If the rough check didn't indicate no insertion needed, find minimising point.
329 d.minimise_below( 325 let res = d.minimise_below(
330 minimise_below, 326 minimise_below,
331 refinement_tolerance, 327 refinement_tolerance,
332 config.refinement.max_steps, 328 config.refinement.max_steps,
333 ) 329 );
334 .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above)) 330
331 res.map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
335 } 332 }
336 } 333 }
337 334
338 fn verify_merge_candidate<G, BT>( 335 fn verify_merge_candidate<M>(
339 &self, 336 &self,
340 d: &mut BTFN<F, G, BT, N>, 337 d: &mut M,
341 μ: &RNDM<F, N>, 338 μ: &RNDM<N, F>,
342 τ: F, 339 τ: F,
343 ε: F, 340 ε: F,
344 config: &FBGenericConfig<F>, 341 config: &InsertionConfig<F>,
345 ) -> bool 342 ) -> bool
346 where 343 where
347 BT: BTSearch<F, N, Agg = Bounds<F>>, 344 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 { 345 {
351 let τα = τ * self.α(); 346 let τα = τ * self.α();
352 let refinement_tolerance = ε * config.refinement.tolerance_mult; 347 let refinement_tolerance = ε * config.refinement.tolerance_mult;
353 let merge_tolerance = config.merge_tolerance_mult * ε; 348 let merge_tolerance = config.merge_tolerance_mult * ε;
354 let keep_above = -τα - merge_tolerance; 349 let keep_above = -τα - merge_tolerance;
365 refinement_tolerance, 360 refinement_tolerance,
366 config.refinement.max_steps, 361 config.refinement.max_steps,
367 )); 362 ));
368 } 363 }
369 364
370 fn verify_merge_candidate_radonsq<G, BT>( 365 fn verify_merge_candidate_radonsq<M>(
371 &self, 366 &self,
372 d: &mut BTFN<F, G, BT, N>, 367 d: &mut M,
373 μ: &RNDM<F, N>, 368 μ: &RNDM<N, F>,
374 τ: F, 369 τ: F,
375 ε: F, 370 ε: F,
376 config: &FBGenericConfig<F>, 371 config: &InsertionConfig<F>,
377 radon_μ: &RNDM<F, N>, 372 radon_μ: &RNDM<N, F>,
378 ) -> bool 373 ) -> bool
379 where 374 where
380 BT: BTSearch<F, N, Agg = Bounds<F>>, 375 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 { 376 {
384 let τα = τ * self.α(); 377 let τα = τ * self.α();
385 let refinement_tolerance = ε * config.refinement.tolerance_mult; 378 let refinement_tolerance = ε * config.refinement.tolerance_mult;
386 let merge_tolerance = config.merge_tolerance_mult * ε; 379 let merge_tolerance = config.merge_tolerance_mult * ε;
387 let slack = radon_μ.norm(Radon); 380 let slack = radon_μ.norm(Radon);
424 self.α() 417 self.α()
425 } 418 }
426 } 419 }
427 420
428 #[replace_float_literals(F::cast_from(literal))] 421 #[replace_float_literals(F::cast_from(literal))]
429 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for NonnegRadonRegTerm<F> 422 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<Loc<N, F>, F>
430 where 423 for NonnegRadonRegTerm<F>
431 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
432 { 424 {
433 fn goodness<G, BT>( 425 fn goodness<M>(
434 &self, 426 &self,
435 d: &mut BTFN<F, G, BT, N>, 427 d: &mut M,
436 _μ: &RNDM<F, N>, 428 _μ: &RNDM<N, F>,
437 y: &Loc<F, N>, 429 y: &Loc<N, F>,
438 z: &Loc<F, N>, 430 z: &Loc<N, F>,
439 τ: F, 431 τ: F,
440 ε: F, 432 ε: F,
441 _config: &FBGenericConfig<F>, 433 _config: &InsertionConfig<F>,
442 ) -> F 434 ) -> F
443 where 435 where
444 BT: BTSearch<F, N, Agg = Bounds<F>>, 436 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 { 437 {
448 let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α())); 438 let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α()));
449 w(z) - w(y) 439 w(z) - w(y)
450 } 440 }
451 441
453 b / self.α() 443 b / self.α()
454 } 444 }
455 } 445 }
456 446
457 #[replace_float_literals(F::cast_from(literal))] 447 #[replace_float_literals(F::cast_from(literal))]
458 impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for RadonRegTerm<F> 448 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( 449 fn solve_findim(
463 &self, 450 &self,
464 mA: &DMatrix<F::MixedType>, 451 mA: &DMatrix<F::MixedType>,
465 g: &DVector<F::MixedType>, 452 g: &DVector<F::MixedType>,
466 τ: F, 453 τ: F,
467 x: &mut DVector<F::MixedType>, 454 x: &mut DVector<F::MixedType>,
468 mA_normest: F, 455 mA_normest: F,
469 ε: F, 456 ε: F,
470 config: &FBGenericConfig<F>, 457 config: &InsertionConfig<F>,
471 ) -> usize { 458 ) -> usize {
472 let inner_tolerance = ε * config.inner.tolerance_mult; 459 let inner_tolerance = ε * config.inner.tolerance_mult;
473 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); 460 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
474 quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) 461 quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
475 } 462 }
479 y: &DVector<F::MixedType>, 466 y: &DVector<F::MixedType>,
480 g: &DVector<F::MixedType>, 467 g: &DVector<F::MixedType>,
481 τ: F, 468 τ: F,
482 x: &mut DVector<F::MixedType>, 469 x: &mut DVector<F::MixedType>,
483 ε: F, 470 ε: F,
484 config: &FBGenericConfig<F>, 471 config: &InsertionConfig<F>,
485 ) -> usize { 472 ) -> usize {
486 let inner_tolerance = ε * config.inner.tolerance_mult; 473 let inner_tolerance = ε * config.inner.tolerance_mult;
487 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); 474 let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
488 l1squared_unconstrained(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it) 475 l1squared_unconstrained(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
489 } 476 }
490 477
491 fn find_tolerance_violation_slack<G, BT>( 478 fn find_tolerance_violation_slack<M>(
492 &self, 479 &self,
493 d: &mut BTFN<F, G, BT, N>, 480 d: &mut M,
494 τ: F, 481 τ: F,
495 ε: F, 482 ε: F,
496 skip_by_rough_check: bool, 483 skip_by_rough_check: bool,
497 config: &FBGenericConfig<F>, 484 config: &InsertionConfig<F>,
498 slack: F, 485 slack: F,
499 ) -> Option<(Loc<F, N>, F, bool)> 486 ) -> Option<(Loc<N, F>, F, bool)>
500 where 487 where
501 BT: BTSearch<F, N, Agg = Bounds<F>>, 488 M: MinMaxMapping<Loc<N, F>, F>,
502 G: SupportGenerator<F, N, Id = BT::Data>,
503 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
504 { 489 {
505 let τα = τ * self.α(); 490 let τα = τ * self.α();
506 let keep_below = τα + slack + ε; 491 let keep_below = τα + slack + ε;
507 let keep_above = -(τα + slack) - ε; 492 let keep_above = -(τα + slack) - ε;
508 let maximise_above = τα + slack + ε * config.insertion_cutoff_factor; 493 let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
539 } 524 }
540 } 525 }
541 } 526 }
542 } 527 }
543 528
544 fn verify_merge_candidate<G, BT>( 529 fn verify_merge_candidate<M>(
545 &self, 530 &self,
546 d: &mut BTFN<F, G, BT, N>, 531 d: &mut M,
547 μ: &RNDM<F, N>, 532 μ: &RNDM<N, F>,
548 τ: F, 533 τ: F,
549 ε: F, 534 ε: F,
550 config: &FBGenericConfig<F>, 535 config: &InsertionConfig<F>,
551 ) -> bool 536 ) -> bool
552 where 537 where
553 BT: BTSearch<F, N, Agg = Bounds<F>>, 538 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 { 539 {
557 let τα = τ * self.α(); 540 let τα = τ * self.α();
558 let refinement_tolerance = ε * config.refinement.tolerance_mult; 541 let refinement_tolerance = ε * config.refinement.tolerance_mult;
559 let merge_tolerance = config.merge_tolerance_mult * ε; 542 let merge_tolerance = config.merge_tolerance_mult * ε;
560 let keep_below = τα + merge_tolerance; 543 let keep_below = τα + merge_tolerance;
583 refinement_tolerance, 566 refinement_tolerance,
584 config.refinement.max_steps, 567 config.refinement.max_steps,
585 )); 568 ));
586 } 569 }
587 570
588 fn verify_merge_candidate_radonsq<G, BT>( 571 fn verify_merge_candidate_radonsq<M>(
589 &self, 572 &self,
590 d: &mut BTFN<F, G, BT, N>, 573 d: &mut M,
591 μ: &RNDM<F, N>, 574 μ: &RNDM<N, F>,
592 τ: F, 575 τ: F,
593 ε: F, 576 ε: F,
594 config: &FBGenericConfig<F>, 577 config: &InsertionConfig<F>,
595 radon_μ: &RNDM<F, N>, 578 radon_μ: &RNDM<N, F>,
596 ) -> bool 579 ) -> bool
597 where 580 where
598 BT: BTSearch<F, N, Agg = Bounds<F>>, 581 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 { 582 {
602 let τα = τ * self.α(); 583 let τα = τ * self.α();
603 let refinement_tolerance = ε * config.refinement.tolerance_mult; 584 let refinement_tolerance = ε * config.refinement.tolerance_mult;
604 let merge_tolerance = config.merge_tolerance_mult * ε; 585 let merge_tolerance = config.merge_tolerance_mult * ε;
605 let slack = radon_μ.norm(Radon); 586 let slack = radon_μ.norm(Radon);
648 self.α() 629 self.α()
649 } 630 }
650 } 631 }
651 632
652 #[replace_float_literals(F::cast_from(literal))] 633 #[replace_float_literals(F::cast_from(literal))]
653 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<F, N> for RadonRegTerm<F> 634 impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<Loc<N, F>, F>
654 where 635 for RadonRegTerm<F>
655 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
656 { 636 {
657 fn goodness<G, BT>( 637 fn goodness<M>(
658 &self, 638 &self,
659 d: &mut BTFN<F, G, BT, N>, 639 d: &mut M,
660 _μ: &RNDM<F, N>, 640 _μ: &RNDM<N, F>,
661 y: &Loc<F, N>, 641 y: &Loc<N, F>,
662 z: &Loc<F, N>, 642 z: &Loc<N, F>,
663 τ: F, 643 τ: F,
664 ε: F, 644 ε: F,
665 _config: &FBGenericConfig<F>, 645 _config: &InsertionConfig<F>,
666 ) -> F 646 ) -> F
667 where 647 where
668 BT: BTSearch<F, N, Agg = Bounds<F>>, 648 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 { 649 {
672 let α = self.α(); 650 let α = self.α();
673 let w = |x| { 651 let w = |x| {
674 let dx = d.apply(x); 652 let dx = d.apply(x);
675 ((-ε + dx) / (τ * α)).max(1.0.min(ε + dx) / (τ * α)) 653 ((-ε + dx) / (τ * α)).max(1.0.min(ε + dx) / (τ * α))

Mercurial Repository: pointsource_algs

mercurial