src/forward_model.rs

/*!
Forward models from discrete measures to observations.
*/

use crate::dataterm::QuadraticDataTerm;
use crate::measures::{Radon, RNDM};
use crate::types::*;
use alg_tools::error::{DynError, DynResult};
use alg_tools::euclidean::{ClosedEuclidean, Euclidean};
pub use alg_tools::linops::*;
use alg_tools::norms::{Norm, NormExponent, L2};
use serde::{Deserialize, Serialize};

pub mod bias;
pub mod sensor_grid;

/// `ForwardeModel`s are bounded preadjointable linear operators  $A ∈ 𝕃(𝒵(Ω); E)$
/// where $𝒵(Ω) ⊂ ℳ(Ω)$ is the space of sums of delta measures, presented by
/// [`crate::measures::DiscreteMeasure`], and $E$ is a [`Euclidean`] space.
pub trait ForwardModel<Domain: Space, F: Float = f64, E: NormExponent = Radon>:
    BoundedLinear<Domain, E, L2, F, Codomain = Self::Observable>
    + GEMV<F, Domain, Self::Observable>
    + Preadjointable<Domain, Self::Observable>
where
    Domain: Norm<E, F>,
{
    /// 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: ClosedEuclidean<F> + Clone;

    /// 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;
}

/// Guess for [`BoundedCurvature`] calculations.
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
pub enum BoundedCurvatureGuess {
    /// No iterate $μ^k$ is worse than $μ=0$.
    BetterThanZero,
}

/// Curvature error control.
pub trait BoundedCurvature<F: Float = f64> {
    /// Returns an estimate of $(ℓ_F, ℓ_∇v, Θ²)$ or individual errors for each.
    fn curvature_bound_components(
        &self,
        guess: BoundedCurvatureGuess,
    ) -> (DynResult<F>, DynResult<F>, DynResult<F>);
}

/// Curvature error control: helper Lipschitz-like bounds for the quadratic dataterms
///  $F(μ) = \frac{1}{2}\|Aμ-b\|^2$.
pub trait BasicCurvatureBoundEstimates<F: Float = f64> {
    /// Returns $(ℓ_F, ℓ_{∇v}^0, Θ²)$ or individual errors for each.
    fn basic_curvature_bound_components(&self) -> (DynResult<F>, DynResult<F>, DynResult<F>);
}

impl<F, A, Z, const N: usize> BoundedCurvature<F> for QuadraticDataTerm<F, RNDM<N, F>, A>
where
    F: Float,
    Z: Clone + Space + Euclidean<F>,
    A: Mapping<RNDM<N, F>, Codomain = Z>,
    A: BasicCurvatureBoundEstimates<F>,
{
    fn curvature_bound_components(
        &self,
        guess: BoundedCurvatureGuess,
    ) -> (DynResult<F>, DynResult<F>, DynResult<F>) {
        match guess {
            BoundedCurvatureGuess::BetterThanZero => {
                let opA = self.operator();
                let b = self.data();
                let (ℓ_F, ℓ_gradv_0, θ2) = opA.basic_curvature_bound_components();
                (ℓ_F, ℓ_gradv_0.map(|l| l * b.norm2()), θ2)
            }
        }
    }
}