Mercurial > repos > pointsource_algs / file comparison

comparison: src/kernels/hat.rs

src/kernels/hat.rs

branch
dev
changeset 61
4f468d35fa29
parent 35
b087e3eab191
equal deleted inserted replaced
60:9738b51d90d7 61:4f468d35fa29
1 //! Implementation of the hat function 1 //! Implementation of the hat function
2 2
3 use crate::types::Lipschitz;
4 use alg_tools::bisection_tree::{Constant, Support};
5 use alg_tools::bounds::{Bounded, Bounds, GlobalAnalysis, LocalAnalysis};
6 use alg_tools::error::DynResult;
7 use alg_tools::loc::Loc;
8 use alg_tools::mapping::{Instance, Mapping};
9 use alg_tools::maputil::array_init;
10 use alg_tools::norms::*;
11 use alg_tools::sets::Cube;
12 use alg_tools::types::*;
3 use numeric_literals::replace_float_literals; 13 use numeric_literals::replace_float_literals;
4 use serde::Serialize; 14 use serde::Serialize;
5 use alg_tools::types::*;
6 use alg_tools::norms::*;
7 use alg_tools::loc::Loc;
8 use alg_tools::sets::Cube;
9 use alg_tools::bisection_tree::{
10 Support,
11 Constant,
12 Bounds,
13 LocalAnalysis,
14 GlobalAnalysis,
15 Bounded,
16 };
17 use alg_tools::mapping::{Mapping, Instance};
18 use alg_tools::maputil::array_init;
19 use crate::types::Lipschitz;
20 15
21 /// Representation of the hat function $f(x)=1-\\|x\\|\_1/ε$ of `width` $ε$ on $ℝ^N$. 16 /// Representation of the hat function $f(x)=1-\\|x\\|\_1/ε$ of `width` $ε$ on $ℝ^N$.
22 #[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)] 17 #[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
23 pub struct Hat<C : Constant, const N : usize> { 18 pub struct Hat<C: Constant, const N: usize> {
24 /// The parameter $ε>0$. 19 /// The parameter $ε>0$.
25 pub width : C, 20 pub width: C,
26 } 21 }
27 22
28 #[replace_float_literals(C::Type::cast_from(literal))] 23 #[replace_float_literals(C::Type::cast_from(literal))]
29 impl<'a, C : Constant, const N : usize> Mapping<Loc<C::Type, N>> for Hat<C, N> { 24 impl<'a, C: Constant, const N: usize> Mapping<Loc<N, C::Type>> for Hat<C, N> {
30 type Codomain = C::Type; 25 type Codomain = C::Type;
31 26
32 #[inline] 27 #[inline]
33 fn apply<I : Instance<Loc<C::Type, N>>>(&self, x : I) -> Self::Codomain { 28 fn apply<I: Instance<Loc<N, C::Type>>>(&self, x: I) -> Self::Codomain {
34 let ε = self.width.value(); 29 let ε = self.width.value();
35 0.0.max(1.0-x.cow().norm(L1)/ε) 30 0.0.max(1.0 - x.decompose().norm(L1) / ε)
36 } 31 }
37 } 32 }
38 33
39 #[replace_float_literals(C::Type::cast_from(literal))] 34 #[replace_float_literals(C::Type::cast_from(literal))]
40 impl<'a, C : Constant, const N : usize> Support<C::Type, N> for Hat<C, N> { 35 impl<'a, C: Constant, const N: usize> Support<N, C::Type> for Hat<C, N> {
41 #[inline] 36 #[inline]
42 fn support_hint(&self) -> Cube<C::Type,N> { 37 fn support_hint(&self) -> Cube<N, C::Type> {
43 let ε = self.width.value(); 38 let ε = self.width.value();
44 array_init(|| [-ε, ε]).into() 39 array_init(|| [-ε, ε]).into()
45 } 40 }
46 41
47 #[inline] 42 #[inline]
48 fn in_support(&self, x : &Loc<C::Type,N>) -> bool { 43 fn in_support(&self, x: &Loc<N, C::Type>) -> bool {
49 x.norm(L1) < self.width.value() 44 x.norm(L1) < self.width.value()
50 } 45 }
51 46
52 /*fn fully_in_support(&self, _cube : &Cube<C::Type,N>) -> bool { 47 /*fn fully_in_support(&self, _cube : &Cube<C::Type,N>) -> bool {
53 todo!("Not implemented, but not used at the moment") 48 todo!("Not implemented, but not used at the moment")
54 }*/ 49 }*/
55 50
56 #[inline] 51 #[inline]
57 fn bisection_hint(&self, cube : &Cube<C::Type,N>) -> [Option<C::Type>; N] { 52 fn bisection_hint(&self, cube: &Cube<N, C::Type>) -> [Option<C::Type>; N] {
58 let ε = self.width.value(); 53 let ε = self.width.value();
59 cube.map(|a, b| { 54 cube.map(|a, b| {
60 if a < 1.0 { 55 if a < 1.0 {
61 if 1.0 < b { 56 if 1.0 < b {
62 Some(1.0) 57 Some(1.0)
63 } else { 58 } else {
64 if a < -ε { 59 if a < -ε {
65 if b > -ε { Some(-ε) } else { None } 60 if b > -ε {
61 Some(-ε)
62 } else {
63 None
64 }
66 } else { 65 } else {
67 None 66 None
68 } 67 }
69 } 68 }
70 } else { 69 } else {
71 if b > ε { Some(ε) } else { None } 70 if b > ε {
71 Some(ε)
72 } else {
73 None
74 }
72 } 75 }
73 }); 76 });
74 todo!("also diagonals") 77 todo!("also diagonals")
75 } 78 }
76 } 79 }
77 80
78
79 #[replace_float_literals(C::Type::cast_from(literal))] 81 #[replace_float_literals(C::Type::cast_from(literal))]
80 impl<'a, C : Constant, const N : usize> 82 impl<'a, C: Constant, const N: usize> GlobalAnalysis<C::Type, Bounds<C::Type>> for Hat<C, N> {
81 GlobalAnalysis<C::Type, Bounds<C::Type>>
82 for Hat<C, N> {
83 #[inline] 83 #[inline]
84 fn global_analysis(&self) -> Bounds<C::Type> { 84 fn global_analysis(&self) -> Bounds<C::Type> {
85 Bounds(0.0, 1.0) 85 Bounds(0.0, 1.0)
86 } 86 }
87 } 87 }
88 88
89 #[replace_float_literals(C::Type::cast_from(literal))] 89 #[replace_float_literals(C::Type::cast_from(literal))]
90 impl<'a, C : Constant, const N : usize> Lipschitz<L1> for Hat<C, N> { 90 impl<'a, C: Constant, const N: usize> Lipschitz<L1> for Hat<C, N> {
91 type FloatType = C::Type; 91 type FloatType = C::Type;
92 92
93 fn lipschitz_factor(&self, _l1 : L1) -> Option<C::Type> { 93 fn lipschitz_factor(&self, _l1: L1) -> DynResult<C::Type> {
94 Some(1.0/self.width.value()) 94 Ok(1.0 / self.width.value())
95 } 95 }
96 } 96 }
97 97
98 #[replace_float_literals(C::Type::cast_from(literal))] 98 #[replace_float_literals(C::Type::cast_from(literal))]
99 impl<'a, C : Constant, const N : usize> Lipschitz<L2> for Hat<C, N> { 99 impl<'a, C: Constant, const N: usize> Lipschitz<L2> for Hat<C, N> {
100 type FloatType = C::Type; 100 type FloatType = C::Type;
101 101
102 fn lipschitz_factor(&self, _l2 : L2) -> Option<C::Type> { 102 fn lipschitz_factor(&self, _l2: L2) -> DynResult<C::Type> {
103 self.lipschitz_factor(L1).map(|l1| 103 self.lipschitz_factor(L1)
104 <L2 as Dominated<C::Type, L1, Loc<C::Type,N>>>::from_norm(&L2, l1, L1) 104 .map(|l1| <L2 as Dominated<C::Type, L1, Loc<N, C::Type>>>::from_norm(&L2, l1, L1))
105 )
106 } 105 }
107 } 106 }
108 107
109 impl<'a, C : Constant, const N : usize> 108 impl<'a, C: Constant, const N: usize> LocalAnalysis<C::Type, Bounds<C::Type>, N> for Hat<C, N> {
110 LocalAnalysis<C::Type, Bounds<C::Type>, N>
111 for Hat<C, N> {
112 #[inline] 109 #[inline]
113 fn local_analysis(&self, cube : &Cube<C::Type, N>) -> Bounds<C::Type> { 110 fn local_analysis(&self, cube: &Cube<N, C::Type>) -> Bounds<C::Type> {
114 // The function is maximised/minimised where the 1-norm is minimised/maximised. 111 // The function is maximised/minimised where the 1-norm is minimised/maximised.
115 let lower = self.apply(cube.maxnorm_point()); 112 let lower = self.apply(cube.maxnorm_point());
116 let upper = self.apply(cube.minnorm_point()); 113 let upper = self.apply(cube.minnorm_point());
117 Bounds(lower, upper) 114 Bounds(lower, upper)
118 } 115 }
119 } 116 }
120 117
121 #[replace_float_literals(C::Type::cast_from(literal))] 118 #[replace_float_literals(C::Type::cast_from(literal))]
122 impl<'a, C : Constant, const N : usize> 119 impl<'a, C: Constant, const N: usize> Norm<Linfinity, C::Type> for Hat<C, N> {
123 Norm<C::Type, Linfinity>
124 for Hat<C, N> {
125 #[inline] 120 #[inline]
126 fn norm(&self, _ : Linfinity) -> C::Type { 121 fn norm(&self, _: Linfinity) -> C::Type {
127 self.bounds().upper() 122 self.bounds().upper()
128 } 123 }
129 } 124 }
130

Mercurial Repository: pointsource_algs

mercurial