Fri, 18 Nov 2022 10:29:50 +0200
Improvements and minor fixes to bisection tree refinement.
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 | ||
29 | impl<SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> | |
30 | where SM: Storage<E,N,M>, SV: Storage<E,M,K>, | |
31 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One, | |
32 | DefaultAllocator : Allocator<E,N,K>, | |
33 | DefaultAllocator : Allocator<E,M,K>, | |
34 | DefaultAllocator : Allocator<E,N,M>, | |
35 | DefaultAllocator : Allocator<E,M,N> { | |
36 | type Codomain = OMatrix<E,N,K>; | |
37 | ||
38 | #[inline] | |
39 | fn apply(&self, x : &Matrix<E,M,K,SV>) -> Self::Codomain { | |
40 | self.mul(x) | |
41 | } | |
42 | } | |
43 | ||
44 | 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> | |
45 | where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: StorageMut<E,N,K>, | |
46 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, | |
47 | DefaultAllocator : Allocator<E,N,K>, | |
48 | DefaultAllocator : Allocator<E,M,K>, | |
49 | DefaultAllocator : Allocator<E,N,M>, | |
50 | DefaultAllocator : Allocator<E,M,N> { | |
51 | ||
52 | #[inline] | |
53 | fn gemv(&self, y : &mut Matrix<E,N,K,SV2>, α : E, x : &Matrix<E,M,K,SV1>, β : E) { | |
54 | Matrix::gemm(y, α, self, x, β) | |
55 | } | |
56 | ||
57 | #[inline] | |
58 | fn apply_mut<'a>(&self, y : &mut Matrix<E,N,K,SV2>, x : &Matrix<E,M,K,SV1>) { | |
59 | self.mul_to(x, y) | |
60 | } | |
61 | } | |
62 | ||
63 | impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM> | |
64 | where SM: StorageMut<E,M>, SV1: Storage<E,M>, | |
65 | M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float, | |
66 | DefaultAllocator : Allocator<E,M> { | |
67 | ||
68 | #[inline] | |
69 | fn axpy(&mut self, α : E, x : &Vector<E,M,SV1>, β : E) { | |
70 | Matrix::axpy(self, α, x, β) | |
71 | } | |
72 | ||
73 | #[inline] | |
74 | fn copy_from(&mut self, y : &Vector<E,M,SV1>) { | |
75 | Matrix::copy_from(self, y) | |
76 | } | |
77 | } | |
78 | ||
79 | impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM> | |
80 | where SM: StorageMut<E,M>, | |
81 | M : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + Float + RealField, | |
82 | DefaultAllocator : Allocator<E,M> { | |
83 | #[inline] | |
84 | fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) { | |
85 | self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) | |
86 | } | |
87 | } | |
88 | ||
89 | impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>,Matrix<E,N,K,SV2>> | |
90 | for Matrix<E,N,M,SM> | |
91 | where SM: Storage<E,N,M>, SV1: Storage<E,M,K>, SV2: Storage<E,N,K>, | |
92 | N : Dim, M : Dim, K : Dim, E : Scalar + ClosedMul + ClosedAdd + Zero + One + SimdComplexField, | |
93 | DefaultAllocator : Allocator<E,N,K>, | |
94 | DefaultAllocator : Allocator<E,M,K>, | |
95 | DefaultAllocator : Allocator<E,N,M>, | |
96 | DefaultAllocator : Allocator<E,M,N> { | |
97 | type AdjointCodomain = OMatrix<E,M,K>; | |
98 | type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a; | |
99 | ||
100 | #[inline] | |
101 | fn adjoint(&self) -> Self::Adjoint<'_> { | |
102 | Matrix::adjoint(self) | |
103 | } | |
104 | } | |
105 | ||
106 | impl<E,M,S,Si> Dot<Vector<E,M,Si>,E> | |
107 | for Vector<E,M,S> | |
108 | where M : Dim, | |
109 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One, | |
110 | S : Storage<E,M>, | |
111 | Si : Storage<E,M>, | |
112 | DefaultAllocator : Allocator<E,M> { | |
113 | ||
114 | #[inline] | |
115 | fn dot(&self, other : &Vector<E,M,Si>) -> E { | |
116 | Vector::<E,M,S>::dot(self, other) | |
117 | } | |
118 | } | |
119 | ||
120 | /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`. | |
121 | #[inline] | |
122 | fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>( | |
123 | /*ed: &EuclideanNorm,*/ | |
124 | m1: &Matrix<T, R1, C1, S1>, | |
125 | m2: &Matrix<T, R2, C2, S2>, | |
126 | ) -> T::SimdRealField | |
127 | where | |
128 | T: SimdComplexField, | |
129 | R1: Dim, | |
130 | C1: Dim, | |
131 | S1: Storage<T, R1, C1>, | |
132 | R2: Dim, | |
133 | C2: Dim, | |
134 | S2: Storage<T, R2, C2>, | |
135 | ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2>, | |
136 | { | |
137 | m1.zip_fold(m2, T::SimdRealField::zero(), |acc, a, b| { | |
138 | let diff = a - b; | |
139 | acc + diff.simd_modulus_squared() | |
140 | }) | |
141 | } | |
142 | ||
143 | // TODO: should allow different input storages in `Euclidean`. | |
144 | ||
145 | impl<E,M,S> Euclidean<E> | |
146 | for Vector<E,M,S> | |
147 | where M : Dim, | |
148 | S : StorageMut<E,M>, | |
149 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
150 | DefaultAllocator : Allocator<E,M> { | |
151 | ||
152 | type Output = OVector<E, M>; | |
153 | ||
154 | #[inline] | |
155 | fn similar_origin(&self) -> OVector<E, M> { | |
156 | OVector::zeros_generic(M::from_usize(self.len()), Const) | |
157 | } | |
158 | ||
159 | #[inline] | |
160 | fn norm2_squared(&self) -> E { | |
161 | Vector::<E,M,S>::norm_squared(self) | |
162 | } | |
163 | ||
164 | #[inline] | |
165 | fn dist2_squared(&self, other : &Self) -> E { | |
166 | metric_distance_squared(self, other) | |
167 | } | |
168 | } | |
169 | ||
170 | impl<E,M,S> StaticEuclidean<E> | |
171 | for Vector<E,M,S> | |
172 | where M : DimName, | |
173 | S : StorageMut<E,M>, | |
174 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
175 | DefaultAllocator : Allocator<E,M> { | |
176 | ||
177 | #[inline] | |
178 | fn origin() -> OVector<E, M> { | |
179 | OVector::zeros() | |
180 | } | |
181 | } | |
182 | ||
183 | impl<E,M,S> Norm<E, L1> | |
184 | for Vector<E,M,S> | |
185 | where M : Dim, | |
186 | S : StorageMut<E,M>, | |
187 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
188 | DefaultAllocator : Allocator<E,M> { | |
189 | ||
190 | #[inline] | |
191 | fn norm(&self, _ : L1) -> E { | |
192 | LpNorm(1).norm(self) | |
193 | } | |
194 | } | |
195 | ||
196 | impl<E,M,S> Dist<E, L1> | |
197 | for Vector<E,M,S> | |
198 | where M : Dim, | |
199 | S : StorageMut<E,M>, | |
200 | E : Float + Scalar + ClosedMul + ClosedAdd + Zero + One + RealField, | |
201 | DefaultAllocator : Allocator<E,M> { | |
202 | #[inline] | |
203 | fn dist(&self, other : &Self, _ : L1) -> E { | |
204 | LpNorm(1).metric_distance(self, other) | |
205 | } | |
206 | } | |
207 | ||
208 | impl<E,M,S> Norm<E, L2> | |
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, _ : L2) -> E { | |
217 | LpNorm(2).norm(self) | |
218 | } | |
219 | } | |
220 | ||
221 | impl<E,M,S> Dist<E, L2> | |
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, _ : L2) -> E { | |
229 | LpNorm(2).metric_distance(self, other) | |
230 | } | |
231 | } | |
232 | ||
233 | impl<E,M,S> Norm<E, Linfinity> | |
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, _ : Linfinity) -> E { | |
242 | UniformNorm.norm(self) | |
243 | } | |
244 | } | |
245 | ||
246 | impl<E,M,S> Dist<E, Linfinity> | |
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, _ : Linfinity) -> E { | |
254 | UniformNorm.metric_distance(self, other) | |
255 | } | |
256 | } | |
257 | ||
5 | 258 | /// Helper trait to hide the symbols of [`nalgebra::RealField`]. |
259 | /// | |
260 | /// By assuming `ToNalgebraRealField` intead of `nalgebra::RealField` as a trait bound, | |
261 | /// functions can piggyback `nalgebra::RealField` without exponsing themselves to it. | |
262 | /// Thus methods from [`num_traits`] can be used directly without similarly named methods | |
263 | /// from [`nalgebra`] conflicting with them. Only when absolutely necessary to work with | |
264 | /// nalgebra, one can convert to the nalgebra view of the same type using the methods of | |
265 | /// this trait. | |
0 | 266 | pub trait ToNalgebraRealField : Float { |
5 | 267 | /// The nalgebra type corresponding to this type. Usually same as `Self`. |
268 | /// | |
269 | /// This type only carries `nalgebra` traits. | |
0 | 270 | type NalgebraType : RealField; |
5 | 271 | /// The “mixed” type corresponding to this type. Usually same as `Self`. |
272 | /// | |
273 | /// This type carries both `num_traits` and `nalgebra` traits. | |
0 | 274 | type MixedType : RealField + Float; |
275 | ||
5 | 276 | /// Convert to the nalgebra view of `self`. |
0 | 277 | fn to_nalgebra(self) -> Self::NalgebraType; |
5 | 278 | |
279 | /// Convert to the mixed (nalgebra and num_traits) view of `self`. | |
0 | 280 | fn to_nalgebra_mixed(self) -> Self::MixedType; |
281 | ||
5 | 282 | /// Convert from the nalgebra view of `self`. |
0 | 283 | fn from_nalgebra(t : Self::NalgebraType) -> Self; |
5 | 284 | |
285 | /// Convert from the mixed (nalgebra and num_traits) view to `self`. | |
0 | 286 | fn from_nalgebra_mixed(t : Self::MixedType) -> Self; |
287 | } | |
288 | ||
289 | impl ToNalgebraRealField for f32 { | |
290 | type NalgebraType = f32; | |
291 | type MixedType = f32; | |
292 | ||
293 | #[inline] | |
294 | fn to_nalgebra(self) -> Self::NalgebraType { self } | |
295 | ||
296 | #[inline] | |
297 | fn to_nalgebra_mixed(self) -> Self::MixedType { self } | |
298 | ||
299 | #[inline] | |
300 | fn from_nalgebra(t : Self::NalgebraType) -> Self { t } | |
301 | ||
302 | #[inline] | |
303 | fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } | |
304 | ||
305 | } | |
306 | ||
307 | impl ToNalgebraRealField for f64 { | |
308 | type NalgebraType = f64; | |
309 | type MixedType = f64; | |
310 | ||
311 | #[inline] | |
312 | fn to_nalgebra(self) -> Self::NalgebraType { self } | |
313 | ||
314 | #[inline] | |
315 | fn to_nalgebra_mixed(self) -> Self::MixedType { self } | |
316 | ||
317 | #[inline] | |
318 | fn from_nalgebra(t : Self::NalgebraType) -> Self { t } | |
319 | ||
320 | #[inline] | |
321 | fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } | |
322 | } | |
323 |