src/norms.rs

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

mercurial