Mercurial > repos > alg_tools / file comparison

comparison: src/discrete_gradient.rs

src/discrete_gradient.rs

changeset 198
3868555d135c
parent 186
afe04e6b4a5b
equal deleted inserted replaced
94:1f19c6bbf07b 198:3868555d135c
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, SimplyAdjointable, 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 { dims, d, d_stride, i: [0; N], k: 0, len }
221 } 211 }
222 } 212 }
223 213
224 impl<'a, const N : usize> Iterator for Iter<'a, N> { 214 impl<'a, const N: usize> Iterator for Iter<'a, N> {
225 type Item = DiscretisationOrInterior; 215 type Item = DiscretisationOrInterior;
226 fn next(&mut self) -> Option<Self::Item> { 216 fn next(&mut self) -> Option<Self::Item> {
227 let res = if self.k >= self.len { 217 let res = if self.k >= self.len {
228 None 218 None
229 } else { 219 } else {
241 }; 231 };
242 self.k += 1; 232 self.k += 1;
243 for j in 0..N { 233 for j in 0..N {
244 if self.i[j] + 1 < self.dims[j] { 234 if self.i[j] + 1 < self.dims[j] {
245 self.i[j] += 1; 235 self.i[j] += 1;
246 break 236 break;
247 } 237 }
248 self.i[j] = 0 238 self.i[j] = 0
249 } 239 }
250 res 240 res
251 } 241 }
252 } 242 }
253 243
254 impl<F, B, const N : usize> Mapping<DVector<F>> 244 impl<F, B, const N: usize> Mapping<DVector<F>> for Grad<F, B, N>
255 for Grad<F, B, N> 245 where
256 where 246 B: Discretisation<F>,
257 B : Discretisation<F>, 247 F: Float + nalgebra::RealField,
258 F : Float + nalgebra::RealField,
259 { 248 {
260 type Codomain = DVector<F>; 249 type Codomain = DVector<F>;
261 fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> { 250 fn apply<I: Instance<DVector<F>>>(&self, i: I) -> DVector<F> {
262 let mut y = DVector::zeros(N * self.len()); 251 let mut y = DVector::zeros(N * self.len());
263 self.apply_add(&mut y, i); 252 self.apply_add(&mut y, i);
264 y 253 y
265 } 254 }
266 } 255 }
267 256
268 #[replace_float_literals(F::cast_from(literal))] 257 #[replace_float_literals(F::cast_from(literal))]
269 impl<F, B, const N : usize> GEMV<F, DVector<F>> 258 impl<F, B, const N: usize> GEMV<F, DVector<F>> for Grad<F, B, N>
270 for Grad<F, B, N> 259 where
271 where 260 B: Discretisation<F>,
272 B : Discretisation<F>, 261 F: Float + nalgebra::RealField,
273 F : Float + nalgebra::RealField, 262 {
274 { 263 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 { 264 if β == 0.0 {
279 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0); 265 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0);
280 } else if β != 1.0 { 266 } else if β != 1.0 {
281 //*y *= β; 267 //*y *= β;
282 y.as_mut_slice().iter_mut().for_each(|x| *x *= β); 268 y.as_mut_slice().iter_mut().for_each(|x| *x *= β);
284 let h = self.h; 270 let h = self.h;
285 let m = self.len(); 271 let m = self.len();
286 i.eval(|x| { 272 i.eval(|x| {
287 assert_eq!(x.len(), m); 273 assert_eq!(x.len(), m);
288 for d in 0..N { 274 for d in 0..N {
289 let mut v = y.generic_view_mut((d*m, 0), (Dyn(m), U1)); 275 let mut v = y.generic_view_mut((d * m, 0), (Dyn(m), U1));
290 for b in Iter::new(&self.dims, d) { 276 for b in Iter::new(&self.dims, d) {
291 self.discretisation.add_diff_mut(&mut v, x, α/h, b) 277 self.discretisation.add_diff_mut(&mut v, x, α / h, b)
292 } 278 }
293 } 279 }
294 }) 280 })
295 } 281 }
296 } 282 }
297 283
298 impl<F, B, const N : usize> Mapping<DVector<F>> 284 impl<F, B, const N: usize> Mapping<DVector<F>> for Div<F, B, N>
299 for Div<F, B, N> 285 where
300 where 286 B: Discretisation<F>,
301 B : Discretisation<F>, 287 F: Float + nalgebra::RealField,
302 F : Float + nalgebra::RealField,
303 { 288 {
304 type Codomain = DVector<F>; 289 type Codomain = DVector<F>;
305 fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> { 290 fn apply<I: Instance<DVector<F>>>(&self, i: I) -> DVector<F> {
306 let mut y = DVector::zeros(self.len()); 291 let mut y = DVector::zeros(self.len());
307 self.apply_add(&mut y, i); 292 self.apply_add(&mut y, i);
308 y 293 y
309 } 294 }
310 } 295 }
311 296
312 #[replace_float_literals(F::cast_from(literal))] 297 #[replace_float_literals(F::cast_from(literal))]
313 impl<F, B, const N : usize> GEMV<F, DVector<F>> 298 impl<F, B, const N: usize> GEMV<F, DVector<F>> for Div<F, B, N>
314 for Div<F, B, N> 299 where
315 where 300 B: Discretisation<F>,
316 B : Discretisation<F>, 301 F: Float + nalgebra::RealField,
317 F : Float + nalgebra::RealField, 302 {
318 { 303 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 { 304 if β == 0.0 {
323 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0); 305 y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0);
324 } else if β != 1.0 { 306 } else if β != 1.0 {
325 //*y *= β; 307 //*y *= β;
326 y.as_mut_slice().iter_mut().for_each(|x| *x *= β); 308 y.as_mut_slice().iter_mut().for_each(|x| *x *= β);
327 } 309 }
328 let h = self.h; 310 let h = self.h;
329 let m = self.len(); 311 let m = self.len();
330 i.eval(|x| { 312 i.eval_decompose(|x| {
331 assert_eq!(x.len(), N * m); 313 assert_eq!(x.len(), N * m);
332 for d in 0..N { 314 for d in 0..N {
333 let v = x.generic_view((d*m, 0), (Dyn(m), U1)); 315 let v = x.generic_view((d * m, 0), (Dyn(m), U1));
334 for b in Iter::new(&self.dims, d) { 316 for b in Iter::new(&self.dims, d) {
335 self.discretisation.add_diff_mut(y, &v, α/h, b) 317 self.discretisation.add_diff_mut(y, &v, α / h, b)
336 } 318 }
337 } 319 }
338 }) 320 })
339 } 321 }
340 } 322 }
341 323
342 impl<F, B, const N : usize> Grad<F, B, N> 324 impl<F, B, const N: usize> Grad<F, B, N>
343 where 325 where
344 B : Discretisation<F>, 326 B: Discretisation<F>,
345 F : Float + nalgebra::RealField 327 F: Float + nalgebra::RealField,
346 { 328 {
347 /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions. 329 /// 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) 330 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>() { 331 if u.len() == dims.iter().product::<usize>() {
352 Some(Grad { dims, h, discretisation } ) 332 Some(Grad { dims, h, discretisation })
353 } else { 333 } else {
354 None 334 None
355 } 335 }
356 } 336 }
357 337
358 fn len(&self) -> usize { 338 fn len(&self) -> usize {
359 self.dims.iter().product::<usize>() 339 self.dims.iter().product::<usize>()
360 } 340 }
361 } 341 }
362 342
363 343 impl<F, B, const N: usize> Div<F, B, N>
364 impl<F, B, const N : usize> Div<F, B, N> 344 where
365 where 345 B: Discretisation<F>,
366 B : Discretisation<F>, 346 F: Float + nalgebra::RealField,
367 F : Float + nalgebra::RealField
368 { 347 {
369 /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions. 348 /// 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) 349 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 { 350 if u.len() == dims.iter().product::<usize>() * N {
374 Some(Div { dims, h, discretisation } ) 351 Some(Div { dims, h, discretisation })
375 } else { 352 } else {
376 None 353 None
377 } 354 }
378 } 355 }
379 356
380 fn len(&self) -> usize { 357 fn len(&self) -> usize {
381 self.dims.iter().product::<usize>() 358 self.dims.iter().product::<usize>()
382 } 359 }
383 } 360 }
384 361
385 impl<F, B, const N : usize> Linear<DVector<F>> 362 impl<F, B, const N: usize> Linear<DVector<F>> for Grad<F, B, N>
386 for Grad<F, B, N> 363 where
387 where 364 B: Discretisation<F>,
388 B : Discretisation<F>, 365 F: Float + nalgebra::RealField,
389 F : Float + nalgebra::RealField, 366 {
390 { 367 }
391 } 368
392 369 impl<F, B, const N: usize> Linear<DVector<F>> for Div<F, B, N>
393 impl<F, B, const N : usize> Linear<DVector<F>> 370 where
394 for Div<F, B, N> 371 B: Discretisation<F>,
395 where 372 F: Float + nalgebra::RealField,
396 B : Discretisation<F>, 373 {
397 F : Float + nalgebra::RealField, 374 }
398 { 375
399 } 376 impl<F, B, const N: usize> BoundedLinear<DVector<F>, L2, L2, F> for Grad<F, B, N>
400 377 where
401 impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F> 378 B: Discretisation<F>,
402 for Grad<F, B, N> 379 F: Float + nalgebra::RealField,
403 where 380 DVector<F>: Norm<L2, F>,
404 B : Discretisation<F>, 381 {
405 F : Float + nalgebra::RealField, 382 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. 383 // Fuck nalgebra.
410 self.discretisation.opnorm_bound(num_traits::Float::abs(self.h)) 384 self.discretisation
411 } 385 .opnorm_bound(num_traits::Float::abs(self.h))
412 } 386 }
413 387 }
414 388
415 impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F> 389 impl<F, B, const N: usize> BoundedLinear<DVector<F>, L2, L2, F> for Div<F, B, N>
416 for Div<F, B, N> 390 where
417 where 391 B: Discretisation<F>,
418 B : Discretisation<F>, 392 F: Float + nalgebra::RealField,
419 F : Float + nalgebra::RealField, 393 DVector<F>: Norm<L2, F>,
420 DVector<F> : Norm<F, L2>, 394 {
421 { 395 fn opnorm_bound(&self, _: L2, _: L2) -> DynResult<F> {
422 fn opnorm_bound(&self, _ : L2, _ : L2) -> F {
423 // Fuck nalgebra. 396 // Fuck nalgebra.
424 self.discretisation.opnorm_bound(num_traits::Float::abs(self.h)) 397 self.discretisation
425 } 398 .opnorm_bound(num_traits::Float::abs(self.h))
426 } 399 }
427 400 }
428 impl<F, B, const N : usize> 401
429 Adjointable<DVector<F>, DVector<F>> 402 impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Grad<F, B, N>
430 for Grad<F, B, N> 403 where
431 where 404 B: Discretisation<F>,
432 B : Discretisation<F>, 405 F: Float + nalgebra::RealField,
433 F : Float + nalgebra::RealField,
434 { 406 {
435 type AdjointCodomain = DVector<F>; 407 type AdjointCodomain = DVector<F>;
436 type Adjoint<'a> = Div<F, B::Opposite, N> where Self : 'a; 408 type Adjoint<'a>
437 409 = Div<F, B::Opposite, N>
438 /// Form the adjoint operator of `self`. 410 where
411 Self: 'a;
412
439 fn adjoint(&self) -> Self::Adjoint<'_> { 413 fn adjoint(&self) -> Self::Adjoint<'_> {
440 Div { 414 Div { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
441 dims : self.dims, 415 }
442 h : -self.h, 416 }
443 discretisation : self.discretisation.opposite(), 417
444 } 418 impl<F, B, const N: usize> SimplyAdjointable<DVector<F>, DVector<F>> for Grad<F, B, N>
445 } 419 where
446 } 420 B: Discretisation<F>,
447 421 F: Float + nalgebra::RealField,
448
449 impl<F, B, const N : usize>
450 Adjointable<DVector<F>, DVector<F>>
451 for Div<F, B, N>
452 where
453 B : Discretisation<F>,
454 F : Float + nalgebra::RealField,
455 { 422 {
456 type AdjointCodomain = DVector<F>; 423 type AdjointCodomain = DVector<F>;
457 type Adjoint<'a> = Grad<F, B::Opposite, N> where Self : 'a; 424 type SimpleAdjoint = Div<F, B::Opposite, N>;
458 425
459 /// Form the adjoint operator of `self`. 426 fn adjoint(&self) -> Self::SimpleAdjoint {
427 Div { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
428 }
429 }
430
431 impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Div<F, B, N>
432 where
433 B: Discretisation<F>,
434 F: Float + nalgebra::RealField,
435 {
436 type AdjointCodomain = DVector<F>;
437 type Adjoint<'a>
438 = Grad<F, B::Opposite, N>
439 where
440 Self: 'a;
441
460 fn adjoint(&self) -> Self::Adjoint<'_> { 442 fn adjoint(&self) -> Self::Adjoint<'_> {
461 Grad { 443 Grad { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
462 dims : self.dims, 444 }
463 h : -self.h, 445 }
464 discretisation : self.discretisation.opposite(), 446
465 } 447 impl<F, B, const N: usize> SimplyAdjointable<DVector<F>, DVector<F>> for Div<F, B, N>
448 where
449 B: Discretisation<F>,
450 F: Float + nalgebra::RealField,
451 {
452 type AdjointCodomain = DVector<F>;
453 type SimpleAdjoint = Grad<F, B::Opposite, N>;
454
455 fn adjoint(&self) -> Self::SimpleAdjoint {
456 Grad { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
466 } 457 }
467 } 458 }
468 459
469 #[cfg(test)] 460 #[cfg(test)]
470 mod tests { 461 mod tests {
471 use super::*; 462 use super::*;
472 463
473 #[test] 464 #[test]
474 fn grad_adjoint() { 465 fn grad_adjoint() {
475 let im = DVector::from( (0..9).map(|t| t as f64).collect::<Vec<_>>()); 466 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<_>>()); 467 let v = DVector::from((0..18).map(|t| t as f64).collect::<Vec<_>>());
477 468
478 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardNeumann).unwrap(); 469 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardNeumann).unwrap();
479 let a = grad.apply(&im).dot(&v); 470 let a = grad.apply(&im).dot(&v);
480 let b = grad.adjoint().apply(&v).dot(&im); 471 let b = grad.adjoint().apply(&v).dot(&im);
481 assert_eq!(a, b); 472 assert_eq!(a, b);
482 473
483 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardDirichlet).unwrap(); 474 let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardDirichlet).unwrap();
484 let a = grad.apply(&im).dot(&v); 475 let a = grad.apply(&im).dot(&v);
485 let b = grad.adjoint().apply(&v).dot(&im); 476 let b = grad.adjoint().apply(&v).dot(&im);
486 assert_eq!(a, b); 477 assert_eq!(a, b);
487 478 }
488 } 479 }
489 }

Mercurial Repository: alg_tools

mercurial