Mon, 24 Oct 2022 09:41:43 +0300
Allow step closure of AlgIterators to indicate succesfull termination or failure.
0 | 1 | |
2 | use numeric_literals::replace_float_literals; | |
3 | use std::iter::Sum; | |
4 | use std::marker::PhantomData; | |
5 | pub use nalgebra::Const; | |
6 | use crate::types::Float; | |
7 | use crate::mapping::Mapping; | |
8 | use crate::linops::Linear; | |
9 | use crate::sets::Set; | |
10 | pub use crate::sets::Cube; | |
11 | pub use crate::loc::Loc; | |
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 | ||
20 | /// Presentation for functions constructed as a sum of components with limited support. | |
21 | /// Lookup of relevant components at a specific point of a [`Loc<F,N>`] domain is done | |
22 | /// with a bisection tree (BT). We do notn impose that the `BT` parameter would be an in | |
23 | /// instance [`BTImpl`] to allow performance improvements via “pre-BTFNs” that do not | |
24 | /// construct a BT until being added to another function. Addition and subtraction always | |
25 | /// use a copy-on-write clone of the BT of the left-hand-side operand, discarding | |
26 | /// the BT of the right-hand-side operand. | |
27 | #[derive(Clone,Debug)] | |
28 | pub struct BTFN< | |
29 | F : Float, | |
30 | G : SupportGenerator<F, N>, | |
31 | BT /*: BTImpl<F, N>*/, | |
32 | const N : usize | |
33 | > /*where G::SupportType : LocalAnalysis<F, A, N>*/ { | |
34 | bt : BT, | |
35 | generator : G, | |
36 | _phantoms : PhantomData<F>, | |
37 | } | |
38 | ||
39 | impl<F : Float, G, BT, const N : usize> | |
40 | BTFN<F, G, BT, N> | |
41 | where G : SupportGenerator<F, N, Id=BT::Data>, | |
42 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
43 | BT : BTImpl<F, N> { | |
44 | ||
45 | /// Returns the [`Support`] corresponding to the identifier `d`. | |
46 | /// Simply wraps [`SupportGenerator::support_for`]. | |
47 | #[inline] | |
48 | fn support_for(&self, d : G::Id) -> G::SupportType { | |
49 | self.generator.support_for(d) | |
50 | } | |
51 | ||
52 | /// Create a new BTFN. The bisection tree `bt` should be pre-initialised to | |
53 | /// correspond to the `generator`. Use [`Self::construct`] if no preinitialised tree | |
54 | /// is available. Use [`Self::new_refresh`] when the aggregator may need updates. | |
55 | pub fn new(bt : BT, generator : G) -> Self { | |
56 | BTFN { | |
57 | bt : bt, | |
58 | generator : generator, | |
59 | _phantoms : std::marker::PhantomData, | |
60 | } | |
61 | } | |
62 | ||
63 | /// Create a new BTFN. The bisection tree `bt` should be pre-initialised to | |
64 | /// correspond to the `generator`, but the aggregator may be out of date. | |
65 | pub fn new_refresh(bt : &BT, generator : G) -> Self { | |
66 | // clone().refresh_aggregator(…) as opposed to convert_aggregator | |
67 | // ensures that type is maintained. Due to Rc-pointer copy-on-write, | |
68 | // the effort is not significantly different. | |
69 | let mut btnew = bt.clone(); | |
70 | btnew.refresh_aggregator(&generator); | |
71 | BTFN::new(btnew, generator) | |
72 | } | |
73 | ||
74 | /// Create a new BTFN. Unlike [`Self::new`], this constructs the bisection tree based on | |
75 | /// the generator. | |
76 | pub fn construct(domain : Cube<F, N>, depth : BT::Depth, generator : G) -> Self { | |
77 | let mut bt = BT::new(domain, depth); | |
78 | for (d, support) in generator.all_data() { | |
79 | bt.insert(d, &support); | |
80 | } | |
81 | Self::new(bt, generator) | |
82 | } | |
83 | ||
84 | /// Construct a new instance for a different aggregator | |
85 | pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N> | |
86 | where ANew : Aggregator, | |
87 | G : SupportGenerator<F, N, Id=BT::Data>, | |
88 | G::SupportType : LocalAnalysis<F, ANew, N> { | |
89 | BTFN::new(self.bt.convert_aggregator(&self.generator), self.generator) | |
90 | } | |
91 | ||
92 | /// Change the [`BTImpl`]. | |
93 | pub fn instantiate< | |
94 | BTNew : BTImpl<F, N, Data=BT::Data>, | |
95 | > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<F, G, BTNew, N> | |
96 | where G::SupportType : LocalAnalysis<F, BTNew::Agg, N> { | |
97 | BTFN::construct(domain, depth, self.generator) | |
98 | } | |
99 | ||
100 | /// Change the generator (after, e.g., a scaling of the latter). | |
101 | pub fn new_generator(&self, generator : G) -> Self { | |
102 | BTFN::new_refresh(&self.bt, generator) | |
103 | } | |
104 | ||
105 | /// Refresh aggregator after updates to generator | |
106 | fn refresh_aggregator(&mut self) { | |
107 | self.bt.refresh_aggregator(&self.generator); | |
108 | } | |
109 | ||
110 | } | |
111 | ||
112 | /// A “pre-BTFN” with no bisection tree. | |
113 | pub type PreBTFN<F, G, const N : usize> = BTFN<F, G, (), N>; | |
114 | ||
115 | impl<F : Float, G, const N : usize> PreBTFN<F, G, N> where G : SupportGenerator<F, N> { | |
116 | ||
117 | /// Create a new “pre-BTFN” with no bisection tree. This can feasibly only be | |
118 | /// multiplied by a scalar or added to or subtracted from a proper [`BTFN`]. | |
119 | pub fn new_pre(generator : G) -> Self { | |
120 | BTFN { | |
121 | bt : (), | |
122 | generator : generator, | |
123 | _phantoms : std::marker::PhantomData, | |
124 | } | |
125 | } | |
126 | } | |
127 | ||
128 | impl<F : Float, G, BT, const N : usize> | |
129 | BTFN<F, G, BT, N> | |
130 | where G : SupportGenerator<F, N, Id=usize>, | |
131 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
132 | BT : BTImpl<F, N, Data=usize> { | |
133 | ||
134 | /// Helper function for implementing [`std::ops::Add`]. | |
135 | fn add_another<G2>(self, other : G2) -> BTFN<F, BothGenerators<G, G2>, BT, N> | |
136 | where G2 : SupportGenerator<F, N, Id=usize>, | |
137 | G2::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
138 | ||
139 | let mut bt = self.bt; | |
140 | let both = BothGenerators(self.generator, other); | |
141 | ||
142 | for (d, support) in both.all_right_data() { | |
143 | bt.insert(d, &support); | |
144 | } | |
145 | ||
146 | BTFN { | |
147 | bt : bt, | |
148 | generator : both, | |
149 | _phantoms : std::marker::PhantomData, | |
150 | } | |
151 | } | |
152 | } | |
153 | ||
154 | macro_rules! make_btfn_add { | |
155 | ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => { | |
156 | impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> | |
157 | std::ops::Add<BTFN<F, G2, BT2, N>> for | |
158 | $lhs | |
159 | where BT1 : BTImpl<F, N, Data=usize>, | |
160 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
161 | G2 : SupportGenerator<F, N, Id=usize>, | |
162 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
163 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { | |
164 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
165 | #[inline] | |
166 | fn add(self, other : BTFN<F, G2, BT2, N>) -> Self::Output { | |
167 | $preprocess(self).add_another(other.generator) | |
168 | } | |
169 | } | |
170 | ||
171 | impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> | |
172 | std::ops::Add<&'b BTFN<F, G2, BT2, N>> for | |
173 | $lhs | |
174 | where BT1 : BTImpl<F, N, Data=usize>, | |
175 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
176 | G2 : SupportGenerator<F, N, Id=usize> + Clone, | |
177 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
178 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { | |
179 | ||
180 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
181 | #[inline] | |
182 | fn add(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output { | |
183 | $preprocess(self).add_another(other.generator.clone()) | |
184 | } | |
185 | } | |
186 | } | |
187 | } | |
188 | ||
189 | make_btfn_add!(BTFN<F, G1, BT1, N>, std::convert::identity, ); | |
190 | make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone, Clone); | |
191 | ||
192 | macro_rules! make_btfn_sub { | |
193 | ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => { | |
194 | impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> | |
195 | std::ops::Sub<BTFN<F, G2, BT2, N>> for | |
196 | $lhs | |
197 | where BT1 : BTImpl<F, N, Data=usize>, | |
198 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
199 | G2 : SupportGenerator<F, N, Id=usize>, | |
200 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
201 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { | |
202 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
203 | #[inline] | |
204 | fn sub(self, other : BTFN<F, G2, BT2, N>) -> Self::Output { | |
205 | $preprocess(self).add_another(other.generator.neg()) | |
206 | } | |
207 | } | |
208 | ||
209 | impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> | |
210 | std::ops::Sub<&'b BTFN<F, G2, BT2, N>> for | |
211 | $lhs | |
212 | where BT1 : BTImpl<F, N, Data=usize>, | |
213 | G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, | |
214 | G2 : SupportGenerator<F, N, Id=usize> + Clone, | |
215 | G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
216 | G2::SupportType : LocalAnalysis<F, BT1::Agg, N>, | |
217 | /*&'b G2 : std::ops::Neg<Output=G2>*/ { | |
218 | ||
219 | type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; | |
220 | #[inline] | |
221 | fn sub(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output { | |
222 | // FIXME: the compiler crashes in an overflow on this. | |
223 | //$preprocess(self).add_another(std::ops::Neg::neg(&other.generator)) | |
224 | $preprocess(self).add_another(other.generator.clone().neg()) | |
225 | } | |
226 | } | |
227 | } | |
228 | } | |
229 | ||
230 | make_btfn_sub!(BTFN<F, G1, BT1, N>, std::convert::identity, ); | |
231 | make_btfn_sub!(&'a BTFN<F, G1, BT1, N>, Clone::clone, Clone); | |
232 | ||
233 | macro_rules! make_btfn_scalarop_rhs { | |
234 | ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { | |
235 | impl<F : Float, G, BT, const N : usize> | |
236 | std::ops::$trait_assign<F> | |
237 | for BTFN<F, G, BT, N> | |
238 | where BT : BTImpl<F, N>, | |
239 | G : SupportGenerator<F, N, Id=BT::Data>, | |
240 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
241 | #[inline] | |
242 | fn $fn_assign(&mut self, t : F) { | |
243 | self.generator.$fn_assign(t); | |
244 | self.refresh_aggregator(); | |
245 | } | |
246 | } | |
247 | ||
248 | impl<F : Float, G, BT, const N : usize> | |
249 | std::ops::$trait<F> | |
250 | for BTFN<F, G, BT, N> | |
251 | where BT : BTImpl<F, N>, | |
252 | G : SupportGenerator<F, N, Id=BT::Data>, | |
253 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
254 | type Output = Self; | |
255 | #[inline] | |
256 | fn $fn(mut self, t : F) -> Self::Output { | |
257 | self.generator.$fn_assign(t); | |
258 | self.refresh_aggregator(); | |
259 | self | |
260 | } | |
261 | } | |
262 | ||
263 | impl<'a, F : Float, G, BT, const N : usize> | |
264 | std::ops::$trait<F> | |
265 | for &'a BTFN<F, G, BT, N> | |
266 | where BT : BTImpl<F, N>, | |
267 | G : SupportGenerator<F, N, Id=BT::Data>, | |
268 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
269 | &'a G : std::ops::$trait<F,Output=G> { | |
270 | type Output = BTFN<F, G, BT, N>; | |
271 | #[inline] | |
272 | fn $fn(self, t : F) -> Self::Output { | |
273 | self.new_generator(self.generator.$fn(t)) | |
274 | } | |
275 | } | |
276 | } | |
277 | } | |
278 | ||
279 | make_btfn_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); | |
280 | make_btfn_scalarop_rhs!(Div, div, DivAssign, div_assign); | |
281 | ||
282 | macro_rules! make_btfn_scalarop_lhs { | |
283 | ($trait:ident, $fn:ident, $fn_assign:ident, $($f:ident)+) => { $( | |
284 | impl<G, BT, const N : usize> | |
285 | std::ops::$trait<BTFN<$f, G, BT, N>> | |
286 | for $f | |
287 | where BT : BTImpl<$f, N>, | |
288 | G : SupportGenerator<$f, N, Id=BT::Data>, | |
289 | G::SupportType : LocalAnalysis<$f, BT::Agg, N> { | |
290 | type Output = BTFN<$f, G, BT, N>; | |
291 | #[inline] | |
292 | fn $fn(self, mut a : BTFN<$f, G, BT, N>) -> Self::Output { | |
293 | a.generator.$fn_assign(self); | |
294 | a.refresh_aggregator(); | |
295 | a | |
296 | } | |
297 | } | |
298 | ||
299 | impl<'a, G, BT, const N : usize> | |
300 | std::ops::$trait<&'a BTFN<$f, G, BT, N>> | |
301 | for $f | |
302 | where BT : BTImpl<$f, N>, | |
303 | G : SupportGenerator<$f, N, Id=BT::Data> + Clone, | |
304 | G::SupportType : LocalAnalysis<$f, BT::Agg, N>, | |
305 | // FIXME: This causes compiler overflow | |
306 | /*&'a G : std::ops::$trait<$f,Output=G>*/ { | |
307 | type Output = BTFN<$f, G, BT, N>; | |
308 | #[inline] | |
309 | fn $fn(self, a : &'a BTFN<$f, G, BT, N>) -> Self::Output { | |
310 | let mut tmp = a.generator.clone(); | |
311 | tmp.$fn_assign(self); | |
312 | a.new_generator(tmp) | |
313 | // FIXME: Prevented by the compiler overflow above. | |
314 | //a.new_generator(a.generator.$fn(a)) | |
315 | } | |
316 | } | |
317 | )+ } | |
318 | } | |
319 | ||
320 | make_btfn_scalarop_lhs!(Mul, mul, mul_assign, f32 f64); | |
321 | make_btfn_scalarop_lhs!(Div, div, div_assign, f32 f64); | |
322 | ||
323 | macro_rules! make_btfn_unaryop { | |
324 | ($trait:ident, $fn:ident) => { | |
325 | impl<F : Float, G, BT, const N : usize> | |
326 | std::ops::$trait | |
327 | for BTFN<F, G, BT, N> | |
328 | where BT : BTImpl<F, N>, | |
329 | G : SupportGenerator<F, N, Id=BT::Data>, | |
330 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
331 | type Output = Self; | |
332 | #[inline] | |
333 | fn $fn(mut self) -> Self::Output { | |
334 | self.generator = self.generator.$fn(); | |
335 | self.refresh_aggregator(); | |
336 | self | |
337 | } | |
338 | } | |
339 | ||
340 | /*impl<'a, F : Float, G, BT, const N : usize> | |
341 | std::ops::$trait | |
342 | for &'a BTFN<F, G, BT, N> | |
343 | where BT : BTImpl<F, N>, | |
344 | G : SupportGenerator<F, N, Id=BT::Data>, | |
345 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
346 | &'a G : std::ops::$trait<Output=G> { | |
347 | type Output = BTFN<F, G, BT, N>; | |
348 | #[inline] | |
349 | fn $fn(self) -> Self::Output { | |
350 | self.new_generator(std::ops::$trait::$fn(&self.generator)) | |
351 | } | |
352 | }*/ | |
353 | } | |
354 | } | |
355 | ||
356 | make_btfn_unaryop!(Neg, neg); | |
357 | ||
358 | ||
359 | ||
360 | // | |
361 | // Mapping | |
362 | // | |
363 | ||
364 | impl<'a, F : Float, G, BT, V, const N : usize> Mapping<&'a Loc<F,N>> | |
365 | for BTFN<F, G, BT, N> | |
366 | where BT : BTImpl<F, N>, | |
367 | G : SupportGenerator<F, N, Id=BT::Data>, | |
368 | G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<&'a Loc<F,N>, Codomain = V>, | |
369 | V : Sum { | |
370 | ||
371 | type Codomain = V; | |
372 | ||
373 | fn value(&self, x : &'a Loc<F,N>) -> Self::Codomain { | |
374 | self.bt.iter_at(x).map(|&d| self.support_for(d).value(x)).sum() | |
375 | } | |
376 | } | |
377 | ||
378 | impl<F : Float, G, BT, V, const N : usize> Mapping<Loc<F,N>> | |
379 | for BTFN<F, G, BT, N> | |
380 | where BT : BTImpl<F, N>, | |
381 | G : SupportGenerator<F, N, Id=BT::Data>, | |
382 | G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F,N>, Codomain = V>, | |
383 | V : Sum { | |
384 | ||
385 | type Codomain = V; | |
386 | ||
387 | fn value(&self, x : Loc<F,N>) -> Self::Codomain { | |
388 | self.bt.iter_at(&x).map(|&d| self.support_for(d).value(x)).sum() | |
389 | } | |
390 | } | |
391 | ||
392 | impl<F : Float, G, BT, const N : usize> GlobalAnalysis<F, BT::Agg> | |
393 | for BTFN<F, G, BT, N> | |
394 | where BT : BTImpl<F, N>, | |
395 | G : SupportGenerator<F, N, Id=BT::Data>, | |
396 | G::SupportType : LocalAnalysis<F, BT::Agg, N> { | |
397 | ||
398 | #[inline] | |
399 | fn global_analysis(&self) -> BT::Agg { | |
400 | self.bt.global_analysis() | |
401 | } | |
402 | } | |
403 | ||
404 | // | |
405 | // Blanket implementation of BTFN as a linear functional over objects | |
406 | // that are linear functionals over BTFN. | |
407 | // | |
408 | ||
409 | impl<X, F : Float, G, BT, const N : usize> Linear<X> | |
410 | for BTFN<F, G, BT, N> | |
411 | where BT : BTImpl<F, N>, | |
412 | G : SupportGenerator<F, N, Id=BT::Data>, | |
413 | G::SupportType : LocalAnalysis<F, BT::Agg, N>, | |
414 | X : Linear<BTFN<F, G, BT, N>, Codomain=F> { | |
415 | type Codomain = F; | |
416 | ||
417 | #[inline] | |
418 | fn apply(&self, x : &X) -> F { | |
419 | x.apply(self) | |
420 | } | |
421 | } | |
422 | ||
423 | /// Helper trait for performing approximate minimisation using P2 elements. | |
424 | pub trait P2Minimise<U, F : Float> : Set<U> { | |
425 | fn p2_minimise<G : Fn(&U) -> F>(&self, g : G) -> (U, F); | |
426 | ||
427 | } | |
428 | ||
429 | impl<F : Float> P2Minimise<Loc<F, 1>, F> for Cube<F, 1> { | |
430 | fn p2_minimise<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> (Loc<F, 1>, F) { | |
431 | let interval = Simplex(self.corners()); | |
432 | interval.p2_model(&g).minimise(&interval) | |
433 | } | |
434 | } | |
435 | ||
436 | #[replace_float_literals(F::cast_from(literal))] | |
437 | impl<F : Float> P2Minimise<Loc<F, 2>, F> for Cube<F, 2> { | |
438 | fn p2_minimise<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> (Loc<F, 2>, F) { | |
439 | if false { | |
440 | // Split into two triangle (simplex) with separate P2 model in each. | |
441 | // The six nodes of each triangle are the corners and the edges. | |
442 | let [a, b, c, d] = self.corners(); | |
443 | let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)]; | |
444 | ||
445 | let ab = midpoint(&a, &b); | |
446 | let bc = midpoint(&b, &c); | |
447 | let ca = midpoint(&c, &a); | |
448 | let cd = midpoint(&c, &d); | |
449 | let da = midpoint(&d, &a); | |
450 | let [vab, vbc, vca, vcd, vda] = [g(&ab), g(&bc), g(&ca), g(&cd), g(&da)]; | |
451 | ||
452 | let s1 = Simplex([a, b, c]); | |
453 | let m1 = P2LocalModel::<F, 2, 3, 3, 3>::new( | |
454 | &[a, b, c, ab, bc, ca], | |
455 | &[va, vb, vc, vab, vbc, vca] | |
456 | ); | |
457 | ||
458 | let r1@(_, v1) = m1.minimise(&s1); | |
459 | ||
460 | let s2 = Simplex([c, d, a]); | |
461 | let m2 = P2LocalModel::<F, 2, 3, 3, 3>::new( | |
462 | &[c, d, a, cd, da, ca], | |
463 | &[vc, vd, va, vcd, vda, vca] | |
464 | ); | |
465 | ||
466 | let r2@(_, v2) = m2.minimise(&s2); | |
467 | ||
468 | if v1 < v2 { r1 } else { r2 } | |
469 | } else { | |
470 | // Single P2 model for the entire cube. | |
471 | let [a, b, c, d] = self.corners(); | |
472 | let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)]; | |
473 | let [e, f] = match 'r' { | |
474 | 'm' => [(&a + &b + &c) / 3.0, (&c + &d + &a) / 3.0], | |
475 | 'c' => [midpoint(&a, &b), midpoint(&a, &d)], | |
476 | 'w' => [(&a + &b * 2.0) / 3.0, (&a + &d * 2.0) / 3.0], | |
477 | 'r' => { | |
478 | // Pseudo-randomise edge midpoints | |
479 | let Loc([x, y]) = a; | |
480 | let tmp : f64 = (x+y).as_(); | |
481 | match tmp.to_bits() % 4 { | |
482 | 0 => [midpoint(&a, &b), midpoint(&a, &d)], | |
483 | 1 => [midpoint(&c, &d), midpoint(&a, &d)], | |
484 | 2 => [midpoint(&a, &b), midpoint(&b, &c)], | |
485 | _ => [midpoint(&c, &d), midpoint(&b, &c)], | |
486 | } | |
487 | }, | |
488 | _ => [self.center(), (&a + &b) / 2.0], | |
489 | }; | |
490 | let [ve, vf] = [g(&e), g(&f)]; | |
491 | ||
492 | let m1 = P2LocalModel::<F, 2, 3, 3, 3>::new( | |
493 | &[a, b, c, d, e, f], | |
494 | &[va, vb, vc, vd, ve, vf], | |
495 | ); | |
496 | ||
497 | m1.minimise(self) | |
498 | } | |
499 | } | |
500 | } | |
501 | ||
502 | struct RefineMax; | |
503 | struct RefineMin; | |
504 | ||
505 | /// TODO: update doc. | |
506 | /// Simple refiner that keeps the difference between upper and lower bounds | |
507 | /// of [`Bounded`] functions within $ε$. The functions have to also be continuous | |
508 | /// for this refiner to ever succeed. | |
509 | /// The type parameter `T` should be either [`RefineMax`] or [`RefineMin`]. | |
510 | struct P2Refiner<F : Float, T> { | |
511 | bound : Option<F>, | |
512 | tolerance : F, | |
513 | max_steps : usize, | |
514 | #[allow(dead_code)] // `how` is just for type system purposes. | |
515 | how : T, | |
516 | } | |
517 | ||
518 | impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> | |
519 | for P2Refiner<F, RefineMax> | |
520 | where Cube<F, N> : P2Minimise<Loc<F, N>, F>, | |
521 | G : SupportGenerator<F, N>, | |
522 | G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F> | |
523 | + LocalAnalysis<F, Bounds<F>, N> { | |
524 | type Result = Option<(Loc<F, N>, F)>; | |
525 | type Sorting = UpperBoundSorting<F>; | |
526 | ||
527 | fn refine( | |
528 | &self, | |
529 | aggregator : &Bounds<F>, | |
530 | cube : &Cube<F, N>, | |
531 | data : &Vec<G::Id>, | |
532 | generator : &G, | |
533 | step : usize | |
534 | ) -> RefinerResult<Bounds<F>, Self::Result> { | |
535 | // g gives the negative of the value of the function presented by `data` and `generator`. | |
536 | let g = move |x : &Loc<F,N>| { | |
537 | let f = move |&d| generator.support_for(d).value(x); | |
538 | -data.iter().map(f).sum::<F>() | |
539 | }; | |
540 | // … so the negative of the minimum is the maximm we want. | |
541 | let (x, neg_v) = cube.p2_minimise(g); | |
542 | let v = -neg_v; | |
543 | ||
544 | if self.bound.map_or(false, |b| aggregator.upper() <= b + self.tolerance) { | |
545 | // The upper bound is below the maximisation threshold. Don't bother with this cube. | |
546 | RefinerResult::Uncertain(*aggregator, None) | |
547 | } else if step < self.max_steps && (aggregator.upper() - v).abs() > self.tolerance { | |
548 | // The function isn't refined enough in `cube`, so return None | |
549 | // to indicate that further subdivision is required. | |
550 | RefinerResult::NeedRefinement | |
551 | } else { | |
552 | // The data is refined enough, so return new hopefully better bounds | |
553 | // and the maximiser. | |
554 | let res = (x, v); | |
555 | let bounds = Bounds(aggregator.lower(), v); | |
556 | RefinerResult::Uncertain(bounds, Some(res)) | |
557 | } | |
558 | } | |
559 | } | |
560 | ||
561 | ||
562 | impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> | |
563 | for P2Refiner<F, RefineMin> | |
564 | where Cube<F, N> : P2Minimise<Loc<F, N>, F>, | |
565 | G : SupportGenerator<F, N>, | |
566 | G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F> | |
567 | + LocalAnalysis<F, Bounds<F>, N> { | |
568 | type Result = Option<(Loc<F, N>, F)>; | |
569 | type Sorting = LowerBoundSorting<F>; | |
570 | ||
571 | fn refine( | |
572 | &self, | |
573 | aggregator : &Bounds<F>, | |
574 | cube : &Cube<F, N>, | |
575 | data : &Vec<G::Id>, | |
576 | generator : &G, | |
577 | step : usize | |
578 | ) -> RefinerResult<Bounds<F>, Self::Result> { | |
579 | // g gives the value of the function presented by `data` and `generator`. | |
580 | let g = move |x : &Loc<F,N>| { | |
581 | let f = move |&d| generator.support_for(d).value(x); | |
582 | data.iter().map(f).sum::<F>() | |
583 | }; | |
584 | // Minimise it. | |
585 | let (x, v) = cube.p2_minimise(g); | |
586 | ||
587 | if self.bound.map_or(false, |b| aggregator.lower() >= b - self.tolerance) { | |
588 | // The lower bound is above the minimisation threshold. Don't bother with this cube. | |
589 | RefinerResult::Uncertain(*aggregator, None) | |
590 | } else if step < self.max_steps && (aggregator.lower() - v).abs() > self.tolerance { | |
591 | // The function isn't refined enough in `cube`, so return None | |
592 | // to indicate that further subdivision is required. | |
593 | RefinerResult::NeedRefinement | |
594 | } else { | |
595 | // The data is refined enough, so return new hopefully better bounds | |
596 | // and the minimiser. | |
597 | let res = (x, v); | |
598 | let bounds = Bounds(v, aggregator.upper()); | |
599 | RefinerResult::Uncertain(bounds, Some(res)) | |
600 | } | |
601 | } | |
602 | } | |
603 | ||
604 | ||
605 | ||
606 | struct BoundRefiner<F : Float, T> { | |
607 | bound : F, | |
608 | tolerance : F, | |
609 | max_steps : usize, | |
610 | #[allow(dead_code)] // `how` is just for type system purposes. | |
611 | how : T, | |
612 | } | |
613 | ||
614 | impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> | |
615 | for BoundRefiner<F, RefineMax> | |
616 | where G : SupportGenerator<F, N> { | |
617 | type Result = bool; | |
618 | type Sorting = UpperBoundSorting<F>; | |
619 | ||
620 | fn refine( | |
621 | &self, | |
622 | aggregator : &Bounds<F>, | |
623 | _cube : &Cube<F, N>, | |
624 | _data : &Vec<G::Id>, | |
625 | _generator : &G, | |
626 | step : usize | |
627 | ) -> RefinerResult<Bounds<F>, Self::Result> { | |
628 | if aggregator.upper() <= self.bound + self.tolerance { | |
629 | // Below upper bound within tolerances. Indicate uncertain success. | |
630 | RefinerResult::Uncertain(*aggregator, true) | |
631 | } else if aggregator.lower() >= self.bound - self.tolerance { | |
632 | // Above upper bound within tolerances. Indicate certain failure. | |
633 | RefinerResult::Certain(false) | |
634 | } else if step < self.max_steps { | |
635 | // No decision possible, but within step bounds - further subdivision is required. | |
636 | RefinerResult::NeedRefinement | |
637 | } else { | |
638 | // No decision possible, but past step bounds | |
639 | RefinerResult::Uncertain(*aggregator, false) | |
640 | } | |
641 | } | |
642 | } | |
643 | ||
644 | ||
645 | impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N> | |
646 | where BT : BTSearch<F, N, Agg=Bounds<F>>, | |
647 | G : SupportGenerator<F, N, Id=BT::Data>, | |
648 | G::SupportType : for<'a> Mapping<&'a Loc<F,N>,Codomain=F> | |
649 | + LocalAnalysis<F, Bounds<F>, N>, | |
650 | Cube<F, N> : P2Minimise<Loc<F, N>, F> { | |
651 | ||
652 | /// Try to maximise the BTFN within stated value `tolerance`, taking at most `max_steps` | |
653 | /// refinement steps. Returns the maximiser and the maximum value. | |
654 | pub fn maximise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) { | |
655 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : None }; | |
656 | self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.").unwrap() | |
657 | } | |
658 | ||
659 | /// Try to maximise the BTFN within stated value `tolerance`, taking at most `max_steps` | |
660 | /// refinement steps. Returns the maximiser and the maximum value when one is found above | |
661 | /// the `bound` threshold, otherwise `None`. | |
662 | pub fn maximise_above(&mut self, bound : F, tolerance : F, max_steps : usize) | |
663 | -> Option<(Loc<F, N>, F)> { | |
664 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : Some(bound) }; | |
665 | self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.") | |
666 | } | |
667 | ||
668 | /// Try to minimise the BTFN within stated value `tolerance`, taking at most `max_steps` | |
669 | /// refinement steps. Returns the minimiser and the minimum value. | |
670 | pub fn minimise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) { | |
671 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : None }; | |
672 | self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.").unwrap() | |
673 | } | |
674 | ||
675 | /// Try to minimise the BTFN within stated value `tolerance`, taking at most `max_steps` | |
676 | /// refinement steps. Returns the minimiser and the minimum value when one is found below | |
677 | /// the `bound` threshold, otherwise `None`. | |
678 | pub fn minimise_below(&mut self, bound : F, tolerance : F, max_steps : usize) | |
679 | -> Option<(Loc<F, N>, F)> { | |
680 | let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : Some(bound) }; | |
681 | self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.") | |
682 | } | |
683 | /// Verify that the BTFN is has stated upper `bound` within stated `tolerance, taking at most | |
684 | /// `max_steps` refinment steps. | |
685 | pub fn has_upper_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool { | |
686 | let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMax }; | |
687 | self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.") | |
688 | } | |
689 | } |