Fri, 18 Nov 2022 10:29:50 +0200
Improvements and minor fixes to bisection tree refinement.
0 | 1 | |
5 | 2 | /*! |
3 | Traits for representing the support of a [`Mapping`], and analysing the mapping on a [`Cube`]. | |
4 | */ | |
0 | 5 | use serde::Serialize; |
6 | use std::ops::{MulAssign,DivAssign,Neg}; | |
5 | 7 | use crate::types::{Float, Num}; |
8 | use crate::maputil::map2; | |
9 | use crate::mapping::Mapping; | |
10 | use crate::sets::Cube; | |
11 | use crate::loc::Loc; | |
0 | 12 | use super::aggregator::Bounds; |
13 | use crate::norms::{Norm, L1, L2, Linfinity}; | |
14 | ||
5 | 15 | /// A trait for encoding constant [`Float`] values |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
16 | pub trait Constant : Copy + Sync + Send + 'static + std::fmt::Debug + Into<Self::Type> { |
5 | 17 | /// The type of the value |
0 | 18 | type Type : Float; |
5 | 19 | /// Returns the value of the constant |
0 | 20 | fn value(&self) -> Self::Type; |
21 | } | |
22 | ||
23 | impl<F : Float> Constant for F { | |
24 | type Type = F; | |
25 | #[inline] | |
26 | fn value(&self) -> F { *self } | |
27 | } | |
28 | ||
29 | ||
30 | /// A trait for working with the supports of [`Mapping`]s. | |
5 | 31 | /// |
0 | 32 | /// Mapping is not a super-trait to allow more general use. |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
33 | pub trait Support<F : Num, const N : usize> : Sized + Sync + Send + 'static { |
5 | 34 | /// Return a cube containing the support of the function represented by `self`. |
35 | /// | |
36 | /// The hint may be larger than the actual support, but must contain it. | |
0 | 37 | fn support_hint(&self) -> Cube<F,N>; |
38 | ||
5 | 39 | /// Indicate whether `x` is in the support of the function represented by `self`. |
0 | 40 | fn in_support(&self, x : &Loc<F,N>) -> bool; |
41 | ||
5 | 42 | // Indicate whether `cube` is fully in the support of the function represented by `self`. |
0 | 43 | //fn fully_in_support(&self, cube : &Cube<F,N>) -> bool; |
44 | ||
5 | 45 | /// Return an optional hint for bisecting the support. |
46 | /// | |
47 | /// The output along each axis a possible coordinate at which to bisect `cube`. | |
48 | /// | |
49 | /// This is useful for nonsmooth functions to make finite element models as used by | |
50 | /// [`BTFN`][super::btfn::BTFN] minimisation/maximisation compatible with points of | |
51 | /// non-differentiability. | |
52 | /// | |
53 | /// The default implementation returns `[None; N]`. | |
0 | 54 | #[inline] |
5 | 55 | #[allow(unused_variables)] |
56 | fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] { | |
57 | [None; N] | |
0 | 58 | } |
59 | ||
5 | 60 | /// Translate `self` by `x`. |
0 | 61 | #[inline] |
62 | fn shift(self, x : Loc<F, N>) -> Shift<Self, F, N> { | |
63 | Shift { shift : x, base_fn : self } | |
64 | } | |
65 | ||
5 | 66 | /// Multiply `self` by the scalar `a`. |
0 | 67 | #[inline] |
68 | fn weigh<C : Constant<Type=F>>(self, a : C) -> Weighted<Self, C> { | |
69 | Weighted { weight : a, base_fn : self } | |
70 | } | |
71 | } | |
72 | ||
5 | 73 | /// Trait for globally analysing a property `A` of a [`Mapping`]. |
74 | /// | |
75 | /// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as | |
76 | /// [`Bounds`][super::aggregator::Bounds]. | |
0 | 77 | pub trait GlobalAnalysis<F : Num, A> { |
5 | 78 | /// Perform global analysis of the property `A` of `Self`. |
79 | /// | |
80 | /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds], | |
81 | /// this function will return global upper and lower bounds for the mapping | |
82 | /// represented by `self`. | |
0 | 83 | fn global_analysis(&self) -> A; |
84 | } | |
85 | ||
86 | // default impl<F, A, N, L> GlobalAnalysis<F, A, N> for L | |
87 | // where L : LocalAnalysis<F, A, N> { | |
88 | // #[inline] | |
89 | // fn global_analysis(&self) -> Bounds<F> { | |
90 | // self.local_analysis(&self.support_hint()) | |
91 | // } | |
92 | // } | |
93 | ||
5 | 94 | /// Trait for locally analysing a property `A` of a [`Mapping`] (implementing [`Support`]) |
95 | /// within a [`Cube`]. | |
96 | /// | |
97 | /// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as | |
98 | /// [`Bounds`][super::aggregator::Bounds]. | |
0 | 99 | pub trait LocalAnalysis<F : Num, A, const N : usize> : GlobalAnalysis<F, A> + Support<F, N> { |
5 | 100 | /// Perform local analysis of the property `A` of `Self`. |
101 | /// | |
102 | /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds], | |
103 | /// this function will return upper and lower bounds within `cube` for the mapping | |
104 | /// represented by `self`. | |
0 | 105 | fn local_analysis(&self, cube : &Cube<F, N>) -> A; |
106 | } | |
107 | ||
5 | 108 | /// Trait for determining the upper and lower bounds of an float-valued [`Mapping`]. |
109 | /// | |
0 | 110 | /// This is a blanket-implemented alias for [`GlobalAnalysis`]`<F, Bounds<F>>` |
111 | /// [`Mapping`] is not a supertrait to allow flexibility in the implementation of either | |
112 | /// reference or non-reference arguments. | |
113 | pub trait Bounded<F : Float> : GlobalAnalysis<F, Bounds<F>> { | |
114 | /// Return lower and upper bounds for the values of of `self`. | |
115 | #[inline] | |
116 | fn bounds(&self) -> Bounds<F> { | |
117 | self.global_analysis() | |
118 | } | |
119 | } | |
120 | ||
121 | impl<F : Float, T : GlobalAnalysis<F, Bounds<F>>> Bounded<F> for T { } | |
122 | ||
5 | 123 | /// Shift of [`Support`] and [`Mapping`]; output of [`Support::shift`]. |
0 | 124 | #[derive(Copy,Clone,Debug,Serialize)] // Serialize! but not implemented by Loc. |
125 | pub struct Shift<T, F, const N : usize> { | |
126 | shift : Loc<F, N>, | |
127 | base_fn : T, | |
128 | } | |
129 | ||
130 | impl<'a, T, V, F : Float, const N : usize> Mapping<&'a Loc<F, N>> for Shift<T,F,N> | |
131 | where T : for<'b> Mapping<&'b Loc<F,N>,Codomain=V> { | |
132 | type Codomain = V; | |
133 | #[inline] | |
134 | fn value(&self, x : &'a Loc<F, N>) -> Self::Codomain { | |
135 | self.base_fn.value(&(x - &self.shift)) | |
136 | } | |
137 | } | |
138 | ||
139 | impl<'a, T, V, F : Float, const N : usize> Mapping<Loc<F, N>> for Shift<T,F,N> | |
140 | where T : for<'b> Mapping<Loc<F,N>,Codomain=V> { | |
141 | type Codomain = V; | |
142 | #[inline] | |
143 | fn value(&self, x : Loc<F, N>) -> Self::Codomain { | |
144 | self.base_fn.value(x - &self.shift) | |
145 | } | |
146 | } | |
147 | ||
148 | impl<'a, T, F : Float, const N : usize> Support<F,N> for Shift<T,F,N> | |
149 | where T : Support<F, N> { | |
150 | #[inline] | |
151 | fn support_hint(&self) -> Cube<F,N> { | |
152 | self.base_fn.support_hint().shift(&self.shift) | |
153 | } | |
154 | ||
155 | #[inline] | |
156 | fn in_support(&self, x : &Loc<F,N>) -> bool { | |
157 | self.base_fn.in_support(&(x - &self.shift)) | |
158 | } | |
159 | ||
160 | // fn fully_in_support(&self, _cube : &Cube<F,N>) -> bool { | |
161 | // //self.base_fn.fully_in_support(cube.shift(&vectorneg(self.shift))) | |
162 | // todo!("Not implemented, but not used at the moment") | |
163 | // } | |
164 | ||
165 | #[inline] | |
166 | fn bisection_hint(&self, cube : &Cube<F,N>) -> [Option<F>; N] { | |
167 | let base_hint = self.base_fn.bisection_hint(cube); | |
168 | map2(base_hint, &self.shift, |h, s| h.map(|z| z + *s)) | |
169 | } | |
170 | ||
171 | } | |
172 | ||
173 | impl<'a, T, F : Float, const N : usize> GlobalAnalysis<F, Bounds<F>> for Shift<T,F,N> | |
174 | where T : LocalAnalysis<F, Bounds<F>, N> { | |
175 | #[inline] | |
176 | fn global_analysis(&self) -> Bounds<F> { | |
177 | self.base_fn.global_analysis() | |
178 | } | |
179 | } | |
180 | ||
181 | impl<'a, T, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Shift<T,F,N> | |
182 | where T : LocalAnalysis<F, Bounds<F>, N> { | |
183 | #[inline] | |
184 | fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> { | |
185 | self.base_fn.local_analysis(&cube.shift(&(-self.shift))) | |
186 | } | |
187 | } | |
188 | ||
189 | macro_rules! impl_shift_norm { | |
190 | ($($norm:ident)*) => { $( | |
191 | impl<'a, T, F : Float, const N : usize> Norm<F, $norm> for Shift<T,F,N> | |
192 | where T : Norm<F, $norm> { | |
193 | #[inline] | |
194 | fn norm(&self, n : $norm) -> F { | |
195 | self.base_fn.norm(n) | |
196 | } | |
197 | } | |
198 | )* } | |
199 | } | |
200 | ||
201 | impl_shift_norm!(L1 L2 Linfinity); | |
202 | ||
5 | 203 | /// Weighting of a [`Support`] and [`Mapping`] by scalar multiplication; |
204 | /// output of [`Support::weigh`]. | |
0 | 205 | #[derive(Copy,Clone,Debug,Serialize)] |
206 | pub struct Weighted<T, C : Constant> { | |
5 | 207 | /// The weight |
0 | 208 | pub weight : C, |
5 | 209 | /// The base [`Support`] or [`Mapping`] being weighted. |
0 | 210 | pub base_fn : T, |
211 | } | |
212 | ||
213 | impl<'a, T, V, F : Float, C, const N : usize> Mapping<&'a Loc<F, N>> for Weighted<T, C> | |
214 | where T : for<'b> Mapping<&'b Loc<F,N>,Codomain=V>, | |
215 | V : std::ops::Mul<F,Output=V>, | |
216 | C : Constant<Type=F> { | |
217 | type Codomain = V; | |
218 | #[inline] | |
219 | fn value(&self, x : &'a Loc<F, N>) -> Self::Codomain { | |
220 | self.base_fn.value(x) * self.weight.value() | |
221 | } | |
222 | } | |
223 | ||
224 | impl<'a, T, V, F : Float, C, const N : usize> Mapping<Loc<F, N>> for Weighted<T, C> | |
225 | where T : for<'b> Mapping<Loc<F,N>,Codomain=V>, | |
226 | V : std::ops::Mul<F,Output=V>, | |
227 | C : Constant<Type=F> { | |
228 | type Codomain = V; | |
229 | #[inline] | |
230 | fn value(&self, x : Loc<F, N>) -> Self::Codomain { | |
231 | self.base_fn.value(x) * self.weight.value() | |
232 | } | |
233 | } | |
234 | ||
235 | impl<'a, T, F : Float, C, const N : usize> Support<F,N> for Weighted<T, C> | |
236 | where T : Support<F, N>, | |
237 | C : Constant<Type=F> { | |
238 | ||
239 | #[inline] | |
240 | fn support_hint(&self) -> Cube<F,N> { | |
241 | self.base_fn.support_hint() | |
242 | } | |
243 | ||
244 | #[inline] | |
245 | fn in_support(&self, x : &Loc<F,N>) -> bool { | |
246 | self.base_fn.in_support(x) | |
247 | } | |
248 | ||
249 | // fn fully_in_support(&self, cube : &Cube<F,N>) -> bool { | |
250 | // self.base_fn.fully_in_support(cube) | |
251 | // } | |
252 | ||
253 | #[inline] | |
254 | fn bisection_hint(&self, cube : &Cube<F,N>) -> [Option<F>; N] { | |
255 | self.base_fn.bisection_hint(cube) | |
256 | } | |
257 | } | |
258 | ||
259 | impl<'a, T, F : Float, C> GlobalAnalysis<F, Bounds<F>> for Weighted<T, C> | |
260 | where T : GlobalAnalysis<F, Bounds<F>>, | |
261 | C : Constant<Type=F> { | |
262 | #[inline] | |
263 | fn global_analysis(&self) -> Bounds<F> { | |
264 | let Bounds(lower, upper) = self.base_fn.global_analysis(); | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
265 | debug_assert!(lower <= upper); |
0 | 266 | match self.weight.value() { |
267 | w if w < F::ZERO => Bounds(w * upper, w * lower), | |
268 | w => Bounds(w * lower, w * upper), | |
269 | } | |
270 | } | |
271 | } | |
272 | ||
273 | impl<'a, T, F : Float, C, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Weighted<T, C> | |
274 | where T : LocalAnalysis<F, Bounds<F>, N>, | |
275 | C : Constant<Type=F> { | |
276 | #[inline] | |
277 | fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> { | |
278 | let Bounds(lower, upper) = self.base_fn.local_analysis(cube); | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
279 | debug_assert!(lower <= upper); |
0 | 280 | match self.weight.value() { |
281 | w if w < F::ZERO => Bounds(w * upper, w * lower), | |
282 | w => Bounds(w * lower, w * upper), | |
283 | } | |
284 | } | |
285 | } | |
286 | ||
287 | macro_rules! make_weighted_scalarop_rhs { | |
288 | ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { | |
289 | impl<F : Float, T> std::ops::$trait_assign<F> for Weighted<T, F> { | |
290 | #[inline] | |
291 | fn $fn_assign(&mut self, t : F) { | |
292 | self.weight.$fn_assign(t); | |
293 | } | |
294 | } | |
295 | ||
296 | impl<'a, F : Float, T> std::ops::$trait<F> for Weighted<T, F> { | |
297 | type Output = Self; | |
298 | #[inline] | |
299 | fn $fn(mut self, t : F) -> Self { | |
300 | self.weight.$fn_assign(t); | |
301 | self | |
302 | } | |
303 | } | |
304 | ||
305 | impl<'a, F : Float, T> std::ops::$trait<F> for &'a Weighted<T, F> | |
306 | where T : Clone { | |
307 | type Output = Weighted<T, F>; | |
308 | #[inline] | |
309 | fn $fn(self, t : F) -> Self::Output { | |
310 | Weighted { weight : self.weight.$fn(t), base_fn : self.base_fn.clone() } | |
311 | } | |
312 | } | |
313 | } | |
314 | } | |
315 | ||
316 | make_weighted_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); | |
317 | make_weighted_scalarop_rhs!(Div, div, DivAssign, div_assign); | |
318 | ||
319 | macro_rules! impl_weighted_norm { | |
320 | ($($norm:ident)*) => { $( | |
321 | impl<'a, T, F : Float> Norm<F, $norm> for Weighted<T,F> | |
322 | where T : Norm<F, $norm> { | |
323 | #[inline] | |
324 | fn norm(&self, n : $norm) -> F { | |
325 | self.base_fn.norm(n) * self.weight.abs() | |
326 | } | |
327 | } | |
328 | )* } | |
329 | } | |
330 | ||
331 | impl_weighted_norm!(L1 L2 Linfinity); | |
332 | ||
333 | ||
5 | 334 | /// Normalisation of [`Support`] and [`Mapping`] to L¹ norm 1. |
335 | /// | |
0 | 336 | /// Currently only scalar-valued functions are supported. |
337 | #[derive(Copy, Clone, Debug, Serialize, PartialEq)] | |
5 | 338 | pub struct Normalised<T>( |
339 | /// The base [`Support`] or [`Mapping`]. | |
340 | pub T | |
341 | ); | |
0 | 342 | |
343 | impl<'a, T, F : Float, const N : usize> Mapping<&'a Loc<F, N>> for Normalised<T> | |
344 | where T : Norm<F, L1> + for<'b> Mapping<&'b Loc<F,N>, Codomain=F> { | |
345 | type Codomain = F; | |
346 | #[inline] | |
347 | fn value(&self, x : &'a Loc<F, N>) -> Self::Codomain { | |
348 | let w = self.0.norm(L1); | |
349 | if w == F::ZERO { F::ZERO } else { self.0.value(x) / w } | |
350 | } | |
351 | } | |
352 | ||
353 | impl<'a, T, F : Float, const N : usize> Mapping<Loc<F, N>> for Normalised<T> | |
354 | where T : Norm<F, L1> + for<'b> Mapping<Loc<F,N>, Codomain=F> { | |
355 | type Codomain = F; | |
356 | #[inline] | |
357 | fn value(&self, x : Loc<F, N>) -> Self::Codomain { | |
358 | let w = self.0.norm(L1); | |
359 | if w == F::ZERO { F::ZERO } else { self.0.value(x) / w } | |
360 | } | |
361 | } | |
362 | ||
363 | impl<'a, T, F : Float, const N : usize> Support<F,N> for Normalised<T> | |
364 | where T : Norm<F, L1> + Support<F, N> { | |
365 | ||
366 | #[inline] | |
367 | fn support_hint(&self) -> Cube<F,N> { | |
368 | self.0.support_hint() | |
369 | } | |
370 | ||
371 | #[inline] | |
372 | fn in_support(&self, x : &Loc<F,N>) -> bool { | |
373 | self.0.in_support(x) | |
374 | } | |
375 | ||
376 | // fn fully_in_support(&self, cube : &Cube<F,N>) -> bool { | |
377 | // self.0.fully_in_support(cube) | |
378 | // } | |
379 | ||
380 | #[inline] | |
381 | fn bisection_hint(&self, cube : &Cube<F,N>) -> [Option<F>; N] { | |
382 | self.0.bisection_hint(cube) | |
383 | } | |
384 | } | |
385 | ||
386 | impl<'a, T, F : Float> GlobalAnalysis<F, Bounds<F>> for Normalised<T> | |
387 | where T : Norm<F, L1> + GlobalAnalysis<F, Bounds<F>> { | |
388 | #[inline] | |
389 | fn global_analysis(&self) -> Bounds<F> { | |
390 | let Bounds(lower, upper) = self.0.global_analysis(); | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
391 | debug_assert!(lower <= upper); |
0 | 392 | let w = self.0.norm(L1); |
393 | debug_assert!(w >= F::ZERO); | |
394 | Bounds(w * lower, w * upper) | |
395 | } | |
396 | } | |
397 | ||
398 | impl<'a, T, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Normalised<T> | |
399 | where T : Norm<F, L1> + LocalAnalysis<F, Bounds<F>, N> { | |
400 | #[inline] | |
401 | fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> { | |
402 | let Bounds(lower, upper) = self.0.local_analysis(cube); | |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
403 | debug_assert!(lower <= upper); |
0 | 404 | let w = self.0.norm(L1); |
405 | debug_assert!(w >= F::ZERO); | |
406 | Bounds(w * lower, w * upper) | |
407 | } | |
408 | } | |
409 | ||
410 | impl<'a, T, F : Float> Norm<F, L1> for Normalised<T> | |
411 | where T : Norm<F, L1> { | |
412 | #[inline] | |
413 | fn norm(&self, _ : L1) -> F { | |
414 | let w = self.0.norm(L1); | |
415 | if w == F::ZERO { F::ZERO } else { F::ONE } | |
416 | } | |
417 | } | |
418 | ||
419 | macro_rules! impl_normalised_norm { | |
420 | ($($norm:ident)*) => { $( | |
421 | impl<'a, T, F : Float> Norm<F, $norm> for Normalised<T> | |
422 | where T : Norm<F, $norm> + Norm<F, L1> { | |
423 | #[inline] | |
424 | fn norm(&self, n : $norm) -> F { | |
425 | let w = self.0.norm(L1); | |
426 | if w == F::ZERO { F::ZERO } else { self.0.norm(n) / w } | |
427 | } | |
428 | } | |
429 | )* } | |
430 | } | |
431 | ||
432 | impl_normalised_norm!(L2 Linfinity); | |
433 | ||
434 | /* | |
435 | impl<F : Num, S : Support<F, N>, const N : usize> LocalAnalysis<F, NullAggregator, N> for S { | |
436 | fn local_analysis(&self, _cube : &Cube<F, N>) -> NullAggregator { NullAggregator } | |
437 | } | |
438 | ||
439 | impl<F : Float, S : Bounded<F>, const N : usize> LocalAnalysis<F, Bounds<F>, N> for S { | |
440 | #[inline] | |
441 | fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> { | |
442 | self.bounds(cube) | |
443 | } | |
444 | }*/ | |
445 | ||
5 | 446 | /// Generator of [`Support`]-implementing component functions based on low storage requirement |
447 | /// [ids][`Self::Id`]. | |
0 | 448 | pub trait SupportGenerator<F : Float, const N : usize> |
8
4e09b7829b51
Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
449 | : MulAssign<F> + DivAssign<F> + Neg<Output=Self> + Clone + Sync + Send + 'static { |
5 | 450 | /// The identification type |
0 | 451 | type Id : 'static + Copy; |
5 | 452 | /// The type of the [`Support`] (often also a [`Mapping`]). |
0 | 453 | type SupportType : 'static + Support<F, N>; |
5 | 454 | /// An iterator over all the [`Support`]s of the generator. |
0 | 455 | type AllDataIter<'a> : Iterator<Item=(Self::Id, Self::SupportType)> where Self : 'a; |
456 | ||
5 | 457 | /// Returns the component identified by `id`. |
458 | /// | |
459 | /// Panics if `id` is an invalid identifier. | |
460 | fn support_for(&self, id : Self::Id) -> Self::SupportType; | |
0 | 461 | |
5 | 462 | /// Returns the number of different components in this generator. |
0 | 463 | fn support_count(&self) -> usize; |
464 | ||
5 | 465 | /// Returns an iterator over all pairs of `(id, support)`. |
0 | 466 | fn all_data(&self) -> Self::AllDataIter<'_>; |
467 | } | |
468 |