Mercurial > repos > pointsource_algs / file revision

/*!
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: helper bounds for (4.2d), (5.2a), (5.2b), (5.15a), and (5.16a).
///
/// Based on Lemma 5.11 and Example 5.12, the helper bound for (5.15a) and (5.16a) is (3.8).
/// Thus, subject to `guess` being correct, returns factor $(ℓ_F, Θ²)$ such that
/// $B_{F'(μ)} dγ ≤ ℓ_F c_2$ and $⟨F'(μ+Δ)-F'(μ)|Δ⟩ ≤ Θ²|γ|(c_2)$, where $Δ=(π_♯^1-π_♯^0)γ$.
///
/// This trait is supposed to be implemented by the data term $F$, in the basic case a
/// [`Mapping`] from [`RNDM<N, F>`] to  a [`Float`] `F`.
/// The generic implementation for operators that satisfy [`BasicCurvatureBoundEstimates`]
/// uses Remark 5.15 and Example 5.16 for (4.2d) and (5.2a), and (5.2b);
/// and Lemma 3.8 for (3.8).
pub trait BoundedCurvature<F: Float = f64> {
    /// Returns $(ℓ_F, Θ²)$ or individual errors for each.
    fn curvature_bound_components(
        &self,
        guess: BoundedCurvatureGuess,
    ) -> (DynResult<F>, DynResult<F>);
}

/// Curvature error control: helper bounds for (4.2d), (5.2a), (5.2b), (5.15a), and (5.16a)
/// for quadratic dataterms $F(μ) = \frac{1}{2}\|Aμ-b\|^2$.
///
/// This trait is to be implemented by the [`Linear`] operator $A$, in the basic from
/// [`RNDM<N, F>`] to a an Euclidean space.
/// It is used by implementations of [`BoundedCurvature`] for $F$.
///
/// Based on Lemma 5.11 and Example 5.12, the helper bound for (5.15a) and (5.16a) is (3.8).
/// This trait provides the factor $θ²$ of (3.8) as determined by Lemma 3.8.
/// To aid in calculating (4.2d), (5.2a), (5.2b), motivated by Example 5.16, it also provides
/// $ℓ_F^0$ such that $∇v^k$ $ℓ_F^0 \|Aμ-b\|$-Lipschitz. Here $v^k := F'(∪^k)$.
pub trait BasicCurvatureBoundEstimates<F: Float = f64> {
    /// Returns $(ℓ_F^0, Θ²)$ or individual errors for each.
    fn basic_curvature_bound_components(&self) -> (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>) {
        match guess {
            BoundedCurvatureGuess::BetterThanZero => {
                let opA = self.operator();
                let b = self.data();
                let (ℓ_F0, θ2) = opA.basic_curvature_bound_components();
                (ℓ_F0.map(|l| l * b.norm2()), θ2)
            }
        }
    }
}