pointsource_algs: comparison src/forward

-:9738b51d90d7
+:4f468d35fa29
 /*!
 Forward models from discrete measures to observations.
 */
-use alg_tools::error::DynError;
+use crate::dataterm::QuadraticDataTerm;
-use alg_tools::euclidean::Euclidean;
+use crate::measures::{Radon, RNDM};
-use alg_tools::instance::Instance;
+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};
-use crate::measures::Radon;
-use crate::types::*;
 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
 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
-for<'a> Self::Observable: Instance<Self::Observable>,
+Domain: Norm<E, F>,
-Domain: Norm<F, E>,
 {
 /// 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<F, Output = Self::Observable> + AXPY<F> + Space + Clone;
+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;
 }
-/// Trait for operators $A$ for which $A_*A$ is bounded by some other operator.
+/// Guess for [`BoundedCurvature`] calculations.
-pub trait AdjointProductBoundedBy<Domain: Space, D>: Linear<Domain> {
+#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
-type FloatType: Float;
+pub enum BoundedCurvatureGuess {
-/// Return $L$ such that $A_*A ≤ LD$.
+/// No iterate $μ^k$ is worse than $μ=0$.
-fn adjoint_product_bound(&self, other: &D) -> Option<Self::FloatType>;
+BetterThanZero,
 }
-/// Trait for operators $A$ for which $A_*A$ is bounded by a diagonal operator.
+/// Curvature error control: helper bounds for (4.2d), (5.2a), (5.2b), (5.15a), and (5.16a).
-pub trait AdjointProductPairBoundedBy<Domain: Space, D1, D2>: Linear<Domain> {
+///
-type FloatType: Float;
+/// Based on Lemma 5.11 and Example 5.12, the helper bound for (5.15a) and (5.16a) is (3.8).
-/// Return $(L, L_z)$ such that $A_*A ≤ (L_1 D_1, L_2 D_2)$.
+/// Thus, subject to `guess` being correct, returns factor $(ℓ_F, Θ²)$ such that
-fn adjoint_product_pair_bound(
+/// $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,
-other1: &D1,
+guess: BoundedCurvatureGuess,
-other_2: &D2,
+) -> (DynResult<F>, DynResult<F>);
-) -> Option<(Self::FloatType, Self::FloatType)>;
 }
-/*
+/// Curvature error control: helper bounds for (4.2d), (5.2a), (5.2b), (5.15a), and (5.16a)
-/// Trait for [`ForwardModel`]s whose preadjoint has Lipschitz values.
+/// for quadratic dataterms $F(μ) = \frac{1}{2}\|Aμ-b\|^2$.
-pub trait LipschitzValues {
+///
-type FloatType : Float;
+/// This trait is to be implemented by the [`Linear`] operator $A$, in the basic from
-/// Return (if one exists) a factor $L$ such that $A_*z$ is $L$-Lipschitz for all
+/// [`RNDM<N, F>`] to a an Euclidean space.
-/// $z$ in the unit ball.
+/// It is used by implementations of [`BoundedCurvature`] for $F$.
-fn value_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
+///
-None
+/// 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>);
+}
-/// Return (if one exists) a factor $L$ such that $∇A_*z$ is $L$-Lipschitz for all
+impl<F, A, Z, const N: usize> BoundedCurvature<F> for QuadraticDataTerm<F, RNDM<N, F>, A>
-/// $z$ in the unit ball.
+where
-fn value_diff_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
+F: Float,
-None
+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)
+}
+}
 }
 }
-*/
-/// Trait for [`ForwardModel`]s that satisfy bounds on curvature.
-pub trait BoundedCurvature {
-type FloatType: Float;
-/// Returns factor $ℓ_F$ and $ℓ_r$ such that
-/// $B_{F'(μ)} dγ ≤ ℓ_F c_2$ and $⟨F'(μ)+F'(μ+Δ)|Δ⟩ ≤ ℓ_r|γ|(c_2)$,
-/// where $Δ=(π_♯^1-π_♯^0)γ$.
-fn curvature_bound_components(&self) -> (Option<Self::FloatType>, Option<Self::FloatType>);
-}

comparison: src/forward_model.rs

src/forward_model.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model.rs

src/forward_model.rs