Mercurial > repos > pointsource_algs / file comparison

comparison: src/regularisation.rs

src/regularisation.rs

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

Mercurial Repository: pointsource_algs

mercurial