Mercurial > repos > alg_tools / file comparison

comparison: src/nalgebra_support.rs

src/nalgebra_support.rs

changeset 90
b3c35d16affe
parent 70
672aec2e1acd
child 82
981069ef919b
equal deleted inserted replaced
25:d14c877e14b7 90:b3c35d16affe
1 /*! 1 /*!
2 Integration with nalgebra. 2 Integration with nalgebra.
3 3
4 This module mainly implements [`Euclidean`], [`Norm`], [`Dot`], [`Linear`], etc. for [`nalgebra`] 4 This module mainly implements [`Euclidean`], [`Norm`], [`Linear`], etc. for [`nalgebra`]
5 matrices and vectors. 5 matrices and vectors.
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 nalgebra::{ 11 use nalgebra::{
12 Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar, 12 Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar,
13 ClosedMul, ClosedAdd, SimdComplexField, Vector, OVector, RealField, 13 ClosedAddAssign, ClosedMulAssign, SimdComplexField, Vector, OVector, RealField,
14 LpNorm, UniformNorm 14 LpNorm, UniformNorm
15 }; 15 };
16 use nalgebra::Norm as NalgebraNorm;
17 use nalgebra::base::constraint::{ 16 use nalgebra::base::constraint::{
18 ShapeConstraint, SameNumberOfRows, SameNumberOfColumns 17 ShapeConstraint, SameNumberOfRows, SameNumberOfColumns
19 }; 18 };
20 use nalgebra::base::dimension::*; 19 use nalgebra::base::dimension::*;
21 use nalgebra::base::allocator::Allocator; 20 use nalgebra::base::allocator::Allocator;
22 use std::ops::Mul; 21 use std::ops::Mul;
23 use num_traits::identities::{Zero, One}; 22 use num_traits::identities::{Zero, One};
24 use crate::linops::*; 23 use crate::linops::*;
25 use crate::euclidean::*; 24 use crate::euclidean::*;
25 use crate::mapping::{Space, BasicDecomposition};
26 use crate::types::Float; 26 use crate::types::Float;
27 use crate::norms::*; 27 use crate::norms::*;
28 28 use crate::instance::Instance;
29 impl<SM,SV,N,M,K,E> Apply<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> 29
30 where SM: Storage<E,N,M>, SV: Storage<E,M,K>, 30 impl<SM,N,M,E> Space for Matrix<E,N,M,SM>
31 N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, 31 where
32 DefaultAllocator : Allocator<E,N,K>, 32 SM: Storage<E,N,M> + Clone,
33 DefaultAllocator : Allocator<E,M,K>, 33 N : Dim, M : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
34 DefaultAllocator : Allocator<E,N,M>, 34 DefaultAllocator : Allocator<N,M>,
35 DefaultAllocator : Allocator<E,M,N> { 35 {
36 type Output = OMatrix<E,N,K>; 36 type Decomp = BasicDecomposition;
37 37 }
38 #[inline] 38
39 fn apply(&self, x : Matrix<E,M,K,SV>) -> Self::Output { 39 impl<SM,SV,N,M,K,E> Mapping<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
40 self.mul(x) 40 where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone,
41 } 41 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
42 } 42 DefaultAllocator : Allocator<N,K>,
43 43 DefaultAllocator : Allocator<M,K>,
44 impl<'a, SM,SV,N,M,K,E> Apply<&'a Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> 44 DefaultAllocator : Allocator<N,M>,
45 where SM: Storage<E,N,M>, SV: Storage<E,M,K>, 45 DefaultAllocator : Allocator<M,N> {
46 N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, 46 type Codomain = OMatrix<E,N,K>;
47 DefaultAllocator : Allocator<E,N,K>, 47
48 DefaultAllocator : Allocator<E,M,K>, 48 #[inline]
49 DefaultAllocator : Allocator<E,N,M>, 49 fn apply<I : Instance<Matrix<E,M,K,SV>>>(
50 DefaultAllocator : Allocator<E,M,N> { 50 &self, x : I
51 type Output = OMatrix<E,N,K>; 51 ) -> Self::Codomain {
52 52 x.either(|owned| self.mul(owned), |refr| self.mul(refr))
53 #[inline] 53 }
54 fn apply(&self, x : &'a Matrix<E,M,K,SV>) -> Self::Output { 54 }
55 self.mul(x) 55
56 }
57 }
58 56
59 impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> 57 impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
60 where SM: Storage<E,N,M>, SV: Storage<E,M,K>, 58 where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone,
61 N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, 59 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
62 DefaultAllocator : Allocator<E,N,K>, 60 DefaultAllocator : Allocator<N,K>,
63 DefaultAllocator : Allocator<E,M,K>, 61 DefaultAllocator : Allocator<M,K>,
64 DefaultAllocator : Allocator<E,N,M>, 62 DefaultAllocator : Allocator<N,M>,
65 DefaultAllocator : Allocator<E,M,N> { 63 DefaultAllocator : Allocator<M,N> {
66 type Codomain = OMatrix<E,N,K>;
67 } 64 }
68 65
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> 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 where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: StorageMut<E,N,K>, 67 where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: StorageMut<E,N,K>,
71 N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, 68 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float,
72 DefaultAllocator : Allocator<E,N,K>, 69 DefaultAllocator : Allocator<N,K>,
73 DefaultAllocator : Allocator<E,M,K>, 70 DefaultAllocator : Allocator<M,K>,
74 DefaultAllocator : Allocator<E,N,M>, 71 DefaultAllocator : Allocator<N,M>,
75 DefaultAllocator : Allocator<E,M,N> { 72 DefaultAllocator : Allocator<M,N> {
76 73
77 #[inline] 74 #[inline]
78 fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : &Matrix<E,M,K,SV1>, β : E) { 75 fn gemv<I : Instance<Matrix<E,M,K,SV1>>>(
79 Matrix::gemm(y, α, self, x, β) 76 &self, y : &mut Matrix<E,N,K,SV2>, α : E, x : I, β : E
80 } 77 ) {
81 78 x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β))
82 #[inline] 79 }
83 fn apply_mut<'a>(&self, y : &mut Matrix<E,N,K,SV2>, x : &Matrix<E,M,K,SV1>) { 80
84 self.mul_to(x, y) 81 #[inline]
82 fn apply_mut<'a, I : Instance<Matrix<E,M,K,SV1>>>(&self, y : &mut Matrix<E,N,K,SV2>, x : I) {
83 x.eval(|x̃| self.mul_to(x̃, y))
85 } 84 }
86 } 85 }
87 86
88 impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM> 87 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>, 88 where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone,
90 M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, 89 M : Dim, E : Scalar + Zero + One + Float,
91 DefaultAllocator : Allocator<E,M> { 90 DefaultAllocator : Allocator<M> {
92 91 type Owned = OVector<E, M>;
93 #[inline] 92
94 fn axpy(&mut self, α : E, x : &Vector<E,M,SV1>, β : E) { 93 #[inline]
95 Matrix::axpy(self, α, x, β) 94 fn axpy<I : Instance<Vector<E,M,SV1>>>(&mut self, α : E, x : I, β : E) {
96 } 95 x.eval(|x̃| Matrix::axpy(self, α, x̃, β))
97 96 }
98 #[inline] 97
99 fn copy_from(&mut self, y : &Vector<E,M,SV1>) { 98 #[inline]
100 Matrix::copy_from(self, y) 99 fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) {
101 } 100 y.eval(|ỹ| Matrix::copy_from(self, ỹ))
102 } 101 }
102
103 #[inline]
104 fn set_zero(&mut self) {
105 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 }
112 }
113
114 /* Implemented automatically as Euclidean.
115 impl<SM,M,E> Projection<E, L2> for Vector<E,M,SM>
116 where SM: StorageMut<E,M> + Clone,
117 M : Dim, E : Scalar + Zero + One + Float + RealField,
118 DefaultAllocator : Allocator<M> {
119 #[inline]
120 fn proj_ball_mut(&mut self, ρ : E, _ : L2) {
121 let n = self.norm(L2);
122 if n > ρ {
123 self.iter_mut().for_each(|v| *v *= ρ/n)
124 }
125 }
126 }*/
103 127
104 impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM> 128 impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM>
105 where SM: StorageMut<E,M>, 129 where SM: StorageMut<E,M> + Clone,
106 M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float + RealField, 130 M : Dim, E : Scalar + Zero + One + Float + RealField,
107 DefaultAllocator : Allocator<E,M> { 131 DefaultAllocator : Allocator<M> {
108 #[inline] 132 #[inline]
109 fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) { 133 fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) {
110 self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) 134 self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ))
111 } 135 }
112 } 136 }
113 137
114 impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>,Matrix<E,N,K,SV2>> 138 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> 139 for Matrix<E,N,M,SM>
116 where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: Storage<E,N,K>, 140 where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: Storage<E,N,K> + Clone,
117 N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + SimdComplexField, 141 N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + SimdComplexField,
118 DefaultAllocator : Allocator<E,N,K>, 142 DefaultAllocator : Allocator<N,K>,
119 DefaultAllocator : Allocator<E,M,K>, 143 DefaultAllocator : Allocator<M,K>,
120 DefaultAllocator : Allocator<E,N,M>, 144 DefaultAllocator : Allocator<N,M>,
121 DefaultAllocator : Allocator<E,M,N> { 145 DefaultAllocator : Allocator<M,N> {
122 type AdjointCodomain = OMatrix<E,M,K>; 146 type AdjointCodomain = OMatrix<E,M,K>;
123 type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a; 147 type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a;
124 148
125 #[inline] 149 #[inline]
126 fn adjoint(&self) -> Self::Adjoint<'_> { 150 fn adjoint(&self) -> Self::Adjoint<'_> {
127 Matrix::adjoint(self) 151 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 } 152 }
143 } 153 }
144 154
145 /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`. 155 /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`.
146 #[inline] 156 #[inline]
168 // TODO: should allow different input storages in `Euclidean`. 178 // TODO: should allow different input storages in `Euclidean`.
169 179
170 impl<E,M,S> Euclidean<E> 180 impl<E,M,S> Euclidean<E>
171 for Vector<E,M,S> 181 for Vector<E,M,S>
172 where M : Dim, 182 where M : Dim,
173 S : StorageMut<E,M>, 183 S : StorageMut<E,M> + Clone,
174 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 184 E : Float + Scalar + Zero + One + RealField,
175 DefaultAllocator : Allocator<E,M> { 185 DefaultAllocator : Allocator<M> {
176 186
177 type Output = OVector<E, M>; 187 type Output = OVector<E, M>;
178 188
179 #[inline] 189 #[inline]
180 fn similar_origin(&self) -> OVector<E, M> { 190 fn dot<I : Instance<Self>>(&self, other : I) -> E {
181 OVector::zeros_generic(M::from_usize(self.len()), Const) 191 Vector::<E,M,S>::dot(self, other.ref_instance())
182 } 192 }
183 193
184 #[inline] 194 #[inline]
185 fn norm2_squared(&self) -> E { 195 fn norm2_squared(&self) -> E {
186 Vector::<E,M,S>::norm_squared(self) 196 Vector::<E,M,S>::norm_squared(self)
187 } 197 }
188 198
189 #[inline] 199 #[inline]
190 fn dist2_squared(&self, other : &Self) -> E { 200 fn dist2_squared<I : Instance<Self>>(&self, other : I) -> E {
191 metric_distance_squared(self, other) 201 metric_distance_squared(self, other.ref_instance())
192 } 202 }
193 } 203 }
194 204
195 impl<E,M,S> StaticEuclidean<E> 205 impl<E,M,S> StaticEuclidean<E>
196 for Vector<E,M,S> 206 for Vector<E,M,S>
197 where M : DimName, 207 where M : DimName,
198 S : StorageMut<E,M>, 208 S : StorageMut<E,M> + Clone,
199 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 209 E : Float + Scalar + Zero + One + RealField,
200 DefaultAllocator : Allocator<E,M> { 210 DefaultAllocator : Allocator<M> {
201 211
202 #[inline] 212 #[inline]
203 fn origin() -> OVector<E, M> { 213 fn origin() -> OVector<E, M> {
204 OVector::zeros() 214 OVector::zeros()
205 } 215 }
206 } 216 }
207 217
218 /// The default norm for `Vector` is [`L2`].
219 impl<E,M,S> Normed<E>
220 for Vector<E,M,S>
221 where M : Dim,
222 S : Storage<E,M> + Clone,
223 E : Float + Scalar + Zero + One + RealField,
224 DefaultAllocator : Allocator<M> {
225
226 type NormExp = L2;
227
228 #[inline]
229 fn norm_exponent(&self) -> Self::NormExp {
230 L2
231 }
232
233 #[inline]
234 fn is_zero(&self) -> bool {
235 Vector::<E,M,S>::norm_squared(self) == E::ZERO
236 }
237 }
238
239 impl<E,M,S> HasDual<E>
240 for Vector<E,M,S>
241 where M : Dim,
242 S : Storage<E,M> + Clone,
243 E : Float + Scalar + Zero + One + RealField,
244 DefaultAllocator : Allocator<M> {
245 // TODO: Doesn't work with different storage formats.
246 type DualSpace = Self;
247 }
248
208 impl<E,M,S> Norm<E, L1> 249 impl<E,M,S> Norm<E, L1>
209 for Vector<E,M,S> 250 for Vector<E,M,S>
210 where M : Dim, 251 where M : Dim,
211 S : StorageMut<E,M>, 252 S : Storage<E,M>,
212 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 253 E : Float + Scalar + Zero + One + RealField,
213 DefaultAllocator : Allocator<E,M> { 254 DefaultAllocator : Allocator<M> {
214 255
215 #[inline] 256 #[inline]
216 fn norm(&self, _ : L1) -> E { 257 fn norm(&self, _ : L1) -> E {
217 LpNorm(1).norm(self) 258 nalgebra::Norm::norm(&LpNorm(1), self)
218 } 259 }
219 } 260 }
220 261
221 impl<E,M,S> Dist<E, L1> 262 impl<E,M,S> Dist<E, L1>
222 for Vector<E,M,S> 263 for Vector<E,M,S>
223 where M : Dim, 264 where M : Dim,
224 S : StorageMut<E,M>, 265 S : Storage<E,M> + Clone,
225 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 266 E : Float + Scalar + Zero + One + RealField,
226 DefaultAllocator : Allocator<E,M> { 267 DefaultAllocator : Allocator<M> {
227 #[inline] 268 #[inline]
228 fn dist(&self, other : &Self, _ : L1) -> E { 269 fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> E {
229 LpNorm(1).metric_distance(self, other) 270 nalgebra::Norm::metric_distance(&LpNorm(1), self, other.ref_instance())
230 } 271 }
231 } 272 }
232 273
233 impl<E,M,S> Norm<E, L2> 274 impl<E,M,S> Norm<E, L2>
234 for Vector<E,M,S> 275 for Vector<E,M,S>
235 where M : Dim, 276 where M : Dim,
236 S : StorageMut<E,M>, 277 S : Storage<E,M>,
237 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 278 E : Float + Scalar + Zero + One + RealField,
238 DefaultAllocator : Allocator<E,M> { 279 DefaultAllocator : Allocator<M> {
239 280
240 #[inline] 281 #[inline]
241 fn norm(&self, _ : L2) -> E { 282 fn norm(&self, _ : L2) -> E {
242 LpNorm(2).norm(self) 283 nalgebra::Norm::norm(&LpNorm(2), self)
243 } 284 }
244 } 285 }
245 286
246 impl<E,M,S> Dist<E, L2> 287 impl<E,M,S> Dist<E, L2>
247 for Vector<E,M,S> 288 for Vector<E,M,S>
248 where M : Dim, 289 where M : Dim,
249 S : StorageMut<E,M>, 290 S : Storage<E,M> + Clone,
250 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 291 E : Float + Scalar + Zero + One + RealField,
251 DefaultAllocator : Allocator<E,M> { 292 DefaultAllocator : Allocator<M> {
252 #[inline] 293 #[inline]
253 fn dist(&self, other : &Self, _ : L2) -> E { 294 fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> E {
254 LpNorm(2).metric_distance(self, other) 295 nalgebra::Norm::metric_distance(&LpNorm(2), self, other.ref_instance())
255 } 296 }
256 } 297 }
257 298
258 impl<E,M,S> Norm<E, Linfinity> 299 impl<E,M,S> Norm<E, Linfinity>
259 for Vector<E,M,S> 300 for Vector<E,M,S>
260 where M : Dim, 301 where M : Dim,
261 S : StorageMut<E,M>, 302 S : Storage<E,M>,
262 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 303 E : Float + Scalar + Zero + One + RealField,
263 DefaultAllocator : Allocator<E,M> { 304 DefaultAllocator : Allocator<M> {
264 305
265 #[inline] 306 #[inline]
266 fn norm(&self, _ : Linfinity) -> E { 307 fn norm(&self, _ : Linfinity) -> E {
267 UniformNorm.norm(self) 308 nalgebra::Norm::norm(&UniformNorm, self)
268 } 309 }
269 } 310 }
270 311
271 impl<E,M,S> Dist<E, Linfinity> 312 impl<E,M,S> Dist<E, Linfinity>
272 for Vector<E,M,S> 313 for Vector<E,M,S>
273 where M : Dim, 314 where M : Dim,
274 S : StorageMut<E,M>, 315 S : Storage<E,M> + Clone,
275 E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, 316 E : Float + Scalar + Zero + One + RealField,
276 DefaultAllocator : Allocator<E,M> { 317 DefaultAllocator : Allocator<M> {
277 #[inline] 318 #[inline]
278 fn dist(&self, other : &Self, _ : Linfinity) -> E { 319 fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> E {
279 UniformNorm.metric_distance(self, other) 320 nalgebra::Norm::metric_distance(&UniformNorm, self, other.ref_instance())
280 } 321 }
281 } 322 }
282 323
283 /// Helper trait to hide the symbols of [`nalgebra::RealField`]. 324 /// Helper trait to hide the symbols of [`nalgebra::RealField`].
284 /// 325 ///

Mercurial Repository: alg_tools

mercurial