--- a/src/discrete_gradient.rs Wed Apr 30 01:06:25 2025 -0500
+++ b/src/discrete_gradient.rs Wed Apr 30 16:39:01 2025 -0500
@@ -1,14 +1,13 @@
/*!
Simple disrete gradient operators
*/
+use crate::error::DynResult;
+use crate::instance::Instance;
+use crate::linops::{Adjointable, BoundedLinear, Linear, Mapping, GEMV};
+use crate::norms::{Norm, L2};
+use crate::types::Float;
+use nalgebra::{DVector, Dyn, Matrix, Storage, StorageMut, U1};
use numeric_literals::replace_float_literals;
-use nalgebra::{
- DVector, Matrix, U1, Storage, StorageMut, Dyn
-};
-use crate::types::Float;
-use crate::instance::Instance;
-use crate::linops::{Mapping, Linear, BoundedLinear, Adjointable, GEMV};
-use crate::norms::{Norm, L2};
#[derive(Copy, Clone, Debug)]
/// Forward differences with Neumann boundary conditions
@@ -27,26 +26,17 @@
pub struct BackwardNeumann;
/// Finite differences gradient
-pub struct Grad<
- F : Float + nalgebra::RealField,
- B : Discretisation<F>,
- const N : usize
-> {
- dims : [usize; N],
- h : F, // may be negative to implement adjoints!
- discretisation : B,
+pub struct Grad<F: Float + nalgebra::RealField, B: Discretisation<F>, const N: usize> {
+ dims: [usize; N],
+ h: F, // may be negative to implement adjoints!
+ discretisation: B,
}
-
/// Finite differences divergence
-pub struct Div<
- F : Float + nalgebra::RealField,
- B : Discretisation<F>,
- const N : usize
-> {
- dims : [usize; N],
- h : F, // may be negative to implement adjoints!
- discretisation : B,
+pub struct Div<F: Float + nalgebra::RealField, B: Discretisation<F>, const N: usize> {
+ dims: [usize; N],
+ h: F, // may be negative to implement adjoints!
+ discretisation: B,
}
/// Internal: classification of a point in a 1D discretisation
@@ -62,53 +52,53 @@
use DiscretisationOrInterior::*;
/// Trait for different discretisations
-pub trait Discretisation<F : Float + nalgebra::RealField> : Copy {
+pub trait Discretisation<F: Float + nalgebra::RealField>: Copy {
/// Opposite discretisation, appropriate for adjoints with negated cell width.
- type Opposite : Discretisation<F>;
+ type Opposite: Discretisation<F>;
/// Add to appropiate index of `v` (as determined by `b`) the appropriate difference
/// of `x` with cell width `h`.
fn add_diff_mut<SMut, S>(
&self,
- v : &mut Matrix<F, Dyn, U1, SMut>,
- x : &Matrix<F, Dyn, U1, S>,
- α : F,
- b : DiscretisationOrInterior,
+ v: &mut Matrix<F, Dyn, U1, SMut>,
+ x: &Matrix<F, Dyn, U1, S>,
+ α: F,
+ b: DiscretisationOrInterior,
) where
- SMut : StorageMut<F, Dyn, U1>,
- S : Storage<F, Dyn, U1>;
+ SMut: StorageMut<F, Dyn, U1>,
+ S: Storage<F, Dyn, U1>;
/// Give the opposite discretisation, appropriate for adjoints with negated `h`.
fn opposite(&self) -> Self::Opposite;
/// Bound for the corresponding operator norm.
#[replace_float_literals(F::cast_from(literal))]
- fn opnorm_bound(&self, h : F) -> F {
+ fn opnorm_bound(&self, h: F) -> DynResult<F> {
// See: Chambolle, “An Algorithm for Total Variation Minimization and Applications”.
// Ok for forward and backward differences.
//
// Fuck nalgebra for polluting everything with its own shit.
- num_traits::Float::sqrt(8.0) / h
+ Ok(num_traits::Float::sqrt(8.0) / h)
}
}
-impl<F : Float + nalgebra::RealField> Discretisation<F> for ForwardNeumann {
+impl<F: Float + nalgebra::RealField> Discretisation<F> for ForwardNeumann {
type Opposite = BackwardDirichlet;
#[inline]
fn add_diff_mut<SMut, S>(
&self,
- v : &mut Matrix<F, Dyn, U1, SMut>,
- x : &Matrix<F, Dyn, U1, S>,
- α : F,
- b : DiscretisationOrInterior,
+ v: &mut Matrix<F, Dyn, U1, SMut>,
+ x: &Matrix<F, Dyn, U1, S>,
+ α: F,
+ b: DiscretisationOrInterior,
) where
- SMut : StorageMut<F, Dyn, U1>,
- S : Storage<F, Dyn, U1>
+ SMut: StorageMut<F, Dyn, U1>,
+ S: Storage<F, Dyn, U1>,
{
match b {
- Interior(c, (_, f)) | LeftBoundary(c, f) => { v[c] += (x[f] - x[c]) * α },
- RightBoundary(_c, _b) => { },
+ Interior(c, (_, f)) | LeftBoundary(c, f) => v[c] += (x[f] - x[c]) * α,
+ RightBoundary(_c, _b) => {}
}
}
@@ -118,23 +108,23 @@
}
}
-impl<F : Float + nalgebra::RealField> Discretisation<F> for BackwardNeumann {
+impl<F: Float + nalgebra::RealField> Discretisation<F> for BackwardNeumann {
type Opposite = ForwardDirichlet;
#[inline]
fn add_diff_mut<SMut, S>(
&self,
- v : &mut Matrix<F, Dyn, U1, SMut>,
- x : &Matrix<F, Dyn, U1, S>,
- α : F,
- b : DiscretisationOrInterior,
+ v: &mut Matrix<F, Dyn, U1, SMut>,
+ x: &Matrix<F, Dyn, U1, S>,
+ α: F,
+ b: DiscretisationOrInterior,
) where
- SMut : StorageMut<F, Dyn, U1>,
- S : Storage<F, Dyn, U1>
+ SMut: StorageMut<F, Dyn, U1>,
+ S: Storage<F, Dyn, U1>,
{
match b {
- Interior(c, (b, _)) | RightBoundary(c, b) => { v[c] += (x[c] - x[b]) * α },
- LeftBoundary(_c, _f) => { },
+ Interior(c, (b, _)) | RightBoundary(c, b) => v[c] += (x[c] - x[b]) * α,
+ LeftBoundary(_c, _f) => {}
}
}
@@ -144,24 +134,24 @@
}
}
-impl<F : Float + nalgebra::RealField> Discretisation<F> for BackwardDirichlet {
+impl<F: Float + nalgebra::RealField> Discretisation<F> for BackwardDirichlet {
type Opposite = ForwardNeumann;
#[inline]
fn add_diff_mut<SMut, S>(
&self,
- v : &mut Matrix<F, Dyn, U1, SMut>,
- x : &Matrix<F, Dyn, U1, S>,
- α : F,
- b : DiscretisationOrInterior,
+ v: &mut Matrix<F, Dyn, U1, SMut>,
+ x: &Matrix<F, Dyn, U1, S>,
+ α: F,
+ b: DiscretisationOrInterior,
) where
- SMut : StorageMut<F, Dyn, U1>,
- S : Storage<F, Dyn, U1>
+ SMut: StorageMut<F, Dyn, U1>,
+ S: Storage<F, Dyn, U1>,
{
match b {
- Interior(c, (b, _f)) => { v[c] += (x[c] - x[b]) * α },
- LeftBoundary(c, _f) => { v[c] += x[c] * α },
- RightBoundary(c, b) => { v[c] -= x[b] * α },
+ Interior(c, (b, _f)) => v[c] += (x[c] - x[b]) * α,
+ LeftBoundary(c, _f) => v[c] += x[c] * α,
+ RightBoundary(c, b) => v[c] -= x[b] * α,
}
}
@@ -171,24 +161,24 @@
}
}
-impl<F : Float + nalgebra::RealField> Discretisation<F> for ForwardDirichlet {
+impl<F: Float + nalgebra::RealField> Discretisation<F> for ForwardDirichlet {
type Opposite = BackwardNeumann;
#[inline]
fn add_diff_mut<SMut, S>(
&self,
- v : &mut Matrix<F, Dyn, U1, SMut>,
- x : &Matrix<F, Dyn, U1, S>,
- α : F,
- b : DiscretisationOrInterior,
+ v: &mut Matrix<F, Dyn, U1, SMut>,
+ x: &Matrix<F, Dyn, U1, S>,
+ α: F,
+ b: DiscretisationOrInterior,
) where
- SMut : StorageMut<F, Dyn, U1>,
- S : Storage<F, Dyn, U1>
+ SMut: StorageMut<F, Dyn, U1>,
+ S: Storage<F, Dyn, U1>,
{
match b {
- Interior(c, (_b, f)) => { v[c] += (x[f] - x[c]) * α },
- LeftBoundary(c, f) => { v[c] += x[f] * α },
- RightBoundary(c, _b) => { v[c] -= x[c] * α },
+ Interior(c, (_b, f)) => v[c] += (x[f] - x[c]) * α,
+ LeftBoundary(c, f) => v[c] += x[f] * α,
+ RightBoundary(c, _b) => v[c] -= x[c] * α,
}
}
@@ -198,30 +188,37 @@
}
}
-struct Iter<'a, const N : usize> {
+struct Iter<'a, const N: usize> {
/// Dimensions
- dims : &'a [usize; N],
+ dims: &'a [usize; N],
/// Dimension along which to calculate differences
- d : usize,
+ d: usize,
/// Stride along coordinate d
- d_stride : usize,
+ d_stride: usize,
/// Cartesian indices
- i : [usize; N],
+ i: [usize; N],
/// Linear index
- k : usize,
+ k: usize,
/// Maximal linear index
- len : usize
+ len: usize,
}
-impl<'a, const N : usize> Iter<'a, N> {
- fn new(dims : &'a [usize; N], d : usize) -> Self {
+impl<'a, const N: usize> Iter<'a, N> {
+ fn new(dims: &'a [usize; N], d: usize) -> Self {
let d_stride = dims[0..d].iter().product::<usize>();
let len = dims.iter().product::<usize>();
- Iter{ dims, d, d_stride, i : [0; N], k : 0, len }
+ Iter {
+ dims,
+ d,
+ d_stride,
+ i: [0; N],
+ k: 0,
+ len,
+ }
}
}
-impl<'a, const N : usize> Iterator for Iter<'a, N> {
+impl<'a, const N: usize> Iterator for Iter<'a, N> {
type Item = DiscretisationOrInterior;
fn next(&mut self) -> Option<Self::Item> {
let res = if self.k >= self.len {
@@ -243,7 +240,7 @@
for j in 0..N {
if self.i[j] + 1 < self.dims[j] {
self.i[j] += 1;
- break
+ break;
}
self.i[j] = 0
}
@@ -251,14 +248,13 @@
}
}
-impl<F, B, const N : usize> Mapping<DVector<F>>
-for Grad<F, B, N>
+impl<F, B, const N: usize> Mapping<DVector<F>> for Grad<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
type Codomain = DVector<F>;
- fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> {
+ fn apply<I: Instance<DVector<F>>>(&self, i: I) -> DVector<F> {
let mut y = DVector::zeros(N * self.len());
self.apply_add(&mut y, i);
y
@@ -266,15 +262,12 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<F, B, const N : usize> GEMV<F, DVector<F>>
-for Grad<F, B, N>
+impl<F, B, const N: usize> GEMV<F, DVector<F>> for Grad<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
- fn gemv<I : Instance<DVector<F>>>(
- &self, y : &mut DVector<F>, α : F, i : I, β : F
- ) {
+ fn gemv<I: Instance<DVector<F>>>(&self, y: &mut DVector<F>, α: F, i: I, β: F) {
if β == 0.0 {
y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0);
} else if β != 1.0 {
@@ -286,23 +279,22 @@
i.eval(|x| {
assert_eq!(x.len(), m);
for d in 0..N {
- let mut v = y.generic_view_mut((d*m, 0), (Dyn(m), U1));
+ let mut v = y.generic_view_mut((d * m, 0), (Dyn(m), U1));
for b in Iter::new(&self.dims, d) {
- self.discretisation.add_diff_mut(&mut v, x, α/h, b)
+ self.discretisation.add_diff_mut(&mut v, x, α / h, b)
}
}
})
}
}
-impl<F, B, const N : usize> Mapping<DVector<F>>
-for Div<F, B, N>
+impl<F, B, const N: usize> Mapping<DVector<F>> for Div<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
type Codomain = DVector<F>;
- fn apply<I : Instance<DVector<F>>>(&self, i : I) -> DVector<F> {
+ fn apply<I: Instance<DVector<F>>>(&self, i: I) -> DVector<F> {
let mut y = DVector::zeros(self.len());
self.apply_add(&mut y, i);
y
@@ -310,15 +302,12 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<F, B, const N : usize> GEMV<F, DVector<F>>
-for Div<F, B, N>
+impl<F, B, const N: usize> GEMV<F, DVector<F>> for Div<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
- fn gemv<I : Instance<DVector<F>>>(
- &self, y : &mut DVector<F>, α : F, i : I, β : F
- ) {
+ fn gemv<I: Instance<DVector<F>>>(&self, y: &mut DVector<F>, α: F, i: I, β: F) {
if β == 0.0 {
y.as_mut_slice().iter_mut().for_each(|x| *x = 0.0);
} else if β != 1.0 {
@@ -330,26 +319,51 @@
i.eval(|x| {
assert_eq!(x.len(), N * m);
for d in 0..N {
- let v = x.generic_view((d*m, 0), (Dyn(m), U1));
+ let v = x.generic_view((d * m, 0), (Dyn(m), U1));
for b in Iter::new(&self.dims, d) {
- self.discretisation.add_diff_mut(y, &v, α/h, b)
+ self.discretisation.add_diff_mut(y, &v, α / h, b)
}
}
})
}
}
-impl<F, B, const N : usize> Grad<F, B, N>
+impl<F, B, const N: usize> Grad<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
/// Creates a new discrete gradient operator for the vector `u`, verifying dimensions.
- pub fn new_for(u : &DVector<F>, h : F, dims : [usize; N], discretisation : B)
- -> Option<Self>
- {
+ pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> {
if u.len() == dims.iter().product::<usize>() {
- Some(Grad { dims, h, discretisation } )
+ Some(Grad {
+ dims,
+ h,
+ discretisation,
+ })
+ } else {
+ None
+ }
+ }
+
+ fn len(&self) -> usize {
+ self.dims.iter().product::<usize>()
+ }
+}
+
+impl<F, B, const N: usize> Div<F, B, N>
+where
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
+{
+ /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions.
+ pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> {
+ if u.len() == dims.iter().product::<usize>() * N {
+ Some(Div {
+ dims,
+ h,
+ discretisation,
+ })
} else {
None
}
@@ -360,108 +374,84 @@
}
}
-
-impl<F, B, const N : usize> Div<F, B, N>
+impl<F, B, const N: usize> Linear<DVector<F>> for Grad<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
- /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions.
- pub fn new_for(u : &DVector<F>, h : F, dims : [usize; N], discretisation : B)
- -> Option<Self>
- {
- if u.len() == dims.iter().product::<usize>() * N {
- Some(Div { dims, h, discretisation } )
- } else {
- None
- }
- }
-
- fn len(&self) -> usize {
- self.dims.iter().product::<usize>()
- }
}
-impl<F, B, const N : usize> Linear<DVector<F>>
-for Grad<F, B, N>
+impl<F, B, const N: usize> Linear<DVector<F>> for Div<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
}
-impl<F, B, const N : usize> Linear<DVector<F>>
-for Div<F, B, N>
+impl<F, B, const N: usize> BoundedLinear<DVector<F>, L2, L2, F> for Grad<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
+ DVector<F>: Norm<F, L2>,
{
-}
-
-impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F>
-for Grad<F, B, N>
-where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
- DVector<F> : Norm<F, L2>,
-{
- fn opnorm_bound(&self, _ : L2, _ : L2) -> F {
+ fn opnorm_bound(&self, _: L2, _: L2) -> DynResult<F> {
// Fuck nalgebra.
- self.discretisation.opnorm_bound(num_traits::Float::abs(self.h))
+ self.discretisation
+ .opnorm_bound(num_traits::Float::abs(self.h))
}
}
-
-impl<F, B, const N : usize> BoundedLinear<DVector<F>, L2, L2, F>
-for Div<F, B, N>
+impl<F, B, const N: usize> BoundedLinear<DVector<F>, L2, L2, F> for Div<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
- DVector<F> : Norm<F, L2>,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
+ DVector<F>: Norm<F, L2>,
{
- fn opnorm_bound(&self, _ : L2, _ : L2) -> F {
+ fn opnorm_bound(&self, _: L2, _: L2) -> DynResult<F> {
// Fuck nalgebra.
- self.discretisation.opnorm_bound(num_traits::Float::abs(self.h))
+ self.discretisation
+ .opnorm_bound(num_traits::Float::abs(self.h))
}
}
-impl<F, B, const N : usize>
-Adjointable<DVector<F>, DVector<F>>
-for Grad<F, B, N>
+impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Grad<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
type AdjointCodomain = DVector<F>;
- type Adjoint<'a> = Div<F, B::Opposite, N> where Self : 'a;
+ type Adjoint<'a>
+ = Div<F, B::Opposite, N>
+ where
+ Self: 'a;
/// Form the adjoint operator of `self`.
fn adjoint(&self) -> Self::Adjoint<'_> {
Div {
- dims : self.dims,
- h : -self.h,
- discretisation : self.discretisation.opposite(),
+ dims: self.dims,
+ h: -self.h,
+ discretisation: self.discretisation.opposite(),
}
}
}
-
-impl<F, B, const N : usize>
-Adjointable<DVector<F>, DVector<F>>
-for Div<F, B, N>
+impl<F, B, const N: usize> Adjointable<DVector<F>, DVector<F>> for Div<F, B, N>
where
- B : Discretisation<F>,
- F : Float + nalgebra::RealField,
+ B: Discretisation<F>,
+ F: Float + nalgebra::RealField,
{
type AdjointCodomain = DVector<F>;
- type Adjoint<'a> = Grad<F, B::Opposite, N> where Self : 'a;
+ type Adjoint<'a>
+ = Grad<F, B::Opposite, N>
+ where
+ Self: 'a;
/// Form the adjoint operator of `self`.
fn adjoint(&self) -> Self::Adjoint<'_> {
Grad {
- dims : self.dims,
- h : -self.h,
- discretisation : self.discretisation.opposite(),
+ dims: self.dims,
+ h: -self.h,
+ discretisation: self.discretisation.opposite(),
}
}
}
@@ -472,8 +462,8 @@
#[test]
fn grad_adjoint() {
- let im = DVector::from( (0..9).map(|t| t as f64).collect::<Vec<_>>());
- let v = DVector::from( (0..18).map(|t| t as f64).collect::<Vec<_>>());
+ let im = DVector::from((0..9).map(|t| t as f64).collect::<Vec<_>>());
+ let v = DVector::from((0..18).map(|t| t as f64).collect::<Vec<_>>());
let grad = Grad::new_for(&im, 1.0, [3, 3], ForwardNeumann).unwrap();
let a = grad.apply(&im).dot(&v);
@@ -484,6 +474,5 @@
let a = grad.apply(&im).dot(&v);
let b = grad.adjoint().apply(&v).dot(&im);
assert_eq!(a, b);
-
}
}