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