Mercurial > repos > alg_tools / file comparison

comparison: src/bisection_tree/btfn.rs

src/bisection_tree/btfn.rs

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

Mercurial Repository: alg_tools

mercurial