Mercurial > repos > pointsource_algs / file diff

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/forward_model/bias.rs	Tue Dec 31 09:25:45 2024 -0500
@@ -0,0 +1,110 @@
+/*!
+Simple parametric forward model.
+ */
+
+use numeric_literals::replace_float_literals;
+use alg_tools::types::{Float, ClosedAdd};
+use alg_tools::mapping::Space;
+use alg_tools::direct_product::Pair;
+use alg_tools::linops::{Linear, RowOp, ColOp, IdOp, ZeroOp, AXPY};
+use alg_tools::error::DynError;
+use alg_tools::norms::{L2, Norm, PairNorm, NormExponent};
+use crate::types::L2Squared;
+use crate::measures::RNDM;
+use super::{ForwardModel, AdjointProductBoundedBy, AdjointProductPairBoundedBy, LipschitzValues};
+use crate::transport::TransportLipschitz;
+
+impl<Domain, F, A, E> ForwardModel<Pair<Domain, A::Observable>, F, PairNorm<E, L2, L2>>
+for RowOp<A, IdOp<A::Observable>>
+where
+    E : NormExponent,
+    Domain : Space + Norm<F, E>,
+    F : Float,
+    A::Observable : ClosedAdd + Norm<F, L2> + 'static,
+    A : ForwardModel<Domain, F, E> + 'static
+{
+    type Observable = A::Observable;
+
+    fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError {
+        self.0.write_observable(b, prefix)
+    }
+
+    /// Returns a zero observable
+    fn zero_observable(&self) -> Self::Observable {
+        self.0.zero_observable()
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<Domain, F, A, D, Z> AdjointProductPairBoundedBy<Pair<Domain, Z>, D, IdOp<Z>>
+for RowOp<A, IdOp<Z>>
+where
+    Domain : Space,
+    F : Float,
+    Z : Clone + Space + ClosedAdd,
+    A : AdjointProductBoundedBy<Domain, D, FloatType=F, Codomain = Z>,
+    D : Linear<Domain>,
+    A::Codomain : ClosedAdd,
+{
+    type FloatType = F;
+
+    fn adjoint_product_pair_bound(&self, d : &D, _ : &IdOp<Z>) -> Option<(F, F)> {
+        self.0.adjoint_product_bound(d).map(|l_0| {
+            // [A_*; B_*][A, B] = [A_*A, A_* B; B_* A, B_* B] ≤ diag(2A_*A, 2B_*B)
+            // ≤ diag(2l_A𝒟_A, 2l_B𝒟_B), where now 𝒟_B=Id and l_B=1.
+            (2.0 * l_0, 2.0)
+        })
+    }
+}
+
+/// This `impl` is bit of an abuse as the codomain of `Apre` is a [`Pair`] of a measure predual,
+/// to which this `impl` applies, and another space.
+impl<F, Apre, Z> LipschitzValues
+for ColOp<Apre, IdOp<Z>>
+where
+    F : Float,
+    Z : Clone + Space + ClosedAdd,
+    Apre : LipschitzValues<FloatType = F>,
+{
+    type FloatType = F;
+    /// Return (if one exists) a factor $L$ such that $A_*z$ is $L$-Lipschitz for all
+    /// $z$ in the unit ball.
+    fn value_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
+        self.0.value_unit_lipschitz_factor()
+    }
+
+    /// Return (if one exists) a factor $L$ such that $∇A_*z$ is $L$-Lipschitz for all
+    /// $z$ in the unit ball.
+    fn value_diff_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
+        self.0.value_diff_unit_lipschitz_factor()
+    }
+}
+
+
+
+impl<'a, F : Float, Y : Space, XD, const N : usize> TransportLipschitz<L2Squared> for
+ZeroOp<'a, RNDM<F, N>, XD, Y, F> {
+    type FloatType = F;
+
+    fn transport_lipschitz_factor(&self, _ : L2Squared) -> Self::FloatType {
+        F::ZERO
+    }
+}
+
+
+/// TODO: should assume `D` to be positive semi-definite and self-adjoint.
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F, D, XD, Y, const N : usize> AdjointProductBoundedBy<RNDM<F, N>, D>
+for ZeroOp<'a, RNDM<F, N>, XD, Y, F>
+where
+    F : Float,
+    Y : AXPY<F> + Clone,
+    D : Linear<RNDM<F, N>>,
+{
+    type FloatType = F;
+    /// Return $L$ such that $A_*A ≤ L𝒟$ is bounded by some `other` operator $𝒟$.
+    fn adjoint_product_bound(&self, _ : &D) -> Option<F> {
+        Some(0.0)
+    }
+}
+