Mercurial > repos > pointsource_algs / file comparison
comparison: src/forward_model/bias.rs
src/forward_model/bias.rs
- branch
- dev
- changeset 35
- b087e3eab191
- child 37
- c5d8bd1a7728
equal
deleted
inserted
replaced
| 34:efa60bc4f743 | 35:b087e3eab191 |
|---|---|
| 1 /*! | |
| 2 Simple parametric forward model. | |
| 3 */ | |
| 4 | |
| 5 use numeric_literals::replace_float_literals; | |
| 6 use alg_tools::types::{Float, ClosedAdd}; | |
| 7 use alg_tools::mapping::Space; | |
| 8 use alg_tools::direct_product::Pair; | |
| 9 use alg_tools::linops::{Linear, RowOp, ColOp, IdOp, ZeroOp, AXPY}; | |
| 10 use alg_tools::error::DynError; | |
| 11 use alg_tools::norms::{L2, Norm, PairNorm, NormExponent}; | |
| 12 use crate::types::L2Squared; | |
| 13 use crate::measures::RNDM; | |
| 14 use super::{ForwardModel, AdjointProductBoundedBy, AdjointProductPairBoundedBy, LipschitzValues}; | |
| 15 use crate::transport::TransportLipschitz; | |
| 16 | |
| 17 impl<Domain, F, A, E> ForwardModel<Pair<Domain, A::Observable>, F, PairNorm<E, L2, L2>> | |
| 18 for RowOp<A, IdOp<A::Observable>> | |
| 19 where | |
| 20 E : NormExponent, | |
| 21 Domain : Space + Norm<F, E>, | |
| 22 F : Float, | |
| 23 A::Observable : ClosedAdd + Norm<F, L2> + 'static, | |
| 24 A : ForwardModel<Domain, F, E> + 'static | |
| 25 { | |
| 26 type Observable = A::Observable; | |
| 27 | |
| 28 fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError { | |
| 29 self.0.write_observable(b, prefix) | |
| 30 } | |
| 31 | |
| 32 /// Returns a zero observable | |
| 33 fn zero_observable(&self) -> Self::Observable { | |
| 34 self.0.zero_observable() | |
| 35 } | |
| 36 } | |
| 37 | |
| 38 #[replace_float_literals(F::cast_from(literal))] | |
| 39 impl<Domain, F, A, D, Z> AdjointProductPairBoundedBy<Pair<Domain, Z>, D, IdOp<Z>> | |
| 40 for RowOp<A, IdOp<Z>> | |
| 41 where | |
| 42 Domain : Space, | |
| 43 F : Float, | |
| 44 Z : Clone + Space + ClosedAdd, | |
| 45 A : AdjointProductBoundedBy<Domain, D, FloatType=F, Codomain = Z>, | |
| 46 D : Linear<Domain>, | |
| 47 A::Codomain : ClosedAdd, | |
| 48 { | |
| 49 type FloatType = F; | |
| 50 | |
| 51 fn adjoint_product_pair_bound(&self, d : &D, _ : &IdOp<Z>) -> Option<(F, F)> { | |
| 52 self.0.adjoint_product_bound(d).map(|l_0| { | |
| 53 // [A_*; B_*][A, B] = [A_*A, A_* B; B_* A, B_* B] ≤ diag(2A_*A, 2B_*B) | |
| 54 // ≤ diag(2l_A𝒟_A, 2l_B𝒟_B), where now 𝒟_B=Id and l_B=1. | |
| 55 (2.0 * l_0, 2.0) | |
| 56 }) | |
| 57 } | |
| 58 } | |
| 59 | |
| 60 /// This `impl` is bit of an abuse as the codomain of `Apre` is a [`Pair`] of a measure predual, | |
| 61 /// to which this `impl` applies, and another space. | |
| 62 impl<F, Apre, Z> LipschitzValues | |
| 63 for ColOp<Apre, IdOp<Z>> | |
| 64 where | |
| 65 F : Float, | |
| 66 Z : Clone + Space + ClosedAdd, | |
| 67 Apre : LipschitzValues<FloatType = F>, | |
| 68 { | |
| 69 type FloatType = F; | |
| 70 /// Return (if one exists) a factor $L$ such that $A_*z$ is $L$-Lipschitz for all | |
| 71 /// $z$ in the unit ball. | |
| 72 fn value_unit_lipschitz_factor(&self) -> Option<Self::FloatType> { | |
| 73 self.0.value_unit_lipschitz_factor() | |
| 74 } | |
| 75 | |
| 76 /// Return (if one exists) a factor $L$ such that $∇A_*z$ is $L$-Lipschitz for all | |
| 77 /// $z$ in the unit ball. | |
| 78 fn value_diff_unit_lipschitz_factor(&self) -> Option<Self::FloatType> { | |
| 79 self.0.value_diff_unit_lipschitz_factor() | |
| 80 } | |
| 81 } | |
| 82 | |
| 83 | |
| 84 | |
| 85 impl<'a, F : Float, Y : Space, XD, const N : usize> TransportLipschitz<L2Squared> for | |
| 86 ZeroOp<'a, RNDM<F, N>, XD, Y, F> { | |
| 87 type FloatType = F; | |
| 88 | |
| 89 fn transport_lipschitz_factor(&self, _ : L2Squared) -> Self::FloatType { | |
| 90 F::ZERO | |
| 91 } | |
| 92 } | |
| 93 | |
| 94 | |
| 95 /// TODO: should assume `D` to be positive semi-definite and self-adjoint. | |
| 96 #[replace_float_literals(F::cast_from(literal))] | |
| 97 impl<'a, F, D, XD, Y, const N : usize> AdjointProductBoundedBy<RNDM<F, N>, D> | |
| 98 for ZeroOp<'a, RNDM<F, N>, XD, Y, F> | |
| 99 where | |
| 100 F : Float, | |
| 101 Y : AXPY<F> + Clone, | |
| 102 D : Linear<RNDM<F, N>>, | |
| 103 { | |
| 104 type FloatType = F; | |
| 105 /// Return $L$ such that $A_*A ≤ L𝒟$ is bounded by some `other` operator $𝒟$. | |
| 106 fn adjoint_product_bound(&self, _ : &D) -> Option<F> { | |
| 107 Some(0.0) | |
| 108 } | |
| 109 } | |
| 110 |
