diff -r 000000000000 -r eb3c7813b67a src/forward_model.rs --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/forward_model.rs Thu Dec 01 23:07:35 2022 +0200 @@ -0,0 +1,678 @@ +/*! +Forward models from discrete measures to observations. +*/ + +use numeric_literals::replace_float_literals; +use nalgebra::base::{ + DMatrix, + DVector +}; +use std::iter::Zip; +use std::ops::RangeFrom; +use std::marker::PhantomData; + +pub use alg_tools::linops::*; +use alg_tools::euclidean::Euclidean; +use alg_tools::norms::{ + L1, Linfinity, Norm +}; +use alg_tools::bisection_tree::*; +use alg_tools::mapping::RealMapping; +use alg_tools::lingrid::*; +use alg_tools::iter::{MapX, Mappable}; +use alg_tools::nalgebra_support::ToNalgebraRealField; +use alg_tools::tabledump::write_csv; +use alg_tools::error::DynError; + +use crate::types::*; +use crate::measures::*; +use crate::seminorms::{ + Lipschitz, + ConvolutionOp, + SimpleConvolutionKernel, +}; +use crate::kernels::{ + Convolution, + AutoConvolution, + BoundedBy, +}; + +pub type RNDM = DiscreteMeasure, F>; + +/// `ForwardeModel`s are bounded preadjointable linear operators $A ∈ 𝕃(𝒵(Ω); E)$ +/// where $𝒵(Ω) ⊂ ℳ(Ω)$ is the space of sums of delta measures, presented by +/// [`DiscreteMeasure`], and $E$ is a [`Euclidean`] space. +pub trait ForwardModel +: BoundedLinear, Codomain=Self::Observable, FloatType=F> ++ GEMV, Self::Observable> ++ Linear, Codomain=Self::Observable> ++ Preadjointable, Self::Observable> { + /// The codomain or value space (of “observables”) for this operator. + /// It is assumed to be a [`Euclidean`] space, and therefore also (identified with) + /// the domain of the preadjoint. + type Observable : Euclidean + + AXPY + + Clone; + + /// Return A_*A and A_* b + fn findim_quadratic_model( + &self, + μ : &DiscreteMeasure, + b : &Self::Observable + ) -> (DMatrix, DVector); + + /// Write an observable into a file. + fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError; + + /// Returns a zero observable + fn zero_observable(&self) -> Self::Observable; + + /// Returns an empty (uninitialised) observable. + /// + /// This is used as a placeholder for temporary [`std::mem::replace`] move operations. + fn empty_observable(&self) -> Self::Observable; +} + +pub type ShiftedSensor = Shift, F, N>; + +/// Trait for physical convolution models. Has blanket implementation for all cases. +pub trait Spread +: 'static + Clone + Support + RealMapping + Bounded {} + +impl Spread for T +where F : Float, + T : 'static + Clone + Support + Bounded + RealMapping {} + +/// Trait for compactly supported sensors. Has blanket implementation for all cases. +pub trait Sensor : Spread + Norm + Norm {} + +impl Sensor for T +where F : Float, + T : Spread + Norm + Norm {} + + +pub trait SensorGridBT : +Clone + BTImpl> +where F : Float, + S : Sensor, + P : Spread {} + +impl +SensorGridBT +for T +where T : Clone + BTImpl>, + F : Float, + S : Sensor, + P : Spread {} + +// We need type alias bounds to access associated types +#[allow(type_alias_bounds)] +type SensorGridBTFN, const N : usize> += BTFN, BT, N>; + +/// Sensor grid forward model +#[derive(Clone)] +pub struct SensorGrid +where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread, + BT : SensorGridBT, { + domain : Cube, + sensor_count : [usize; N], + sensor : S, + spread : P, + base_sensor : Convolution, + bt : BT, +} + +impl SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + + pub fn new( + domain : Cube, + sensor_count : [usize; N], + sensor : S, + spread : P, + depth : BT::Depth + ) -> Self { + let base_sensor = Convolution(sensor.clone(), spread.clone()); + let bt = BT::new(domain, depth); + let mut sensorgrid = SensorGrid { + domain, + sensor_count, + sensor, + spread, + base_sensor, + bt, + }; + + for (x, id) in sensorgrid.grid().into_iter().zip(0usize..) { + let s = sensorgrid.shifted_sensor(x); + sensorgrid.bt.insert(id, &s); + } + + sensorgrid + } + + pub fn grid(&self) -> LinGrid { + lingrid_centered(&self.domain, &self.sensor_count) + } + + pub fn n_sensors(&self) -> usize { + self.sensor_count.iter().product() + } + + #[inline] + fn shifted_sensor(&self, x : Loc) -> ShiftedSensor { + self.base_sensor.clone().shift(x) + } + + #[inline] + fn _zero_observable(&self) -> DVector { + DVector::zeros(self.n_sensors()) + } +} + +impl Apply> for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + + type Output = DVector; + + #[inline] + fn apply(&self, μ : RNDM) -> DVector { + self.apply(&μ) + } +} + +impl<'a, F, S, P, BT, const N : usize> Apply<&'a RNDM> for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + + type Output = DVector; + + fn apply(&self, μ : &'a RNDM) -> DVector { + let mut res = self._zero_observable(); + self.apply_add(&mut res, μ); + res + } +} + +impl Linear> for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + type Codomain = DVector; +} + + +#[replace_float_literals(F::cast_from(literal))] +impl GEMV, DVector> for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + + fn gemv(&self, y : &mut DVector, α : F, μ : &RNDM, β : F) { + let grid = self.grid(); + if β == 0.0 { + y.fill(0.0) + } else if β != 1.0 { + *y *= β; // Need to multiply first, as we have to be able to add to y. + } + if α == 1.0 { + self.apply_add(y, μ) + } else { + for δ in μ.iter_spikes() { + for &d in self.bt.iter_at(&δ.x) { + let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); + y[d] += sensor.apply(&δ.x) * (α * δ.α); + } + } + } + } + + fn apply_add(&self, y : &mut DVector, μ : &RNDM) { + let grid = self.grid(); + for δ in μ.iter_spikes() { + for &d in self.bt.iter_at(&δ.x) { + let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); + y[d] += sensor.apply(&δ.x) * δ.α; + } + } + } + +} + +impl Apply, F>> +for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + + type Output = DVector; + + #[inline] + fn apply(&self, δ : DeltaMeasure, F>) -> DVector { + self.apply(&δ) + } +} + +impl<'a, F, S, P, BT, const N : usize> Apply<&'a DeltaMeasure, F>> +for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + + type Output = DVector; + + fn apply(&self, δ : &DeltaMeasure, F>) -> DVector { + let mut res = DVector::zeros(self.n_sensors()); + let grid = self.grid(); + for &d in self.bt.iter_at(&δ.x) { + let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); + res[d] += sensor.apply(&δ.x) * δ.α; + } + res + } +} + +impl Linear, F>> for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + type Codomain = DVector; +} + +impl BoundedLinear> for SensorGrid +where F : Float, + BT : SensorGridBT>, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis { + type FloatType = F; + + /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped + /// with the Radon norm, and $ℝ^n$ with the Euclidean norm. + fn opnorm_bound(&self) -> F { + // With {x_i}_{i=1}^n the grid centres and φ the kernel, we have + // |Aμ|_2 = sup_{|z|_2 ≤ 1} ⟨z,Αμ⟩ = sup_{|z|_2 ≤ 1} ⟨A^*z|μ⟩ + // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ + // = sup_{|z|_2 ≤ 1} |∑ φ(· - x_i)z_i|_∞ |μ|_ℳ + // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ ∑ |z_i| |μ|_ℳ + // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ √n |z|_2 |μ|_ℳ + // = |φ|_∞ √n |μ|_ℳ. + // Hence + let n = F::cast_from(self.n_sensors()); + self.base_sensor.bounds().uniform() * n.sqrt() + } +} + +type SensorGridPreadjoint<'a, A, F, const N : usize> = PreadjointHelper<'a, A, RNDM>; + + +impl +Preadjointable, DVector> +for SensorGrid +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis, + Weighted, F> : LocalAnalysis { + type PreadjointCodomain = BTFN, BT, N>; + type Preadjoint<'a> = SensorGridPreadjoint<'a, Self, F, N> where Self : 'a; + + fn preadjoint(&self) -> Self::Preadjoint<'_> { + PreadjointHelper::new(self) + } +} + +#[derive(Clone,Debug)] +pub struct SensorGridSupportGenerator +where F : Float, + S : Sensor, + P : Spread { + base_sensor : Convolution, + grid : LinGrid, + weights : DVector +} + +impl SensorGridSupportGenerator +where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread { + + #[inline] + fn construct_sensor(&self, id : usize, w : F) -> Weighted, F> { + let x = self.grid.entry_linear_unchecked(id); + self.base_sensor.clone().shift(x).weigh(w) + } + + #[inline] + fn construct_sensor_and_id<'a>(&'a self, (id, w) : (usize, &'a F)) + -> (usize, Weighted, F>) { + (id.into(), self.construct_sensor(id, *w)) + } +} + +impl SupportGenerator +for SensorGridSupportGenerator +where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread { + type Id = usize; + type SupportType = Weighted, F>; + type AllDataIter<'a> = MapX<'a, Zip, + std::slice::Iter<'a, F>>, + Self, + (Self::Id, Self::SupportType)> + where Self : 'a; + + #[inline] + fn support_for(&self, d : Self::Id) -> Self::SupportType { + self.construct_sensor(d, self.weights[d]) + } + + #[inline] + fn support_count(&self) -> usize { + self.weights.len() + } + + #[inline] + fn all_data(&self) -> Self::AllDataIter<'_> { + (0..).zip(self.weights.as_slice().iter()).mapX(self, Self::construct_sensor_and_id) + } +} + +/// Helper structure for constructing preadjoints of `S` where `S : Linear`. +/// [`Linear`] needs to be implemented for each instance, but [`Adjointable`] +/// and [`BoundedLinear`] have blanket implementations. +#[derive(Clone,Debug)] +pub struct PreadjointHelper<'a, S : 'a, X> { + forward_op : &'a S, + _domain : PhantomData +} + +impl<'a, S : 'a, X> PreadjointHelper<'a, S, X> { + pub fn new(forward_op : &'a S) -> Self { + PreadjointHelper { forward_op, _domain: PhantomData } + } +} + +impl<'a, X, Ypre, S> Adjointable +for PreadjointHelper<'a, S, X> +where Self : Linear, + S : Clone + Linear { + type AdjointCodomain = S::Codomain; + type Adjoint<'b> = S where Self : 'b; + fn adjoint(&self) -> Self::Adjoint<'_> { + self.forward_op.clone() + } +} + +impl<'a, X, Ypre, S> BoundedLinear +for PreadjointHelper<'a, S, X> +where Self : Linear, + S : 'a + Clone + BoundedLinear { + type FloatType = S::FloatType; + fn opnorm_bound(&self) -> Self::FloatType { + self.forward_op.opnorm_bound() + } +} + + +impl<'a, 'b, F, S, P, BT, const N : usize> Apply<&'b DVector> +for PreadjointHelper<'a, SensorGrid, RNDM> +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis, + Weighted, F> : LocalAnalysis { + + type Output = SensorGridBTFN; + + fn apply(&self, x : &'b DVector) -> Self::Output { + self.apply(x.clone()) + } +} + +impl<'a, F, S, P, BT, const N : usize> Apply> +for PreadjointHelper<'a, SensorGrid, RNDM> +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis, + Weighted, F> : LocalAnalysis { + + type Output = SensorGridBTFN; + + fn apply(&self, x : DVector) -> Self::Output { + let fwd = &self.forward_op; + let generator = SensorGridSupportGenerator{ + base_sensor : fwd.base_sensor.clone(), + grid : fwd.grid(), + weights : x + }; + BTFN::new_refresh(&fwd.bt, generator) + } +} + +impl<'a, F, S, P, BT, const N : usize> Linear> +for PreadjointHelper<'a, SensorGrid, RNDM> +where F : Float, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis, + Weighted, F> : LocalAnalysis { + + type Codomain = SensorGridBTFN; +} + +impl ForwardModel, F> +for SensorGrid +where F : Float + ToNalgebraRealField + nalgebra::RealField, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + ShiftedSensor : LocalAnalysis, + Weighted, F> : LocalAnalysis { + type Observable = DVector; + + fn findim_quadratic_model( + &self, + μ : &DiscreteMeasure, F>, + b : &Self::Observable + ) -> (DMatrix, DVector) { + assert_eq!(b.len(), self.n_sensors()); + let mut mA = DMatrix::zeros(self.n_sensors(), μ.len()); + let grid = self.grid(); + for (mut mAcol, δ) in mA.column_iter_mut().zip(μ.iter_spikes()) { + for &d in self.bt.iter_at(&δ.x) { + let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); + mAcol[d] += sensor.apply(&δ.x); + } + } + let mAt = mA.transpose(); + (&mAt * mA, &mAt * b) + } + + fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError { + let it = self.grid().into_iter().zip(b.iter()).map(|(x, &v)| (x, v)); + write_csv(it, prefix + ".txt") + } + + #[inline] + fn zero_observable(&self) -> Self::Observable { + self._zero_observable() + } + + #[inline] + fn empty_observable(&self) -> Self::Observable { + DVector::zeros(0) + } + +} + +/// Implements the calculation a factor $L$ such that $A_*A ≤ L 𝒟$ for $A$ the forward model +/// and $𝒟$ a seminorm of suitable form. +/// +/// **This assumes (but does not check) that the sensors are not overlapping.** +#[replace_float_literals(F::cast_from(literal))] +impl Lipschitz> +for SensorGrid +where F : Float + nalgebra::RealField + ToNalgebraRealField, + BT : SensorGridBT, + S : Sensor, + P : Spread, + Convolution : Spread, + K : SimpleConvolutionKernel, + AutoConvolution

