Mercurial > repos > alg_tools / file comparison

comparison: src/discrete_gradient.rs

src/discrete_gradient.rs

branch
dev
changeset 110
a1278320be26
parent 74
2c76df38d02b
child 124
6aa955ad8122
equal deleted inserted replaced
109:943c6b3b9414 110:a1278320be26
1 /*! 1 /*!
2 Simple disrete gradient operators 2 Simple disrete gradient operators
3 */ 3 */
4 use crate::error::DynResult;
5 use crate::instance::Instance;
6 use crate::linops::{Adjointable, BoundedLinear, Linear, Mapping, GEMV};
7 use crate::norms::{Norm, L2};
8 use crate::types::Float;
9 use nalgebra::{DVector, Dyn, Matrix, Storage, StorageMut, U1};
4 use numeric_literals::replace_float_literals; 10 use numeric_literals::replace_float_literals;
5 use nalgebra::{
6 DVector, Matrix, U1, Storage, StorageMut, Dyn
7 };
8 use crate::types::Float;
9 use crate::instance::Instance;
10 use crate::linops::{Mapping, Linear, BoundedLinear, Adjointable, GEMV};
11 use crate::norms::{Norm, L2};
12 11
13 #[derive(Copy, Clone, Debug)] 12 #[derive(Copy, Clone, Debug)]
14 /// Forward differences with Neumann boundary conditions 13 /// Forward differences with Neumann boundary conditions
15 pub struct ForwardNeumann; 14 pub struct ForwardNeumann;
16 15
25 #[derive(Copy, Clone, Debug)] 24 #[derive(Copy, Clone, Debug)]
26 /// Backward differences with Neumann boundary conditions 25 /// Backward differences with Neumann boundary conditions
27 pub struct BackwardNeumann; 26 pub struct BackwardNeumann;
28 27
29 /// Finite differences gradient 28 /// Finite differences gradient
30 pub struct Grad< 29 pub struct Grad<F: Float + nalgebra::RealField, B: Discretisation<F>, const N: usize> {
31 F : Float + nalgebra::RealField, 30 dims: [usize; N],
32 B : Discretisation<F>, 31 h: F, // may be negative to implement adjoints!
33 const N : usize 32 discretisation: B,
34 > { 33 }
35 dims : [usize; N],
36 h : F, // may be negative to implement adjoints!
37 discretisation : B,
38 }
39
40 34
41 /// Finite differences divergence 35 /// Finite differences divergence
42 pub struct Div< 36 pub struct Div<F: Float + nalgebra::RealField, B: Discretisation<F>, const N: usize> {
43 F : Float + nalgebra::RealField, 37 dims: [usize; N],
44 B : Discretisation<F>, 38 h: F, // may be negative to implement adjoints!
45 const N : usize 39 discretisation: B,
46 > {
47 dims : [usize; N],
48 h : F, // may be negative to implement adjoints!
49 discretisation : B,
50 } 40 }
51 41
52 /// Internal: classification of a point in a 1D discretisation 42 /// Internal: classification of a point in a 1D discretisation
53 pub enum DiscretisationOrInterior { 43 pub enum DiscretisationOrInterior {
54 /// center, forward 44 /// center, forward
60 } 50 }
61 51
62 use DiscretisationOrInterior::*; 52 use DiscretisationOrInterior::*;
63 53
64 /// Trait for different discretisations 54 /// Trait for different discretisations
65 pub trait Discretisation<F : Float + nalgebra::RealField> : Copy { 55 pub trait Discretisation<F: Float + nalgebra::RealField>: Copy {
66 /// Opposite discretisation, appropriate for adjoints with negated cell width. 56 /// Opposite discretisation, appropriate for adjoints with negated cell width.
67 type Opposite : Discretisation<F>; 57 type Opposite: Discretisation<F>;
68 58
69 /// Add to appropiate index of `v` (as determined by `b`) the appropriate difference 59 /// Add to appropiate index of `v` (as determined by `b`) the appropriate difference
70 /// of `x` with cell width `h`. 60 /// of `x` with cell width `h`.
71 fn add_diff_mut<SMut, S>( 61 fn add_diff_mut<SMut, S>(
72 &self, 62 &self,
73 v : &mut Matrix<F, Dyn, U1, SMut>, 63 v: &mut Matrix<F, Dyn, U1, SMut>,
74 x : &Matrix<F, Dyn, U1, S>, 64 x: &Matrix<F, Dyn, U1, S>,
75 α : F, 65 α: F,
76 b : DiscretisationOrInterior, 66 b: DiscretisationOrInterior,
77 ) where 67 ) where
78 SMut : StorageMut<F, Dyn, U1>, 68 SMut: StorageMut<F, Dyn, U1>,
79 S : Storage<F, Dyn, U1>; 69 S: Storage<F, Dyn, U1>;
80 70
81 /// Give the opposite discretisation, appropriate for adjoints with negated `h`. 71 /// Give the opposite discretisation, appropriate for adjoints with negated `h`.
82 fn opposite(&self) -> Self::Opposite; 72 fn opposite(&self) -> Self::Opposite;
83 73
84 /// Bound for the corresponding operator norm. 74 /// Bound for the corresponding operator norm.
85 #[replace_float_literals(F::cast_from(literal))] 75 #[replace_float_literals(F::cast_from(literal))]
86 fn opnorm_bound(&self, h : F) -> F { 76 fn opnorm_bound(&self, h: F) -> DynResult<F> {
87 // See: Chambolle, “An Algorithm for Total Variation Minimization and Applications”. 77 // See: Chambolle, “An Algorithm for Total Variation Minimization and Applications”.
88 // Ok for forward and backward differences. 78 // Ok for forward and backward differences.
89 // 79 //
90 // Fuck nalgebra for polluting everything with its own shit. 80 // Fuck nalgebra for polluting everything with its own shit.
91 num_traits::Float::sqrt(8.0) / h 81 Ok(num_traits::Float::sqrt(8.0) / h)
92 } 82 }
93 } 83 }
94 84
95 impl<F : Float + nalgebra::RealField> Discretisation<F> for ForwardNeumann { 85 impl<F: Float + nalgebra::RealField> Discretisation<F> for ForwardNeumann {
96 type Opposite = BackwardDirichlet; 86 type Opposite = BackwardDirichlet;
97 87
98 #[inline] 88 #[inline]
99 fn add_diff_mut<SMut, S>( 89 fn add_diff_mut<SMut, S>(
100 &self, 90 &self,
101 v : &mut Matrix<F, Dyn, U1, SMut>, 91 v: &mut Matrix<F, Dyn, U1, SMut>,
102 x : &Matrix<F, Dyn, U1, S>, 92 x: &Matrix<F, Dyn, U1, S>,
103 α : F, 93 α: F,
104 b : DiscretisationOrInterior, 94 b: DiscretisationOrInterior,
105 ) where 95 ) where
106 SMut : StorageMut<F, Dyn, U1>, 96 SMut: StorageMut<F, Dyn, U1>,
107 S : Storage<F, Dyn, U1> 97 S: Storage<F, Dyn, U1>,
108 { 98 {
109 match b { 99 match b {
110 Interior(c, (_, f)) | LeftBoundary(c, f) => { v[c] += (x[f] - x[c]) * α }, 100 Interior(c, (_, f)) | LeftBoundary(c, f) => v[c] += (x[f] - x[c]) * α,
111 RightBoundary(_c, _b) => { }, 101 RightBoundary(_c, _b) => {}
112 } 102 }
113 } 103 }
114 104
115 #[inline] 105 #[inline]
116 fn opposite(&self) -> Self::Opposite { 106 fn opposite(&self) -> Self::Opposite {
117 BackwardDirichlet 107 BackwardDirichlet
118 } 108 }
119 } 109 }
120 110
121 impl<F : Float + nalgebra::RealField> Discretisation<F> for BackwardNeumann { 111 impl<F: Float + nalgebra::RealField> Discretisation<F> for BackwardNeumann {
122 type Opposite = ForwardDirichlet; 112 type Opposite = ForwardDirichlet;
123 113
124 #[inline] 114 #[inline]
125 fn add_diff_mut<SMut, S>( 115 fn add_diff_mut<SMut, S>(
126 &self, 116 &self,
127 v : &mut Matrix<F, Dyn, U1, SMut>, 117 v: &mut Matrix<F, Dyn, U1, SMut>,
128 x : &Matrix<F, Dyn, U1, S>, 118 x: &Matrix<F, Dyn, U1, S>,
129 α : F, 119 α: F,
130 b : DiscretisationOrInterior, 120 b: DiscretisationOrInterior,
131 ) where 121 ) where
132 SMut : StorageMut<F, Dyn, U1>, 122 SMut: StorageMut<F, Dyn, U1>,
133 S : Storage<F, Dyn, U1> 123 S: Storage<F, Dyn, U1>,
134 { 124 {
135 match b { 125 match b {
136 Interior(c, (b, _)) | RightBoundary(c, b) => { v[c] += (x[c] - x[b]) * α }, 126 Interior(c, (b, _)) | RightBoundary(c, b) => v[c] += (x[c] - x[b]) * α,
137 LeftBoundary(_c, _f) => { }, 127 LeftBoundary(_c, _f) => {}
138 } 128 }
139 } 129 }
140 130
141 #[inline] 131 #[inline]
142 fn opposite(&self) -> Self::Opposite { 132 fn opposite(&self) -> Self::Opposite {
143 ForwardDirichlet 133 ForwardDirichlet
144 } 134 }
145 } 135 }
146 136
147 impl<F : Float + nalgebra::RealField> Discretisation<F> for BackwardDirichlet { 137 impl<F: Float + nalgebra::RealField> Discretisation<F> for BackwardDirichlet {
148 type Opposite = ForwardNeumann; 138 type Opposite = ForwardNeumann;
149 139
150 #[inline] 140 #[inline]
151 fn add_diff_mut<SMut, S>( 141 fn add_diff_mut<SMut, S>(
152 &self, 142 &self,
153 v : &mut Matrix<F, Dyn, U1, SMut>, 143 v: &mut Matrix<F, Dyn, U1, SMut>,
154 x : &Matrix<F, Dyn, U1, S>, 144 x: &Matrix<F, Dyn, U1, S>,
155 α : F, 145 α: F,
156 b : DiscretisationOrInterior, 146 b: DiscretisationOrInterior,
157 ) where 147 ) where
158 SMut : StorageMut<F, Dyn, U1>, 148 SMut: StorageMut<F, Dyn, U1>,
159 S : Storage<F, Dyn, U1> 149 S: Storage<F, Dyn, U1>,
160 { 150 {
161 match b { 151 match b {
162 Interior(c, (b, _f)) => { v[c] += (x[c] - x[b]) * α }, 152 Interior(c, (b, _f)) => v[c] += (x[c] - x[b]) * α,
163 LeftBoundary(c, _f) => { v[c] += x[c] * α }, 153 LeftBoundary(c, _f) => v[c] += x[c] * α,
164 RightBoundary(c, b) => { v[c] -= x[b] * α }, 154 RightBoundary(c, b) => v[c] -= x[b] * α,
165 } 155 }
166 } 156 }
167 157
168 #[inline] 158 #[inline]
169 fn opposite(&self) -> Self::Opposite { 159 fn opposite(&self) -> Self::Opposite {
170 ForwardNeumann 160 ForwardNeumann
171 } 161 }
172 } 162 }
173 163
174 impl<F : Float + nalgebra::RealField> Discretisation<F> for ForwardDirichlet { 164 impl<F: Float + nalgebra::RealField> Discretisation<F> for ForwardDirichlet {
175 type Opposite = BackwardNeumann; 165 type Opposite = BackwardNeumann;
176 166
177 #[inline] 167 #[inline]
178 fn add_diff_mut<SMut, S>( 168 fn add_diff_mut<SMut, S>(
179 &self, 169 &self,
180 v : &mut Matrix<F, Dyn, U1, SMut>, 170 v: &mut Matrix<F, Dyn, U1, SMut>,
181 x : &Matrix<F, Dyn, U1, S>, 171 x: &Matrix<F, Dyn, U1, S>,
182 α : F, 172 α: F,
183 b : DiscretisationOrInterior, 173 b: DiscretisationOrInterior,
184 ) where 174 ) where
185 SMut : StorageMut<F, Dyn, U1>, 175 SMut: StorageMut<F, Dyn, U1>,
186 S : Storage<F, Dyn, U1> 176 S: Storage<F, Dyn, U1>,
187 { 177 {
188 match b { 178 match b {
189 Interior(c, (_b, f)) => { v[c] += (x[f] - x[c]) * α }, 179 Interior(c, (_b, f)) => v[c] += (x[f] - x[c]) * α,
190 LeftBoundary(c, f) => { v[c] += x[f] * α }, 180 LeftBoundary(c, f) => v[c] += x[f] * α,
191 RightBoundary(c, _b) => { v[c] -= x[c] * α }, 181 RightBoundary(c, _b) => v[c] -= x[c] * α,
192 } 182 }
193 } 183 }
194 184
195 #[inline] 185 #[inline]
196 fn opposite(&self) -> Self::Opposite { 186 fn opposite(&self) -> Self::Opposite {
197 BackwardNeumann 187 BackwardNeumann
198 } 188 }
199 } 189 }
200 190
201 struct Iter<'a, const N : usize> { 191 struct Iter<'a, const N: usize> {
202 /// Dimensions 192 /// Dimensions
203 dims : &'a [usize; N], 193 dims: &'a [usize; N],
204 /// Dimension along which to calculate differences 194 /// Dimension along which to calculate differences
205 d : usize, 195 d: usize,
206 /// Stride along coordinate d 196 /// Stride along coordinate d
207 d_stride : usize, 197 d_stride: usize,
208 /// Cartesian indices 198 /// Cartesian indices
209 i : [usize; N], 199 i: [usize; N],
210 /// Linear index 200 /// Linear index
211 k : usize, 201 k: usize,
212 /// Maximal linear index 202 /// Maximal linear index
213 len : usize 203 len: usize,
214 } 204 }
215 205
216 impl<'a, const N : usize> Iter<'a, N> { 206 impl<'a, const N: usize> Iter<'a, N> {
217 fn new(dims : &'a [usize; N], d : usize) -> Self { 207 fn new(dims: &'a [usize; N], d: usize) -> Self {
218 let d_stride = dims[0..d].iter().product::<usize>(); 208 let d_stride = dims[0..d].iter().product::<usize>();
219 let len = dims.iter().product::<usize>(); 209 let len = dims.iter().product::<usize>();
220 Iter{ dims, d, d_stride, i : [0; N], k : 0, len } 210 Iter {
221 } 211 dims,
222 } 212 d,
223 213 d_stride,
224 impl<'a, const N : usize> Iterator for Iter<'a, N> { 214 i: [0; N],
215 k: 0,
216 len,
217 }
218 }
219 }
220
221 impl<'a, const N: usize> Iterator for Iter<'a, N> {
225 type Item = DiscretisationOrInterior; 222 type Item = DiscretisationOrInterior;
226 fn next(&mut self) -> Option<Self::Item> { 223 fn next(&mut self) -> Option<Self::Item> {
227 let res = if self.k >= self.len { 224 let res = if self.k >= self.len {
228 None 225 None
229 } else { 226 } else {
241 }; 238 };
242 self.k += 1; 239 self.k += 1;
243 for j in 0..N { 240 for j in 0..N {
244 if self.i[j] + 1 < self.dims[j] { 241 if self.i[j] + 1 < self.dims[j] {
245 self.i[j] += 1; 242 self.i[j] += 1;
246 break 243 break;
247 } 244 }
248 self.i[j] = 0 245 self.i[j] = 0
249 } 246 }
250 res 247 res
251 } 248 }
252 } 249 }
253 250
254 impl<F, B, const N : usize> Mapping<DVector<F>> 251 impl<F, B, const N: usize> Mapping<DVector<F>> for Grad<F, B, N>
255 for Grad<F, B, N> 252 where
256 where 253 B: Discretisation<F>,
257 B : Discretisation<F>, 254 F: Float + nalgebra::RealField,
258 F : Float + nalgebra::RealField,
259 { 255 {
260 type Codomain = DVector<F>; 256 type Codomain = DVector<F>;
261 fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> { 257 fn apply<I: Instance<DVector<F>>>(&self, i: I) -> DVector<F> {
262 let mut y = DVector::zeros(N * self.len()); 258 let mut y = DVector::zeros(N * self.len());
263 self.apply_add(&mut y, i); 259 self.apply_add(&mut y, i);
264 y 260 y
265 } 261 }
266 } 262 }
267 263
268 #[replace_float_literals(F::cast_from(literal))] 264 #[replace_float_literals(F::cast_from(literal))]
269 impl<F, B, const N : usize> GEMV<F, DVector<F>> 265 impl<F, B, const N: usize> GEMV<F, DVector<F>> for Grad<F, B, N>
270 for Grad<F, B, N> 266 where
271 where 267 B: Discretisation<F>,
272 B : Discretisation<F>, 268 F: Float + nalgebra::RealField,
273 F : Float + nalgebra::RealField, 269 {
274 { 270 fn gemv<I: Instance<DVector<F>>>(&self, y: &mut DVector<F>, α: F, i: I, β: F) {
275 fn gemv<I : Instance<DVector<F>>>(
276 &self, y : &mut DVector<F>, α : F, i : I, β : F
277 ) {
278 if β == 0.0 { 271 if β == 0.0 {
279 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0); 272 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0);
280 } else if β != 1.0 { 273 } else if β != 1.0 {
281 //*y *= β; 274 //*y *= β;
282 y.as_mut_slice().iter_mut().for_each(|x| *x *= β); 275 y.as_mut_slice().iter_mut().for_each(|x| *x *= β);
284 let h = self.h; 277 let h = self.h;
285 let m = self.len(); 278 let m = self.len();
286 i.eval(|x| { 279 i.eval(|x| {
287 assert_eq!(x.len(), m); 280 assert_eq!(x.len(), m);
288 for d in 0..N { 281 for d in 0..N {
289 let mut v = y.generic_view_mut((d*m, 0), (Dyn(m), U1)); 282 let mut v = y.generic_view_mut((d * m, 0), (Dyn(m), U1));
290 for b in Iter::new(&self.dims, d) { 283 for b in Iter::new(&self.dims, d) {
291 self.discretisation.add_diff_mut(&mut v, x, α/h, b) 284 self.discretisation.add_diff_mut(&mut v, x, α / h, b)
292 } 285 }
293 } 286 }
294 }) 287 })
295 } 288 }
296 } 289 }
297 290
298 impl<F, B, const N : usize> Mapping<DVector<F>> 291 impl<F, B, const N: usize> Mapping<DVector<F>> for Div<F, B, N>
299 for Div<F, B, N> 292 where
300 where 293 B: Discretisation<F>,
301 B : Discretisation<F>, 294 F: Float + nalgebra::RealField,
302 F : Float + nalgebra::RealField,
303 { 295 {
304 type Codomain = DVector<F>; 296 type Codomain = DVector<F>;
305 fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> { 297 fn apply<I: Instance<DVector<F>>>(&self, i: I) -> DVector<F> {
306 let mut y = DVector::zeros(self.len()); 298 let mut y = DVector::zeros(self.len());
307 self.apply_add(&mut y, i); 299 self.apply_add(&mut y, i);
308 y 300 y
309 } 301 }
310 } 302 }
311 303
312 #[replace_float_literals(F::cast_from(literal))] 304 #[replace_float_literals(F::cast_from(literal))]
313 impl<F, B, const N : usize> GEMV<F, DVector<F>> 305 impl<F, B, const N: usize> GEMV<F, DVector<F>> for Div<F, B, N>
314 for Div<F, B, N> 306 where
315 where 307 B: Discretisation<F>,
316 B : Discretisation<F>, 308 F: Float + nalgebra::RealField,
317 F : Float + nalgebra::RealField, 309 {
318 { 310 fn gemv<I: Instance<DVector<F>>>(&self, y: &mut DVector<F>, α: F, i: I, β: F) {
319 fn gemv<I : Instance<DVector<F>>>(
320 &self, y : &mut DVector<F>, α : F, i : I, β : F
321 ) {
322 if β == 0.0 { 311 if β == 0.0 {
323 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0); 312 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0);
324 } else if β != 1.0 { 313 } else if β != 1.0 {
325 //*y *= β; 314 //*y *= β;
326 y.as_mut_slice().iter_mut().for_each(|x| *x *= β); 315 y.as_mut_slice().iter_mut().for_each(|x| *x *= β);
328 let h = self.h; 317 let h = self.h;
329 let m = self.len(); 318 let m = self.len();
330 i.eval(|x| { 319 i.eval(|x| {
331 assert_eq!(x.len(), N * m); 320 assert_eq!(x.len(), N * m);
332 for d in 0..N { 321 for d in 0..N {
333 let v = x.generic_view((d*m, 0), (Dyn(m), U1)); 322 let v = x.generic_view((d * m, 0), (Dyn(m), U1));
334 for b in Iter::new(&self.dims, d) { 323 for b in Iter::new(&self.dims, d) {
335 self.discretisation.add_diff_mut(y, &v, α/h, b) 324 self.discretisation.add_diff_mut(y, &v, α / h, b)
336 } 325 }
337 } 326 }
338 }) 327 })
339 } 328 }
340 } 329 }
341 330
342 impl<F, B, const N : usize> Grad<F, B, N> 331 impl<F, B, const N: usize> Grad<F, B, N>
343 where 332 where
344 B : Discretisation<F>, 333 B: Discretisation<F>,
345 F : Float + nalgebra::RealField 334 F: Float + nalgebra::RealField,
346 { 335 {
347 /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions. 336 /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions.
348 pub fn new_for(u : &DVector<F>, h : F, dims : [usize; N], discretisation : B) 337 pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> {
349 -> Option<Self>
350 {
351 if u.len() == dims.iter().product::<usize>() { 338 if u.len() == dims.iter().product::<usize>() {
352 Some(Grad { dims, h, discretisation } ) 339 Some(Grad {
340 dims,
341 h,
342 discretisation,
343 })
353 } else { 344 } else {
354 None 345 None
355 } 346 }
356 } 347 }
357 348
358 fn len(&self) -> usize { 349 fn len(&self) -> usize {
359 self.dims.iter().product::<usize>() 350 self.dims.iter().product::<usize>()
360 } 351 }
361 } 352 }
362 353
363 354 impl<F, B, const N: usize> Div<F, B, N>
364 impl<F, B, const N : usize> Div<F, B, N> 355 where
365 where 356 B: Discretisation<F>,
366 B : Discretisation<F>, 357 F: Float + nalgebra::RealField,
367 F : Float + nalgebra::RealField
368 { 358 {
369 /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions. 359 /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions.
370 pub fn new_for(u : &DVector<F>, h : F, dims : [usize; N], discretisation : B) 360 pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> {
371 -> Option<Self>
372 {
373 if u.len() == dims.iter().product::<usize>() * N { 361 if u.len() == dims.iter().product::<usize>() * N {
374 Some(Div { dims, h, discretisation } ) 362 Some(Div {
363 dims,
364 h,
365 discretisation,
366 })
375 } else { 367 } else {
376 None 368 None
377 } 369 }
378 } 370 }
379 371
380 fn len(&self) -> usize { 372 fn len(&self) -> usize {
381 self.dims.iter().product::<usize>() 373 self.dims.iter().product::<usize>()
382 } 374 }
383 } 375 }
384 376
385 impl<F, B, const N : usize> Linear<DVector<F>> 377 impl<F, B, const N: usize> Linear<DVector<F>> for Grad<F, B, N>
386 for Grad<F, B, N> 378 where
387 where 379 B: Discretisation<F>,
388 B : Discretisation<F>, 380 F: Float + nalgebra::RealField,
389 F : Float + nalgebra::RealField, 381 {
390 { 382 }
391 } 383
392 384 impl<F, B, const N: usize> Linear<DVector<F>> for Div<F, B, N>
393 impl<F, B, const N : usize> Linear<DVector<F>> 385 where
394 for Div<F, B, N> 386 B: Discretisation<F>,
395 where 387 F: Float + nalgebra::RealField,
396 B : Discretisation<F>, 388 {
397 F : Float + nalgebra::RealField, 389 }
398 { 390
399 } 391 impl<F, B, const N: usize> BoundedLinear<DVector<F>, L2, L2, F> for Grad<F, B, N>
400 392 where
401 impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F> 393 B: Discretisation<F>,
402 for Grad<F, B, N> 394 F: Float + nalgebra::RealField,
403 where 395 DVector<F>: Norm<F, L2>,
404 B : Discretisation<F>, 396 {
405 F : Float + nalgebra::RealField, 397 fn opnorm_bound(&self, _: L2, _: L2) -> DynResult<F> {
406 DVector<F> : Norm<F, L2>,
407 {
408 fn opnorm_bound(&self, _ : L2, _ : L2) -> F {
409 // Fuck nalgebra. 398 // Fuck nalgebra.
410 self.discretisation.opnorm_bound(num_traits::Float::abs(self.h)) 399 self.discretisation
411 } 400 .opnorm_bound(num_traits::Float::abs(self.h))
412 } 401 }
413 402 }
414 403
415 impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F> 404 impl<F, B, const N: usize> BoundedLinear<DVector<F>, L2, L2, F> for Div<F, B, N>
416 for Div<F, B, N> 405 where
417 where 406 B: Discretisation<F>,
418 B : Discretisation<F>, 407 F: Float + nalgebra::RealField,
419 F : Float + nalgebra::RealField, 408 DVector<F>: Norm<F, L2>,
420 DVector<F> : Norm<F, L2>, 409 {
421 { 410 fn opnorm_bound(&self, _: L2, _: L2) -> DynResult<F> {
422 fn opnorm_bound(&self, _ : L2, _ : L2) -> F {
423 // Fuck nalgebra. 411 // Fuck nalgebra.
424 self.discretisation.opnorm_bound(num_traits::Float::abs(self.h)) 412 self.discretisation
425 } 413 .opnorm_bound(num_traits::Float::abs(self.h))
426 } 414 }
427 415 }
428 impl<F, B, const N : usize> 416
429 Adjointable<DVector<F>, DVector<F>> 417 impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Grad<F, B, N>
430 for Grad<F, B, N> 418 where
431 where 419 B: Discretisation<F>,
432 B : Discretisation<F>, 420 F: Float + nalgebra::RealField,
433 F : Float + nalgebra::RealField,
434 { 421 {
435 type AdjointCodomain = DVector<F>; 422 type AdjointCodomain = DVector<F>;
436 type Adjoint<'a> = Div<F, B::Opposite, N> where Self : 'a; 423 type Adjoint<'a>
424 = Div<F, B::Opposite, N>
425 where
426 Self: 'a;
437 427
438 /// Form the adjoint operator of `self`. 428 /// Form the adjoint operator of `self`.
439 fn adjoint(&self) -> Self::Adjoint<'_> { 429 fn adjoint(&self) -> Self::Adjoint<'_> {
440 Div { 430 Div {
441 dims : self.dims, 431 dims: self.dims,
442 h : -self.h, 432 h: -self.h,
443 discretisation : self.discretisation.opposite(), 433 discretisation: self.discretisation.opposite(),
444 } 434 }
445 } 435 }
446 } 436 }
447 437
448 438 impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Div<F, B, N>
449 impl<F, B, const N : usize> 439 where
450 Adjointable<DVector<F>, DVector<F>> 440 B: Discretisation<F>,
451 for Div<F, B, N> 441 F: Float + nalgebra::RealField,
452 where
453 B : Discretisation<F>,
454 F : Float + nalgebra::RealField,
455 { 442 {
456 type AdjointCodomain = DVector<F>; 443 type AdjointCodomain = DVector<F>;
457 type Adjoint<'a> = Grad<F, B::Opposite, N> where Self : 'a; 444 type Adjoint<'a>
445 = Grad<F, B::Opposite, N>
446 where
447 Self: 'a;
458 448
459 /// Form the adjoint operator of `self`. 449 /// Form the adjoint operator of `self`.
460 fn adjoint(&self) -> Self::Adjoint<'_> { 450 fn adjoint(&self) -> Self::Adjoint<'_> {
461 Grad { 451 Grad {
462 dims : self.dims, 452 dims: self.dims,
463 h : -self.h, 453 h: -self.h,
464 discretisation : self.discretisation.opposite(), 454 discretisation: self.discretisation.opposite(),
465 } 455 }
466 } 456 }
467 } 457 }
468 458
469 #[cfg(test)] 459 #[cfg(test)]
470 mod tests { 460 mod tests {
471 use super::*; 461 use super::*;
472 462
473 #[test] 463 #[test]
474 fn grad_adjoint() { 464 fn grad_adjoint() {
475 let im = DVector::from( (0..9).map(|t| t as f64).collect::<Vec<_>>()); 465 let im = DVector::from((0..9).map(|t| t as f64).collect::<Vec<_>>());
476 let v = DVector::from( (0..18).map(|t| t as f64).collect::<Vec<_>>()); 466 let v = DVector::from((0..18).map(|t| t as f64).collect::<Vec<_>>());
477 467
478 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardNeumann).unwrap(); 468 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardNeumann).unwrap();
479 let a = grad.apply(&im).dot(&v); 469 let a = grad.apply(&im).dot(&v);
480 let b = grad.adjoint().apply(&v).dot(&im); 470 let b = grad.adjoint().apply(&v).dot(&im);
481 assert_eq!(a, b); 471 assert_eq!(a, b);
482 472
483 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardDirichlet).unwrap(); 473 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardDirichlet).unwrap();
484 let a = grad.apply(&im).dot(&v); 474 let a = grad.apply(&im).dot(&v);
485 let b = grad.adjoint().apply(&v).dot(&im); 475 let b = grad.adjoint().apply(&v).dot(&im);
486 assert_eq!(a, b); 476 assert_eq!(a, b);
487 477 }
488 } 478 }
489 }

Mercurial Repository: alg_tools

mercurial