Fri, 07 Oct 2022 13:11:08 +0300
Taylor 2 model attempts
0 | 1 | //! Implementation of the gaussian kernel. |
2 | ||
3 | use float_extras::f64::erf; | |
4 | use numeric_literals::replace_float_literals; | |
5 | use serde::Serialize; | |
6 | use alg_tools::types::*; | |
7 | use alg_tools::euclidean::Euclidean; | |
8 | use alg_tools::norms::*; | |
9 | use alg_tools::loc::Loc; | |
10 | use alg_tools::sets::Cube; | |
11 | use alg_tools::bisection_tree::{ | |
12 | Support, | |
13 | Constant, | |
14 | Bounds, | |
15 | LocalAnalysis, | |
16 | GlobalAnalysis, | |
17 | Weighted, | |
18 | Bounded, | |
3 | 19 | Taylor2Model, |
20 | Taylor2ModelParams, | |
0 | 21 | }; |
22 | use alg_tools::mapping::Apply; | |
23 | use alg_tools::maputil::array_init; | |
24 | ||
25 | use crate::fourier::Fourier; | |
26 | use super::base::*; | |
27 | use super::ball_indicator::CubeIndicator; | |
28 | ||
29 | /// Storage presentation of the the anisotropic gaussian kernel of `variance` $σ^2$. | |
30 | /// | |
31 | /// This is the function $f(x) = C e^{-\\|x\\|\_2^2/(2σ^2)}$ for $x ∈ ℝ^N$ | |
32 | /// with $C=1/(2πσ^2)^{N/2}$. | |
33 | #[derive(Copy,Clone,Debug,Serialize,Eq)] | |
34 | pub struct Gaussian<S : Constant, const N : usize> { | |
35 | /// The variance $σ^2$. | |
36 | pub variance : S, | |
37 | } | |
38 | ||
39 | impl<S1, S2, const N : usize> PartialEq<Gaussian<S2, N>> for Gaussian<S1, N> | |
40 | where S1 : Constant, | |
41 | S2 : Constant<Type=S1::Type> { | |
42 | fn eq(&self, other : &Gaussian<S2, N>) -> bool { | |
43 | self.variance.value() == other.variance.value() | |
44 | } | |
45 | } | |
46 | ||
47 | impl<S1, S2, const N : usize> PartialOrd<Gaussian<S2, N>> for Gaussian<S1, N> | |
48 | where S1 : Constant, | |
49 | S2 : Constant<Type=S1::Type> { | |
50 | ||
51 | fn partial_cmp(&self, other : &Gaussian<S2, N>) -> Option<std::cmp::Ordering> { | |
52 | // A gaussian is ≤ another gaussian if the Fourier transforms satisfy the | |
53 | // corresponding inequality. That in turns holds if and only if the variances | |
54 | // satisfy the opposite inequality. | |
55 | let σ1sq = self.variance.value(); | |
56 | let σ2sq = other.variance.value(); | |
57 | σ2sq.partial_cmp(&σ1sq) | |
58 | } | |
59 | } | |
60 | ||
61 | ||
62 | #[replace_float_literals(S::Type::cast_from(literal))] | |
63 | impl<'a, S, const N : usize> Apply<&'a Loc<S::Type, N>> for Gaussian<S, N> | |
64 | where S : Constant { | |
65 | type Output = S::Type; | |
66 | // This is not normalised to neither to have value 1 at zero or integral 1 | |
67 | // (unless the cut-off ε=0). | |
68 | #[inline] | |
69 | fn apply(&self, x : &'a Loc<S::Type, N>) -> Self::Output { | |
70 | let d_squared = x.norm2_squared(); | |
71 | let σ2 = self.variance.value(); | |
72 | let scale = self.scale(); | |
73 | (-d_squared / (2.0 * σ2)).exp() / scale | |
74 | } | |
75 | } | |
76 | ||
77 | impl<S, const N : usize> Apply<Loc<S::Type, N>> for Gaussian<S, N> | |
78 | where S : Constant { | |
79 | type Output = S::Type; | |
80 | // This is not normalised to neither to have value 1 at zero or integral 1 | |
81 | // (unless the cut-off ε=0). | |
82 | #[inline] | |
83 | fn apply(&self, x : Loc<S::Type, N>) -> Self::Output { | |
84 | self.apply(&x) | |
85 | } | |
86 | } | |
87 | ||
88 | ||
89 | #[replace_float_literals(S::Type::cast_from(literal))] | |
90 | impl<'a, S, const N : usize> Gaussian<S, N> | |
91 | where S : Constant { | |
92 | ||
93 | /// Returns the (reciprocal) scaling constant $1/C=(2πσ^2)^{N/2}$. | |
94 | #[inline] | |
95 | pub fn scale(&self) -> S::Type { | |
96 | let π = S::Type::PI; | |
97 | let σ2 = self.variance.value(); | |
98 | (2.0*π*σ2).powi(N as i32).sqrt() | |
99 | } | |
100 | } | |
101 | ||
102 | impl<'a, S, const N : usize> Support<S::Type, N> for Gaussian<S, N> | |
103 | where S : Constant { | |
104 | #[inline] | |
105 | fn support_hint(&self) -> Cube<S::Type,N> { | |
106 | array_init(|| [S::Type::NEG_INFINITY, S::Type::INFINITY]).into() | |
107 | } | |
108 | ||
109 | #[inline] | |
110 | fn in_support(&self, _x : &Loc<S::Type,N>) -> bool { | |
111 | true | |
112 | } | |
113 | } | |
114 | ||
115 | #[replace_float_literals(S::Type::cast_from(literal))] | |
116 | impl<S, const N : usize> GlobalAnalysis<S::Type, Bounds<S::Type>> for Gaussian<S, N> | |
117 | where S : Constant { | |
118 | #[inline] | |
119 | fn global_analysis(&self) -> Bounds<S::Type> { | |
120 | Bounds(0.0, 1.0/self.scale()) | |
121 | } | |
122 | } | |
123 | ||
124 | impl<S, const N : usize> LocalAnalysis<S::Type, Bounds<S::Type>, N> for Gaussian<S, N> | |
125 | where S : Constant { | |
126 | #[inline] | |
127 | fn local_analysis(&self, cube : &Cube<S::Type, N>) -> Bounds<S::Type> { | |
128 | // The function is maximised/minimised where the 2-norm is minimised/maximised. | |
129 | let lower = self.apply(cube.maxnorm_point()); | |
130 | let upper = self.apply(cube.minnorm_point()); | |
131 | Bounds(lower, upper) | |
132 | } | |
133 | } | |
134 | ||
135 | #[replace_float_literals(C::Type::cast_from(literal))] | |
136 | impl<'a, C : Constant, const N : usize> Norm<C::Type, L1> | |
137 | for Gaussian<C, N> { | |
138 | #[inline] | |
139 | fn norm(&self, _ : L1) -> C::Type { | |
140 | 1.0 | |
141 | } | |
142 | } | |
143 | ||
144 | #[replace_float_literals(C::Type::cast_from(literal))] | |
145 | impl<'a, C : Constant, const N : usize> Norm<C::Type, Linfinity> | |
146 | for Gaussian<C, N> { | |
147 | #[inline] | |
148 | fn norm(&self, _ : Linfinity) -> C::Type { | |
149 | self.bounds().upper() | |
150 | } | |
151 | } | |
152 | ||
153 | #[replace_float_literals(C::Type::cast_from(literal))] | |
154 | impl<'a, C : Constant, const N : usize> Fourier<C::Type> | |
155 | for Gaussian<C, N> { | |
156 | type Domain = Loc<C::Type, N>; | |
157 | type Transformed = Weighted<Gaussian<C::Type, N>, C::Type>; | |
158 | ||
159 | #[inline] | |
160 | fn fourier(&self) -> Self::Transformed { | |
161 | let π = C::Type::PI; | |
162 | let σ2 = self.variance.value(); | |
163 | let g = Gaussian { variance : 1.0 / (4.0*π*π*σ2) }; | |
164 | g.weigh(g.scale()) | |
165 | } | |
166 | } | |
167 | ||
168 | /// Representation of the “cut” gaussian $f χ\_{[-a, a]^n}$ | |
169 | /// where $a>0$ and $f$ is a gaussian kernel on $ℝ^n$. | |
170 | pub type BasicCutGaussian<C, S, const N : usize> = SupportProductFirst<CubeIndicator<C, N>, | |
171 | Gaussian<S, N>>; | |
172 | ||
173 | ||
174 | /// This implements $χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$ | |
175 | /// where $a,b>0$ and $f$ is a gaussian kernel on $ℝ^n$. | |
176 | #[replace_float_literals(F::cast_from(literal))] | |
177 | impl<'a, F : Float, R, C, S, const N : usize> Apply<&'a Loc<F, N>> | |
178 | for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
179 | where R : Constant<Type=F>, | |
180 | C : Constant<Type=F>, | |
181 | S : Constant<Type=F> { | |
182 | ||
183 | type Output = F; | |
184 | ||
185 | #[inline] | |
186 | fn apply(&self, y : &'a Loc<F, N>) -> F { | |
187 | let Convolution(ref ind, | |
188 | SupportProductFirst(ref cut, | |
189 | ref gaussian)) = self; | |
190 | let a = cut.r.value(); | |
191 | let b = ind.r.value(); | |
192 | let σ = gaussian.variance.value().sqrt(); | |
193 | let π = F::PI; | |
194 | let t = F::SQRT_2 * σ; | |
195 | let c = σ * (8.0/π).sqrt(); | |
196 | ||
197 | // This is just a product of one-dimensional versions | |
198 | let unscaled = y.product_map(|x| { | |
199 | let c1 = -(a.min(b + x)); //(-a).max(-x-b); | |
200 | let c2 = a.min(b - x); | |
201 | if c1 >= c2 { | |
202 | 0.0 | |
203 | } else { | |
204 | let e1 = F::cast_from(erf((c1 / t).as_())); | |
205 | let e2 = F::cast_from(erf((c2 / t).as_())); | |
206 | debug_assert!(e2 >= e1); | |
207 | c * (e2 - e1) | |
208 | } | |
209 | }); | |
210 | ||
211 | unscaled / gaussian.scale() | |
212 | } | |
213 | } | |
214 | ||
215 | impl<F : Float, R, C, S, const N : usize> Apply<Loc<F, N>> | |
216 | for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
217 | where R : Constant<Type=F>, | |
218 | C : Constant<Type=F>, | |
219 | S : Constant<Type=F> { | |
220 | ||
221 | type Output = F; | |
222 | ||
223 | #[inline] | |
224 | fn apply(&self, y : Loc<F, N>) -> F { | |
225 | self.apply(&y) | |
226 | } | |
227 | } | |
228 | ||
3 | 229 | #[replace_float_literals(F::cast_from(literal))] |
230 | impl<'a, F : Float, R, C, S, const N : usize> Taylor2Model<F, N> | |
231 | for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
232 | where R : Constant<Type=F>, | |
233 | C : Constant<Type=F>, | |
234 | S : Constant<Type=F>, | |
235 | Loc<Loc<F, 3>, N> : Product2Taylor<F, N> { | |
236 | ||
237 | fn taylor2model(&self, y : Loc<F, N>) -> Taylor2ModelParams<F, N> { | |
238 | let Convolution(ref ind, | |
239 | SupportProductFirst(ref cut, | |
240 | ref gaussian)) = self; | |
241 | let a = cut.r.value(); | |
242 | let b = ind.r.value(); | |
243 | let σ = gaussian.variance.value().sqrt(); | |
244 | let π = F::PI; | |
245 | let t = F::SQRT_2 * σ; | |
246 | let c = σ * (8.0/π).sqrt(); | |
247 | let sc = gaussian.scale(); | |
248 | ||
249 | // This is just a product of one-dimensional versions | |
250 | y.map(|x| { | |
251 | let c1 = -(a.min(b + x)); //(-a).max(-x-b); | |
252 | let c2 = a.min(b - x); | |
253 | if c1 >= c2 { | |
254 | Loc([0.0, 0.0, 0.0]) | |
255 | } else { | |
256 | let s1 = c1 / t; | |
257 | let s2 = c2 / t; | |
258 | let e1 = F::cast_from(erf(s1.as_())); | |
259 | let e2 = F::cast_from(erf(s2.as_())); | |
260 | let d1 = F::FRAC_2_SQRT_PI * (-s1*s1).exp(); | |
261 | let d2 = F::FRAC_2_SQRT_PI * (-s2*s2).exp(); | |
262 | let dd1 = -2.0 * s1 * d1; | |
263 | let dd2 = -2.0 * s2 * s2; | |
264 | debug_assert!(e2 >= e1); | |
265 | Loc([c * (e2 - e1), c * (d2 - d1), c * (dd2 - dd1)]) | |
266 | } | |
267 | }).product2taylor_scaled(sc) | |
268 | } | |
269 | ||
270 | } | |
271 | ||
0 | 272 | impl<F : Float, R, C, S, const N : usize> |
273 | Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
274 | where R : Constant<Type=F>, | |
275 | C : Constant<Type=F>, | |
276 | S : Constant<Type=F> { | |
277 | ||
278 | #[inline] | |
279 | fn get_r(&self) -> F { | |
280 | let Convolution(ref ind, | |
281 | SupportProductFirst(ref cut, ..)) = self; | |
282 | ind.r.value() + cut.r.value() | |
283 | } | |
284 | } | |
285 | ||
286 | impl<F : Float, R, C, S, const N : usize> Support<F, N> | |
287 | for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
288 | where R : Constant<Type=F>, | |
289 | C : Constant<Type=F>, | |
290 | S : Constant<Type=F> { | |
291 | #[inline] | |
292 | fn support_hint(&self) -> Cube<F, N> { | |
293 | let r = self.get_r(); | |
294 | array_init(|| [-r, r]).into() | |
295 | } | |
296 | ||
297 | #[inline] | |
298 | fn in_support(&self, y : &Loc<F, N>) -> bool { | |
299 | let r = self.get_r(); | |
300 | y.iter().all(|x| x.abs() <= r) | |
301 | } | |
302 | ||
303 | #[inline] | |
304 | fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] { | |
305 | let r = self.get_r(); | |
306 | // From c1 = -a.min(b + x) and c2 = a.min(b - x) with c_1 < c_2, | |
307 | // solve bounds for x. that is 0 ≤ a.min(b + x) + a.min(b - x). | |
308 | // If b + x ≤ a and b - x ≤ a, the sum is 2b ≥ 0. | |
309 | // If b + x ≥ a and b - x ≥ a, the sum is 2a ≥ 0. | |
310 | // If b + x ≤ a and b - x ≥ a, the sum is b + x + a ⟹ need x ≥ -a - b = -r. | |
311 | // If b + x ≥ a and b - x ≤ a, the sum is a + b - x ⟹ need x ≤ a + b = r. | |
312 | cube.map(|c, d| symmetric_peak_hint(r, c, d)) | |
313 | } | |
314 | } | |
315 | ||
316 | impl<F : Float, R, C, S, const N : usize> GlobalAnalysis<F, Bounds<F>> | |
317 | for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
318 | where R : Constant<Type=F>, | |
319 | C : Constant<Type=F>, | |
320 | S : Constant<Type=F> { | |
321 | #[inline] | |
322 | fn global_analysis(&self) -> Bounds<F> { | |
323 | Bounds(F::ZERO, self.apply(Loc::ORIGIN)) | |
324 | } | |
325 | } | |
326 | ||
327 | impl<F : Float, R, C, S, const N : usize> LocalAnalysis<F, Bounds<F>, N> | |
328 | for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
329 | where R : Constant<Type=F>, | |
330 | C : Constant<Type=F>, | |
331 | S : Constant<Type=F> { | |
332 | #[inline] | |
333 | fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> { | |
334 | // The function is maximised/minimised where the absolute value is minimised/maximised. | |
335 | let lower = self.apply(cube.maxnorm_point()); | |
336 | let upper = self.apply(cube.minnorm_point()); | |
337 | Bounds(lower, upper) | |
338 | } | |
339 | } | |
340 |