Fri, 18 Nov 2022 10:28:47 +0200
Add package metadata
0 | 1 | |
2 | use numeric_literals::replace_float_literals; | |
3 | use std::iter::Sum; | |
4 | use std::marker::PhantomData; | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
5 | use std::sync::Arc; |
0 | 6 | use crate::types::Float; |
7 | use crate::mapping::Mapping; | |
8 | use crate::linops::Linear; | |
9 | use crate::sets::Set; | |
5 | 10 | use crate::sets::Cube; |
11 | use crate::loc::Loc; | |
0 | 12 | use super::support::*; |
13 | use super::bt::*; | |
14 | use super::refine::*; | |
15 | use super::aggregator::*; | |
16 | use super::either::*; | |
17 | use crate::fe_model::base::RealLocalModel; | |
18 | use crate::fe_model::p2_local_model::*; | |
19 | ||
5 | 20 | /// Presentation for (mathematical) functions constructed as a sum of components functions with |
21 | /// typically small support. | |
22 | /// | |
23 | /// The domain of the function is [`Loc`]`<F, N>`, where `F` is the type of floating point numbers, | |
24 | /// and `N` the dimension. | |
25 | /// | |
26 | /// The `generator` lists the component functions that have to implement [`Support`]. | |
27 | /// Identifiers of the components ([`SupportGenerator::Id`], usually `usize`) are stored stored | |
28 | /// in a [bisection tree][BTImpl], when one is provided as `bt`. However `bt` may also be `()` | |
29 | /// for a [`PreBTFN`] that is only useful for vector space operations with a full [`BTFN`]. | |
0 | 30 | #[derive(Clone,Debug)] |
31 | pub struct BTFN< | |
32 | F : Float, | |
33 | G : SupportGenerator<F, N>, | |
34 | BT /*: BTImpl<F, N>*/, | |
35 | const N : usize | |
36 | > /*where G::SupportType : LocalAnalysis<F, A, N>*/ { | |
37 | bt : BT, | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
38 | generator : Arc<G>, |
0 | 39 | _phantoms : PhantomData<F>, |
40 | } | |
41 | ||
42 | impl<F : Float, G, BT, const N : usize> | |
43 | BTFN<F, G, BT, N> | |
44 | where G : SupportGenerator<F, N, Id=BT::Data>, | |
45 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
46 | BT : BTImpl<F, N> { | |
47 | ||
5 | 48 | /// Create a new BTFN from a support generator and a pre-initialised bisection tree. |
49 | /// | |
50 | /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`. | |
51 | /// Use [`Self::construct`] if no preinitialised tree is available. Use [`Self::new_refresh`] | |
52 | /// when the aggregators of the tree may need updates. | |
53 | /// | |
54 | /// See the documentation for [`BTFN`] on the role of the `generator`. | |
0 | 55 | pub fn new(bt : BT, generator : G) -> Self { |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
56 | Self::new_arc(bt, Arc::new(generator)) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
57 | } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
58 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
59 | fn new_arc(bt : BT, generator : Arc<G>) -> Self { |
0 | 60 | BTFN { |
61 | bt : bt, | |
62 | generator : generator, | |
63 | _phantoms : std::marker::PhantomData, | |
64 | } | |
65 | } | |
66 | ||
5 | 67 | /// Create a new BTFN support generator and a pre-initialised bisection tree, |
68 | /// cloning the tree and refreshing aggregators. | |
69 | /// | |
70 | /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`, but | |
71 | /// the aggregator may be out of date. | |
72 | /// | |
73 | /// See the documentation for [`BTFN`] on the role of the `generator`. | |
0 | 74 | pub fn new_refresh(bt : &BT, generator : G) -> Self { |
75 | // clone().refresh_aggregator(…) as opposed to convert_aggregator | |
76 | // ensures that type is maintained. Due to Rc-pointer copy-on-write, | |
77 | // the effort is not significantly different. | |
78 | let mut btnew = bt.clone(); | |
79 | btnew.refresh_aggregator(&generator); | |
80 | BTFN::new(btnew, generator) | |
81 | } | |
82 | ||
5 | 83 | /// Create a new BTFN from a support generator, domain, and depth for a new [`BT`]. |
84 | /// | |
85 | /// The top node of the created [`BT`] will have the given `domain`. | |
86 | /// | |
87 | /// See the documentation for [`BTFN`] on the role of the `generator`. | |
0 | 88 | pub fn construct(domain : Cube<F, N>, depth : BT::Depth, generator : G) -> Self { |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
89 | Self::construct_arc(domain, depth, Arc::new(generator)) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
90 | } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
91 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
92 | fn construct_arc(domain : Cube<F, N>, depth : BT::Depth, generator : Arc<G>) -> Self { |
0 | 93 | let mut bt = BT::new(domain, depth); |
94 | for (d, support) in generator.all_data() { | |
95 | bt.insert(d, &support); | |
96 | } | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
97 | Self::new_arc(bt, generator) |
0 | 98 | } |
99 | ||
5 | 100 | /// Convert the aggregator of the [`BTFN`] to a different one. |
101 | /// | |
102 | /// This will construct a [`BTFN`] with the same components and generator as the (consumed) | |
103 | /// `self`, but a new `BT` with [`Aggregator`]s of type `ANew`. | |
0 | 104 | pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N> |
105 | where ANew : Aggregator, | |
106 | G : SupportGenerator<F, N, Id=BT::Data>, | |
107 | G::SupportType : LocalAnalysis<F, ANew, N> { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
108 | BTFN::new_arc(self.bt.convert_aggregator(&*self.generator), self.generator) |
0 | 109 | } |
110 | ||
111 | /// Change the generator (after, e.g., a scaling of the latter). | |
5 | 112 | fn new_generator(&self, generator : G) -> Self { |
0 | 113 | BTFN::new_refresh(&self.bt, generator) |
114 | } | |
115 | ||
116 | /// Refresh aggregator after updates to generator | |
117 | fn refresh_aggregator(&mut self) { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
118 | self.bt.refresh_aggregator(&*self.generator); |
0 | 119 | } |
120 | ||
121 | } | |
122 | ||
5 | 123 | impl<F : Float, G, BT, const N : usize> |
124 | BTFN<F, G, BT, N> | |
125 | where G : SupportGenerator<F, N> { | |
126 | /// Change the [bisection tree][BTImpl] of the [`BTFN`] to a different one. | |
127 | /// | |
128 | /// This can be used to convert a [`PreBTFN`] to a full [`BTFN`], or the change | |
129 | /// the aggreagator; see also [`self.convert_aggregator`]. | |
130 | pub fn instantiate< | |
131 | BTNew : BTImpl<F, N, Data=G::Id>, | |
132 | > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<F, G, BTNew, N> | |
133 | where G::SupportType : LocalAnalysis<F, BTNew::Agg, N> { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
134 | BTFN::construct_arc(domain, depth, self.generator) |
5 | 135 | } |
136 | } | |
137 | ||
138 | /// A BTFN with no bisection tree. | |
139 | /// | |
140 | /// Most BTFN methods are not available, but if a BTFN is going to be summed with another | |
141 | /// before other use, it will be more efficient to not construct an unnecessary bisection tree | |
142 | /// that would be shortly dropped. | |
0 | 143 | pub type PreBTFN<F, G, const N : usize> = BTFN<F, G, (), N>; |
144 | ||
145 | impl<F : Float, G, const N : usize> PreBTFN<F, G, N> where G : SupportGenerator<F, N> { | |
146 | ||
5 | 147 | /// Create a new [`PreBTFN`] with no bisection tree. |
0 | 148 | pub fn new_pre(generator : G) -> Self { |
149 | BTFN { | |
150 | bt : (), | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
151 | generator : Arc::new(generator), |
0 | 152 | _phantoms : std::marker::PhantomData, |
153 | } | |
154 | } | |
155 | } | |
156 | ||
157 | impl<F : Float, G, BT, const N : usize> | |
158 | BTFN<F, G, BT, N> | |
159 | where G : SupportGenerator<F, N, Id=usize>, | |
160 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
161 | BT : BTImpl<F, N, Data=usize> { | |
162 | ||
163 | /// Helper function for implementing [`std::ops::Add`]. | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
164 | fn add_another<G2>(&self, g2 : Arc<G2>) -> BTFN<F, BothGenerators<G, G2>, BT, N> |
0 | 165 | where G2 : SupportGenerator<F, N, Id=usize>, |
166 | G2::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
167 | ||
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
168 | let mut bt = self.bt.clone(); |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
169 | let both = BothGenerators(Arc::clone(&self.generator), g2); |
0 | 170 | |
171 | for (d, support) in both.all_right_data() { | |
172 | bt.insert(d, &support); | |
173 | } | |
174 | ||
175 | BTFN { | |
176 | bt : bt, | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
177 | generator : Arc::new(both), |
0 | 178 | _phantoms : std::marker::PhantomData, |
179 | } | |
180 | } | |
181 | } | |
182 | ||
183 | macro_rules! make_btfn_add { | |
184 | ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => { | |
185 | impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> | |
186 | std::ops::Add<BTFN<F, G2, BT2, N>> for | |
187 | $lhs | |
188 | where BT1 : BTImpl<F, N, Data=usize>, | |
189 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
190 | G2 : SupportGenerator<F, N, Id=usize>, | |
191 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
192 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { | |
193 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
194 | #[inline] | |
195 | fn add(self, other : BTFN<F, G2, BT2, N>) -> Self::Output { | |
196 | $preprocess(self).add_another(other.generator) | |
197 | } | |
198 | } | |
199 | ||
200 | impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> | |
201 | std::ops::Add<&'b BTFN<F, G2, BT2, N>> for | |
202 | $lhs | |
203 | where BT1 : BTImpl<F, N, Data=usize>, | |
204 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
205 | G2 : SupportGenerator<F, N, Id=usize>, |
0 | 206 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, |
207 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { | |
208 | ||
209 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
210 | #[inline] | |
211 | fn add(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output { | |
212 | $preprocess(self).add_another(other.generator.clone()) | |
213 | } | |
214 | } | |
215 | } | |
216 | } | |
217 | ||
218 | make_btfn_add!(BTFN<F, G1, BT1, N>, std::convert::identity, ); | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
219 | make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone, ); |
0 | 220 | |
221 | macro_rules! make_btfn_sub { | |
222 | ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => { | |
223 | impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> | |
224 | std::ops::Sub<BTFN<F, G2, BT2, N>> for | |
225 | $lhs | |
226 | where BT1 : BTImpl<F, N, Data=usize>, | |
227 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
228 | G2 : SupportGenerator<F, N, Id=usize>, | |
229 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
230 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { | |
231 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
232 | #[inline] | |
233 | fn sub(self, other : BTFN<F, G2, BT2, N>) -> Self::Output { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
234 | $preprocess(self).add_another(Arc::new( |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
235 | Arc::try_unwrap(other.generator) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
236 | .unwrap_or_else(|arc| (*arc).clone()) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
237 | .neg() |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
238 | )) |
0 | 239 | } |
240 | } | |
241 | ||
242 | impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> | |
243 | std::ops::Sub<&'b BTFN<F, G2, BT2, N>> for | |
244 | $lhs | |
245 | where BT1 : BTImpl<F, N, Data=usize>, | |
246 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
247 | G2 : SupportGenerator<F, N, Id=usize> + Clone, | |
248 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
249 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
5 | 250 | &'b G2 : std::ops::Neg<Output=G2> { |
0 | 251 | |
252 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
253 | #[inline] | |
254 | fn sub(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
255 | $preprocess(self).add_another(Arc::new((*other.generator).clone().neg())) |
0 | 256 | } |
257 | } | |
258 | } | |
259 | } | |
260 | ||
261 | make_btfn_sub!(BTFN<F, G1, BT1, N>, std::convert::identity, ); | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
262 | make_btfn_sub!(&'a BTFN<F, G1, BT1, N>, std::convert::identity, ); |
0 | 263 | |
264 | macro_rules! make_btfn_scalarop_rhs { | |
265 | ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { | |
266 | impl<F : Float, G, BT, const N : usize> | |
267 | std::ops::$trait_assign<F> | |
268 | for BTFN<F, G, BT, N> | |
269 | where BT : BTImpl<F, N>, | |
270 | G : SupportGenerator<F, N, Id=BT::Data>, | |
271 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
272 | #[inline] | |
273 | fn $fn_assign(&mut self, t : F) { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
274 | Arc::make_mut(&mut self.generator).$fn_assign(t); |
0 | 275 | self.refresh_aggregator(); |
276 | } | |
277 | } | |
278 | ||
279 | impl<F : Float, G, BT, const N : usize> | |
280 | std::ops::$trait<F> | |
281 | for BTFN<F, G, BT, N> | |
282 | where BT : BTImpl<F, N>, | |
283 | G : SupportGenerator<F, N, Id=BT::Data>, | |
284 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
285 | type Output = Self; | |
286 | #[inline] | |
287 | fn $fn(mut self, t : F) -> Self::Output { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
288 | Arc::make_mut(&mut self.generator).$fn_assign(t); |
0 | 289 | self.refresh_aggregator(); |
290 | self | |
291 | } | |
292 | } | |
293 | ||
294 | impl<'a, F : Float, G, BT, const N : usize> | |
295 | std::ops::$trait<F> | |
296 | for &'a BTFN<F, G, BT, N> | |
297 | where BT : BTImpl<F, N>, | |
298 | G : SupportGenerator<F, N, Id=BT::Data>, | |
299 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
300 | &'a G : std::ops::$trait<F,Output=G> { | |
301 | type Output = BTFN<F, G, BT, N>; | |
302 | #[inline] | |
303 | fn $fn(self, t : F) -> Self::Output { | |
304 | self.new_generator(self.generator.$fn(t)) | |
305 | } | |
306 | } | |
307 | } | |
308 | } | |
309 | ||
310 | make_btfn_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); | |
311 | make_btfn_scalarop_rhs!(Div, div, DivAssign, div_assign); | |
312 | ||
313 | macro_rules! make_btfn_scalarop_lhs { | |
314 | ($trait:ident, $fn:ident, $fn_assign:ident, $($f:ident)+) => { $( | |
315 | impl<G, BT, const N : usize> | |
316 | std::ops::$trait<BTFN<$f, G, BT, N>> | |
317 | for $f | |
318 | where BT : BTImpl<$f, N>, | |
319 | G : SupportGenerator<$f, N, Id=BT::Data>, | |
320 | G::SupportType : LocalAnalysis<$f, BT::Agg, N> { | |
321 | type Output = BTFN<$f, G, BT, N>; | |
322 | #[inline] | |
323 | fn $fn(self, mut a : BTFN<$f, G, BT, N>) -> Self::Output { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
324 | Arc::make_mut(&mut a.generator).$fn_assign(self); |
0 | 325 | a.refresh_aggregator(); |
326 | a | |
327 | } | |
328 | } | |
329 | ||
330 | impl<'a, G, BT, const N : usize> | |
331 | std::ops::$trait<&'a BTFN<$f, G, BT, N>> | |
332 | for $f | |
333 | where BT : BTImpl<$f, N>, | |
334 | G : SupportGenerator<$f, N, Id=BT::Data> + Clone, | |
335 | G::SupportType : LocalAnalysis<$f, BT::Agg, N>, | |
336 | // FIXME: This causes compiler overflow | |
337 | /*&'a G : std::ops::$trait<$f,Output=G>*/ { | |
338 | type Output = BTFN<$f, G, BT, N>; | |
339 | #[inline] | |
340 | fn $fn(self, a : &'a BTFN<$f, G, BT, N>) -> Self::Output { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
341 | let mut tmp = (*a.generator).clone(); |
0 | 342 | tmp.$fn_assign(self); |
343 | a.new_generator(tmp) | |
344 | // FIXME: Prevented by the compiler overflow above. | |
345 | //a.new_generator(a.generator.$fn(a)) | |
346 | } | |
347 | } | |
348 | )+ } | |
349 | } | |
350 | ||
351 | make_btfn_scalarop_lhs!(Mul, mul, mul_assign, f32 f64); | |
352 | make_btfn_scalarop_lhs!(Div, div, div_assign, f32 f64); | |
353 | ||
354 | macro_rules! make_btfn_unaryop { | |
355 | ($trait:ident, $fn:ident) => { | |
356 | impl<F : Float, G, BT, const N : usize> | |
357 | std::ops::$trait | |
358 | for BTFN<F, G, BT, N> | |
359 | where BT : BTImpl<F, N>, | |
360 | G : SupportGenerator<F, N, Id=BT::Data>, | |
361 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
362 | type Output = Self; | |
363 | #[inline] | |
364 | fn $fn(mut self) -> Self::Output { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
365 | self.generator = Arc::new(Arc::unwrap_or_clone(self.generator).$fn()); |
0 | 366 | self.refresh_aggregator(); |
367 | self | |
368 | } | |
369 | } | |
370 | ||
371 | /*impl<'a, F : Float, G, BT, const N : usize> | |
372 | std::ops::$trait | |
373 | for &'a BTFN<F, G, BT, N> | |
374 | where BT : BTImpl<F, N>, | |
375 | G : SupportGenerator<F, N, Id=BT::Data>, | |
376 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
377 | &'a G : std::ops::$trait<Output=G> { | |
378 | type Output = BTFN<F, G, BT, N>; | |
379 | #[inline] | |
380 | fn $fn(self) -> Self::Output { | |
381 | self.new_generator(std::ops::$trait::$fn(&self.generator)) | |
382 | } | |
383 | }*/ | |
384 | } | |
385 | } | |
386 | ||
387 | make_btfn_unaryop!(Neg, neg); | |
388 | ||
389 | ||
390 | ||
391 | // | |
392 | // Mapping | |
393 | // | |
394 | ||
395 | impl<'a, F : Float, G, BT, V, const N : usize> Mapping<&'a Loc<F,N>> | |
396 | for BTFN<F, G, BT, N> | |
397 | where BT : BTImpl<F, N>, | |
398 | G : SupportGenerator<F, N, Id=BT::Data>, | |
399 | G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<&'a Loc<F,N>, Codomain = V>, | |
400 | V : Sum { | |
401 | ||
402 | type Codomain = V; | |
403 | ||
404 | fn value(&self, x : &'a Loc<F,N>) -> Self::Codomain { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
405 | self.bt.iter_at(x) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
406 | .map(|&d| self.generator.support_for(d).value(x)).sum() |
0 | 407 | } |
408 | } | |
409 | ||
410 | impl<F : Float, G, BT, V, const N : usize> Mapping<Loc<F,N>> | |
411 | for BTFN<F, G, BT, N> | |
412 | where BT : BTImpl<F, N>, | |
413 | G : SupportGenerator<F, N, Id=BT::Data>, | |
414 | G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F,N>, Codomain = V>, | |
415 | V : Sum { | |
416 | ||
417 | type Codomain = V; | |
418 | ||
419 | fn value(&self, x : Loc<F,N>) -> Self::Codomain { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
420 | self.bt.iter_at(&x) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
421 | .map(|&d| self.generator.support_for(d).value(x)).sum() |
0 | 422 | } |
423 | } | |
424 | ||
425 | impl<F : Float, G, BT, const N : usize> GlobalAnalysis<F, BT::Agg> | |
426 | for BTFN<F, G, BT, N> | |
427 | where BT : BTImpl<F, N>, | |
428 | G : SupportGenerator<F, N, Id=BT::Data>, | |
429 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
430 | ||
431 | #[inline] | |
432 | fn global_analysis(&self) -> BT::Agg { | |
433 | self.bt.global_analysis() | |
434 | } | |
435 | } | |
436 | ||
437 | // | |
438 | // Blanket implementation of BTFN as a linear functional over objects | |
439 | // that are linear functionals over BTFN. | |
440 | // | |
441 | ||
442 | impl<X, F : Float, G, BT, const N : usize> Linear<X> | |
443 | for BTFN<F, G, BT, N> | |
444 | where BT : BTImpl<F, N>, | |
445 | G : SupportGenerator<F, N, Id=BT::Data>, | |
446 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
447 | X : Linear<BTFN<F, G, BT, N>, Codomain=F> { | |
448 | type Codomain = F; | |
449 | ||
450 | #[inline] | |
451 | fn apply(&self, x : &X) -> F { | |
452 | x.apply(self) | |
453 | } | |
454 | } | |
455 | ||
456 | /// Helper trait for performing approximate minimisation using P2 elements. | |
5 | 457 | /// |
458 | /// `U` is the domain, generally [`Loc`]`<F, N>`, and `F` the type of floating point numbers. | |
459 | /// `Self` is generally a set of `U`, for example, [`Cube`]`<F, N>`. | |
0 | 460 | pub trait P2Minimise<U, F : Float> : Set<U> { |
5 | 461 | /// Minimise `g` over the set presented by `Self`. |
462 | /// | |
463 | /// The function returns `(x, v)` where `x` is the minimiser `v` an approximation of `g(x)`. | |
0 | 464 | fn p2_minimise<G : Fn(&U) -> F>(&self, g : G) -> (U, F); |
465 | ||
466 | } | |
467 | ||
468 | impl<F : Float> P2Minimise<Loc<F, 1>, F> for Cube<F, 1> { | |
469 | fn p2_minimise<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> (Loc<F, 1>, F) { | |
470 | let interval = Simplex(self.corners()); | |
471 | interval.p2_model(&g).minimise(&interval) | |
472 | } | |
473 | } | |
474 | ||
475 | #[replace_float_literals(F::cast_from(literal))] | |
476 | impl<F : Float> P2Minimise<Loc<F, 2>, F> for Cube<F, 2> { | |
477 | fn p2_minimise<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> (Loc<F, 2>, F) { | |
478 | if false { | |
479 | // Split into two triangle (simplex) with separate P2 model in each. | |
480 | // The six nodes of each triangle are the corners and the edges. | |
481 | let [a, b, c, d] = self.corners(); | |
482 | let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)]; | |
483 | ||
484 | let ab = midpoint(&a, &b); | |
485 | let bc = midpoint(&b, &c); | |
486 | let ca = midpoint(&c, &a); | |
487 | let cd = midpoint(&c, &d); | |
488 | let da = midpoint(&d, &a); | |
489 | let [vab, vbc, vca, vcd, vda] = [g(&ab), g(&bc), g(&ca), g(&cd), g(&da)]; | |
490 | ||
491 | let s1 = Simplex([a, b, c]); | |
5 | 492 | let m1 = P2LocalModel::<F, 2, 3>::new( |
0 | 493 | &[a, b, c, ab, bc, ca], |
494 | &[va, vb, vc, vab, vbc, vca] | |
495 | ); | |
496 | ||
497 | let r1@(_, v1) = m1.minimise(&s1); | |
498 | ||
499 | let s2 = Simplex([c, d, a]); | |
5 | 500 | let m2 = P2LocalModel::<F, 2, 3>::new( |
0 | 501 | &[c, d, a, cd, da, ca], |
502 | &[vc, vd, va, vcd, vda, vca] | |
503 | ); | |
504 | ||
505 | let r2@(_, v2) = m2.minimise(&s2); | |
506 | ||
507 | if v1 < v2 { r1 } else { r2 } | |
508 | } else { | |
509 | // Single P2 model for the entire cube. | |
510 | let [a, b, c, d] = self.corners(); | |
511 | let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)]; | |
512 | let [e, f] = match 'r' { | |
513 | 'm' => [(&a + &b + &c) / 3.0, (&c + &d + &a) / 3.0], | |
514 | 'c' => [midpoint(&a, &b), midpoint(&a, &d)], | |
515 | 'w' => [(&a + &b * 2.0) / 3.0, (&a + &d * 2.0) / 3.0], | |
516 | 'r' => { | |
517 | // Pseudo-randomise edge midpoints | |
518 | let Loc([x, y]) = a; | |
519 | let tmp : f64 = (x+y).as_(); | |
520 | match tmp.to_bits() % 4 { | |
521 | 0 => [midpoint(&a, &b), midpoint(&a, &d)], | |
522 | 1 => [midpoint(&c, &d), midpoint(&a, &d)], | |
523 | 2 => [midpoint(&a, &b), midpoint(&b, &c)], | |
524 | _ => [midpoint(&c, &d), midpoint(&b, &c)], | |
525 | } | |
526 | }, | |
527 | _ => [self.center(), (&a + &b) / 2.0], | |
528 | }; | |
529 | let [ve, vf] = [g(&e), g(&f)]; | |
530 | ||
5 | 531 | let m1 = P2LocalModel::<F, 2, 3>::new( |
0 | 532 | &[a, b, c, d, e, f], |
533 | &[va, vb, vc, vd, ve, vf], | |
534 | ); | |
535 | ||
536 | m1.minimise(self) | |
537 | } | |
538 | } | |
539 | } | |
540 | ||
5 | 541 | /// Helper type to use [`P2Refiner`] for maximisation. |
0 | 542 | struct RefineMax; |
5 | 543 | |
544 | /// Helper type to use [`P2Refiner`] for minimisation. | |
0 | 545 | struct RefineMin; |
546 | ||
5 | 547 | /// A bisection tree [`Refiner`] for maximising or minimising a [`BTFN`]. |
548 | /// | |
0 | 549 | /// The type parameter `T` should be either [`RefineMax`] or [`RefineMin`]. |
550 | struct P2Refiner<F : Float, T> { | |
5 | 551 | /// The maximum / minimum should be above / below this threshold. |
552 | /// If the threshold cannot be satisfied, the refiner will return `None`. | |
0 | 553 | bound : Option<F>, |
5 | 554 | /// Tolerance for function value estimation. |
0 | 555 | tolerance : F, |
5 | 556 | /// Maximum number of steps to execute the refiner for |
0 | 557 | max_steps : usize, |
5 | 558 | /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes. |
0 | 559 | #[allow(dead_code)] // `how` is just for type system purposes. |
560 | how : T, | |
561 | } | |
562 | ||
563 | impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> | |
564 | for P2Refiner<F, RefineMax> | |
565 | where Cube<F, N> : P2Minimise<Loc<F, N>, F>, | |
566 | G : SupportGenerator<F, N>, | |
567 | G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F> | |
568 | + LocalAnalysis<F, Bounds<F>, N> { | |
569 | type Result = Option<(Loc<F, N>, F)>; | |
570 | type Sorting = UpperBoundSorting<F>; | |
571 | ||
572 | fn refine( | |
573 | &self, | |
574 | aggregator : &Bounds<F>, | |
575 | cube : &Cube<F, N>, | |
9
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
576 | data : &[G::Id], |
0 | 577 | generator : &G, |
578 | step : usize | |
579 | ) -> RefinerResult<Bounds<F>, Self::Result> { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
580 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
581 | if self.bound.map_or(false, |b| aggregator.upper() <= b + self.tolerance) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
582 | // The upper bound is below the maximisation threshold. Don't bother with this cube. |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
583 | return RefinerResult::Uncertain(*aggregator, None) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
584 | } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
585 | |
0 | 586 | // g gives the negative of the value of the function presented by `data` and `generator`. |
587 | let g = move |x : &Loc<F,N>| { | |
588 | let f = move |&d| generator.support_for(d).value(x); | |
589 | -data.iter().map(f).sum::<F>() | |
590 | }; | |
591 | // … so the negative of the minimum is the maximm we want. | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
592 | let (x, _neg_v) = cube.p2_minimise(g); |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
593 | //let v = -neg_v; |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
594 | let v = -g(&x); |
0 | 595 | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
596 | if step < self.max_steps && (aggregator.upper() > v + self.tolerance |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
597 | /*|| aggregator.lower() > v - self.tolerance*/) { |
0 | 598 | // The function isn't refined enough in `cube`, so return None |
599 | // to indicate that further subdivision is required. | |
600 | RefinerResult::NeedRefinement | |
601 | } else { | |
602 | // The data is refined enough, so return new hopefully better bounds | |
603 | // and the maximiser. | |
604 | let res = (x, v); | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
605 | let bounds = Bounds(v, v); |
0 | 606 | RefinerResult::Uncertain(bounds, Some(res)) |
607 | } | |
608 | } | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
609 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
610 | fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
611 | match (*r1, r2) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
612 | (Some((_, v1)), Some((_, v2))) => if v1 < v2 { *r1 = r2 } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
613 | (None, Some(_)) => *r1 = r2, |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
614 | (_, _) => {}, |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
615 | } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
616 | } |
0 | 617 | } |
618 | ||
619 | ||
620 | impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> | |
621 | for P2Refiner<F, RefineMin> | |
622 | where Cube<F, N> : P2Minimise<Loc<F, N>, F>, | |
623 | G : SupportGenerator<F, N>, | |
624 | G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F> | |
625 | + LocalAnalysis<F, Bounds<F>, N> { | |
626 | type Result = Option<(Loc<F, N>, F)>; | |
627 | type Sorting = LowerBoundSorting<F>; | |
628 | ||
629 | fn refine( | |
630 | &self, | |
631 | aggregator : &Bounds<F>, | |
632 | cube : &Cube<F, N>, | |
9
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
633 | data : &[G::Id], |
0 | 634 | generator : &G, |
635 | step : usize | |
636 | ) -> RefinerResult<Bounds<F>, Self::Result> { | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
637 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
638 | if self.bound.map_or(false, |b| aggregator.lower() >= b - self.tolerance) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
639 | // The lower bound is above the minimisation threshold. Don't bother with this cube. |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
640 | return RefinerResult::Uncertain(*aggregator, None) |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
641 | } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
642 | |
0 | 643 | // g gives the value of the function presented by `data` and `generator`. |
644 | let g = move |x : &Loc<F,N>| { | |
645 | let f = move |&d| generator.support_for(d).value(x); | |
646 | data.iter().map(f).sum::<F>() | |
647 | }; | |
648 | // Minimise it. | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
649 | let (x, _v) = cube.p2_minimise(g); |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
650 | let v = g(&x); |
0 | 651 | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
652 | if step < self.max_steps && (aggregator.lower() < v - self.tolerance |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
653 | /*|| aggregator.upper() < v + self.tolerance*/) { |
0 | 654 | // The function isn't refined enough in `cube`, so return None |
655 | // to indicate that further subdivision is required. | |
656 | RefinerResult::NeedRefinement | |
657 | } else { | |
658 | // The data is refined enough, so return new hopefully better bounds | |
659 | // and the minimiser. | |
660 | let res = (x, v); | |
9
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
661 | let l = aggregator.lower(); |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
662 | let bounds = if l > v { |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
663 | eprintln!("imprecision!"); |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
664 | Bounds(l, l) |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
665 | } else { |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
666 | Bounds(v, v) |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
667 | }; |
0 | 668 | RefinerResult::Uncertain(bounds, Some(res)) |
669 | } | |
670 | } | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
671 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
672 | fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
673 | match (*r1, r2) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
674 | (Some((_, v1)), Some((_, v2))) => if v1 > v2 { *r1 = r2 } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
675 | (_, Some(_)) => *r1 = r2, |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
676 | (_, _) => {}, |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
677 | } |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
678 | } |
0 | 679 | } |
680 | ||
681 | ||
5 | 682 | /// A bisection tree [`Refiner`] for checking that a [`BTFN`] is within a stated |
683 | //// upper or lower bound. | |
684 | /// | |
685 | /// The type parameter `T` should be either [`RefineMax`] for upper bound or [`RefineMin`] | |
686 | /// for lower bound. | |
0 | 687 | |
688 | struct BoundRefiner<F : Float, T> { | |
5 | 689 | /// The upper/lower bound to check for |
0 | 690 | bound : F, |
5 | 691 | /// Tolerance for function value estimation. |
0 | 692 | tolerance : F, |
5 | 693 | /// Maximum number of steps to execute the refiner for |
0 | 694 | max_steps : usize, |
695 | #[allow(dead_code)] // `how` is just for type system purposes. | |
5 | 696 | /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes. |
0 | 697 | how : T, |
698 | } | |
699 | ||
700 | impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> | |
701 | for BoundRefiner<F, RefineMax> | |
702 | where G : SupportGenerator<F, N> { | |
703 | type Result = bool; | |
704 | type Sorting = UpperBoundSorting<F>; | |
705 | ||
706 | fn refine( | |
707 | &self, | |
708 | aggregator : &Bounds<F>, | |
709 | _cube : &Cube<F, N>, | |
9
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
710 | _data : &[G::Id], |
0 | 711 | _generator : &G, |
712 | step : usize | |
713 | ) -> RefinerResult<Bounds<F>, Self::Result> { | |
714 | if aggregator.upper() <= self.bound + self.tolerance { | |
715 | // Below upper bound within tolerances. Indicate uncertain success. | |
716 | RefinerResult::Uncertain(*aggregator, true) | |
717 | } else if aggregator.lower() >= self.bound - self.tolerance { | |
718 | // Above upper bound within tolerances. Indicate certain failure. | |
719 | RefinerResult::Certain(false) | |
720 | } else if step < self.max_steps { | |
721 | // No decision possible, but within step bounds - further subdivision is required. | |
722 | RefinerResult::NeedRefinement | |
723 | } else { | |
724 | // No decision possible, but past step bounds | |
725 | RefinerResult::Uncertain(*aggregator, false) | |
726 | } | |
727 | } | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
728 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
729 | fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
730 | *r1 = *r1 && r2; |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
731 | } |
0 | 732 | } |
733 | ||
5 | 734 | impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> |
735 | for BoundRefiner<F, RefineMin> | |
736 | where G : SupportGenerator<F, N> { | |
737 | type Result = bool; | |
738 | type Sorting = UpperBoundSorting<F>; | |
739 | ||
740 | fn refine( | |
741 | &self, | |
742 | aggregator : &Bounds<F>, | |
743 | _cube : &Cube<F, N>, | |
9
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
744 | _data : &[G::Id], |
5 | 745 | _generator : &G, |
746 | step : usize | |
747 | ) -> RefinerResult<Bounds<F>, Self::Result> { | |
748 | if aggregator.lower() >= self.bound - self.tolerance { | |
749 | // Above lower bound within tolerances. Indicate uncertain success. | |
750 | RefinerResult::Uncertain(*aggregator, true) | |
751 | } else if aggregator.upper() <= self.bound + self.tolerance { | |
752 | // Below lower bound within tolerances. Indicate certain failure. | |
753 | RefinerResult::Certain(false) | |
754 | } else if step < self.max_steps { | |
755 | // No decision possible, but within step bounds - further subdivision is required. | |
756 | RefinerResult::NeedRefinement | |
757 | } else { | |
758 | // No decision possible, but past step bounds | |
759 | RefinerResult::Uncertain(*aggregator, false) | |
760 | } | |
761 | } | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
762 | |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
763 | fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
764 | *r1 = *r1 && r2; |
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
765 | } |
5 | 766 | } |
0 | 767 | |
9
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
768 | // FIXME: The most likely reason for the “Refiner failure” expectation in the methods below |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
769 | // is numerical inaccuracy: the `glb` maintained in `HeapContainer` (`refine.rs`) becomes bigger |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
770 | // than the *upper bound* of nodes attempted to be inserted into the `heap` in the container. |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
771 | // But the `glb` is there exactly to prevent that. Due to numerical inaccuracy, however, a |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
772 | // newly subdivided node may have lower upper bound than the original lower bound that should |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
773 | // have been above the `glb` since the node was picked from the queue. Due to the subdivision |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
774 | // process, if a node whose lower bound is at the `glb` is picked, all of its refined subnodes |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
775 | // should have lower bound at least the old `glb`, so in a single-threaded situation there should |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
776 | // always be nodes above the `glb` in the queue. In a multi-threaded situation a node below the |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
777 | // `glb` may be picked by some thread. When that happens, that thread inserts no new nodes into |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
778 | // the queue. If the queue empties as a result of that, the thread goes to wait for other threads |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
779 | // to produce results. Since some node had a node whose lower bound was above the `glb`, eventually |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
780 | // there should be a result, or new nodes above the `glb` inserted into the queue. Then the waiting |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
781 | // threads can also continue processing. If, however, numerical inaccuracy destroyes the `glb`, |
f40dfaf2166d
Improvements and minor fixes to bisection tree refinement.
Tuomo Valkonen <tuomov@iki.fi>
parents:
8
diff
changeset
|
782 | // the queue may run out, and we get “Refiner failure”. |
0 | 783 | impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N> |
784 | where BT : BTSearch<F, N, Agg=Bounds<F>>, | |
785 | G : SupportGenerator<F, N, Id=BT::Data>, | |
786 | G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F> | |
787 | + LocalAnalysis<F, Bounds<F>, N>, | |
788 | Cube<F, N> : P2Minimise<Loc<F, N>, F> { | |
789 | ||
5 | 790 | /// Maximise the `BTFN` within stated value `tolerance`. |
791 | /// | |
792 | /// At most `max_steps` refinement steps are taken. | |
793 | /// Returns the approximate maximiser and the corresponding function value. | |
0 | 794 | pub fn maximise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) { |
795 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : None }; | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
796 | self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.").unwrap() |
0 | 797 | } |
798 | ||
5 | 799 | /// Maximise the `BTFN` within stated value `tolerance` subject to a lower bound. |
800 | /// | |
801 | /// At most `max_steps` refinement steps are taken. | |
802 | /// Returns the approximate maximiser and the corresponding function value when one is found | |
803 | /// above the `bound` threshold, otherwise `None`. | |
0 | 804 | pub fn maximise_above(&mut self, bound : F, tolerance : F, max_steps : usize) |
805 | -> Option<(Loc<F, N>, F)> { | |
806 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : Some(bound) }; | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
807 | self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") |
0 | 808 | } |
809 | ||
5 | 810 | /// Minimise the `BTFN` within stated value `tolerance`. |
811 | /// | |
812 | /// At most `max_steps` refinement steps are taken. | |
813 | /// Returns the approximate minimiser and the corresponding function value. | |
0 | 814 | pub fn minimise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) { |
815 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : None }; | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
816 | self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.").unwrap() |
0 | 817 | } |
818 | ||
5 | 819 | /// Minimise the `BTFN` within stated value `tolerance` subject to a lower bound. |
820 | /// | |
821 | /// At most `max_steps` refinement steps are taken. | |
822 | /// Returns the approximate minimiser and the corresponding function value when one is found | |
823 | /// above the `bound` threshold, otherwise `None`. | |
0 | 824 | pub fn minimise_below(&mut self, bound : F, tolerance : F, max_steps : usize) |
825 | -> Option<(Loc<F, N>, F)> { | |
826 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : Some(bound) }; | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
827 | self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") |
0 | 828 | } |
5 | 829 | |
830 | /// Verify that the `BTFN` has a given upper `bound` within indicated `tolerance`. | |
831 | /// | |
832 | /// At most `max_steps` refinement steps are taken. | |
0 | 833 | pub fn has_upper_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool { |
834 | let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMax }; | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
835 | self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") |
0 | 836 | } |
5 | 837 | |
838 | /// Verify that the `BTFN` has a given lower `bound` within indicated `tolerance`. | |
839 | /// | |
840 | /// At most `max_steps` refinement steps are taken. | |
841 | pub fn has_lower_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool { | |
842 | let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMin }; | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
843 | self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") |
5 | 844 | } |
0 | 845 | } |