Mercurial > repos > alg_tools / file comparison

comparison: src/bisection_tree/btfn.rs

src/bisection_tree/btfn.rs

branch
dev
changeset 96
962c8e346ab9
parent 63
f7b87d84864d
child 97
4e80fb049dca
equal deleted inserted replaced
95:9cb8225a3a41 96:962c8e346ab9
1 1 use crate::mapping::{
2 BasicDecomposition, DifferentiableImpl, DifferentiableMapping, Instance, Mapping, Space,
3 };
4 use crate::types::Float;
2 use numeric_literals::replace_float_literals; 5 use numeric_literals::replace_float_literals;
3 use std::iter::Sum; 6 use std::iter::Sum;
4 use std::marker::PhantomData; 7 use std::marker::PhantomData;
5 use std::sync::Arc; 8 use std::sync::Arc;
6 use crate::types::Float;
7 use crate::mapping::{
8 Instance, Mapping, DifferentiableImpl, DifferentiableMapping, Space,
9 BasicDecomposition,
10 };
11 //use crate::linops::{Apply, Linear}; 9 //use crate::linops::{Apply, Linear};
12 use crate::sets::Set; 10 use super::aggregator::*;
13 use crate::sets::Cube; 11 use super::bt::*;
14 use crate::loc::Loc; 12 use super::either::*;
13 use super::refine::*;
15 use super::support::*; 14 use super::support::*;
16 use super::bt::*;
17 use super::refine::*;
18 use super::aggregator::*;
19 use super::either::*;
20 use crate::fe_model::base::RealLocalModel; 15 use crate::fe_model::base::RealLocalModel;
21 use crate::fe_model::p2_local_model::*; 16 use crate::fe_model::p2_local_model::*;
22 17 use crate::loc::Loc;
23 /// Presentation for (mathematical) functions constructed as a sum of components functions with 18 use crate::sets::Cube;
19 use crate::sets::Set;
20
21 /// Presentation for (mathematical) functions constructed as a sum of components functions with
24 /// typically small support. 22 /// typically small support.
25 /// 23 ///
26 /// The domain of the function is [`Loc`]`<F, N>`, where `F` is the type of floating point numbers, 24 /// The domain of the function is [`Loc`]`<F, N>`, where `F` is the type of floating point numbers,
27 /// and `N` the dimension. 25 /// and `N` the dimension.
28 /// 26 ///
29 /// The `generator` lists the component functions that have to implement [`Support`]. 27 /// The `generator` lists the component functions that have to implement [`Support`].
30 /// Identifiers of the components ([`SupportGenerator::Id`], usually `usize`) are stored stored 28 /// Identifiers of the components ([`SupportGenerator::Id`], usually `usize`) are stored stored
31 /// in a [bisection tree][BTImpl], when one is provided as `bt`. However `bt` may also be `()` 29 /// in a [bisection tree][BTImpl], when one is provided as `bt`. However `bt` may also be `()`
32 /// for a [`PreBTFN`] that is only useful for vector space operations with a full [`BTFN`]. 30 /// for a [`PreBTFN`] that is only useful for vector space operations with a full [`BTFN`].
33 #[derive(Clone,Debug)] 31 #[derive(Clone, Debug)]
34 pub struct BTFN< 32 pub struct BTFN<F: Float, G: SupportGenerator<F, N>, BT /*: BTImpl<F, N>*/, const N: usize> /*where G::SupportType : LocalAnalysis<F, A, N>*/
35 F : Float, 33 {
36 G : SupportGenerator<F, N>, 34 bt: BT,
37 BT /*: BTImpl<F, N>*/, 35 generator: Arc<G>,
38 const N : usize 36 _phantoms: PhantomData<F>,
39 > /*where G::SupportType : LocalAnalysis<F, A, N>*/ { 37 }
40 bt : BT, 38
41 generator : Arc<G>, 39 impl<F: Float, G, BT, const N: usize> Space for BTFN<F, G, BT, N>
42 _phantoms : PhantomData<F>, 40 where
43 } 41 G: SupportGenerator<F, N, Id = BT::Data>,
44 42 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
45 impl<F : Float, G, BT, const N : usize> 43 BT: BTImpl<F, N>,
46 Space for BTFN<F, G, BT, N>
47 where
48 G : SupportGenerator<F, N, Id=BT::Data>,
49 G::SupportType : LocalAnalysis<F, BT::Agg, N>,
50 BT : BTImpl<F, N>
51 { 44 {
52 type Decomp = BasicDecomposition; 45 type Decomp = BasicDecomposition;
53 } 46 }
54 47
55 impl<F : Float, G, BT, const N : usize> 48 impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N>
56 BTFN<F, G, BT, N> 49 where
57 where 50 G: SupportGenerator<F, N, Id = BT::Data>,
58 G : SupportGenerator<F, N, Id=BT::Data>, 51 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
59 G::SupportType : LocalAnalysis<F, BT::Agg, N>, 52 BT: BTImpl<F, N>,
60 BT : BTImpl<F, N> 53 {
61 {
62
63 /// Create a new BTFN from a support generator and a pre-initialised bisection tree. 54 /// Create a new BTFN from a support generator and a pre-initialised bisection tree.
64 /// 55 ///
65 /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`. 56 /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`.
66 /// Use [`Self::construct`] if no preinitialised tree is available. Use [`Self::new_refresh`] 57 /// Use [`Self::construct`] if no preinitialised tree is available. Use [`Self::new_refresh`]
67 /// when the aggregators of the tree may need updates. 58 /// when the aggregators of the tree may need updates.
68 /// 59 ///
69 /// See the documentation for [`BTFN`] on the role of the `generator`. 60 /// See the documentation for [`BTFN`] on the role of the `generator`.
70 pub fn new(bt : BT, generator : G) -> Self { 61 pub fn new(bt: BT, generator: G) -> Self {
71 Self::new_arc(bt, Arc::new(generator)) 62 Self::new_arc(bt, Arc::new(generator))
72 } 63 }
73 64
74 fn new_arc(bt : BT, generator : Arc<G>) -> Self { 65 fn new_arc(bt: BT, generator: Arc<G>) -> Self {
75 BTFN { 66 BTFN {
76 bt : bt, 67 bt: bt,
77 generator : generator, 68 generator: generator,
78 _phantoms : std::marker::PhantomData, 69 _phantoms: std::marker::PhantomData,
79 } 70 }
80 } 71 }
81 72
82 /// Create a new BTFN support generator and a pre-initialised bisection tree, 73 /// Create a new BTFN support generator and a pre-initialised bisection tree,
83 /// cloning the tree and refreshing aggregators. 74 /// cloning the tree and refreshing aggregators.
84 /// 75 ///
85 /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`, but 76 /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`, but
86 /// the aggregator may be out of date. 77 /// the aggregator may be out of date.
87 /// 78 ///
88 /// See the documentation for [`BTFN`] on the role of the `generator`. 79 /// See the documentation for [`BTFN`] on the role of the `generator`.
89 pub fn new_refresh(bt : &BT, generator : G) -> Self { 80 pub fn new_refresh(bt: &BT, generator: G) -> Self {
90 // clone().refresh_aggregator(…) as opposed to convert_aggregator 81 // clone().refresh_aggregator(…) as opposed to convert_aggregator
91 // ensures that type is maintained. Due to Rc-pointer copy-on-write, 82 // ensures that type is maintained. Due to Rc-pointer copy-on-write,
92 // the effort is not significantly different. 83 // the effort is not significantly different.
93 let mut btnew = bt.clone(); 84 let mut btnew = bt.clone();
94 btnew.refresh_aggregator(&generator); 85 btnew.refresh_aggregator(&generator);
98 /// Create a new BTFN from a support generator, domain, and depth for a new [`BT`]. 89 /// Create a new BTFN from a support generator, domain, and depth for a new [`BT`].
99 /// 90 ///
100 /// The top node of the created [`BT`] will have the given `domain`. 91 /// The top node of the created [`BT`] will have the given `domain`.
101 /// 92 ///
102 /// See the documentation for [`BTFN`] on the role of the `generator`. 93 /// See the documentation for [`BTFN`] on the role of the `generator`.
103 pub fn construct(domain : Cube<F, N>, depth : BT::Depth, generator : G) -> Self { 94 pub fn construct(domain: Cube<F, N>, depth: BT::Depth, generator: G) -> Self {
104 Self::construct_arc(domain, depth, Arc::new(generator)) 95 Self::construct_arc(domain, depth, Arc::new(generator))
105 } 96 }
106 97
107 fn construct_arc(domain : Cube<F, N>, depth : BT::Depth, generator : Arc<G>) -> Self { 98 fn construct_arc(domain: Cube<F, N>, depth: BT::Depth, generator: Arc<G>) -> Self {
108 let mut bt = BT::new(domain, depth); 99 let mut bt = BT::new(domain, depth);
109 for (d, support) in generator.all_data() { 100 for (d, support) in generator.all_data() {
110 bt.insert(d, &support); 101 bt.insert(d, &support);
111 } 102 }
112 Self::new_arc(bt, generator) 103 Self::new_arc(bt, generator)
115 /// Convert the aggregator of the [`BTFN`] to a different one. 106 /// Convert the aggregator of the [`BTFN`] to a different one.
116 /// 107 ///
117 /// This will construct a [`BTFN`] with the same components and generator as the (consumed) 108 /// This will construct a [`BTFN`] with the same components and generator as the (consumed)
118 /// `self`, but a new `BT` with [`Aggregator`]s of type `ANew`. 109 /// `self`, but a new `BT` with [`Aggregator`]s of type `ANew`.
119 pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N> 110 pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N>
120 where ANew : Aggregator, 111 where
121 G : SupportGenerator<F, N, Id=BT::Data>, 112 ANew: Aggregator,
122 G::SupportType : LocalAnalysis<F, ANew, N> { 113 G: SupportGenerator<F, N, Id = BT::Data>,
114 G::SupportType: LocalAnalysis<F, ANew, N>,
115 {
123 BTFN::new_arc(self.bt.convert_aggregator(&*self.generator), self.generator) 116 BTFN::new_arc(self.bt.convert_aggregator(&*self.generator), self.generator)
124 } 117 }
125 118
126 /// Change the generator (after, e.g., a scaling of the latter). 119 /// Change the generator (after, e.g., a scaling of the latter).
127 fn new_generator(&self, generator : G) -> Self { 120 fn new_generator(&self, generator: G) -> Self {
128 BTFN::new_refresh(&self.bt, generator) 121 BTFN::new_refresh(&self.bt, generator)
129 } 122 }
130 123
131 /// Refresh aggregator after updates to generator 124 /// Refresh aggregator after updates to generator
132 fn refresh_aggregator(&mut self) { 125 fn refresh_aggregator(&mut self) {
133 self.bt.refresh_aggregator(&*self.generator); 126 self.bt.refresh_aggregator(&*self.generator);
134 } 127 }
135 128 }
136 } 129
137 130 impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N>
138 impl<F : Float, G, BT, const N : usize> 131 where
139 BTFN<F, G, BT, N> 132 G: SupportGenerator<F, N>,
140 where G : SupportGenerator<F, N> { 133 {
141 /// Change the [bisection tree][BTImpl] of the [`BTFN`] to a different one. 134 /// Change the [bisection tree][BTImpl] of the [`BTFN`] to a different one.
142 /// 135 ///
143 /// This can be used to convert a [`PreBTFN`] to a full [`BTFN`], or the change 136 /// This can be used to convert a [`PreBTFN`] to a full [`BTFN`], or the change
144 /// the aggreagator; see also [`self.convert_aggregator`]. 137 /// the aggreagator; see also [`self.convert_aggregator`].
145 pub fn instantiate< 138 pub fn instantiate<BTNew: BTImpl<F, N, Data = G::Id>>(
146 BTNew : BTImpl<F, N, Data=G::Id>, 139 self,
147 > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<F, G, BTNew, N> 140 domain: Cube<F, N>,
148 where G::SupportType : LocalAnalysis<F, BTNew::Agg, N> { 141 depth: BTNew::Depth,
142 ) -> BTFN<F, G, BTNew, N>
143 where
144 G::SupportType: LocalAnalysis<F, BTNew::Agg, N>,
145 {
149 BTFN::construct_arc(domain, depth, self.generator) 146 BTFN::construct_arc(domain, depth, self.generator)
150 } 147 }
151 } 148 }
152 149
153 /// A BTFN with no bisection tree. 150 /// A BTFN with no bisection tree.
154 /// 151 ///
155 /// Most BTFN methods are not available, but if a BTFN is going to be summed with another 152 /// Most BTFN methods are not available, but if a BTFN is going to be summed with another
156 /// before other use, it will be more efficient to not construct an unnecessary bisection tree 153 /// before other use, it will be more efficient to not construct an unnecessary bisection tree
157 /// that would be shortly dropped. 154 /// that would be shortly dropped.
158 pub type PreBTFN<F, G, const N : usize> = BTFN<F, G, (), N>; 155 pub type PreBTFN<F, G, const N: usize> = BTFN<F, G, (), N>;
159 156
160 impl<F : Float, G, const N : usize> PreBTFN<F, G, N> where G : SupportGenerator<F, N> { 157 impl<F: Float, G, const N: usize> PreBTFN<F, G, N>
161 158 where
159 G: SupportGenerator<F, N>,
160 {
162 /// Create a new [`PreBTFN`] with no bisection tree. 161 /// Create a new [`PreBTFN`] with no bisection tree.
163 pub fn new_pre(generator : G) -> Self { 162 pub fn new_pre(generator: G) -> Self {
164 BTFN { 163 BTFN {
165 bt : (), 164 bt: (),
166 generator : Arc::new(generator), 165 generator: Arc::new(generator),
167 _phantoms : std::marker::PhantomData, 166 _phantoms: std::marker::PhantomData,
168 } 167 }
169 } 168 }
170 } 169 }
171 170
172 impl<F : Float, G, BT, const N : usize> 171 impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N>
173 BTFN<F, G, BT, N> 172 where
174 where G : SupportGenerator<F, N, Id=usize>, 173 G: SupportGenerator<F, N, Id = usize>,
175 G::SupportType : LocalAnalysis<F, BT::Agg, N>, 174 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
176 BT : BTImpl<F, N, Data=usize> { 175 BT: BTImpl<F, N, Data = usize>,
177 176 {
178 /// Helper function for implementing [`std::ops::Add`]. 177 /// Helper function for implementing [`std::ops::Add`].
179 fn add_another<G2>(&self, g2 : Arc<G2>) -> BTFN<F, BothGenerators<G, G2>, BT, N> 178 fn add_another<G2>(&self, g2: Arc<G2>) -> BTFN<F, BothGenerators<G, G2>, BT, N>
180 where G2 : SupportGenerator<F, N, Id=usize>, 179 where
181 G2::SupportType : LocalAnalysis<F, BT::Agg, N> { 180 G2: SupportGenerator<F, N, Id = usize>,
182 181 G2::SupportType: LocalAnalysis<F, BT::Agg, N>,
182 {
183 let mut bt = self.bt.clone(); 183 let mut bt = self.bt.clone();
184 let both = BothGenerators(Arc::clone(&self.generator), g2); 184 let both = BothGenerators(Arc::clone(&self.generator), g2);
185 185
186 for (d, support) in both.all_right_data() { 186 for (d, support) in both.all_right_data() {
187 bt.insert(d, &support); 187 bt.insert(d, &support);
188 } 188 }
189 189
190 BTFN { 190 BTFN {
191 bt : bt, 191 bt: bt,
192 generator : Arc::new(both), 192 generator: Arc::new(both),
193 _phantoms : std::marker::PhantomData, 193 _phantoms: std::marker::PhantomData,
194 } 194 }
195 } 195 }
196 } 196 }
197 197
198 macro_rules! make_btfn_add { 198 macro_rules! make_btfn_add {
229 } 229 }
230 } 230 }
231 } 231 }
232 232
233 make_btfn_add!(BTFN<F, G1, BT1, N>, std::convert::identity, ); 233 make_btfn_add!(BTFN<F, G1, BT1, N>, std::convert::identity, );
234 make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone, ); 234 make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone,);
235 235
236 macro_rules! make_btfn_sub { 236 macro_rules! make_btfn_sub {
237 ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => { 237 ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => {
238 impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> 238 impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize>
239 std::ops::Sub<BTFN<F, G2, BT2, N>> for 239 std::ops::Sub<BTFN<F, G2, BT2, N>> for
272 } 272 }
273 } 273 }
274 } 274 }
275 275
276 make_btfn_sub!(BTFN<F, G1, BT1, N>, std::convert::identity, ); 276 make_btfn_sub!(BTFN<F, G1, BT1, N>, std::convert::identity, );
277 make_btfn_sub!(&'a BTFN<F, G1, BT1, N>, std::convert::identity, ); 277 make_btfn_sub!(&'a BTFN<F, G1, BT1, N>, std::convert::identity,);
278 278
279 macro_rules! make_btfn_scalarop_rhs { 279 macro_rules! make_btfn_scalarop_rhs {
280 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { 280 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
281 impl<F : Float, G, BT, const N : usize> 281 impl<F: Float, G, BT, const N: usize> std::ops::$trait_assign<F> for BTFN<F, G, BT, N>
282 std::ops::$trait_assign<F> 282 where
283 for BTFN<F, G, BT, N> 283 BT: BTImpl<F, N>,
284 where BT : BTImpl<F, N>, 284 G: SupportGenerator<F, N, Id = BT::Data>,
285 G : SupportGenerator<F, N, Id=BT::Data>, 285 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
286 G::SupportType : LocalAnalysis<F, BT::Agg, N> { 286 {
287 #[inline] 287 #[inline]
288 fn $fn_assign(&mut self, t : F) { 288 fn $fn_assign(&mut self, t: F) {
289 Arc::make_mut(&mut self.generator).$fn_assign(t); 289 Arc::make_mut(&mut self.generator).$fn_assign(t);
290 self.refresh_aggregator(); 290 self.refresh_aggregator();
291 } 291 }
292 } 292 }
293 293
294 impl<F : Float, G, BT, const N : usize> 294 impl<F: Float, G, BT, const N: usize> std::ops::$trait<F> for BTFN<F, G, BT, N>
295 std::ops::$trait<F> 295 where
296 for BTFN<F, G, BT, N> 296 BT: BTImpl<F, N>,
297 where BT : BTImpl<F, N>, 297 G: SupportGenerator<F, N, Id = BT::Data>,
298 G : SupportGenerator<F, N, Id=BT::Data>, 298 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
299 G::SupportType : LocalAnalysis<F, BT::Agg, N> { 299 {
300 type Output = Self; 300 type Output = Self;
301 #[inline] 301 #[inline]
302 fn $fn(mut self, t : F) -> Self::Output { 302 fn $fn(mut self, t: F) -> Self::Output {
303 Arc::make_mut(&mut self.generator).$fn_assign(t); 303 Arc::make_mut(&mut self.generator).$fn_assign(t);
304 self.refresh_aggregator(); 304 self.refresh_aggregator();
305 self 305 self
306 } 306 }
307 } 307 }
308 308
309 impl<'a, F : Float, G, BT, const N : usize> 309 impl<'a, F: Float, G, BT, const N: usize> std::ops::$trait<F> for &'a BTFN<F, G, BT, N>
310 std::ops::$trait<F> 310 where
311 for &'a BTFN<F, G, BT, N> 311 BT: BTImpl<F, N>,
312 where BT : BTImpl<F, N>, 312 G: SupportGenerator<F, N, Id = BT::Data>,
313 G : SupportGenerator<F, N, Id=BT::Data>, 313 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
314 G::SupportType : LocalAnalysis<F, BT::Agg, N>, 314 &'a G: std::ops::$trait<F, Output = G>,
315 &'a G : std::ops::$trait<F,Output=G> { 315 {
316 type Output = BTFN<F, G, BT, N>; 316 type Output = BTFN<F, G, BT, N>;
317 #[inline] 317 #[inline]
318 fn $fn(self, t : F) -> Self::Output { 318 fn $fn(self, t: F) -> Self::Output {
319 self.new_generator(self.generator.$fn(t)) 319 self.new_generator(self.generator.$fn(t))
320 } 320 }
321 } 321 }
322 } 322 };
323 } 323 }
324 324
325 make_btfn_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); 325 make_btfn_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
326 make_btfn_scalarop_rhs!(Div, div, DivAssign, div_assign); 326 make_btfn_scalarop_rhs!(Div, div, DivAssign, div_assign);
327 327
366 make_btfn_scalarop_lhs!(Mul, mul, mul_assign, f32 f64); 366 make_btfn_scalarop_lhs!(Mul, mul, mul_assign, f32 f64);
367 make_btfn_scalarop_lhs!(Div, div, div_assign, f32 f64); 367 make_btfn_scalarop_lhs!(Div, div, div_assign, f32 f64);
368 368
369 macro_rules! make_btfn_unaryop { 369 macro_rules! make_btfn_unaryop {
370 ($trait:ident, $fn:ident) => { 370 ($trait:ident, $fn:ident) => {
371 impl<F : Float, G, BT, const N : usize> 371 impl<F: Float, G, BT, const N: usize> std::ops::$trait for BTFN<F, G, BT, N>
372 std::ops::$trait 372 where
373 for BTFN<F, G, BT, N> 373 BT: BTImpl<F, N>,
374 where BT : BTImpl<F, N>, 374 G: SupportGenerator<F, N, Id = BT::Data>,
375 G : SupportGenerator<F, N, Id=BT::Data>, 375 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
376 G::SupportType : LocalAnalysis<F, BT::Agg, N> { 376 {
377 type Output = Self; 377 type Output = Self;
378 #[inline] 378 #[inline]
379 fn $fn(mut self) -> Self::Output { 379 fn $fn(mut self) -> Self::Output {
380 self.generator = Arc::new(Arc::unwrap_or_clone(self.generator).$fn()); 380 self.generator = Arc::new(Arc::unwrap_or_clone(self.generator).$fn());
381 self.refresh_aggregator(); 381 self.refresh_aggregator();
394 #[inline] 394 #[inline]
395 fn $fn(self) -> Self::Output { 395 fn $fn(self) -> Self::Output {
396 self.new_generator(std::ops::$trait::$fn(&self.generator)) 396 self.new_generator(std::ops::$trait::$fn(&self.generator))
397 } 397 }
398 }*/ 398 }*/
399 } 399 };
400 } 400 }
401 401
402 make_btfn_unaryop!(Neg, neg); 402 make_btfn_unaryop!(Neg, neg);
403 403
404 // 404 //
405 // Apply, Mapping, Differentiate 405 // Apply, Mapping, Differentiate
406 // 406 //
407 407
408 impl<F : Float, G, BT, V, const N : usize> Mapping<Loc<F, N>> 408 impl<F: Float, G, BT, V, const N: usize> Mapping<Loc<F, N>> for BTFN<F, G, BT, N>
409 for BTFN<F, G, BT, N> 409 where
410 where 410 BT: BTImpl<F, N>,
411 BT : BTImpl<F, N>, 411 G: SupportGenerator<F, N, Id = BT::Data>,
412 G : SupportGenerator<F, N, Id=BT::Data>, 412 G::SupportType: LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F, N>, Codomain = V>,
413 G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F, N>, Codomain = V>, 413 V: Sum + Space,
414 V : Sum + Space, 414 {
415 {
416
417 type Codomain = V; 415 type Codomain = V;
418 416
419 fn apply<I : Instance<Loc<F,N>>>(&self, x : I) -> Self::Codomain { 417 fn apply<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Codomain {
420 let xc = x.cow(); 418 let xc = x.cow();
421 self.bt.iter_at(&*xc) 419 self.bt
422 .map(|&d| self.generator.support_for(d).apply(&*xc)).sum() 420 .iter_at(&*xc)
423 } 421 .map(|&d| self.generator.support_for(d).apply(&*xc))
424 } 422 .sum()
425 423 }
426 impl<F : Float, G, BT, V, const N : usize> DifferentiableImpl<Loc<F, N>> 424 }
427 for BTFN<F, G, BT, N> 425
428 where 426 impl<F: Float, G, BT, V, const N: usize> DifferentiableImpl<Loc<F, N>> for BTFN<F, G, BT, N>
429 BT : BTImpl<F, N>, 427 where
430 G : SupportGenerator<F, N, Id=BT::Data>, 428 BT: BTImpl<F, N>,
431 G::SupportType : LocalAnalysis<F, BT::Agg, N> 429 G: SupportGenerator<F, N, Id = BT::Data>,
432 + DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>, 430 G::SupportType:
433 V : Sum + Space, 431 LocalAnalysis<F, BT::Agg, N> + DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>,
434 { 432 V: Sum + Space,
435 433 {
436 type Derivative = V; 434 type Derivative = V;
437 435
438 fn differential_impl<I : Instance<Loc<F, N>>>(&self, x :I) -> Self::Derivative { 436 fn differential_impl<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Derivative {
439 let xc = x.cow(); 437 let xc = x.cow();
440 self.bt.iter_at(&*xc) 438 self.bt
439 .iter_at(&*xc)
441 .map(|&d| self.generator.support_for(d).differential(&*xc)) 440 .map(|&d| self.generator.support_for(d).differential(&*xc))
442 .sum() 441 .sum()
443 } 442 }
444 } 443 }
445 444
446 // 445 //
447 // GlobalAnalysis 446 // GlobalAnalysis
448 // 447 //
449 448
450 impl<F : Float, G, BT, const N : usize> GlobalAnalysis<F, BT::Agg> 449 impl<F: Float, G, BT, const N: usize> GlobalAnalysis<F, BT::Agg> for BTFN<F, G, BT, N>
451 for BTFN<F, G, BT, N> 450 where
452 where BT : BTImpl<F, N>, 451 BT: BTImpl<F, N>,
453 G : SupportGenerator<F, N, Id=BT::Data>, 452 G: SupportGenerator<F, N, Id = BT::Data>,
454 G::SupportType : LocalAnalysis<F, BT::Agg, N> { 453 G::SupportType: LocalAnalysis<F, BT::Agg, N>,
455 454 {
456 #[inline] 455 #[inline]
457 fn global_analysis(&self) -> BT::Agg { 456 fn global_analysis(&self) -> BT::Agg {
458 self.bt.global_analysis() 457 self.bt.global_analysis()
459 } 458 }
460 } 459 }
503 502
504 /// Helper trait for performing approximate minimisation using P2 elements. 503 /// Helper trait for performing approximate minimisation using P2 elements.
505 /// 504 ///
506 /// `U` is the domain, generally [`Loc`]`<F, N>`, and `F` the type of floating point numbers. 505 /// `U` is the domain, generally [`Loc`]`<F, N>`, and `F` the type of floating point numbers.
507 /// `Self` is generally a set of `U`, for example, [`Cube`]`<F, N>`. 506 /// `Self` is generally a set of `U`, for example, [`Cube`]`<F, N>`.
508 pub trait P2Minimise<U : Space, F : Float> : Set<U> { 507 pub trait P2Minimise<U: Space, F: Float>: Set<U> {
509 /// Minimise `g` over the set presented by `Self`. 508 /// Minimise `g` over the set presented by `Self`.
510 /// 509 ///
511 /// The function returns `(x, v)` where `x` is the minimiser `v` an approximation of `g(x)`. 510 /// The function returns `(x, v)` where `x` is the minimiser `v` an approximation of `g(x)`.
512 fn p2_minimise<G : Fn(&U) -> F>(&self, g : G) -> (U, F); 511 fn p2_minimise<G: Fn(&U) -> F>(&self, g: G) -> (U, F);
513 512 }
514 } 513
515 514 impl<F: Float> P2Minimise<Loc<F, 1>, F> for Cube<F, 1> {
516 impl<F : Float> P2Minimise<Loc<F, 1>, F> for Cube<F, 1> { 515 fn p2_minimise<G: Fn(&Loc<F, 1>) -> F>(&self, g: G) -> (Loc<F, 1>, F) {
517 fn p2_minimise<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> (Loc<F, 1>, F) {
518 let interval = Simplex(self.corners()); 516 let interval = Simplex(self.corners());
519 interval.p2_model(&g).minimise(&interval) 517 interval.p2_model(&g).minimise(&interval)
520 } 518 }
521 } 519 }
522 520
523 #[replace_float_literals(F::cast_from(literal))] 521 #[replace_float_literals(F::cast_from(literal))]
524 impl<F : Float> P2Minimise<Loc<F, 2>, F> for Cube<F, 2> { 522 impl<F: Float> P2Minimise<Loc<F, 2>, F> for Cube<F, 2> {
525 fn p2_minimise<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> (Loc<F, 2>, F) { 523 fn p2_minimise<G: Fn(&Loc<F, 2>) -> F>(&self, g: G) -> (Loc<F, 2>, F) {
526 if false { 524 if false {
527 // Split into two triangle (simplex) with separate P2 model in each. 525 // Split into two triangle (simplex) with separate P2 model in each.
528 // The six nodes of each triangle are the corners and the edges. 526 // The six nodes of each triangle are the corners and the edges.
529 let [a, b, c, d] = self.corners(); 527 let [a, b, c, d] = self.corners();
530 let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)]; 528 let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)];
535 let cd = midpoint(&c, &d); 533 let cd = midpoint(&c, &d);
536 let da = midpoint(&d, &a); 534 let da = midpoint(&d, &a);
537 let [vab, vbc, vca, vcd, vda] = [g(&ab), g(&bc), g(&ca), g(&cd), g(&da)]; 535 let [vab, vbc, vca, vcd, vda] = [g(&ab), g(&bc), g(&ca), g(&cd), g(&da)];
538 536
539 let s1 = Simplex([a, b, c]); 537 let s1 = Simplex([a, b, c]);
540 let m1 = P2LocalModel::<F, 2, 3>::new( 538 let m1 =
541 &[a, b, c, ab, bc, ca], 539 P2LocalModel::<F, 2, 3>::new(&[a, b, c, ab, bc, ca], &[va, vb, vc, vab, vbc, vca]);
542 &[va, vb, vc, vab, vbc, vca] 540
543 ); 541 let r1 @ (_, v1) = m1.minimise(&s1);
544
545 let r1@(_, v1) = m1.minimise(&s1);
546 542
547 let s2 = Simplex([c, d, a]); 543 let s2 = Simplex([c, d, a]);
548 let m2 = P2LocalModel::<F, 2, 3>::new( 544 let m2 =
549 &[c, d, a, cd, da, ca], 545 P2LocalModel::<F, 2, 3>::new(&[c, d, a, cd, da, ca], &[vc, vd, va, vcd, vda, vca]);
550 &[vc, vd, va, vcd, vda, vca] 546
551 ); 547 let r2 @ (_, v2) = m2.minimise(&s2);
552 548
553 let r2@(_, v2) = m2.minimise(&s2); 549 if v1 < v2 {
554 550 r1
555 if v1 < v2 { r1 } else { r2 } 551 } else {
552 r2
553 }
556 } else { 554 } else {
557 // Single P2 model for the entire cube. 555 // Single P2 model for the entire cube.
558 let [a, b, c, d] = self.corners(); 556 let [a, b, c, d] = self.corners();
559 let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)]; 557 let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)];
560 let [e, f] = match 'r' { 558 let [e, f] = match 'r' {
561 'm' => [(&a + &b + &c) / 3.0, (&c + &d + &a) / 3.0], 559 'm' => [(&a + &b + &c) / 3.0, (&c + &d + &a) / 3.0],
562 'c' => [midpoint(&a, &b), midpoint(&a, &d)], 560 'c' => [midpoint(&a, &b), midpoint(&a, &d)],
563 'w' => [(&a + &b * 2.0) / 3.0, (&a + &d * 2.0) / 3.0], 561 'w' => [(&a + &b * 2.0) / 3.0, (&a + &d * 2.0) / 3.0],
564 'r' => { 562 'r' => {
565 // Pseudo-randomise edge midpoints 563 // Pseudo-randomise edge midpoints
566 let Loc([x, y]) = a; 564 let Loc([x, y]) = a;
567 let tmp : f64 = (x+y).as_(); 565 let tmp: f64 = (x + y).as_();
568 match tmp.to_bits() % 4 { 566 match tmp.to_bits() % 4 {
569 0 => [midpoint(&a, &b), midpoint(&a, &d)], 567 0 => [midpoint(&a, &b), midpoint(&a, &d)],
570 1 => [midpoint(&c, &d), midpoint(&a, &d)], 568 1 => [midpoint(&c, &d), midpoint(&a, &d)],
571 2 => [midpoint(&a, &b), midpoint(&b, &c)], 569 2 => [midpoint(&a, &b), midpoint(&b, &c)],
572 _ => [midpoint(&c, &d), midpoint(&b, &c)], 570 _ => [midpoint(&c, &d), midpoint(&b, &c)],
573 } 571 }
574 }, 572 }
575 _ => [self.center(), (&a + &b) / 2.0], 573 _ => [self.center(), (&a + &b) / 2.0],
576 }; 574 };
577 let [ve, vf] = [g(&e), g(&f)]; 575 let [ve, vf] = [g(&e), g(&f)];
578 576
579 let m1 = P2LocalModel::<F, 2, 3>::new( 577 let m1 = P2LocalModel::<F, 2, 3>::new(&[a, b, c, d, e, f], &[va, vb, vc, vd, ve, vf]);
580 &[a, b, c, d, e, f],
581 &[va, vb, vc, vd, ve, vf],
582 );
583 578
584 m1.minimise(self) 579 m1.minimise(self)
585 } 580 }
586 } 581 }
587 } 582 }
593 struct RefineMin; 588 struct RefineMin;
594 589
595 /// A bisection tree [`Refiner`] for maximising or minimising a [`BTFN`]. 590 /// A bisection tree [`Refiner`] for maximising or minimising a [`BTFN`].
596 /// 591 ///
597 /// The type parameter `T` should be either [`RefineMax`] or [`RefineMin`]. 592 /// The type parameter `T` should be either [`RefineMax`] or [`RefineMin`].
598 struct P2Refiner<F : Float, T> { 593 struct P2Refiner<F: Float, T> {
599 /// The maximum / minimum should be above / below this threshold. 594 /// The maximum / minimum should be above / below this threshold.
600 /// If the threshold cannot be satisfied, the refiner will return `None`. 595 /// If the threshold cannot be satisfied, the refiner will return `None`.
601 bound : Option<F>, 596 bound: Option<F>,
602 /// Tolerance for function value estimation. 597 /// Tolerance for function value estimation.
603 tolerance : F, 598 tolerance: F,
604 /// Maximum number of steps to execute the refiner for 599 /// Maximum number of steps to execute the refiner for
605 max_steps : usize, 600 max_steps: usize,
606 /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes. 601 /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes.
607 #[allow(dead_code)] // `how` is just for type system purposes. 602 #[allow(dead_code)] // `how` is just for type system purposes.
608 how : T, 603 how: T,
609 } 604 }
610 605
611 impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> 606 impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for P2Refiner<F, RefineMax>
612 for P2Refiner<F, RefineMax> 607 where
613 where Cube<F, N> : P2Minimise<Loc<F, N>, F>, 608 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
614 G : SupportGenerator<F, N>, 609 G: SupportGenerator<F, N>,
615 G::SupportType : Mapping<Loc<F, N>, Codomain=F> 610 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
616 + LocalAnalysis<F, Bounds<F>, N> { 611 {
617 type Result = Option<(Loc<F, N>, F)>; 612 type Result = Option<(Loc<F, N>, F)>;
618 type Sorting = UpperBoundSorting<F>; 613 type Sorting = UpperBoundSorting<F>;
619 614
620 fn refine( 615 fn refine(
621 &self, 616 &self,
622 aggregator : &Bounds<F>, 617 aggregator: &Bounds<F>,
623 cube : &Cube<F, N>, 618 cube: &Cube<F, N>,
624 data : &[G::Id], 619 data: &[G::Id],
625 generator : &G, 620 generator: &G,
626 step : usize 621 step: usize,
627 ) -> RefinerResult<Bounds<F>, Self::Result> { 622 ) -> RefinerResult<Bounds<F>, Self::Result> {
628 623 if self
629 if self.bound.map_or(false, |b| aggregator.upper() <= b + self.tolerance) { 624 .bound
625 .map_or(false, |b| aggregator.upper() <= b + self.tolerance)
626 {
630 // The upper bound is below the maximisation threshold. Don't bother with this cube. 627 // The upper bound is below the maximisation threshold. Don't bother with this cube.
631 return RefinerResult::Uncertain(*aggregator, None) 628 return RefinerResult::Uncertain(*aggregator, None);
632 } 629 }
633 630
634 // g gives the negative of the value of the function presented by `data` and `generator`. 631 // g gives the negative of the value of the function presented by `data` and `generator`.
635 let g = move |x : &Loc<F, N>| { 632 let g = move |x: &Loc<F, N>| {
636 let f = move |&d| generator.support_for(d).apply(x); 633 let f = move |&d| generator.support_for(d).apply(x);
637 -data.iter().map(f).sum::<F>() 634 -data.iter().map(f).sum::<F>()
638 }; 635 };
639 // … so the negative of the minimum is the maximm we want. 636 // … so the negative of the minimum is the maximm we want.
640 let (x, _neg_v) = cube.p2_minimise(g); 637 let (x, _neg_v) = cube.p2_minimise(g);
641 //let v = -neg_v; 638 //let v = -neg_v;
642 let v = -g(&x); 639 let v = -g(&x);
643 640
644 if step < self.max_steps && (aggregator.upper() > v + self.tolerance 641 if step < self.max_steps
645 /*|| aggregator.lower() > v - self.tolerance*/) { 642 && (aggregator.upper() > v + self.tolerance/*|| aggregator.lower() > v - self.tolerance*/)
643 {
646 // The function isn't refined enough in `cube`, so return None 644 // The function isn't refined enough in `cube`, so return None
647 // to indicate that further subdivision is required. 645 // to indicate that further subdivision is required.
648 RefinerResult::NeedRefinement 646 RefinerResult::NeedRefinement
649 } else { 647 } else {
650 // The data is refined enough, so return new hopefully better bounds 648 // The data is refined enough, so return new hopefully better bounds
653 let bounds = Bounds(v, v); 651 let bounds = Bounds(v, v);
654 RefinerResult::Uncertain(bounds, Some(res)) 652 RefinerResult::Uncertain(bounds, Some(res))
655 } 653 }
656 } 654 }
657 655
658 fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { 656 fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
659 match (*r1, r2) { 657 match (*r1, r2) {
660 (Some((_, v1)), Some((_, v2))) => if v1 < v2 { *r1 = r2 } 658 (Some((_, v1)), Some((_, v2))) => {
659 if v1 < v2 {
660 *r1 = r2
661 }
662 }
661 (None, Some(_)) => *r1 = r2, 663 (None, Some(_)) => *r1 = r2,
662 (_, _) => {}, 664 (_, _) => {}
663 } 665 }
664 } 666 }
665 } 667 }
666 668
667 669 impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for P2Refiner<F, RefineMin>
668 impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> 670 where
669 for P2Refiner<F, RefineMin> 671 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
670 where Cube<F, N> : P2Minimise<Loc<F, N>, F>, 672 G: SupportGenerator<F, N>,
671 G : SupportGenerator<F, N>, 673 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
672 G::SupportType : Mapping<Loc<F, N>, Codomain=F> 674 {
673 + LocalAnalysis<F, Bounds<F>, N> {
674 type Result = Option<(Loc<F, N>, F)>; 675 type Result = Option<(Loc<F, N>, F)>;
675 type Sorting = LowerBoundSorting<F>; 676 type Sorting = LowerBoundSorting<F>;
676 677
677 fn refine( 678 fn refine(
678 &self, 679 &self,
679 aggregator : &Bounds<F>, 680 aggregator: &Bounds<F>,
680 cube : &Cube<F, N>, 681 cube: &Cube<F, N>,
681 data : &[G::Id], 682 data: &[G::Id],
682 generator : &G, 683 generator: &G,
683 step : usize 684 step: usize,
684 ) -> RefinerResult<Bounds<F>, Self::Result> { 685 ) -> RefinerResult<Bounds<F>, Self::Result> {
685 686 if self
686 if self.bound.map_or(false, |b| aggregator.lower() >= b - self.tolerance) { 687 .bound
688 .map_or(false, |b| aggregator.lower() >= b - self.tolerance)
689 {
687 // The lower bound is above the minimisation threshold. Don't bother with this cube. 690 // The lower bound is above the minimisation threshold. Don't bother with this cube.
688 return RefinerResult::Uncertain(*aggregator, None) 691 return RefinerResult::Uncertain(*aggregator, None);
689 } 692 }
690 693
691 // g gives the value of the function presented by `data` and `generator`. 694 // g gives the value of the function presented by `data` and `generator`.
692 let g = move |x : &Loc<F, N>| { 695 let g = move |x: &Loc<F, N>| {
693 let f = move |&d| generator.support_for(d).apply(x); 696 let f = move |&d| generator.support_for(d).apply(x);
694 data.iter().map(f).sum::<F>() 697 data.iter().map(f).sum::<F>()
695 }; 698 };
696 // Minimise it. 699 // Minimise it.
697 let (x, _v) = cube.p2_minimise(g); 700 let (x, _v) = cube.p2_minimise(g);
698 let v = g(&x); 701 let v = g(&x);
699 702
700 if step < self.max_steps && (aggregator.lower() < v - self.tolerance 703 if step < self.max_steps
701 /*|| aggregator.upper() < v + self.tolerance*/) { 704 && (aggregator.lower() < v - self.tolerance/*|| aggregator.upper() < v + self.tolerance*/)
705 {
702 // The function isn't refined enough in `cube`, so return None 706 // The function isn't refined enough in `cube`, so return None
703 // to indicate that further subdivision is required. 707 // to indicate that further subdivision is required.
704 RefinerResult::NeedRefinement 708 RefinerResult::NeedRefinement
705 } else { 709 } else {
706 // The data is refined enough, so return new hopefully better bounds 710 // The data is refined enough, so return new hopefully better bounds
715 }; 719 };
716 RefinerResult::Uncertain(bounds, Some(res)) 720 RefinerResult::Uncertain(bounds, Some(res))
717 } 721 }
718 } 722 }
719 723
720 fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { 724 fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
721 match (*r1, r2) { 725 match (*r1, r2) {
722 (Some((_, v1)), Some((_, v2))) => if v1 > v2 { *r1 = r2 } 726 (Some((_, v1)), Some((_, v2))) => {
727 if v1 > v2 {
728 *r1 = r2
729 }
730 }
723 (_, Some(_)) => *r1 = r2, 731 (_, Some(_)) => *r1 = r2,
724 (_, _) => {}, 732 (_, _) => {}
725 } 733 }
726 } 734 }
727 } 735 }
728
729 736
730 /// A bisection tree [`Refiner`] for checking that a [`BTFN`] is within a stated 737 /// A bisection tree [`Refiner`] for checking that a [`BTFN`] is within a stated
731 //// upper or lower bound. 738 //// upper or lower bound.
732 /// 739 ///
733 /// The type parameter `T` should be either [`RefineMax`] for upper bound or [`RefineMin`] 740 /// The type parameter `T` should be either [`RefineMax`] for upper bound or [`RefineMin`]
734 /// for lower bound. 741 /// for lower bound.
735 742
736 struct BoundRefiner<F : Float, T> { 743 struct BoundRefiner<F: Float, T> {
737 /// The upper/lower bound to check for 744 /// The upper/lower bound to check for
738 bound : F, 745 bound: F,
739 /// Tolerance for function value estimation. 746 /// Tolerance for function value estimation.
740 tolerance : F, 747 tolerance: F,
741 /// Maximum number of steps to execute the refiner for 748 /// Maximum number of steps to execute the refiner for
742 max_steps : usize, 749 max_steps: usize,
743 #[allow(dead_code)] // `how` is just for type system purposes. 750 #[allow(dead_code)] // `how` is just for type system purposes.
744 /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes. 751 /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes.
745 how : T, 752 how: T,
746 } 753 }
747 754
748 impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> 755 impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for BoundRefiner<F, RefineMax>
749 for BoundRefiner<F, RefineMax> 756 where
750 where G : SupportGenerator<F, N> { 757 G: SupportGenerator<F, N>,
758 {
751 type Result = bool; 759 type Result = bool;
752 type Sorting = UpperBoundSorting<F>; 760 type Sorting = UpperBoundSorting<F>;
753 761
754 fn refine( 762 fn refine(
755 &self, 763 &self,
756 aggregator : &Bounds<F>, 764 aggregator: &Bounds<F>,
757 _cube : &Cube<F, N>, 765 _cube: &Cube<F, N>,
758 _data : &[G::Id], 766 _data: &[G::Id],
759 _generator : &G, 767 _generator: &G,
760 step : usize 768 step: usize,
761 ) -> RefinerResult<Bounds<F>, Self::Result> { 769 ) -> RefinerResult<Bounds<F>, Self::Result> {
762 if aggregator.upper() <= self.bound + self.tolerance { 770 if aggregator.upper() <= self.bound + self.tolerance {
763 // Below upper bound within tolerances. Indicate uncertain success. 771 // Below upper bound within tolerances. Indicate uncertain success.
764 RefinerResult::Uncertain(*aggregator, true) 772 RefinerResult::Uncertain(*aggregator, true)
765 } else if aggregator.lower() >= self.bound - self.tolerance { 773 } else if aggregator.lower() >= self.bound - self.tolerance {
772 // No decision possible, but past step bounds 780 // No decision possible, but past step bounds
773 RefinerResult::Uncertain(*aggregator, false) 781 RefinerResult::Uncertain(*aggregator, false)
774 } 782 }
775 } 783 }
776 784
777 fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { 785 fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
778 *r1 = *r1 && r2; 786 *r1 = *r1 && r2;
779 } 787 }
780 } 788 }
781 789
782 impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N> 790 impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for BoundRefiner<F, RefineMin>
783 for BoundRefiner<F, RefineMin> 791 where
784 where G : SupportGenerator<F, N> { 792 G: SupportGenerator<F, N>,
793 {
785 type Result = bool; 794 type Result = bool;
786 type Sorting = UpperBoundSorting<F>; 795 type Sorting = UpperBoundSorting<F>;
787 796
788 fn refine( 797 fn refine(
789 &self, 798 &self,
790 aggregator : &Bounds<F>, 799 aggregator: &Bounds<F>,
791 _cube : &Cube<F, N>, 800 _cube: &Cube<F, N>,
792 _data : &[G::Id], 801 _data: &[G::Id],
793 _generator : &G, 802 _generator: &G,
794 step : usize 803 step: usize,
795 ) -> RefinerResult<Bounds<F>, Self::Result> { 804 ) -> RefinerResult<Bounds<F>, Self::Result> {
796 if aggregator.lower() >= self.bound - self.tolerance { 805 if aggregator.lower() >= self.bound - self.tolerance {
797 // Above lower bound within tolerances. Indicate uncertain success. 806 // Above lower bound within tolerances. Indicate uncertain success.
798 RefinerResult::Uncertain(*aggregator, true) 807 RefinerResult::Uncertain(*aggregator, true)
799 } else if aggregator.upper() <= self.bound + self.tolerance { 808 } else if aggregator.upper() <= self.bound + self.tolerance {
806 // No decision possible, but past step bounds 815 // No decision possible, but past step bounds
807 RefinerResult::Uncertain(*aggregator, false) 816 RefinerResult::Uncertain(*aggregator, false)
808 } 817 }
809 } 818 }
810 819
811 fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) { 820 fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
812 *r1 = *r1 && r2; 821 *r1 = *r1 && r2;
813 } 822 }
814 } 823 }
815 824
816 // FIXME: The most likely reason for the “Refiner failure” expectation in the methods below 825 // FIXME: The most likely reason for the “Refiner failure” expectation in the methods below
826 // the queue. If the queue empties as a result of that, the thread goes to wait for other threads 835 // the queue. If the queue empties as a result of that, the thread goes to wait for other threads
827 // to produce results. Since some node had a node whose lower bound was above the `glb`, eventually 836 // to produce results. Since some node had a node whose lower bound was above the `glb`, eventually
828 // there should be a result, or new nodes above the `glb` inserted into the queue. Then the waiting 837 // there should be a result, or new nodes above the `glb` inserted into the queue. Then the waiting
829 // threads can also continue processing. If, however, numerical inaccuracy destroyes the `glb`, 838 // threads can also continue processing. If, however, numerical inaccuracy destroyes the `glb`,
830 // the queue may run out, and we get “Refiner failure”. 839 // the queue may run out, and we get “Refiner failure”.
831 impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N> 840 impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N>
832 where BT : BTSearch<F, N, Agg=Bounds<F>>, 841 where
833 G : SupportGenerator<F, N, Id=BT::Data>, 842 BT: BTSearch<F, N, Agg = Bounds<F>>,
834 G::SupportType : Mapping<Loc<F, N>, Codomain=F> 843 G: SupportGenerator<F, N, Id = BT::Data>,
835 + LocalAnalysis<F, Bounds<F>, N>, 844 G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
836 Cube<F, N> : P2Minimise<Loc<F, N>, F> { 845 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
837 846 {
838 /// Maximise the `BTFN` within stated value `tolerance`. 847 /// Maximise the `BTFN` within stated value `tolerance`.
839 /// 848 ///
840 /// At most `max_steps` refinement steps are taken. 849 /// At most `max_steps` refinement steps are taken.
841 /// Returns the approximate maximiser and the corresponding function value. 850 /// Returns the approximate maximiser and the corresponding function value.
842 pub fn maximise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) { 851 pub fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F) {
843 let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : None }; 852 let refiner = P2Refiner {
844 self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.").unwrap() 853 tolerance,
854 max_steps,
855 how: RefineMax,
856 bound: None,
857 };
858 self.bt
859 .search_and_refine(refiner, &self.generator)
860 .expect("Refiner failure.")
861 .unwrap()
845 } 862 }
846 863
847 /// Maximise the `BTFN` within stated value `tolerance` subject to a lower bound. 864 /// Maximise the `BTFN` within stated value `tolerance` subject to a lower bound.
848 /// 865 ///
849 /// At most `max_steps` refinement steps are taken. 866 /// At most `max_steps` refinement steps are taken.
850 /// Returns the approximate maximiser and the corresponding function value when one is found 867 /// Returns the approximate maximiser and the corresponding function value when one is found
851 /// above the `bound` threshold, otherwise `None`. 868 /// above the `bound` threshold, otherwise `None`.
852 pub fn maximise_above(&mut self, bound : F, tolerance : F, max_steps : usize) 869 pub fn maximise_above(
853 -> Option<(Loc<F, N>, F)> { 870 &mut self,
854 let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : Some(bound) }; 871 bound: F,
855 self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") 872 tolerance: F,
873 max_steps: usize,
874 ) -> Option<(Loc<F, N>, F)> {
875 let refiner = P2Refiner {
876 tolerance,
877 max_steps,
878 how: RefineMax,
879 bound: Some(bound),
880 };
881 self.bt
882 .search_and_refine(refiner, &self.generator)
883 .expect("Refiner failure.")
856 } 884 }
857 885
858 /// Minimise the `BTFN` within stated value `tolerance`. 886 /// Minimise the `BTFN` within stated value `tolerance`.
859 /// 887 ///
860 /// At most `max_steps` refinement steps are taken. 888 /// At most `max_steps` refinement steps are taken.
861 /// Returns the approximate minimiser and the corresponding function value. 889 /// Returns the approximate minimiser and the corresponding function value.
862 pub fn minimise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) { 890 pub fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F) {
863 let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : None }; 891 let refiner = P2Refiner {
864 self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.").unwrap() 892 tolerance,
893 max_steps,
894 how: RefineMin,
895 bound: None,
896 };
897 self.bt
898 .search_and_refine(refiner, &self.generator)
899 .expect("Refiner failure.")
900 .unwrap()
865 } 901 }
866 902
867 /// Minimise the `BTFN` within stated value `tolerance` subject to a lower bound. 903 /// Minimise the `BTFN` within stated value `tolerance` subject to a lower bound.
868 /// 904 ///
869 /// At most `max_steps` refinement steps are taken. 905 /// At most `max_steps` refinement steps are taken.
870 /// Returns the approximate minimiser and the corresponding function value when one is found 906 /// Returns the approximate minimiser and the corresponding function value when one is found
871 /// above the `bound` threshold, otherwise `None`. 907 /// above the `bound` threshold, otherwise `None`.
872 pub fn minimise_below(&mut self, bound : F, tolerance : F, max_steps : usize) 908 pub fn minimise_below(
873 -> Option<(Loc<F, N>, F)> { 909 &mut self,
874 let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : Some(bound) }; 910 bound: F,
875 self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") 911 tolerance: F,
912 max_steps: usize,
913 ) -> Option<(Loc<F, N>, F)> {
914 let refiner = P2Refiner {
915 tolerance,
916 max_steps,
917 how: RefineMin,
918 bound: Some(bound),
919 };
920 self.bt
921 .search_and_refine(refiner, &self.generator)
922 .expect("Refiner failure.")
876 } 923 }
877 924
878 /// Verify that the `BTFN` has a given upper `bound` within indicated `tolerance`. 925 /// Verify that the `BTFN` has a given upper `bound` within indicated `tolerance`.
879 /// 926 ///
880 /// At most `max_steps` refinement steps are taken. 927 /// At most `max_steps` refinement steps are taken.
881 pub fn has_upper_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool { 928 pub fn has_upper_bound(&mut self, bound: F, tolerance: F, max_steps: usize) -> bool {
882 let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMax }; 929 let refiner = BoundRefiner {
883 self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") 930 bound,
931 tolerance,
932 max_steps,
933 how: RefineMax,
934 };
935 self.bt
936 .search_and_refine(refiner, &self.generator)
937 .expect("Refiner failure.")
884 } 938 }
885 939
886 /// Verify that the `BTFN` has a given lower `bound` within indicated `tolerance`. 940 /// Verify that the `BTFN` has a given lower `bound` within indicated `tolerance`.
887 /// 941 ///
888 /// At most `max_steps` refinement steps are taken. 942 /// At most `max_steps` refinement steps are taken.
889 pub fn has_lower_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool { 943 pub fn has_lower_bound(&mut self, bound: F, tolerance: F, max_steps: usize) -> bool {
890 let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMin }; 944 let refiner = BoundRefiner {
891 self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.") 945 bound,
892 } 946 tolerance,
893 } 947 max_steps,
948 how: RefineMin,
949 };
950 self.bt
951 .search_and_refine(refiner, &self.generator)
952 .expect("Refiner failure.")
953 }
954 }

Mercurial Repository: alg_tools

mercurial