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