Mercurial > repos > pointsource_algs / file comparison

comparison: src/kernels/mollifier.rs

src/kernels/mollifier.rs

branch
dev
changeset 61
4f468d35fa29
parent 35
b087e3eab191
equal deleted inserted replaced
60:9738b51d90d7 61:4f468d35fa29
1
2 //! Implementation of the standard mollifier 1 //! Implementation of the standard mollifier
3 2
4 use rgsl::hypergeometric::hyperg_U; 3 use alg_tools::bisection_tree::{Bounds, Constant, GlobalAnalysis, LocalAnalysis, Support};
4 use alg_tools::euclidean::Euclidean;
5 use alg_tools::loc::Loc;
6 use alg_tools::mapping::{Instance, Mapping};
7 use alg_tools::maputil::array_init;
8 use alg_tools::norms::*;
9 use alg_tools::sets::Cube;
10 use alg_tools::types::*;
5 use float_extras::f64::tgamma as gamma; 11 use float_extras::f64::tgamma as gamma;
6 use numeric_literals::replace_float_literals; 12 use numeric_literals::replace_float_literals;
13 use rgsl::hypergeometric::hyperg_U;
7 use serde::Serialize; 14 use serde::Serialize;
8 use alg_tools::types::*;
9 use alg_tools::euclidean::Euclidean;
10 use alg_tools::norms::*;
11 use alg_tools::loc::Loc;
12 use alg_tools::sets::Cube;
13 use alg_tools::bisection_tree::{
14 Support,
15 Constant,
16 Bounds,
17 LocalAnalysis,
18 GlobalAnalysis
19 };
20 use alg_tools::mapping::{Mapping, Instance};
21 use alg_tools::maputil::array_init;
22 15
23 /// Reresentation of the (unnormalised) standard mollifier. 16 /// Reresentation of the (unnormalised) standard mollifier.
24 /// 17 ///
25 /// For the `width` parameter $ε>0$, this is 18 /// For the `width` parameter $ε>0$, this is
26 /// <div>$$ 19 /// <div>$$
27 /// f(x)=\begin{cases} 20 /// f(x)=\begin{cases}
28 /// e^{\frac{ε^2}{\|x\|_2^2-ε^2}}, & \|x\|_2 < ε, \\ 21 /// e^{\frac{ε^2}{\|x\|_2^2-ε^2}}, & \|x\|_2 < ε, \\
29 /// 0, & \text{otherwise}. 22 /// 0, & \text{otherwise}.
30 /// \end{cases} 23 /// \end{cases}
31 /// $$</div> 24 /// $$</div>
32 #[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)] 25 #[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
33 pub struct Mollifier<C : Constant, const N : usize> { 26 pub struct Mollifier<C: Constant, const N: usize> {
34 /// The parameter $ε$ of the mollifier. 27 /// The parameter $ε$ of the mollifier.
35 pub width : C, 28 pub width: C,
36 } 29 }
37 30
38 #[replace_float_literals(C::Type::cast_from(literal))] 31 #[replace_float_literals(C::Type::cast_from(literal))]
39 impl<C : Constant, const N : usize> Mapping<Loc<C::Type, N>> for Mollifier<C, N> { 32 impl<C: Constant, const N: usize> Mapping<Loc<N, C::Type>> for Mollifier<C, N> {
40 type Codomain = C::Type; 33 type Codomain = C::Type;
41 34
42 #[inline] 35 #[inline]
43 fn apply<I : Instance<Loc<C::Type, N>>>(&self, x : I) -> Self::Codomain { 36 fn apply<I: Instance<Loc<N, C::Type>>>(&self, x: I) -> Self::Codomain {
44 let ε = self.width.value(); 37 let ε = self.width.value();
45 let ε2 = ε*ε; 38 let ε2 = ε * ε;
46 let n2 = x.eval(|x| x.norm2_squared()); 39 let n2 = x.eval(|x| x.norm2_squared());
47 if n2 < ε2 { 40 if n2 < ε2 {
48 (n2 / (n2 - ε2)).exp() 41 (n2 / (n2 - ε2)).exp()
49 } else { 42 } else {
50 0.0 43 0.0
51 } 44 }
52 } 45 }
53 } 46 }
54 47
55 48 impl<'a, C: Constant, const N: usize> Support<N, C::Type> for Mollifier<C, N> {
56 impl<'a, C : Constant, const N : usize> Support<C::Type, N> for Mollifier<C, N> {
57 #[inline] 49 #[inline]
58 fn support_hint(&self) -> Cube<C::Type,N> { 50 fn support_hint(&self) -> Cube<N, C::Type> {
59 let ε = self.width.value(); 51 let ε = self.width.value();
60 array_init(|| [-ε, ε]).into() 52 array_init(|| [-ε, ε]).into()
61 } 53 }
62 54
63 #[inline] 55 #[inline]
64 fn in_support(&self, x : &Loc<C::Type,N>) -> bool { 56 fn in_support(&self, x: &Loc<N, C::Type>) -> bool {
65 x.norm2() < self.width.value() 57 x.norm2() < self.width.value()
66 } 58 }
67 59
68 /*fn fully_in_support(&self, _cube : &Cube<C::Type,N>) -> bool { 60 /*fn fully_in_support(&self, _cube : &Cube<C::Type,N>) -> bool {
69 todo!("Not implemented, but not used at the moment") 61 todo!("Not implemented, but not used at the moment")
70 }*/ 62 }*/
71 } 63 }
72 64
73 #[replace_float_literals(C::Type::cast_from(literal))] 65 #[replace_float_literals(C::Type::cast_from(literal))]
74 impl<'a, C : Constant, const N : usize> GlobalAnalysis<C::Type, Bounds<C::Type>> 66 impl<'a, C: Constant, const N: usize> GlobalAnalysis<C::Type, Bounds<C::Type>> for Mollifier<C, N> {
75 for Mollifier<C, N> {
76 #[inline] 67 #[inline]
77 fn global_analysis(&self) -> Bounds<C::Type> { 68 fn global_analysis(&self) -> Bounds<C::Type> {
78 // The function is maximised/minimised where the 2-norm is minimised/maximised. 69 // The function is maximised/minimised where the 2-norm is minimised/maximised.
79 Bounds(0.0, 1.0) 70 Bounds(0.0, 1.0)
80 } 71 }
81 } 72 }
82 73
83 impl<'a, C : Constant, const N : usize> LocalAnalysis<C::Type, Bounds<C::Type>, N> 74 impl<'a, C: Constant, const N: usize> LocalAnalysis<C::Type, Bounds<C::Type>, N>
84 for Mollifier<C, N> { 75 for Mollifier<C, N>
76 {
85 #[inline] 77 #[inline]
86 fn local_analysis(&self, cube : &Cube<C::Type, N>) -> Bounds<C::Type> { 78 fn local_analysis(&self, cube: &Cube<N, C::Type>) -> Bounds<C::Type> {
87 // The function is maximised/minimised where the 2-norm is minimised/maximised. 79 // The function is maximised/minimised where the 2-norm is minimised/maximised.
88 let lower = self.apply(cube.maxnorm_point()); 80 let lower = self.apply(cube.maxnorm_point());
89 let upper = self.apply(cube.minnorm_point()); 81 let upper = self.apply(cube.minnorm_point());
90 Bounds(lower, upper) 82 Bounds(lower, upper)
91 } 83 }
97 /// [https://math.stackexchange.com/questions/4359683/integral-of-the-usual-mollifier-function-finding-its-necessary-constant](). 89 /// [https://math.stackexchange.com/questions/4359683/integral-of-the-usual-mollifier-function-finding-its-necessary-constant]().
98 /// 90 ///
99 /// If `rescaled` is `true`, return the integral of the scaled mollifier that has value one at the 91 /// If `rescaled` is `true`, return the integral of the scaled mollifier that has value one at the
100 /// origin. 92 /// origin.
101 #[inline] 93 #[inline]
102 pub fn mollifier_norm1(n_ : usize, rescaled : bool) -> f64 { 94 pub fn mollifier_norm1(n_: usize, rescaled: bool) -> f64 {
103 assert!(n_ > 0); 95 assert!(n_ > 0);
104 let n = n_ as f64; 96 let n = n_ as f64;
105 let q = 2.0; 97 let q = 2.0;
106 let p = 2.0; 98 let p = 2.0;
107 let base = (2.0*gamma(1.0 + 1.0/p)).powi(n_ as i32) 99 let base = (2.0*gamma(1.0 + 1.0/p)).powi(n_ as i32)
108 /*/ gamma(1.0 + n / p) 100 /*/ gamma(1.0 + n / p)
109 * gamma(1.0 + n / q)*/ 101 * gamma(1.0 + n / q)*/
110 * hyperg_U(1.0 + n / q, 2.0, 1.0); 102 * hyperg_U(1.0 + n / q, 2.0, 1.0);
111 if rescaled { base } else { base / f64::E } 103 if rescaled {
104 base
105 } else {
106 base / f64::E
107 }
112 } 108 }
113 109
114 impl<'a, C : Constant, const N : usize> Norm<C::Type, L1> 110 impl<'a, C: Constant, const N: usize> Norm<L1, C::Type> for Mollifier<C, N> {
115 for Mollifier<C, N> {
116 #[inline] 111 #[inline]
117 fn norm(&self, _ : L1) -> C::Type { 112 fn norm(&self, _: L1) -> C::Type {
118 let ε = self.width.value(); 113 let ε = self.width.value();
119 C::Type::cast_from(mollifier_norm1(N, true)) * ε.powi(N as i32) 114 C::Type::cast_from(mollifier_norm1(N, true)) * ε.powi(N as i32)
120 } 115 }
121 } 116 }
122 117
123 #[replace_float_literals(C::Type::cast_from(literal))] 118 #[replace_float_literals(C::Type::cast_from(literal))]
124 impl<'a, C : Constant, const N : usize> Norm<C::Type, Linfinity> 119 impl<'a, C: Constant, const N: usize> Norm<Linfinity, C::Type> for Mollifier<C, N> {
125 for Mollifier<C, N> {
126 #[inline] 120 #[inline]
127 fn norm(&self, _ : Linfinity) -> C::Type { 121 fn norm(&self, _: Linfinity) -> C::Type {
128 1.0 122 1.0
129 } 123 }
130 } 124 }

Mercurial Repository: pointsource_algs

mercurial