: BoundedBy { + + type FloatType = F; + + fn lipschitz_factor(&self, seminorm : &ConvolutionOp) -> Option { + // Sensors should not take on negative values to allow + // A_*A to be upper bounded by a simple convolution of `spread`. + if self.sensor.bounds().lower() < 0.0 { + return None + } + + // Calculate the factor $L_1$ for betwee $ℱ[ψ * ψ] ≤ L_1 ℱ[ρ]$ for $ψ$ the base spread + // and $ρ$ the kernel of the seminorm. + let l1 = AutoConvolution(self.spread.clone()).bounding_factor(seminorm.kernel())?; + + // Calculate the factor for transitioning from $A_*A$ to `AutoConvolution

`, where A + // consists of several `Convolution` for the physical model `P` and the sensor `S`. + let l0 = self.sensor.norm(Linfinity) * self.sensor.norm(L1); + + // The final transition factor is: + Some(l0 * l1) + } +} + +macro_rules! make_sensorgridsupportgenerator_scalarop_rhs { + ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { + impl + std::ops::$trait_assign + for SensorGridSupportGenerator + where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread { + fn $fn_assign(&mut self, t : F) { + self.weights.$fn_assign(t); + } + } + + impl + std::ops::$trait + for SensorGridSupportGenerator + where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread { + type Output = SensorGridSupportGenerator; + fn $fn(mut self, t : F) -> Self::Output { + std::ops::$trait_assign::$fn_assign(&mut self.weights, t); + self + } + } + + impl<'a, F, S, P, const N : usize> + std::ops::$trait + for &'a SensorGridSupportGenerator + where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread { + type Output = SensorGridSupportGenerator; + fn $fn(self, t : F) -> Self::Output { + SensorGridSupportGenerator{ + base_sensor : self.base_sensor.clone(), + grid : self.grid, + weights : (&self.weights).$fn(t) + } + } + } + } +} + +make_sensorgridsupportgenerator_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); +make_sensorgridsupportgenerator_scalarop_rhs!(Div, div, DivAssign, div_assign); + +macro_rules! make_sensorgridsupportgenerator_unaryop { + ($trait:ident, $fn:ident) => { + impl + std::ops::$trait + for SensorGridSupportGenerator + where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread { + type Output = SensorGridSupportGenerator; + fn $fn(mut self) -> Self::Output { + self.weights = self.weights.$fn(); + self + } + } + + impl<'a, F, S, P, const N : usize> + std::ops::$trait + for &'a SensorGridSupportGenerator + where F : Float, + S : Sensor, + P : Spread, + Convolution : Spread { + type Output = SensorGridSupportGenerator; + fn $fn(self) -> Self::Output { + SensorGridSupportGenerator{ + base_sensor : self.base_sensor.clone(), + grid : self.grid, + weights : (&self.weights).$fn() + } + } + } + } +} + +make_sensorgridsupportgenerator_unaryop!(Neg, neg);