Mercurial > repos > pointsource_algs / file comparison
comparison: src/forward_pdps.rs
src/forward_pdps.rs
- branch
- dev
- changeset 35
- b087e3eab191
- child 37
- c5d8bd1a7728
equal
deleted
inserted
replaced
| 34:efa60bc4f743 | 35:b087e3eab191 |
|---|---|
| 1 /*! | |
| 2 Solver for the point source localisation problem using a | |
| 3 primal-dual proximal splitting with a forward step. | |
| 4 */ | |
| 5 | |
| 6 use numeric_literals::replace_float_literals; | |
| 7 use serde::{Serialize, Deserialize}; | |
| 8 | |
| 9 use alg_tools::iterate::AlgIteratorFactory; | |
| 10 use alg_tools::euclidean::Euclidean; | |
| 11 use alg_tools::sets::Cube; | |
| 12 use alg_tools::loc::Loc; | |
| 13 use alg_tools::mapping::{Mapping, Instance}; | |
| 14 use alg_tools::norms::Norm; | |
| 15 use alg_tools::direct_product::Pair; | |
| 16 use alg_tools::bisection_tree::{ | |
| 17 BTFN, | |
| 18 PreBTFN, | |
| 19 Bounds, | |
| 20 BTNodeLookup, | |
| 21 BTNode, | |
| 22 BTSearch, | |
| 23 P2Minimise, | |
| 24 SupportGenerator, | |
| 25 LocalAnalysis, | |
| 26 //Bounded, | |
| 27 }; | |
| 28 use alg_tools::mapping::RealMapping; | |
| 29 use alg_tools::nalgebra_support::ToNalgebraRealField; | |
| 30 use alg_tools::linops::{ | |
| 31 BoundedLinear, AXPY, GEMV, Adjointable, IdOp, | |
| 32 }; | |
| 33 use alg_tools::convex::{Conjugable, Prox}; | |
| 34 use alg_tools::norms::{L2, Linfinity, PairNorm}; | |
| 35 | |
| 36 use crate::types::*; | |
| 37 use crate::measures::{DiscreteMeasure, Radon, RNDM}; | |
| 38 use crate::measures::merging::SpikeMerging; | |
| 39 use crate::forward_model::{ | |
| 40 ForwardModel, | |
| 41 AdjointProductPairBoundedBy, | |
| 42 }; | |
| 43 use crate::seminorms::DiscreteMeasureOp; | |
| 44 use crate::plot::{ | |
| 45 SeqPlotter, | |
| 46 Plotting, | |
| 47 PlotLookup | |
| 48 }; | |
| 49 use crate::fb::*; | |
| 50 use crate::regularisation::RegTerm; | |
| 51 use crate::dataterm::calculate_residual; | |
| 52 | |
| 53 /// Settings for [`pointsource_forward_pdps_pair`]. | |
| 54 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] | |
| 55 #[serde(default)] | |
| 56 pub struct ForwardPDPSConfig<F : Float> { | |
| 57 /// Primal step length scaling. | |
| 58 pub τ0 : F, | |
| 59 /// Primal step length scaling. | |
| 60 pub σp0 : F, | |
| 61 /// Dual step length scaling. | |
| 62 pub σd0 : F, | |
| 63 /// Generic parameters | |
| 64 pub insertion : FBGenericConfig<F>, | |
| 65 } | |
| 66 | |
| 67 #[replace_float_literals(F::cast_from(literal))] | |
| 68 impl<F : Float> Default for ForwardPDPSConfig<F> { | |
| 69 fn default() -> Self { | |
| 70 let τ0 = 0.99; | |
| 71 ForwardPDPSConfig { | |
| 72 τ0, | |
| 73 σd0 : 0.1, | |
| 74 σp0 : 0.99, | |
| 75 insertion : Default::default() | |
| 76 } | |
| 77 } | |
| 78 } | |
| 79 | |
| 80 type MeasureZ<F, Z, const N : usize> = Pair<RNDM<F, N>, Z>; | |
| 81 | |
| 82 /// Iteratively solve the pointsource localisation with an additional variable | |
| 83 /// using primal-dual proximal splitting with a forward step. | |
| 84 #[replace_float_literals(F::cast_from(literal))] | |
| 85 pub fn pointsource_forward_pdps_pair< | |
| 86 'a, F, I, A, GA, 𝒟, BTA, BT𝒟, G𝒟, S, K, Reg, Z, R, Y, /*KOpM, */ KOpZ, H, const N : usize | |
| 87 >( | |
| 88 opA : &'a A, | |
| 89 b : &A::Observable, | |
| 90 reg : Reg, | |
| 91 op𝒟 : &'a 𝒟, | |
| 92 config : &ForwardPDPSConfig<F>, | |
| 93 iterator : I, | |
| 94 mut plotter : SeqPlotter<F, N>, | |
| 95 //opKμ : KOpM, | |
| 96 opKz : &KOpZ, | |
| 97 fnR : &R, | |
| 98 fnH : &H, | |
| 99 mut z : Z, | |
| 100 mut y : Y, | |
| 101 ) -> MeasureZ<F, Z, N> | |
| 102 where | |
| 103 F : Float + ToNalgebraRealField, | |
| 104 I : AlgIteratorFactory<IterInfo<F, N>>, | |
| 105 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable> + Instance<A::Observable>, | |
| 106 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone, | |
| 107 A : ForwardModel< | |
| 108 MeasureZ<F, Z, N>, | |
| 109 F, | |
| 110 PairNorm<Radon, L2, L2>, | |
| 111 PreadjointCodomain = Pair<BTFN<F, GA, BTA, N>, Z>, | |
| 112 > | |
| 113 + AdjointProductPairBoundedBy<MeasureZ<F, Z, N>, 𝒟, IdOp<Z>, FloatType=F>, | |
| 114 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>, | |
| 115 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone, | |
| 116 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>, | |
| 117 Codomain = BTFN<F, G𝒟, BT𝒟, N>>, | |
| 118 BT𝒟 : BTSearch<F, N, Data=usize, Agg=Bounds<F>>, | |
| 119 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, | |
| 120 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, | |
| 121 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>, | |
| 122 Cube<F, N>: P2Minimise<Loc<F, N>, F>, | |
| 123 PlotLookup : Plotting<N>, | |
| 124 RNDM<F, N> : SpikeMerging<F>, | |
| 125 Reg : RegTerm<F, N>, | |
| 126 KOpZ : BoundedLinear<Z, L2, L2, F, Codomain=Y> | |
| 127 + GEMV<F, Z> | |
| 128 + Adjointable<Z, Y, AdjointCodomain = Z>, | |
| 129 for<'b> KOpZ::Adjoint<'b> : GEMV<F, Y>, | |
| 130 Y : AXPY<F> + Euclidean<F, Output=Y> + Clone + ClosedAdd, | |
| 131 for<'b> &'b Y : Instance<Y>, | |
| 132 Z : AXPY<F, Owned=Z> + Euclidean<F, Output=Z> + Clone + Norm<F, L2>, | |
| 133 for<'b> &'b Z : Instance<Z>, | |
| 134 R : Prox<Z, Codomain=F>, | |
| 135 H : Conjugable<Y, F, Codomain=F>, | |
| 136 for<'b> H::Conjugate<'b> : Prox<Y>, | |
| 137 { | |
| 138 | |
| 139 // Check parameters | |
| 140 assert!(config.τ0 > 0.0 && | |
| 141 config.τ0 < 1.0 && | |
| 142 config.σp0 > 0.0 && | |
| 143 config.σp0 < 1.0 && | |
| 144 config.σd0 > 0.0 && | |
| 145 config.σp0 * config.σd0 <= 1.0, | |
| 146 "Invalid step length parameters"); | |
| 147 | |
| 148 // Initialise iterates | |
| 149 let mut μ = DiscreteMeasure::new(); | |
| 150 let mut residual = calculate_residual(Pair(&μ, &z), opA, b); | |
| 151 | |
| 152 // Set up parameters | |
| 153 let op𝒟norm = op𝒟.opnorm_bound(Radon, Linfinity); | |
| 154 let bigM = 0.0; //opKμ.adjoint_product_bound(&op𝒟).unwrap().sqrt(); | |
| 155 let nKz = opKz.opnorm_bound(L2, L2); | |
| 156 let opIdZ = IdOp::new(); | |
| 157 let (l, l_z) = opA.adjoint_product_pair_bound(&op𝒟, &opIdZ).unwrap(); | |
| 158 // We need to satisfy | |
| 159 // | |
| 160 // τσ_dM(1-σ_p L_z)/(1 - τ L) + [σ_p L_z + σ_pσ_d‖K_z‖^2] < 1 | |
| 161 // ^^^^^^^^^^^^^^^^^^^^^^^^^ | |
| 162 // with 1 > σ_p L_z and 1 > τ L. | |
| 163 // | |
| 164 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition | |
| 165 // ^^^^^ < 1. then we solve τ from the rest. | |
| 166 let σ_d = config.σd0 / nKz; | |
| 167 let σ_p = config.σp0 / (l_z + config.σd0 * nKz); | |
| 168 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0} | |
| 169 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L) | |
| 170 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0}) | |
| 171 let φ = 1.0 - config.σp0; | |
| 172 let a = 1.0 - σ_p * l_z; | |
| 173 let τ = config.τ0 * φ / ( σ_d * bigM * a + φ * l ); | |
| 174 // Acceleration is not currently supported | |
| 175 // let γ = dataterm.factor_of_strong_convexity(); | |
| 176 let ω = 1.0; | |
| 177 | |
| 178 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled | |
| 179 // by τ compared to the conditional gradient approach. | |
| 180 let tolerance = config.insertion.tolerance * τ * reg.tolerance_scaling(); | |
| 181 let mut ε = tolerance.initial(); | |
| 182 | |
| 183 let starH = fnH.conjugate(); | |
| 184 | |
| 185 // Statistics | |
| 186 let full_stats = |residual : &A::Observable, μ : &RNDM<F, N>, z : &Z, ε, stats| IterInfo { | |
| 187 value : residual.norm2_squared_div2() + fnR.apply(z) | |
| 188 + reg.apply(μ) + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)), | |
| 189 n_spikes : μ.len(), | |
| 190 ε, | |
| 191 // postprocessing: config.insertion.postprocessing.then(|| μ.clone()), | |
| 192 .. stats | |
| 193 }; | |
| 194 let mut stats = IterInfo::new(); | |
| 195 | |
| 196 // Run the algorithm | |
| 197 for state in iterator.iter_init(|| full_stats(&residual, &μ, &z, ε, stats.clone())) { | |
| 198 // Calculate initial transport | |
| 199 let Pair(τv, τz) = opA.preadjoint().apply(residual * τ); | |
| 200 let z_base = z.clone(); | |
| 201 let μ_base = μ.clone(); | |
| 202 | |
| 203 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. | |
| 204 let (d, _within_tolerances) = insert_and_reweigh( | |
| 205 &mut μ, &τv, &μ_base, None, | |
| 206 op𝒟, op𝒟norm, | |
| 207 τ, ε, &config.insertion, | |
| 208 ®, &state, &mut stats, | |
| 209 ); | |
| 210 | |
| 211 // // Merge spikes. | |
| 212 // // This expects the prune below to prune γ. | |
| 213 // // TODO: This may not work correctly in all cases. | |
| 214 // let ins = &config.insertion; | |
| 215 // if ins.merge_now(&state) { | |
| 216 // if let SpikeMergingMethod::None = ins.merging { | |
| 217 // } else { | |
| 218 // stats.merged += μ.merge_spikes(ins.merging, |μ_candidate| { | |
| 219 // let ν = μ_candidate.sub_matching(&γ1)-&μ_base_minus_γ0; | |
| 220 // let mut d = &τv̆ + op𝒟.preapply(ν); | |
| 221 // reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, ins) | |
| 222 // }); | |
| 223 // } | |
| 224 // } | |
| 225 | |
| 226 // Prune spikes with zero weight. | |
| 227 stats.pruned += prune_with_stats(&mut μ); | |
| 228 | |
| 229 // Do z variable primal update | |
| 230 z.axpy(-σ_p/τ, τz, 1.0); // TODO: simplify nasty factors | |
| 231 opKz.adjoint().gemv(&mut z, -σ_p, &y, 1.0); | |
| 232 z = fnR.prox(σ_p, z); | |
| 233 // Do dual update | |
| 234 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] | |
| 235 opKz.gemv(&mut y, σ_d*(1.0 + ω), &z, 1.0); | |
| 236 // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b | |
| 237 opKz.gemv(&mut y, -σ_d*ω, z_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b | |
| 238 y = starH.prox(σ_d, y); | |
| 239 | |
| 240 // Update residual | |
| 241 residual = calculate_residual(Pair(&μ, &z), opA, b); | |
| 242 | |
| 243 // Update step length parameters | |
| 244 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); | |
| 245 | |
| 246 // Give statistics if requested | |
| 247 let iter = state.iteration(); | |
| 248 stats.this_iters += 1; | |
| 249 | |
| 250 state.if_verbose(|| { | |
| 251 plotter.plot_spikes(iter, Some(&d), Some(&τv), &μ); | |
| 252 full_stats(&residual, &μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) | |
| 253 }); | |
| 254 | |
| 255 // Update main tolerance for next iteration | |
| 256 ε = tolerance.update(ε, iter); | |
| 257 } | |
| 258 | |
| 259 let fit = |μ̃ : &RNDM<F, N>| { | |
| 260 (opA.apply(Pair(μ̃, &z))-b).norm2_squared_div2() | |
| 261 //+ fnR.apply(z) + reg.apply(μ) | |
| 262 + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z)) | |
| 263 }; | |
| 264 | |
| 265 μ.merge_spikes_fitness(config.insertion.merging, fit, |&v| v); | |
| 266 μ.prune(); | |
| 267 Pair(μ, z) | |
| 268 } |
