Mercurial > repos > alg_tools / file comparison

comparison: src/nalgebra_support.rs

src/nalgebra_support.rs

changeset 198
3868555d135c
parent 184
b7b60b3b3eff
equal deleted inserted replaced
94:1f19c6bbf07b 198:3868555d135c
6 It also provides [`ToNalgebraRealField`] as a vomit-inducingly ugly workaround to nalgebra 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 7 force-feeding its own versions of the same basic mathematical methods on `f32` and `f64` as
8 [`num_traits`] does. 8 [`num_traits`] does.
9 */ 9 */
10 10
11 use crate::euclidean::*;
12 use crate::instance::{Decomposition, Instance, MyCow, Ownable, Space};
13 use crate::linops::*;
14 use crate::norms::*;
15 use crate::types::Float;
16 use nalgebra::base::allocator::Allocator;
17 use nalgebra::base::constraint::{DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint};
18 use nalgebra::base::dimension::*;
11 use nalgebra::{ 19 use nalgebra::{
12 Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar, 20 ArrayStorage, ClosedAddAssign, ClosedMulAssign, DefaultAllocator, Dim, LpNorm, Matrix,
13 ClosedAddAssign, ClosedMulAssign, SimdComplexField, Vector, OVector, RealField, 21 MatrixView, OMatrix, OVector, RawStorage, RealField, SimdComplexField, Storage, StorageMut,
14 LpNorm, UniformNorm 22 UniformNorm, VecStorage, Vector, ViewStorage, ViewStorageMut, U1,
15 }; 23 };
16 use nalgebra::base::constraint::{ 24 use num_traits::identities::Zero;
17 ShapeConstraint, SameNumberOfRows, SameNumberOfColumns
18 };
19 use nalgebra::base::dimension::*;
20 use nalgebra::base::allocator::Allocator;
21 use std::ops::Mul; 25 use std::ops::Mul;
22 use num_traits::identities::{Zero, One}; 26
23 use crate::linops::*; 27 impl<S, M, N, E> Ownable for Matrix<E, M, N, S>
24 use crate::euclidean::*; 28 where
25 use crate::mapping::{Space, BasicDecomposition}; 29 S: Storage<E, M, N>,
26 use crate::types::Float; 30 M: Dim,
27 use crate::norms::*; 31 N: Dim,
28 use crate::instance::Instance; 32 E: Float,
29 33 DefaultAllocator: Allocator<M, N>,
30 impl<SM,N,M,E> Space for Matrix<E,N,M,SM> 34 {
31 where 35 type OwnedVariant = OMatrix<E, M, N>;
32 SM: Storage<E,N,M> + Clone, 36
33 N : Dim, M : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, 37 #[inline]
34 DefaultAllocator : Allocator<N,M>, 38 fn into_owned(self) -> Self::OwnedVariant {
35 { 39 Matrix::into_owned(self)
36 type Decomp = BasicDecomposition; 40 }
37 } 41
38 42 /// Returns an owned instance of a reference.
39 impl<SM,SV,N,M,K,E> Mapping<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> 43 fn clone_owned(&self) -> Self::OwnedVariant {
40 where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone, 44 Matrix::clone_owned(self)
41 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, 45 }
42 DefaultAllocator : Allocator<N,K>, 46
43 DefaultAllocator : Allocator<M,K>, 47 fn cow_owned<'b>(self) -> MyCow<'b, Self::OwnedVariant>
44 DefaultAllocator : Allocator<N,M>, 48 where
45 DefaultAllocator : Allocator<M,N> { 49 Self: 'b,
46 type Codomain = OMatrix<E,N,K>; 50 {
47 51 MyCow::Owned(self.into_owned())
48 #[inline] 52 }
49 fn apply<I : Instance<Matrix<E,M,K,SV>>>( 53
50 &self, x : I 54 fn ref_cow_owned<'b>(&'b self) -> MyCow<'b, Self::OwnedVariant>
51 ) -> Self::Codomain { 55 where
52 x.either(|owned| self.mul(owned), |refr| self.mul(refr)) 56 Self: 'b,
53 } 57 {
54 } 58 MyCow::Owned(self.clone_owned())
55 59 }
56 60 }
57 impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> 61
58 where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone, 62 trait StridesOk<E, N, M = U1, S = <DefaultAllocator as Allocator<N, M>>::Buffer<E>>:
59 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, 63 DimEq<Dyn, S::RStride> + DimEq<Dyn, S::CStride>
60 DefaultAllocator : Allocator<N,K>, 64 where
61 DefaultAllocator : Allocator<M,K>, 65 S: RawStorage<E, N, M>,
62 DefaultAllocator : Allocator<N,M>, 66 E: Float,
63 DefaultAllocator : Allocator<M,N> { 67 N: Dim,
64 } 68 M: Dim,
65 69 DefaultAllocator: Allocator<N, M>,
66 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 {
67 where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: StorageMut<E,N,K>, 71 }
68 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float, 72
69 DefaultAllocator : Allocator<N,K>, 73 impl<E, M> StridesOk<E, Dyn, M, VecStorage<E, Dyn, M>> for ShapeConstraint
70 DefaultAllocator : Allocator<M,K>, 74 where
71 DefaultAllocator : Allocator<N,M>, 75 M: Dim,
72 DefaultAllocator : Allocator<M,N> { 76 E: Float,
73 77 DefaultAllocator: Allocator<Dyn, M>,
74 #[inline] 78 {
75 fn gemv<I : Instance<Matrix<E,M,K,SV1>>>( 79 }
76 &self, y : &mut Matrix<E,N,K,SV2>, α : E, x : I, β : E 80
81 impl<E, const N: usize, const M: usize> StridesOk<E, Const<N>, Const<M>, ArrayStorage<E, N, M>>
82 for ShapeConstraint
83 where
84 E: Float,
85 {
86 }
87
88 macro_rules! strides_ok {
89 ($R:ty, $C:ty where $($qual:tt)*) => {
90 impl<'a, E, N, M, $($qual)*> StridesOk<E, N, M, ViewStorage<'a, E, N, M, $R, $C>> for ShapeConstraint
91 where
92 N: Dim,
93 M: Dim,
94 E: Float,
95 DefaultAllocator: Allocator<N, M>,
96 {
97 }
98 impl<'a, E, N, M, $($qual)*> StridesOk<E, N, M, ViewStorageMut<'a, E, N, M, $R, $C>> for ShapeConstraint
99 where
100 N: Dim,
101 M: Dim,
102 E: Float,
103 DefaultAllocator: Allocator<N, M>,
104 {
105 }
106 };
107 }
108
109 strides_ok!(Dyn, Dyn where );
110 strides_ok!(Dyn, Const<C> where const C : usize);
111 strides_ok!(Const<R>, Dyn where const R : usize);
112 strides_ok!(Const<R>, Const<C> where const R : usize, const C : usize);
113
114 impl<SM, N, M, E> Space for Matrix<E, N, M, SM>
115 where
116 SM: Storage<E, N, M>,
117 N: Dim,
118 M: Dim,
119 E: Float,
120 DefaultAllocator: Allocator<N, M>,
121 ShapeConstraint: StridesOk<E, N, M>,
122 {
123 type Principal = OMatrix<E, N, M>;
124 type Decomp = MatrixDecomposition;
125 }
126
127 #[derive(Copy, Clone, Debug)]
128 pub struct MatrixDecomposition;
129
130 impl<E, M, K, S> Decomposition<Matrix<E, M, K, S>> for MatrixDecomposition
131 where
132 S: Storage<E, M, K>,
133 M: Dim,
134 K: Dim,
135 E: Float,
136 DefaultAllocator: Allocator<M, K>,
137 ShapeConstraint: StridesOk<E, M, K>,
138 {
139 type Decomposition<'b>
140 = MyCow<'b, OMatrix<E, M, K>>
141 where
142 Matrix<E, M, K, S>: 'b;
143 type Reference<'b>
144 = MatrixView<'b, E, M, K, Dyn, Dyn>
145 where
146 Matrix<E, M, K, S>: 'b;
147
148 #[inline]
149 fn lift<'b>(r: Self::Reference<'b>) -> Self::Decomposition<'b>
150 where
151 S: 'b,
152 {
153 MyCow::Owned(r.into_owned())
154 }
155 }
156
157 impl<S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition> for Matrix<E, M, K, S2>
158 where
159 S1: Storage<E, M, K>,
160 S2: Storage<E, M, K>,
161 M: Dim,
162 K: Dim,
163 E: Float,
164 DefaultAllocator: Allocator<M, K>,
165 ShapeConstraint: StridesOk<E, M, K, S2> + StridesOk<E, M, K>,
166 {
167 #[inline]
168 fn eval_ref<'b, R>(
169 &'b self,
170 f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
171 ) -> R
172 where
173 Self: 'b,
174 Matrix<E, M, K, S1>: 'b,
175 {
176 f(self.as_view::<M, K, Dyn, Dyn>())
177 }
178
179 #[inline]
180 fn own(self) -> OMatrix<E, M, K> {
181 self.into_owned()
182 }
183
184 #[inline]
185 fn cow<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
186 where
187 Self: 'b,
188 {
189 self.cow_owned()
190 }
191
192 #[inline]
193 fn decompose<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
194 where
195 Self: 'b,
196 {
197 self.cow_owned()
198 }
199 }
200
201 impl<'a, S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition>
202 for &'a Matrix<E, M, K, S2>
203 where
204 S1: Storage<E, M, K>,
205 S2: Storage<E, M, K>,
206 M: Dim,
207 K: Dim,
208 E: Float,
209 DefaultAllocator: Allocator<M, K>,
210 ShapeConstraint: StridesOk<E, M, K, S2> + StridesOk<E, M, K>,
211 {
212 fn eval_ref<'b, R>(
213 &'b self,
214 f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
215 ) -> R
216 where
217 Self: 'b,
218 Matrix<E, M, K, S1>: 'b,
219 {
220 f((*self).as_view::<M, K, Dyn, Dyn>())
221 }
222
223 #[inline]
224 fn own(self) -> OMatrix<E, M, K> {
225 self.into_owned()
226 }
227
228 #[inline]
229 fn cow<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
230 where
231 Self: 'b,
232 {
233 self.cow_owned()
234 }
235
236 #[inline]
237 fn decompose<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
238 where
239 Self: 'b,
240 {
241 self.cow_owned()
242 }
243 }
244
245 impl<'a, S1, M, SM, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition>
246 for MyCow<'a, Matrix<E, M, K, SM>>
247 where
248 S1: Storage<E, M, K>,
249 SM: Storage<E, M, K>,
250 M: Dim,
251 K: Dim,
252 E: Float,
253 DefaultAllocator: Allocator<M, K>,
254 ShapeConstraint: StridesOk<E, M, K, SM> + StridesOk<E, M, K>,
255 {
256 #[inline]
257 fn eval_ref<'b, R>(
258 &'b self,
259 f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
260 ) -> R
261 where
262 Self: 'b,
263 Matrix<E, M, K, S1>: 'b,
264 {
265 f(self.as_view::<M, K, Dyn, Dyn>())
266 }
267
268 #[inline]
269 fn own(self) -> OMatrix<E, M, K> {
270 self.into_owned()
271 }
272
273 #[inline]
274 fn cow<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
275 where
276 Self: 'b,
277 {
278 self.cow_owned()
279 }
280
281 #[inline]
282 fn decompose<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
283 where
284 Self: 'b,
285 {
286 self.cow_owned()
287 }
288 }
289
290 impl<SM, N, M, K, E> Mapping<OMatrix<E, M, K>> for Matrix<E, N, M, SM>
291 where
292 SM: Storage<E, N, M>,
293 N: Dim,
294 M: Dim,
295 K: Dim,
296 E: Float,
297 DefaultAllocator: Allocator<M, K> + Allocator<N, K>,
298 ShapeConstraint: StridesOk<E, M, K> + StridesOk<E, N, K>,
299 {
300 type Codomain = OMatrix<E, N, K>;
301
302 #[inline]
303 fn apply<I: Instance<OMatrix<E, M, K>>>(&self, x: I) -> Self::Codomain {
304 x.eval_ref(|refr| self.mul(refr))
305 }
306 }
307
308 impl<'a, SM, N, M, K, E> Linear<OMatrix<E, M, K>> for Matrix<E, N, M, SM>
309 where
310 SM: Storage<E, N, M>,
311 N: Dim,
312 M: Dim,
313 K: Dim,
314 E: Float + ClosedMulAssign + ClosedAddAssign,
315 DefaultAllocator: Allocator<N, K> + Allocator<M, K>,
316 ShapeConstraint: StridesOk<E, M, K> + StridesOk<E, N, K>,
317 {
318 }
319
320 impl<SM, SV2, N, M, K, E> GEMV<E, OMatrix<E, M, K>, Matrix<E, N, K, SV2>> for Matrix<E, N, M, SM>
321 where
322 SM: Storage<E, N, M>,
323 SV2: StorageMut<E, N, K>,
324 N: Dim,
325 M: Dim,
326 K: Dim,
327 E: Float,
328 DefaultAllocator: Allocator<N, K> + Allocator<M, K>,
329 ShapeConstraint: StridesOk<E, M, K> + StridesOk<E, N, K>,
330 {
331 #[inline]
332 fn gemv<I: Instance<OMatrix<E, M, K>>>(
333 &self, y: &mut Matrix<E, N, K, SV2>, α: E, x: I, β: E
77 ) { 334 ) {
78 x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β)) 335 x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β))
79 } 336 }
80 337
81 #[inline] 338 #[inline]
82 fn apply_mut<'a, I : Instance<Matrix<E,M,K,SV1>>>(&self, y : &mut Matrix<E,N,K,SV2>, x : I) { 339 fn apply_mut<'a, I: Instance<OMatrix<E, M, K>>>(&self, y: &mut Matrix<E, N, K, SV2>, x: I) {
83 x.eval(|x̃| self.mul_to(x̃, y)) 340 x.eval(|x̃| self.mul_to(x̃, y))
84 } 341 }
85 } 342 }
86 343
87 impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM> 344 impl<S, M, N, E> VectorSpace for Matrix<E, M, N, S>
88 where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone, 345 where
89 M : Dim, E : Scalar + Zero + One + Float, 346 S: Storage<E, M, N>,
90 DefaultAllocator : Allocator<M> { 347 M: Dim,
91 type Owned = OVector<E, M>; 348 N: Dim,
92 349 E: Float,
93 #[inline] 350 DefaultAllocator: Allocator<M, N>,
94 fn axpy<I : Instance<Vector<E,M,SV1>>>(&mut self, α : E, x : I, β : E) { 351 ShapeConstraint: StridesOk<E, M, N>,
95 x.eval(|x̃| Matrix::axpy(self, α, x̃, β)) 352 {
96 } 353 type Field = E;
97 354 type PrincipalV = OMatrix<E, M, N>;
98 #[inline] 355
99 fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) { 356 #[inline]
100 y.eval(|ỹ| Matrix::copy_from(self, ỹ)) 357 fn similar_origin(&self) -> Self::PrincipalV {
358 let (n, m) = self.shape_generic();
359 OMatrix::zeros_generic(n, m)
360 }
361 }
362
363 // This can only be implemented for the “principal” OMatrix as parameter, as otherwise
364 // we run into problems of multiple implementations when calling the methods.
365 impl<M, N, E, S> AXPY<OMatrix<E, M, N>> for Matrix<E, M, N, S>
366 where
367 S: StorageMut<E, M, N>,
368 M: Dim,
369 N: Dim,
370 E: Float,
371 DefaultAllocator: Allocator<M, N>,
372 ShapeConstraint: StridesOk<E, M, N>,
373 {
374 #[inline]
375 fn axpy<I: Instance<OMatrix<E, M, N>>>(&mut self, α: E, x: I, β: E) {
376 x.eval(|x̃| {
377 assert_eq!(self.ncols(), x̃.ncols());
378 // nalgebra does not implement axpy for matrices, and flattenining
379 // also seems difficult, so loop over columns.
380 for (mut y, ỹ) in self.column_iter_mut().zip(x̃.column_iter()) {
381 Vector::axpy(&mut y, α, &ỹ, β)
382 }
383 })
384 }
385
386 #[inline]
387 fn copy_from<I: Instance<OMatrix<E, M, N>>>(&mut self, y: I) {
388 y.eval_ref(|ỹ| Matrix::copy_from(self, &ỹ))
101 } 389 }
102 390
103 #[inline] 391 #[inline]
104 fn set_zero(&mut self) { 392 fn set_zero(&mut self) {
105 self.iter_mut().for_each(|e| *e = E::ZERO); 393 self.iter_mut().for_each(|e| *e = E::ZERO);
106 }
107
108 #[inline]
109 fn similar_origin(&self) -> Self::Owned {
110 OVector::zeros_generic(M::from_usize(self.len()), Const)
111 } 394 }
112 } 395 }
113 396
114 /* Implemented automatically as Euclidean. 397 /* Implemented automatically as Euclidean.
115 impl<SM,M,E> Projection<E, L2> for Vector<E,M,SM> 398 impl<SM,M,E> Projection<E, L2> for Vector<E,M,SM>
123 self.iter_mut().for_each(|v| *v *= ρ/n) 406 self.iter_mut().for_each(|v| *v *= ρ/n)
124 } 407 }
125 } 408 }
126 }*/ 409 }*/
127 410
128 impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM> 411 impl<SM, M, E> Projection<E, Linfinity> for Vector<E, M, SM>
129 where SM: StorageMut<E,M> + Clone, 412 where
130 M : Dim, E : Scalar + Zero + One + Float + RealField, 413 SM: StorageMut<E, M>,
131 DefaultAllocator : Allocator<M> { 414 M: Dim,
132 #[inline] 415 E: Float + RealField,
133 fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) { 416 DefaultAllocator: Allocator<M>,
134 self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) 417 ShapeConstraint: StridesOk<E, M>,
135 } 418 {
136 } 419 #[inline]
137 420 fn proj_ball(self, ρ: E, exp: Linfinity) -> <Self as Space>::Principal {
138 impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> 421 let mut owned = self.into_owned();
139 for Matrix<E,N,M,SM> 422 owned.proj_ball_mut(ρ, exp);
140 where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: Storage<E,N,K> + Clone, 423 owned
141 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + SimdComplexField, 424 }
142 DefaultAllocator : Allocator<N,K>, 425 }
143 DefaultAllocator : Allocator<M,K>, 426
144 DefaultAllocator : Allocator<N,M>, 427 impl<SM, M, E> ProjectionMut<E, Linfinity> for Vector<E, M, SM>
145 DefaultAllocator : Allocator<M,N> { 428 where
146 type AdjointCodomain = OMatrix<E,M,K>; 429 SM: StorageMut<E, M>,
147 type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a; 430 M: Dim,
431 E: Float + RealField,
432 DefaultAllocator: Allocator<M>,
433 ShapeConstraint: StridesOk<E, M>,
434 {
435 #[inline]
436 fn proj_ball_mut(&mut self, ρ: E, _: Linfinity) {
437 self.iter_mut()
438 .for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ))
439 }
440 }
441
442 impl<'own, SM, N, M, K, E> Adjointable<OMatrix<E, M, K>, OMatrix<E, N, K>> for Matrix<E, N, M, SM>
443 where
444 SM: Storage<E, N, M>,
445 N: Dim,
446 M: Dim,
447 K: Dim,
448 E: Float + RealField,
449 DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<M, N>,
450 ShapeConstraint: StridesOk<E, N, K> + StridesOk<E, M, K>,
451 {
452 type AdjointCodomain = OMatrix<E, M, K>;
453 type Adjoint<'a>
454 = OMatrix<E, M, N>
455 where
456 SM: 'a;
148 457
149 #[inline] 458 #[inline]
150 fn adjoint(&self) -> Self::Adjoint<'_> { 459 fn adjoint(&self) -> Self::Adjoint<'_> {
151 Matrix::adjoint(self) 460 Matrix::adjoint(self)
152 } 461 }
153 } 462 }
154 463
464 impl<'own, SM, N, M, K, E> SimplyAdjointable<OMatrix<E, M, K>, OMatrix<E, N, K>>
465 for Matrix<E, N, M, SM>
466 where
467 SM: Storage<E, N, M>,
468 N: Dim,
469 M: Dim,
470 K: Dim,
471 E: Float + RealField,
472 DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<M, N>,
473 ShapeConstraint: StridesOk<E, N, K> + StridesOk<E, M, K>,
474 {
475 type AdjointCodomain = OMatrix<E, M, K>;
476 type SimpleAdjoint = OMatrix<E, M, N>;
477
478 #[inline]
479 fn adjoint(&self) -> Self::SimpleAdjoint {
480 Matrix::adjoint(self)
481 }
482 }
155 /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`. 483 /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`.
156 #[inline] 484 #[inline]
157 fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>( 485 fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>(
158 /*ed: &EuclideanNorm,*/ 486 /*ed: &EuclideanNorm,*/
159 m1: &Matrix<T, R1, C1, S1>, 487 m1: &Matrix<T, R1, C1, S1>,
160 m2: &Matrix<T, R2, C2, S2>, 488 m2: &Matrix<T, R2, C2, S2>,
161 ) -> T::SimdRealField 489 ) -> T::SimdRealField
162 where 490 where
163 T: SimdComplexField, 491 T: SimdComplexField,
164 R1: Dim, 492 R1: Dim,
165 C1: Dim, 493 C1: Dim,
166 S1: Storage<T, R1, C1>, 494 S1: Storage<T, R1, C1>,
167 R2: Dim, 495 R2: Dim,
168 C2: Dim, 496 C2: Dim,
173 let diff = a - b; 501 let diff = a - b;
174 acc + diff.simd_modulus_squared() 502 acc + diff.simd_modulus_squared()
175 }) 503 })
176 } 504 }
177 505
178 // TODO: should allow different input storages in `Euclidean`. 506 impl<E, M, N, S> Euclidean<E> for Matrix<E, M, N, S>
179 507 where
180 impl<E,M,S> Euclidean<E> 508 M: Dim,
181 for Vector<E,M,S> 509 N: Dim,
182 where M : Dim, 510 S: Storage<E, M, N>,
183 S : StorageMut<E,M> + Clone, 511 E: Float + RealField,
184 E : Float + Scalar + Zero + One + RealField, 512 DefaultAllocator: Allocator<M, N>,
185 DefaultAllocator : Allocator<M> { 513 ShapeConstraint: StridesOk<E, M, N>,
186 514 {
187 type Output = OVector<E, M>; 515 type PrincipalE = OMatrix<E, M, N>;
188 516
189 #[inline] 517 #[inline]
190 fn dot<I : Instance<Self>>(&self, other : I) -> E { 518 fn dot<I: Instance<Self>>(&self, other: I) -> E {
191 Vector::<E,M,S>::dot(self, other.ref_instance()) 519 other.eval_ref(|ref r| Matrix::<E, M, N, S>::dot(self, r))
192 } 520 }
193 521
194 #[inline] 522 #[inline]
195 fn norm2_squared(&self) -> E { 523 fn norm2_squared(&self) -> E {
196 Vector::<E,M,S>::norm_squared(self) 524 Matrix::<E, M, N, S>::norm_squared(self)
197 } 525 }
198 526
199 #[inline] 527 #[inline]
200 fn dist2_squared<I : Instance<Self>>(&self, other : I) -> E { 528 fn dist2_squared<I: Instance<Self>>(&self, other: I) -> E {
201 metric_distance_squared(self, other.ref_instance()) 529 other.eval_ref(|ref r| metric_distance_squared(self, r))
202 } 530 }
203 } 531 }
204 532
205 impl<E,M,S> StaticEuclidean<E> 533 impl<E, M, S> StaticEuclidean<E> for Vector<E, M, S>
206 for Vector<E,M,S> 534 where
207 where M : DimName, 535 M: DimName,
208 S : StorageMut<E,M> + Clone, 536 S: Storage<E, M>,
209 E : Float + Scalar + Zero + One + RealField, 537 E: Float + RealField,
210 DefaultAllocator : Allocator<M> { 538 DefaultAllocator: Allocator<M>,
211 539 ShapeConstraint: StridesOk<E, M>,
540 {
212 #[inline] 541 #[inline]
213 fn origin() -> OVector<E, M> { 542 fn origin() -> OVector<E, M> {
214 OVector::zeros() 543 OVector::zeros()
215 } 544 }
216 } 545 }
217 546
218 /// The default norm for `Vector` is [`L2`]. 547 /// The default norm for `Vector` is [`L2`].
219 impl<E,M,S> Normed<E> 548 impl<E, M, N, S> Normed<E> for Matrix<E, M, N, S>
220 for Vector<E,M,S> 549 where
221 where M : Dim, 550 M: Dim,
222 S : Storage<E,M> + Clone, 551 N: Dim,
223 E : Float + Scalar + Zero + One + RealField, 552 S: Storage<E, M, N>,
224 DefaultAllocator : Allocator<M> { 553 E: Float + RealField,
225 554 DefaultAllocator: Allocator<M, N>,
555 ShapeConstraint: StridesOk<E, M, N>,
556 {
226 type NormExp = L2; 557 type NormExp = L2;
227 558
228 #[inline] 559 #[inline]
229 fn norm_exponent(&self) -> Self::NormExp { 560 fn norm_exponent(&self) -> Self::NormExp {
230 L2 561 L2
231 } 562 }
232 563
233 #[inline] 564 #[inline]
234 fn is_zero(&self) -> bool { 565 fn is_zero(&self) -> bool {
235 Vector::<E,M,S>::norm_squared(self) == E::ZERO 566 Matrix::<E, M, N, S>::norm_squared(self) == E::ZERO
236 } 567 }
237 } 568 }
238 569
239 impl<E,M,S> HasDual<E> 570 impl<E, M, N, S> HasDual<E> for Matrix<E, M, N, S>
240 for Vector<E,M,S> 571 where
241 where M : Dim, 572 M: Dim,
242 S : Storage<E,M> + Clone, 573 N: Dim,
243 E : Float + Scalar + Zero + One + RealField, 574 S: Storage<E, M, N>,
244 DefaultAllocator : Allocator<M> { 575 E: Float + RealField,
245 // TODO: Doesn't work with different storage formats. 576 DefaultAllocator: Allocator<M, N>,
246 type DualSpace = Self; 577 ShapeConstraint: StridesOk<E, M, N>,
247 } 578 {
248 579 type DualSpace = OMatrix<E, M, N>;
249 impl<E,M,S> Norm<E, L1> 580
250 for Vector<E,M,S> 581 fn dual_origin(&self) -> OMatrix<E, M, N> {
251 where M : Dim, 582 let (m, n) = self.shape_generic();
252 S : Storage<E,M>, 583 OMatrix::zeros_generic(m, n)
253 E : Float + Scalar + Zero + One + RealField, 584 }
254 DefaultAllocator : Allocator<M> { 585 }
255 586
256 #[inline] 587 impl<E, M, S> Norm<L1, E> for Vector<E, M, S>
257 fn norm(&self, _ : L1) -> E { 588 where
589 M: Dim,
590 S: Storage<E, M>,
591 E: Float + RealField,
592 {
593 #[inline]
594 fn norm(&self, _: L1) -> E {
258 nalgebra::Norm::norm(&LpNorm(1), self) 595 nalgebra::Norm::norm(&LpNorm(1), self)
259 } 596 }
260 } 597 }
261 598
262 impl<E,M,S> Dist<E, L1> 599 impl<E, M, S> Dist<L1, E> for Vector<E, M, S>
263 for Vector<E,M,S> 600 where
264 where M : Dim, 601 M: Dim,
265 S : Storage<E,M> + Clone, 602 S: Storage<E, M>,
266 E : Float + Scalar + Zero + One + RealField, 603 E: Float + RealField,
267 DefaultAllocator : Allocator<M> { 604 DefaultAllocator: Allocator<M>,
268 #[inline] 605 ShapeConstraint: StridesOk<E, M>,
269 fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> E { 606 {
270 nalgebra::Norm::metric_distance(&LpNorm(1), self, other.ref_instance()) 607 #[inline]
271 } 608 fn dist<I: Instance<Self>>(&self, other: I, _: L1) -> E {
272 } 609 other.eval_ref(|ref r| nalgebra::Norm::metric_distance(&LpNorm(1), self, r))
273 610 }
274 impl<E,M,S> Norm<E, L2> 611 }
275 for Vector<E,M,S> 612
276 where M : Dim, 613 impl<E, M, N, S> Norm<L2, E> for Matrix<E, M, N, S>
277 S : Storage<E,M>, 614 where
278 E : Float + Scalar + Zero + One + RealField, 615 M: Dim,
279 DefaultAllocator : Allocator<M> { 616 N: Dim,
280 617 S: Storage<E, M, N>,
281 #[inline] 618 E: Float + RealField,
282 fn norm(&self, _ : L2) -> E { 619 {
620 #[inline]
621 fn norm(&self, _: L2) -> E {
283 nalgebra::Norm::norm(&LpNorm(2), self) 622 nalgebra::Norm::norm(&LpNorm(2), self)
284 } 623 }
285 } 624 }
286 625
287 impl<E,M,S> Dist<E, L2> 626 impl<E, M, N, S> Dist<L2, E> for Matrix<E, M, N, S>
288 for Vector<E,M,S> 627 where
289 where M : Dim, 628 M: Dim,
290 S : Storage<E,M> + Clone, 629 N: Dim,
291 E : Float + Scalar + Zero + One + RealField, 630 S: Storage<E, M, N>,
292 DefaultAllocator : Allocator<M> { 631 E: Float + RealField,
293 #[inline] 632 DefaultAllocator: Allocator<M, N>,
294 fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> E { 633 ShapeConstraint: StridesOk<E, M, N>,
295 nalgebra::Norm::metric_distance(&LpNorm(2), self, other.ref_instance()) 634 {
296 } 635 #[inline]
297 } 636 fn dist<I: Instance<Self>>(&self, other: I, _: L2) -> E {
298 637 other.eval_ref(|ref r| nalgebra::Norm::metric_distance(&LpNorm(2), self, r))
299 impl<E,M,S> Norm<E, Linfinity> 638 }
300 for Vector<E,M,S> 639 }
301 where M : Dim, 640
302 S : Storage<E,M>, 641 impl<E, M, N, S> Norm<Linfinity, E> for Matrix<E, M, N, S>
303 E : Float + Scalar + Zero + One + RealField, 642 where
304 DefaultAllocator : Allocator<M> { 643 M: Dim,
305 644 N: Dim,
306 #[inline] 645 S: Storage<E, M, N>,
307 fn norm(&self, _ : Linfinity) -> E { 646 E: Float + RealField,
647 {
648 #[inline]
649 fn norm(&self, _: Linfinity) -> E {
308 nalgebra::Norm::norm(&UniformNorm, self) 650 nalgebra::Norm::norm(&UniformNorm, self)
309 } 651 }
310 } 652 }
311 653
312 impl<E,M,S> Dist<E, Linfinity> 654 impl<E, M, N, S> Dist<Linfinity, E> for Matrix<E, M, N, S>
313 for Vector<E,M,S> 655 where
314 where M : Dim, 656 M: Dim,
315 S : Storage<E,M> + Clone, 657 N: Dim,
316 E : Float + Scalar + Zero + One + RealField, 658 S: Storage<E, M, N>,
317 DefaultAllocator : Allocator<M> { 659 E: Float + RealField,
318 #[inline] 660 DefaultAllocator: Allocator<M, N>,
319 fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> E { 661 ShapeConstraint: StridesOk<E, M, N>,
320 nalgebra::Norm::metric_distance(&UniformNorm, self, other.ref_instance()) 662 {
663 #[inline]
664 fn dist<I: Instance<Self>>(&self, other: I, _: Linfinity) -> E {
665 other.eval_ref(|ref r| nalgebra::Norm::metric_distance(&UniformNorm, self, r))
321 } 666 }
322 } 667 }
323 668
324 /// Helper trait to hide the symbols of [`nalgebra::RealField`]. 669 /// Helper trait to hide the symbols of [`nalgebra::RealField`].
325 /// 670 ///
327 /// functions can piggyback `nalgebra::RealField` without exponsing themselves to it. 672 /// functions can piggyback `nalgebra::RealField` without exponsing themselves to it.
328 /// Thus methods from [`num_traits`] can be used directly without similarly named methods 673 /// Thus methods from [`num_traits`] can be used directly without similarly named methods
329 /// from [`nalgebra`] conflicting with them. Only when absolutely necessary to work with 674 /// from [`nalgebra`] conflicting with them. Only when absolutely necessary to work with
330 /// nalgebra, one can convert to the nalgebra view of the same type using the methods of 675 /// nalgebra, one can convert to the nalgebra view of the same type using the methods of
331 /// this trait. 676 /// this trait.
332 pub trait ToNalgebraRealField : Float { 677 pub trait ToNalgebraRealField: Float {
333 /// The nalgebra type corresponding to this type. Usually same as `Self`. 678 /// The nalgebra type corresponding to this type. Usually same as `Self`.
334 /// 679 ///
335 /// This type only carries `nalgebra` traits. 680 /// This type only carries `nalgebra` traits.
336 type NalgebraType : RealField; 681 type NalgebraType: RealField;
337 /// The “mixed” type corresponding to this type. Usually same as `Self`. 682 /// The “mixed” type corresponding to this type. Usually same as `Self`.
338 /// 683 ///
339 /// This type carries both `num_traits` and `nalgebra` traits. 684 /// This type carries both `num_traits` and `nalgebra` traits.
340 type MixedType : RealField + Float; 685 type MixedType: RealField + Float;
341 686
342 /// Convert to the nalgebra view of `self`. 687 /// Convert to the nalgebra view of `self`.
343 fn to_nalgebra(self) -> Self::NalgebraType; 688 fn to_nalgebra(self) -> Self::NalgebraType;
344 689
345 /// Convert to the mixed (nalgebra and num_traits) view of `self`. 690 /// Convert to the mixed (nalgebra and num_traits) view of `self`.
346 fn to_nalgebra_mixed(self) -> Self::MixedType; 691 fn to_nalgebra_mixed(self) -> Self::MixedType;
347 692
348 /// Convert from the nalgebra view of `self`. 693 /// Convert from the nalgebra view of `self`.
349 fn from_nalgebra(t : Self::NalgebraType) -> Self; 694 fn from_nalgebra(t: Self::NalgebraType) -> Self;
350 695
351 /// Convert from the mixed (nalgebra and num_traits) view to `self`. 696 /// Convert from the mixed (nalgebra and num_traits) view to `self`.
352 fn from_nalgebra_mixed(t : Self::MixedType) -> Self; 697 fn from_nalgebra_mixed(t: Self::MixedType) -> Self;
353 } 698 }
354 699
355 impl ToNalgebraRealField for f32 { 700 impl ToNalgebraRealField for f32 {
356 type NalgebraType = f32; 701 type NalgebraType = f32;
357 type MixedType = f32; 702 type MixedType = f32;
358 703
359 #[inline] 704 #[inline]
360 fn to_nalgebra(self) -> Self::NalgebraType { self } 705 fn to_nalgebra(self) -> Self::NalgebraType {
361 706 self
362 #[inline] 707 }
363 fn to_nalgebra_mixed(self) -> Self::MixedType { self } 708
364 709 #[inline]
365 #[inline] 710 fn to_nalgebra_mixed(self) -> Self::MixedType {
366 fn from_nalgebra(t : Self::NalgebraType) -> Self { t } 711 self
367 712 }
368 #[inline] 713
369 fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } 714 #[inline]
370 715 fn from_nalgebra(t: Self::NalgebraType) -> Self {
716 t
717 }
718
719 #[inline]
720 fn from_nalgebra_mixed(t: Self::MixedType) -> Self {
721 t
722 }
371 } 723 }
372 724
373 impl ToNalgebraRealField for f64 { 725 impl ToNalgebraRealField for f64 {
374 type NalgebraType = f64; 726 type NalgebraType = f64;
375 type MixedType = f64; 727 type MixedType = f64;
376 728
377 #[inline] 729 #[inline]
378 fn to_nalgebra(self) -> Self::NalgebraType { self } 730 fn to_nalgebra(self) -> Self::NalgebraType {
379 731 self
380 #[inline] 732 }
381 fn to_nalgebra_mixed(self) -> Self::MixedType { self } 733
382 734 #[inline]
383 #[inline] 735 fn to_nalgebra_mixed(self) -> Self::MixedType {
384 fn from_nalgebra(t : Self::NalgebraType) -> Self { t } 736 self
385 737 }
386 #[inline] 738
387 fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } 739 #[inline]
388 } 740 fn from_nalgebra(t: Self::NalgebraType) -> Self {
389 741 t
742 }
743
744 #[inline]
745 fn from_nalgebra_mixed(t: Self::MixedType) -> Self {
746 t
747 }
748 }

Mercurial Repository: alg_tools

mercurial