comparison: src/discrete_gradient.rs
src/discrete_gradient.rs
- branch
- dev
- changeset 186
- afe04e6b4a5b
- parent 133
- 2b13f8a0c8ba
equal
deleted
inserted
replaced
| 185:e6829fbe2737 | 186:afe04e6b4a5b |
|---|---|
| 1 /*! | 1 /*! |
| 2 Simple disrete gradient operators | 2 Simple disrete gradient operators |
| 3 */ | 3 */ |
| 4 use crate::error::DynResult; | 4 use crate::error::DynResult; |
| 5 use crate::instance::Instance; | 5 use crate::instance::Instance; |
| 6 use crate::linops::{Adjointable, BoundedLinear, Linear, Mapping, GEMV}; | 6 use crate::linops::{Adjointable, BoundedLinear, Linear, Mapping, SimplyAdjointable, GEMV}; |
| 7 use crate::norms::{Norm, L2}; | 7 use crate::norms::{Norm, L2}; |
| 8 use crate::types::Float; | 8 use crate::types::Float; |
| 9 use nalgebra::{DVector, Dyn, Matrix, Storage, StorageMut, U1}; | 9 use nalgebra::{DVector, Dyn, Matrix, Storage, StorageMut, U1}; |
| 10 use numeric_literals::replace_float_literals; | 10 use numeric_literals::replace_float_literals; |
| 11 | 11 |
| 205 | 205 |
| 206 impl<'a, const N: usize> Iter<'a, N> { | 206 impl<'a, const N: usize> Iter<'a, N> { |
| 207 fn new(dims: &'a [usize; N], d: usize) -> Self { | 207 fn new(dims: &'a [usize; N], d: usize) -> Self { |
| 208 let d_stride = dims[0..d].iter().product::<usize>(); | 208 let d_stride = dims[0..d].iter().product::<usize>(); |
| 209 let len = dims.iter().product::<usize>(); | 209 let len = dims.iter().product::<usize>(); |
| 210 Iter { | 210 Iter { dims, d, d_stride, i: [0; N], k: 0, len } |
| 211 dims, | |
| 212 d, | |
| 213 d_stride, | |
| 214 i: [0; N], | |
| 215 k: 0, | |
| 216 len, | |
| 217 } | |
| 218 } | 211 } |
| 219 } | 212 } |
| 220 | 213 |
| 221 impl<'a, const N: usize> Iterator for Iter<'a, N> { | 214 impl<'a, const N: usize> Iterator for Iter<'a, N> { |
| 222 type Item = DiscretisationOrInterior; | 215 type Item = DiscretisationOrInterior; |
| 334 F: Float + nalgebra::RealField, | 327 F: Float + nalgebra::RealField, |
| 335 { | 328 { |
| 336 /// 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. |
| 337 pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> { | 330 pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> { |
| 338 if u.len() == dims.iter().product::<usize>() { | 331 if u.len() == dims.iter().product::<usize>() { |
| 339 Some(Grad { | 332 Some(Grad { dims, h, discretisation }) |
| 340 dims, | |
| 341 h, | |
| 342 discretisation, | |
| 343 }) | |
| 344 } else { | 333 } else { |
| 345 None | 334 None |
| 346 } | 335 } |
| 347 } | 336 } |
| 348 | 337 |
| 357 F: Float + nalgebra::RealField, | 346 F: Float + nalgebra::RealField, |
| 358 { | 347 { |
| 359 /// 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. |
| 360 pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> { | 349 pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> { |
| 361 if u.len() == dims.iter().product::<usize>() * N { | 350 if u.len() == dims.iter().product::<usize>() * N { |
| 362 Some(Div { | 351 Some(Div { dims, h, discretisation }) |
| 363 dims, | |
| 364 h, | |
| 365 discretisation, | |
| 366 }) | |
| 367 } else { | 352 } else { |
| 368 None | 353 None |
| 369 } | 354 } |
| 370 } | 355 } |
| 371 | 356 |
| 423 type Adjoint<'a> | 408 type Adjoint<'a> |
| 424 = Div<F, B::Opposite, N> | 409 = Div<F, B::Opposite, N> |
| 425 where | 410 where |
| 426 Self: 'a; | 411 Self: 'a; |
| 427 | 412 |
| 428 /// Form the adjoint operator of `self`. | |
| 429 fn adjoint(&self) -> Self::Adjoint<'_> { | 413 fn adjoint(&self) -> Self::Adjoint<'_> { |
| 430 Div { | 414 Div { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() } |
| 431 dims: self.dims, | 415 } |
| 432 h: -self.h, | 416 } |
| 433 discretisation: self.discretisation.opposite(), | 417 |
| 434 } | 418 impl<F, B, const N: usize> SimplyAdjointable<DVector<F>, DVector<F>> for Grad<F, B, N> |
| 419 where | |
| 420 B: Discretisation<F>, | |
| 421 F: Float + nalgebra::RealField, | |
| 422 { | |
| 423 type AdjointCodomain = DVector<F>; | |
| 424 type SimpleAdjoint = Div<F, B::Opposite, N>; | |
| 425 | |
| 426 fn adjoint(&self) -> Self::SimpleAdjoint { | |
| 427 Div { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() } | |
| 435 } | 428 } |
| 436 } | 429 } |
| 437 | 430 |
| 438 impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Div<F, B, N> | 431 impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Div<F, B, N> |
| 439 where | 432 where |
| 444 type Adjoint<'a> | 437 type Adjoint<'a> |
| 445 = Grad<F, B::Opposite, N> | 438 = Grad<F, B::Opposite, N> |
| 446 where | 439 where |
| 447 Self: 'a; | 440 Self: 'a; |
| 448 | 441 |
| 449 /// Form the adjoint operator of `self`. | |
| 450 fn adjoint(&self) -> Self::Adjoint<'_> { | 442 fn adjoint(&self) -> Self::Adjoint<'_> { |
| 451 Grad { | 443 Grad { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() } |
| 452 dims: self.dims, | 444 } |
| 453 h: -self.h, | 445 } |
| 454 discretisation: self.discretisation.opposite(), | 446 |
| 455 } | 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() } | |
| 456 } | 457 } |
| 457 } | 458 } |
| 458 | 459 |
| 459 #[cfg(test)] | 460 #[cfg(test)] |
| 460 mod tests { | 461 mod tests { |
