Mercurial > repos > pointsource_algs / file comparison
comparison: src/sliding_pdps.rs
src/sliding_pdps.rs
- branch
- dev
- changeset 44
- 03251c546744
- parent 41
- b6bdb6cb4d44
- child 45
- 5200e7090e06
equal
deleted
inserted
replaced
| 43:aacd9af21b3a | 44:03251c546744 |
|---|---|
| 25 use crate::measures::{DiscreteMeasure, Radon, RNDM}; | 25 use crate::measures::{DiscreteMeasure, Radon, RNDM}; |
| 26 use crate::measures::merging::SpikeMerging; | 26 use crate::measures::merging::SpikeMerging; |
| 27 use crate::forward_model::{ | 27 use crate::forward_model::{ |
| 28 ForwardModel, | 28 ForwardModel, |
| 29 AdjointProductPairBoundedBy, | 29 AdjointProductPairBoundedBy, |
| 30 LipschitzValues, | 30 BoundedCurvature, |
| 31 }; | 31 }; |
| 32 // use crate::transport::TransportLipschitz; | 32 // use crate::transport::TransportLipschitz; |
| 33 //use crate::tolerance::Tolerance; | 33 //use crate::tolerance::Tolerance; |
| 34 use crate::plot::{ | 34 use crate::plot::{ |
| 35 SeqPlotter, | 35 SeqPlotter, |
| 111 MeasureZ<F, Z, N>, | 111 MeasureZ<F, Z, N>, |
| 112 F, | 112 F, |
| 113 PairNorm<Radon, L2, L2>, | 113 PairNorm<Radon, L2, L2>, |
| 114 PreadjointCodomain = Pair<S, Z>, | 114 PreadjointCodomain = Pair<S, Z>, |
| 115 > | 115 > |
| 116 + AdjointProductPairBoundedBy<MeasureZ<F, Z, N>, P, IdOp<Z>, FloatType=F>, | 116 + AdjointProductPairBoundedBy<MeasureZ<F, Z, N>, P, IdOp<Z>, FloatType=F> |
| 117 + BoundedCurvature<FloatType=F>, | |
| 117 S : DifferentiableRealMapping<F, N>, | 118 S : DifferentiableRealMapping<F, N>, |
| 118 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable> + Instance<A::Observable>, | 119 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable> + Instance<A::Observable>, |
| 119 for<'b> A::Preadjoint<'b> : LipschitzValues<FloatType=F>, | |
| 120 PlotLookup : Plotting<N>, | 120 PlotLookup : Plotting<N>, |
| 121 RNDM<F, N> : SpikeMerging<F>, | 121 RNDM<F, N> : SpikeMerging<F>, |
| 122 Reg : SlidingRegTerm<F, N>, | 122 Reg : SlidingRegTerm<F, N>, |
| 123 P : ProxPenalty<F, S, Reg, N>, | 123 P : ProxPenalty<F, S, Reg, N>, |
| 124 // KOpM : Linear<RNDM<F, N>, Codomain=Y> | 124 // KOpM : Linear<RNDM<F, N>, Codomain=Y> |
| 186 let a = 1.0 - σ_p * l_z; | 186 let a = 1.0 - σ_p * l_z; |
| 187 let τ = config.τ0 * φ / ( σ_d * bigM * a + φ * l ); | 187 let τ = config.τ0 * φ / ( σ_d * bigM * a + φ * l ); |
| 188 let ψ = 1.0 - τ * l; | 188 let ψ = 1.0 - τ * l; |
| 189 let β = σ_p * config.σd0 * nKz / a; // σ_p * σ_d * (nKz * nK_z) / a; | 189 let β = σ_p * config.σd0 * nKz / a; // σ_p * σ_d * (nKz * nK_z) / a; |
| 190 assert!(β < 1.0); | 190 assert!(β < 1.0); |
| 191 // Now we need κ‖K_μ(π_♯^1 - π_♯^0)γ‖^2 ≤ (1/θ - τ[ℓ_v + ℓ]) ∫ c_2 dγ for κ defined as: | 191 // Now we need κ‖K_μ(π_♯^1 - π_♯^0)γ‖^2 ≤ (1/θ - τ[ℓ_F + ℓ]) ∫ c_2 dγ for κ defined as: |
| 192 let κ = τ * σ_d * ψ / ((1.0 - β) * ψ - τ * σ_d * bigM); | 192 let κ = τ * σ_d * ψ / ((1.0 - β) * ψ - τ * σ_d * bigM); |
| 193 // The factor two in the manuscript disappears due to the definition of 𝚹 being | 193 // The factor two in the manuscript disappears due to the definition of 𝚹 being |
| 194 // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2. | 194 // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2. |
| 195 let (maybe_ℓ_v0, maybe_transport_lip) = opA.curvature_bound_components(); | |
| 196 let transport_lip = maybe_transport_lip.unwrap(); | |
| 195 let calculate_θ = |ℓ_v, max_transport| { | 197 let calculate_θ = |ℓ_v, max_transport| { |
| 196 config.transport.θ0 / (τ*(ℓ + ℓ_v) + κ * bigθ * max_transport) | 198 let ℓ_F = ℓ_v + transport_lip * max_transport; |
| 199 config.transport.θ0 / (τ*(ℓ + ℓ_F) + κ * bigθ * max_transport) | |
| 197 }; | 200 }; |
| 198 let mut θ_or_adaptive = match opA.preadjoint().value_diff_unit_lipschitz_factor() { | 201 let mut θ_or_adaptive = match maybe_ℓ_v0 { |
| 199 // We only estimate w (the uniform Lipschitz for of v), if we also estimate ℓ_v | |
| 200 // (the uniform Lipschitz factor of ∇v). | |
| 201 // We assume that the residual is decreasing. | 202 // We assume that the residual is decreasing. |
| 202 Some(ℓ_v0) => TransportStepLength::AdaptiveMax{ | 203 Some(ℓ_v0) => TransportStepLength::AdaptiveMax { |
| 203 l: ℓ_v0 * b.norm2(), | 204 l: ℓ_v0 * b.norm2(), // TODO: could estimate computing the real reesidual |
| 204 max_transport : 0.0, | 205 max_transport : 0.0, |
| 205 g : calculate_θ | 206 g : calculate_θ |
| 206 }, | 207 }, |
| 207 None => TransportStepLength::FullyAdaptive{ | 208 None => TransportStepLength::FullyAdaptive { |
| 208 l : F::EPSILON, | 209 l : F::EPSILON, |
| 209 max_transport : 0.0, | 210 max_transport : 0.0, |
| 210 g : calculate_θ | 211 g : calculate_θ |
| 211 }, | 212 }, |
| 212 }; | 213 }; |
