Mercurial > repos > alg_tools / file comparison

comparison: src/norms.rs

src/norms.rs

branch
dev
changeset 100
411c6be29fe5
parent 72
44a4f258a1ff
child 106
1256e7f7f7ad
equal deleted inserted replaced
99:9e5b9fc81c52 100:411c6be29fe5
1 /*! 1 /*!
2 Norms, projections, etc. 2 Norms, projections, etc.
3 */ 3 */
4 4
5 use crate::euclidean::*;
6 use crate::mapping::{Instance, Mapping, Space};
7 use crate::types::*;
5 use serde::Serialize; 8 use serde::Serialize;
6 use std::marker::PhantomData; 9 use std::marker::PhantomData;
7 use crate::types::*;
8 use crate::euclidean::*;
9 use crate::mapping::{Mapping, Space, Instance};
10 10
11 // 11 //
12 // Abstract norms 12 // Abstract norms
13 // 13 //
14 14
15 #[derive(Copy,Clone,Debug)] 15 #[derive(Copy, Clone, Debug)]
16 /// Helper structure to convert a [`NormExponent`] into a [`Mapping`] 16 /// Helper structure to convert a [`NormExponent`] into a [`Mapping`]
17 pub struct NormMapping<F : Float, E : NormExponent>{ 17 pub struct NormMapping<F: Float, E: NormExponent> {
18 pub(crate) exponent : E, 18 pub(crate) exponent: E,
19 _phantoms : PhantomData<F> 19 _phantoms: PhantomData<F>,
20 } 20 }
21 21
22 /// An exponent for norms. 22 /// An exponent for norms.
23 /// 23 ///
24 // Just a collection of desirable attributes for a marker type 24 // Just a collection of desirable attributes for a marker type
25 pub trait NormExponent : Copy + Send + Sync + 'static { 25 pub trait NormExponent: Copy {
26 /// Return the norm as a mappin 26 /// Return the norm as a mappin
27 fn as_mapping<F : Float>(self) -> NormMapping<F, Self> { 27 fn as_mapping<F: Float>(self) -> NormMapping<F, Self> {
28 NormMapping{ exponent : self, _phantoms : PhantomData } 28 NormMapping {
29 exponent: self,
30 _phantoms: PhantomData,
31 }
29 } 32 }
30 } 33 }
31 34
32 /// Exponent type for the 1-[`Norm`]. 35 /// Exponent type for the 1-[`Norm`].
33 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] 36 #[derive(Copy, Debug, Clone, Serialize, Eq, PartialEq)]
34 pub struct L1; 37 pub struct L1;
35 impl NormExponent for L1 {} 38 impl NormExponent for L1 {}
36 39
37 /// Exponent type for the 2-[`Norm`]. 40 /// Exponent type for the 2-[`Norm`].
38 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] 41 #[derive(Copy, Debug, Clone, Serialize, Eq, PartialEq)]
39 pub struct L2; 42 pub struct L2;
40 impl NormExponent for L2 {} 43 impl NormExponent for L2 {}
41 44
42 /// Exponent type for the ∞-[`Norm`]. 45 /// Exponent type for the ∞-[`Norm`].
43 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] 46 #[derive(Copy, Debug, Clone, Serialize, Eq, PartialEq)]
44 pub struct Linfinity; 47 pub struct Linfinity;
45 impl NormExponent for Linfinity {} 48 impl NormExponent for Linfinity {}
46 49
47 /// Exponent type for 2,1-[`Norm`]. 50 /// Exponent type for 2,1-[`Norm`].
48 /// (1-norm over a domain Ω, 2-norm of a vector at each point of the domain.) 51 /// (1-norm over a domain Ω, 2-norm of a vector at each point of the domain.)
49 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] 52 #[derive(Copy, Debug, Clone, Serialize, Eq, PartialEq)]
50 pub struct L21; 53 pub struct L21;
51 impl NormExponent for L21 {} 54 impl NormExponent for L21 {}
52 55
53 /// Norms for pairs (a, b). ‖(a,b)‖ = ‖(‖a‖_A, ‖b‖_B)‖_J 56 /// Norms for pairs (a, b). ‖(a,b)‖ = ‖(‖a‖_A, ‖b‖_B)‖_J
54 /// For use with [`crate::direct_product::Pair`] 57 /// For use with [`crate::direct_product::Pair`]
55 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] 58 #[derive(Copy, Debug, Clone, Serialize, Eq, PartialEq)]
56 pub struct PairNorm<A, B, J>(pub A, pub B, pub J); 59 pub struct PairNorm<A, B, J>(pub A, pub B, pub J);
57 60
58 impl<A, B, J> NormExponent for PairNorm<A, B, J> 61 impl<A, B, J> NormExponent for PairNorm<A, B, J>
59 where A : NormExponent, B : NormExponent, J : NormExponent {} 62 where
60 63 A: NormExponent,
64 B: NormExponent,
65 J: NormExponent,
66 {
67 }
61 68
62 /// A Huber/Moreau–Yosida smoothed [`L1`] norm. (Not a norm itself.) 69 /// A Huber/Moreau–Yosida smoothed [`L1`] norm. (Not a norm itself.)
63 /// 70 ///
64 /// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher 71 /// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher
65 /// values more smoothing. Behaviour with γ < 0 is undefined. 72 /// values more smoothing. Behaviour with γ < 0 is undefined.
66 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] 73 #[derive(Copy, Debug, Clone, Serialize, Eq, PartialEq)]
67 pub struct HuberL1<F : Float>(pub F); 74 pub struct HuberL1<F: Float>(pub F);
68 impl<F : Float> NormExponent for HuberL1<F> {} 75 impl<F: Float> NormExponent for HuberL1<F> {}
69 76
70 /// A Huber/Moreau–Yosida smoothed [`L21`] norm. (Not a norm itself.) 77 /// A Huber/Moreau–Yosida smoothed [`L21`] norm. (Not a norm itself.)
71 /// 78 ///
72 /// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher 79 /// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher
73 /// values more smoothing. Behaviour with γ < 0 is undefined. 80 /// values more smoothing. Behaviour with γ < 0 is undefined.
74 #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] 81 #[derive(Copy, Debug, Clone, Serialize, Eq, PartialEq)]
75 pub struct HuberL21<F : Float>(pub F); 82 pub struct HuberL21<F: Float>(pub F);
76 impl<F : Float> NormExponent for HuberL21<F> {} 83 impl<F: Float> NormExponent for HuberL21<F> {}
77
78 84
79 /// A normed space (type) with exponent or other type `Exponent` for the norm. 85 /// A normed space (type) with exponent or other type `Exponent` for the norm.
80 /// 86 ///
81 /// Use as 87 /// Use as
82 /// ``` 88 /// ```
84 /// # use alg_tools::loc::Loc; 90 /// # use alg_tools::loc::Loc;
85 /// let x = Loc([1.0, 2.0, 3.0]); 91 /// let x = Loc([1.0, 2.0, 3.0]);
86 /// 92 ///
87 /// println!("{}, {} {}", x.norm(L1), x.norm(L2), x.norm(Linfinity)) 93 /// println!("{}, {} {}", x.norm(L1), x.norm(L2), x.norm(Linfinity))
88 /// ``` 94 /// ```
89 pub trait Norm<F : Num, Exponent : NormExponent> { 95 pub trait Norm<F: Num, Exponent: NormExponent> {
90 /// Calculate the norm. 96 /// Calculate the norm.
91 fn norm(&self, _p : Exponent) -> F; 97 fn norm(&self, _p: Exponent) -> F;
92 } 98 }
93 99
94 /// Indicates that the `Self`-[`Norm`] is dominated by the `Exponent`-`Norm` on the space 100 /// Indicates that the `Self`-[`Norm`] is dominated by the `Exponent`-`Norm` on the space
95 /// `Elem` with the corresponding field `F`. 101 /// `Elem` with the corresponding field `F`.
96 pub trait Dominated<F : Num, Exponent : NormExponent, Elem> { 102 pub trait Dominated<F: Num, Exponent: NormExponent, Elem> {
97 /// Indicates the factor $c$ for the inequality $‖x‖ ≤ C ‖x‖_p$. 103 /// Indicates the factor $c$ for the inequality $‖x‖ ≤ C ‖x‖_p$.
98 fn norm_factor(&self, p : Exponent) -> F; 104 fn norm_factor(&self, p: Exponent) -> F;
99 /// Given a norm-value $‖x‖_p$, calculates $C‖x‖_p$ such that $‖x‖ ≤ C‖x‖_p$ 105 /// Given a norm-value $‖x‖_p$, calculates $C‖x‖_p$ such that $‖x‖ ≤ C‖x‖_p$
100 #[inline] 106 #[inline]
101 fn from_norm(&self, p_norm : F, p : Exponent) -> F { 107 fn from_norm(&self, p_norm: F, p: Exponent) -> F {
102 p_norm * self.norm_factor(p) 108 p_norm * self.norm_factor(p)
103 } 109 }
104 } 110 }
105 111
106 /// Trait for distances with respect to a norm. 112 /// Trait for distances with respect to a norm.
107 pub trait Dist<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Space { 113 pub trait Dist<F: Num, Exponent: NormExponent>: Norm<F, Exponent> + Space {
108 /// Calculate the distance 114 /// Calculate the distance
109 fn dist<I : Instance<Self>>(&self, other : I, _p : Exponent) -> F; 115 fn dist<I: Instance<Self>>(&self, other: I, _p: Exponent) -> F;
110 } 116 }
111 117
112 /// Trait for Euclidean projections to the `Exponent`-[`Norm`]-ball. 118 /// Trait for Euclidean projections to the `Exponent`-[`Norm`]-ball.
113 /// 119 ///
114 /// Use as 120 /// Use as
117 /// # use alg_tools::loc::Loc; 123 /// # use alg_tools::loc::Loc;
118 /// let x = Loc([1.0, 2.0, 3.0]); 124 /// let x = Loc([1.0, 2.0, 3.0]);
119 /// 125 ///
120 /// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity)); 126 /// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity));
121 /// ``` 127 /// ```
122 pub trait Projection<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Sized 128 pub trait Projection<F: Num, Exponent: NormExponent>: Norm<F, Exponent> + Sized
123 where F : Float { 129 where
130 F: Float,
131 {
124 /// Projection of `self` to the `q`-norm-ball of radius ρ. 132 /// Projection of `self` to the `q`-norm-ball of radius ρ.
125 fn proj_ball(mut self, ρ : F, q : Exponent) -> Self { 133 fn proj_ball(mut self, ρ: F, q: Exponent) -> Self {
126 self.proj_ball_mut(ρ, q); 134 self.proj_ball_mut(ρ, q);
127 self 135 self
128 } 136 }
129 137
130 /// In-place projection of `self` to the `q`-norm-ball of radius ρ. 138 /// In-place projection of `self` to the `q`-norm-ball of radius ρ.
131 fn proj_ball_mut(&mut self, ρ : F, q : Exponent); 139 fn proj_ball_mut(&mut self, ρ: F, q: Exponent);
132 } 140 }
133 141
134 /*impl<F : Float, E : Euclidean<F>> Norm<F, L2> for E { 142 /*impl<F : Float, E : Euclidean<F>> Norm<F, L2> for E {
135 #[inline] 143 #[inline]
136 fn norm(&self, _p : L2) -> F { self.norm2() } 144 fn norm(&self, _p : L2) -> F { self.norm2() }
137 145
138 fn dist(&self, other : &Self, _p : L2) -> F { self.dist2(other) } 146 fn dist(&self, other : &Self, _p : L2) -> F { self.dist2(other) }
139 }*/ 147 }*/
140 148
141 impl<F : Float, E : Euclidean<F> + Norm<F, L2>> Projection<F, L2> for E { 149 impl<F: Float, E: Euclidean<F> + Norm<F, L2>> Projection<F, L2> for E {
142 #[inline] 150 #[inline]
143 fn proj_ball(self, ρ : F, _p : L2) -> Self { self.proj_ball2(ρ) } 151 fn proj_ball(self, ρ: F, _p: L2) -> Self {
144 152 self.proj_ball2(ρ)
145 #[inline] 153 }
146 fn proj_ball_mut(&mut self, ρ : F, _p : L2) { self.proj_ball2_mut(ρ) } 154
147 } 155 #[inline]
148 156 fn proj_ball_mut(&mut self, ρ: F, _p: L2) {
149 impl<F : Float> HuberL1<F> { 157 self.proj_ball2_mut(ρ)
150 fn apply(self, xnsq : F) -> F { 158 }
159 }
160
161 impl<F: Float> HuberL1<F> {
162 fn apply(self, xnsq: F) -> F {
151 let HuberL1(γ) = self; 163 let HuberL1(γ) = self;
152 let xn = xnsq.sqrt(); 164 let xn = xnsq.sqrt();
153 if γ == F::ZERO { 165 if γ == F::ZERO {
154 xn 166 xn
155 } else { 167 } else {
156 if xn > γ { 168 if xn > γ {
157 xn-γ / F::TWO 169 xn - γ / F::TWO
158 } else if xn<(-γ) { 170 } else if xn < (-γ) {
159 -xn-γ / F::TWO 171 -xn - γ / F::TWO
160 } else { 172 } else {
161 xnsq / (F::TWO * γ) 173 xnsq / (F::TWO * γ)
162 } 174 }
163 } 175 }
164 } 176 }
165 } 177 }
166 178
167 impl<F : Float, E : Euclidean<F>> Norm<F, HuberL1<F>> for E { 179 impl<F: Float, E: Euclidean<F>> Norm<F, HuberL1<F>> for E {
168 fn norm(&self, huber : HuberL1<F>) -> F { 180 fn norm(&self, huber: HuberL1<F>) -> F {
169 huber.apply(self.norm2_squared()) 181 huber.apply(self.norm2_squared())
170 } 182 }
171 } 183 }
172 184
173 impl<F : Float, E : Euclidean<F>> Dist<F, HuberL1<F>> for E { 185 impl<F: Float, E: Euclidean<F>> Dist<F, HuberL1<F>> for E {
174 fn dist<I : Instance<Self>>(&self, other : I, huber : HuberL1<F>) -> F { 186 fn dist<I: Instance<Self>>(&self, other: I, huber: HuberL1<F>) -> F {
175 huber.apply(self.dist2_squared(other)) 187 huber.apply(self.dist2_squared(other))
176 } 188 }
177 } 189 }
178 190
179 // impl<F : Float, E : Norm<F, L2>> Norm<F, L21> for Vec<E> { 191 // impl<F : Float, E : Norm<F, L2>> Norm<F, L21> for Vec<E> {
188 // } 200 // }
189 // } 201 // }
190 202
191 impl<E, F, Domain> Mapping<Domain> for NormMapping<F, E> 203 impl<E, F, Domain> Mapping<Domain> for NormMapping<F, E>
192 where 204 where
193 F : Float, 205 F: Float,
194 E : NormExponent, 206 E: NormExponent,
195 Domain : Space + Norm<F, E>, 207 Domain: Space + Norm<F, E>,
196 { 208 {
197 type Codomain = F; 209 type Codomain = F;
198 210
199 #[inline] 211 #[inline]
200 fn apply<I : Instance<Domain>>(&self, x : I) -> F { 212 fn apply<I: Instance<Domain>>(&self, x: I) -> F {
201 x.eval(|r| r.norm(self.exponent)) 213 x.eval(|r| r.norm(self.exponent))
202 } 214 }
203 } 215 }
204 216
205 pub trait Normed<F : Num = f64> : Space + Norm<F, Self::NormExp> { 217 pub trait Normed<F: Num = f64>: Space + Norm<F, Self::NormExp> {
206 type NormExp : NormExponent; 218 type NormExp: NormExponent;
207 219
208 fn norm_exponent(&self) -> Self::NormExp; 220 fn norm_exponent(&self) -> Self::NormExp;
209 221
210 #[inline] 222 #[inline]
211 fn norm_(&self) -> F { 223 fn norm_(&self) -> F {
215 // fn similar_origin(&self) -> Self; 227 // fn similar_origin(&self) -> Self;
216 228
217 fn is_zero(&self) -> bool; 229 fn is_zero(&self) -> bool;
218 } 230 }
219 231
220 pub trait HasDual<F : Num = f64> : Normed<F> { 232 pub trait HasDual<F: Num = f64>: Normed<F> {
221 type DualSpace : Normed<F>; 233 type DualSpace: Normed<F>;
222 } 234 }
223 235
224 /// Automatically implemented trait for reflexive spaces 236 /// Automatically implemented trait for reflexive spaces
225 pub trait Reflexive<F : Num = f64> : HasDual<F> 237 pub trait Reflexive<F: Num = f64>: HasDual<F>
226 where 238 where
227 Self::DualSpace : HasDual<F, DualSpace = Self> 239 Self::DualSpace: HasDual<F, DualSpace = Self>,
228 { } 240 {
229 241 }
230 impl<F : Num, X : HasDual<F>> Reflexive<F> for X 242
231 where 243 impl<F: Num, X: HasDual<F>> Reflexive<F> for X where X::DualSpace: HasDual<F, DualSpace = X> {}
232 X::DualSpace : HasDual<F, DualSpace = X> 244
233 { } 245 pub trait HasDualExponent: NormExponent {
234 246 type DualExp: NormExponent;
235 pub trait HasDualExponent : NormExponent {
236 type DualExp : NormExponent;
237 247
238 fn dual_exponent(&self) -> Self::DualExp; 248 fn dual_exponent(&self) -> Self::DualExp;
239 } 249 }
240 250
241 impl HasDualExponent for L2 { 251 impl HasDualExponent for L2 {
242 type DualExp = L2; 252 type DualExp = L2;
243 253
244 #[inline] 254 #[inline]
245 fn dual_exponent(&self) -> Self::DualExp { 255 fn dual_exponent(&self) -> Self::DualExp {
246 L2 256 L2
247 } 257 }
248 } 258 }
249 259
250 impl HasDualExponent for L1 { 260 impl HasDualExponent for L1 {
251 type DualExp = Linfinity; 261 type DualExp = Linfinity;
252 262
253 #[inline] 263 #[inline]
254 fn dual_exponent(&self) -> Self::DualExp { 264 fn dual_exponent(&self) -> Self::DualExp {
255 Linfinity 265 Linfinity
256 } 266 }
257 } 267 }
258 268
259
260 impl HasDualExponent for Linfinity { 269 impl HasDualExponent for Linfinity {
261 type DualExp = L1; 270 type DualExp = L1;
262 271
263 #[inline] 272 #[inline]
264 fn dual_exponent(&self) -> Self::DualExp { 273 fn dual_exponent(&self) -> Self::DualExp {
265 L1 274 L1
266 } 275 }
267 } 276 }
269 #[macro_export] 278 #[macro_export]
270 macro_rules! impl_weighted_norm { 279 macro_rules! impl_weighted_norm {
271 ($exponent : ty) => { 280 ($exponent : ty) => {
272 impl<C, F, D> Norm<F, Weighted<$exponent, C>> for D 281 impl<C, F, D> Norm<F, Weighted<$exponent, C>> for D
273 where 282 where
274 F : Float, 283 F: Float,
275 D : Norm<F, $exponent>, 284 D: Norm<F, $exponent>,
276 C : Constant<Type = F>, 285 C: Constant<Type = F>,
277 { 286 {
278 fn norm(&self, e : Weighted<$exponent, C>) -> F { 287 fn norm(&self, e: Weighted<$exponent, C>) -> F {
279 let v = e.weight.value(); 288 let v = e.weight.value();
280 assert!(v > F::ZERO); 289 assert!(v > F::ZERO);
281 v * self.norm(e.base_fn) 290 v * self.norm(e.base_fn)
282 } 291 }
283 } 292 }
284 293
285 impl<C : Constant> NormExponent for Weighted<$exponent, C> {} 294 impl<C: Constant> NormExponent for Weighted<$exponent, C> {}
286 295
287 impl<C : Constant> HasDualExponent for Weighted<$exponent, C> 296 impl<C: Constant> HasDualExponent for Weighted<$exponent, C>
288 where $exponent : HasDualExponent { 297 where
298 $exponent: HasDualExponent,
299 {
289 type DualExp = Weighted<<$exponent as HasDualExponent>::DualExp, C::Type>; 300 type DualExp = Weighted<<$exponent as HasDualExponent>::DualExp, C::Type>;
290 301
291 fn dual_exponent(&self) -> Self::DualExp { 302 fn dual_exponent(&self) -> Self::DualExp {
292 Weighted { 303 Weighted {
293 weight : C::Type::ONE / self.weight.value(), 304 weight: C::Type::ONE / self.weight.value(),
294 base_fn : self.base_fn.dual_exponent() 305 base_fn: self.base_fn.dual_exponent(),
295 } 306 }
296 } 307 }
297 } 308 }
298 309
299 impl<C, F, T> Projection<F, Weighted<$exponent , C>> for T 310 impl<C, F, T> Projection<F, Weighted<$exponent, C>> for T
300 where 311 where
301 T : Projection<F, $exponent >, 312 T: Projection<F, $exponent>,
302 F : Float, 313 F: Float,
303 C : Constant<Type = F>, 314 C: Constant<Type = F>,
304 { 315 {
305 fn proj_ball(self, ρ : F, q : Weighted<$exponent , C>) -> Self { 316 fn proj_ball(self, ρ: F, q: Weighted<$exponent, C>) -> Self {
306 self.proj_ball(ρ / q.weight.value(), q.base_fn) 317 self.proj_ball(ρ / q.weight.value(), q.base_fn)
307 } 318 }
308 319
309 fn proj_ball_mut(&mut self, ρ : F, q : Weighted<$exponent , C>) { 320 fn proj_ball_mut(&mut self, ρ: F, q: Weighted<$exponent, C>) {
310 self.proj_ball_mut(ρ / q.weight.value(), q.base_fn) 321 self.proj_ball_mut(ρ / q.weight.value(), q.base_fn)
311 } 322 }
312 } 323 }
313 } 324 };
314 } 325 }
315 326
316 //impl_weighted_norm!(L1); 327 //impl_weighted_norm!(L1);
317 //impl_weighted_norm!(L2); 328 //impl_weighted_norm!(L2);
318 //impl_weighted_norm!(Linfinity); 329 //impl_weighted_norm!(Linfinity);
319

Mercurial Repository: alg_tools

mercurial