Wed, 07 Dec 2022 06:57:26 +0200
README beautification
5 | 1 | /*! |
2 | Integration with nalgebra. | |
3 | ||
4 | This module mainly implements [`Euclidean`], [`Norm`], [`Dot`], [`Linear`], etc. for [`nalgebra`] | |
5 | matrices and vectors. | |
6 | It also provides [`ToNalgebraRealField`] as a vomit-inducingly ugly workaround to nalgebra | |
7 | force-feeding its own versions of the same basic mathematical methods on `f32` and `f64` as | |
8 | [`num_traits`] does. | |
9 | */ | |
0 | 10 | |
11 | use nalgebra::{ | |
12 | Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar, | |
13 | ClosedMul, ClosedAdd, SimdComplexField, Vector, OVector, RealField, | |
14 | LpNorm, UniformNorm | |
15 | }; | |
16 | use nalgebra::Norm as NalgebraNorm; | |
17 | use nalgebra::base::constraint::{ | |
18 | ShapeConstraint, SameNumberOfRows, SameNumberOfColumns | |
19 | }; | |
20 | use nalgebra::base::dimension::*; | |
21 | use nalgebra::base::allocator::Allocator; | |
22 | use std::ops::Mul; | |
23 | use num_traits::identities::{Zero, One}; | |
24 | use crate::linops::*; | |
5 | 25 | use crate::euclidean::*; |
0 | 26 | use crate::types::Float; |
27 | use crate::norms::*; | |
28 | ||
13
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
29 | impl<SM,SV,N,M,K,E> Apply<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
30 | where SM: Storage<E,N,M>, SV: Storage<E,M,K>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
31 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
32 | DefaultAllocator : Allocator<E,N,K>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
33 | DefaultAllocator : Allocator<E,M,K>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
34 | DefaultAllocator : Allocator<E,N,M>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
35 | DefaultAllocator : Allocator<E,M,N> { |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
36 | type Output = OMatrix<E,N,K>; |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
37 | |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
38 | #[inline] |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
39 | fn apply(&self, x : Matrix<E,M,K,SV>) -> Self::Output { |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
40 | self.mul(x) |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
41 | } |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
42 | } |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
43 | |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
44 | impl<'a, SM,SV,N,M,K,E> Apply<&'a Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
45 | where SM: Storage<E,N,M>, SV: Storage<E,M,K>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
46 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
47 | DefaultAllocator : Allocator<E,N,K>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
48 | DefaultAllocator : Allocator<E,M,K>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
49 | DefaultAllocator : Allocator<E,N,M>, |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
50 | DefaultAllocator : Allocator<E,M,N> { |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
51 | type Output = OMatrix<E,N,K>; |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
52 | |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
53 | #[inline] |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
54 | fn apply(&self, x : &'a Matrix<E,M,K,SV>) -> Self::Output { |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
55 | self.mul(x) |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
56 | } |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
57 | } |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
58 | |
465fa2121ccb
Better Linear and Mapping structure that can provide consuming and reference `apply`.
Tuomo Valkonen <tuomov@iki.fi>
parents:
5
diff
changeset
|
59 | impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> |
0 | 60 | where SM: Storage<E,N,M>, SV: Storage<E,M,K>, |
61 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, | |
62 | DefaultAllocator : Allocator<E,N,K>, | |
63 | DefaultAllocator : Allocator<E,M,K>, | |
64 | DefaultAllocator : Allocator<E,N,M>, | |
65 | DefaultAllocator : Allocator<E,M,N> { | |
66 | type Codomain = OMatrix<E,N,K>; | |
67 | } | |
68 | ||
69 | impl<SM,SV1,SV2,N,M,K,E> GEMV<E, Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM> | |
70 | where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: StorageMut<E,N,K>, | |
71 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, | |
72 | DefaultAllocator : Allocator<E,N,K>, | |
73 | DefaultAllocator : Allocator<E,M,K>, | |
74 | DefaultAllocator : Allocator<E,N,M>, | |
75 | DefaultAllocator : Allocator<E,M,N> { | |
76 | ||
77 | #[inline] | |
78 | fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : &Matrix<E,M,K,SV1>, β : E) { | |
79 | Matrix::gemm(y, α, self, x, β) | |
80 | } | |
81 | ||
82 | #[inline] | |
83 | fn apply_mut<'a>(&self, y : &mut Matrix<E,N,K,SV2>, x : &Matrix<E,M,K,SV1>) { | |
84 | self.mul_to(x, y) | |
85 | } | |
86 | } | |
87 | ||
88 | impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM> | |
89 | where SM: StorageMut<E,M>, SV1: Storage<E,M>, | |
90 | M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, | |
91 | DefaultAllocator : Allocator<E,M> { | |
92 | ||
93 | #[inline] | |
94 | fn axpy(&mut self, α : E, x : &Vector<E,M,SV1>, β : E) { | |
95 | Matrix::axpy(self, α, x, β) | |
96 | } | |
97 | ||
98 | #[inline] | |
99 | fn copy_from(&mut self, y : &Vector<E,M,SV1>) { | |
100 | Matrix::copy_from(self, y) | |
101 | } | |
102 | } | |
103 | ||
104 | impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM> | |
105 | where SM: StorageMut<E,M>, | |
106 | M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float + RealField, | |
107 | DefaultAllocator : Allocator<E,M> { | |
108 | #[inline] | |
109 | fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) { | |
110 | self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) | |
111 | } | |
112 | } | |
113 | ||
114 | impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>,Matrix<E,N,K,SV2>> | |
115 | for Matrix<E,N,M,SM> | |
116 | where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: Storage<E,N,K>, | |
117 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + SimdComplexField, | |
118 | DefaultAllocator : Allocator<E,N,K>, | |
119 | DefaultAllocator : Allocator<E,M,K>, | |
120 | DefaultAllocator : Allocator<E,N,M>, | |
121 | DefaultAllocator : Allocator<E,M,N> { | |
122 | type AdjointCodomain = OMatrix<E,M,K>; | |
123 | type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a; | |
124 | ||
125 | #[inline] | |
126 | fn adjoint(&self) -> Self::Adjoint<'_> { | |
127 | Matrix::adjoint(self) | |
128 | } | |
129 | } | |
130 | ||
131 | impl<E,M,S,Si> Dot<Vector<E,M,Si>,E> | |
132 | for Vector<E,M,S> | |
133 | where M : Dim, | |
134 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One, | |
135 | S : Storage<E,M>, | |
136 | Si : Storage<E,M>, | |
137 | DefaultAllocator : Allocator<E,M> { | |
138 | ||
139 | #[inline] | |
140 | fn dot(&self, other : &Vector<E,M,Si>) -> E { | |
141 | Vector::<E,M,S>::dot(self, other) | |
142 | } | |
143 | } | |
144 | ||
145 | /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`. | |
146 | #[inline] | |
147 | fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>( | |
148 | /*ed: &EuclideanNorm,*/ | |
149 | m1: &Matrix<T, R1, C1, S1>, | |
150 | m2: &Matrix<T, R2, C2, S2>, | |
151 | ) -> T::SimdRealField | |
152 | where | |
153 | T: SimdComplexField, | |
154 | R1: Dim, | |
155 | C1: Dim, | |
156 | S1: Storage<T, R1, C1>, | |
157 | R2: Dim, | |
158 | C2: Dim, | |
159 | S2: Storage<T, R2, C2>, | |
160 | ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2>, | |
161 | { | |
162 | m1.zip_fold(m2, T::SimdRealField::zero(), |acc, a, b| { | |
163 | let diff = a - b; | |
164 | acc + diff.simd_modulus_squared() | |
165 | }) | |
166 | } | |
167 | ||
168 | // TODO: should allow different input storages in `Euclidean`. | |
169 | ||
170 | impl<E,M,S> Euclidean<E> | |
171 | for Vector<E,M,S> | |
172 | where M : Dim, | |
173 | S : StorageMut<E,M>, | |
174 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
175 | DefaultAllocator : Allocator<E,M> { | |
176 | ||
177 | type Output = OVector<E, M>; | |
178 | ||
179 | #[inline] | |
180 | fn similar_origin(&self) -> OVector<E, M> { | |
181 | OVector::zeros_generic(M::from_usize(self.len()), Const) | |
182 | } | |
183 | ||
184 | #[inline] | |
185 | fn norm2_squared(&self) -> E { | |
186 | Vector::<E,M,S>::norm_squared(self) | |
187 | } | |
188 | ||
189 | #[inline] | |
190 | fn dist2_squared(&self, other : &Self) -> E { | |
191 | metric_distance_squared(self, other) | |
192 | } | |
193 | } | |
194 | ||
195 | impl<E,M,S> StaticEuclidean<E> | |
196 | for Vector<E,M,S> | |
197 | where M : DimName, | |
198 | S : StorageMut<E,M>, | |
199 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
200 | DefaultAllocator : Allocator<E,M> { | |
201 | ||
202 | #[inline] | |
203 | fn origin() -> OVector<E, M> { | |
204 | OVector::zeros() | |
205 | } | |
206 | } | |
207 | ||
208 | impl<E,M,S> Norm<E, L1> | |
209 | for Vector<E,M,S> | |
210 | where M : Dim, | |
211 | S : StorageMut<E,M>, | |
212 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
213 | DefaultAllocator : Allocator<E,M> { | |
214 | ||
215 | #[inline] | |
216 | fn norm(&self, _ : L1) -> E { | |
217 | LpNorm(1).norm(self) | |
218 | } | |
219 | } | |
220 | ||
221 | impl<E,M,S> Dist<E, L1> | |
222 | for Vector<E,M,S> | |
223 | where M : Dim, | |
224 | S : StorageMut<E,M>, | |
225 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
226 | DefaultAllocator : Allocator<E,M> { | |
227 | #[inline] | |
228 | fn dist(&self, other : &Self, _ : L1) -> E { | |
229 | LpNorm(1).metric_distance(self, other) | |
230 | } | |
231 | } | |
232 | ||
233 | impl<E,M,S> Norm<E, L2> | |
234 | for Vector<E,M,S> | |
235 | where M : Dim, | |
236 | S : StorageMut<E,M>, | |
237 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
238 | DefaultAllocator : Allocator<E,M> { | |
239 | ||
240 | #[inline] | |
241 | fn norm(&self, _ : L2) -> E { | |
242 | LpNorm(2).norm(self) | |
243 | } | |
244 | } | |
245 | ||
246 | impl<E,M,S> Dist<E, L2> | |
247 | for Vector<E,M,S> | |
248 | where M : Dim, | |
249 | S : StorageMut<E,M>, | |
250 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
251 | DefaultAllocator : Allocator<E,M> { | |
252 | #[inline] | |
253 | fn dist(&self, other : &Self, _ : L2) -> E { | |
254 | LpNorm(2).metric_distance(self, other) | |
255 | } | |
256 | } | |
257 | ||
258 | impl<E,M,S> Norm<E, Linfinity> | |
259 | for Vector<E,M,S> | |
260 | where M : Dim, | |
261 | S : StorageMut<E,M>, | |
262 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
263 | DefaultAllocator : Allocator<E,M> { | |
264 | ||
265 | #[inline] | |
266 | fn norm(&self, _ : Linfinity) -> E { | |
267 | UniformNorm.norm(self) | |
268 | } | |
269 | } | |
270 | ||
271 | impl<E,M,S> Dist<E, Linfinity> | |
272 | for Vector<E,M,S> | |
273 | where M : Dim, | |
274 | S : StorageMut<E,M>, | |
275 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
276 | DefaultAllocator : Allocator<E,M> { | |
277 | #[inline] | |
278 | fn dist(&self, other : &Self, _ : Linfinity) -> E { | |
279 | UniformNorm.metric_distance(self, other) | |
280 | } | |
281 | } | |
282 | ||
5 | 283 | /// Helper trait to hide the symbols of [`nalgebra::RealField`]. |
284 | /// | |
285 | /// By assuming `ToNalgebraRealField` intead of `nalgebra::RealField` as a trait bound, | |
286 | /// functions can piggyback `nalgebra::RealField` without exponsing themselves to it. | |
287 | /// Thus methods from [`num_traits`] can be used directly without similarly named methods | |
288 | /// from [`nalgebra`] conflicting with them. Only when absolutely necessary to work with | |
289 | /// nalgebra, one can convert to the nalgebra view of the same type using the methods of | |
290 | /// this trait. | |
0 | 291 | pub trait ToNalgebraRealField : Float { |
5 | 292 | /// The nalgebra type corresponding to this type. Usually same as `Self`. |
293 | /// | |
294 | /// This type only carries `nalgebra` traits. | |
0 | 295 | type NalgebraType : RealField; |
5 | 296 | /// The “mixed” type corresponding to this type. Usually same as `Self`. |
297 | /// | |
298 | /// This type carries both `num_traits` and `nalgebra` traits. | |
0 | 299 | type MixedType : RealField + Float; |
300 | ||
5 | 301 | /// Convert to the nalgebra view of `self`. |
0 | 302 | fn to_nalgebra(self) -> Self::NalgebraType; |
5 | 303 | |
304 | /// Convert to the mixed (nalgebra and num_traits) view of `self`. | |
0 | 305 | fn to_nalgebra_mixed(self) -> Self::MixedType; |
306 | ||
5 | 307 | /// Convert from the nalgebra view of `self`. |
0 | 308 | fn from_nalgebra(t : Self::NalgebraType) -> Self; |
5 | 309 | |
310 | /// Convert from the mixed (nalgebra and num_traits) view to `self`. | |
0 | 311 | fn from_nalgebra_mixed(t : Self::MixedType) -> Self; |
312 | } | |
313 | ||
314 | impl ToNalgebraRealField for f32 { | |
315 | type NalgebraType = f32; | |
316 | type MixedType = f32; | |
317 | ||
318 | #[inline] | |
319 | fn to_nalgebra(self) -> Self::NalgebraType { self } | |
320 | ||
321 | #[inline] | |
322 | fn to_nalgebra_mixed(self) -> Self::MixedType { self } | |
323 | ||
324 | #[inline] | |
325 | fn from_nalgebra(t : Self::NalgebraType) -> Self { t } | |
326 | ||
327 | #[inline] | |
328 | fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } | |
329 | ||
330 | } | |
331 | ||
332 | impl ToNalgebraRealField for f64 { | |
333 | type NalgebraType = f64; | |
334 | type MixedType = f64; | |
335 | ||
336 | #[inline] | |
337 | fn to_nalgebra(self) -> Self::NalgebraType { self } | |
338 | ||
339 | #[inline] | |
340 | fn to_nalgebra_mixed(self) -> Self::MixedType { self } | |
341 | ||
342 | #[inline] | |
343 | fn from_nalgebra(t : Self::NalgebraType) -> Self { t } | |
344 | ||
345 | #[inline] | |
346 | fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } | |
347 | } | |
348 |