Mercurial > repos > pointsource_algs / file comparison
comparison: src/forward_pdps.rs
src/forward_pdps.rs
- changeset 72
- e9a460a0e638
- parent 66
- fe47ad484deb
- child 75
- 677a5fd1b014
equal
deleted
inserted
replaced
| 71:e2953ffd4e0b | 72:e9a460a0e638 |
|---|---|
| 126 + GEMV<F, Z, Y> | 126 + GEMV<F, Z, Y> |
| 127 + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>, | 127 + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>, |
| 128 KOpZ::SimpleAdjoint: GEMV<F, Y, Z>, | 128 KOpZ::SimpleAdjoint: GEMV<F, Y, Z>, |
| 129 Y: ClosedEuclidean<F>, | 129 Y: ClosedEuclidean<F>, |
| 130 for<'b> &'b Y: Instance<Y>, | 130 for<'b> &'b Y: Instance<Y>, |
| 131 Z: ClosedEuclidean<F>, | 131 Z: ClosedEuclidean<F> + AXPY, |
| 132 for<'b> &'b Z: Instance<Z>, | 132 for<'b> &'b Z: Instance<Z>, |
| 133 R: Prox<Z, Codomain = F>, | 133 R: Prox<Z, Codomain = F>, |
| 134 H: Conjugable<Y, F, Codomain = F>, | 134 H: Conjugable<Y, F, Codomain = F>, |
| 135 for<'b> H::Conjugate<'b>: Prox<Y>, | 135 for<'b> H::Conjugate<'b>: Prox<Y>, |
| 136 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>, | 136 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>, |
| 198 let mut stats = IterInfo::new(); | 198 let mut stats = IterInfo::new(); |
| 199 | 199 |
| 200 // Run the algorithm | 200 // Run the algorithm |
| 201 for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) { | 201 for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) { |
| 202 // Calculate initial transport | 202 // Calculate initial transport |
| 203 let Pair(mut τv, τz) = f.differential(Pair(&μ, &z)) * τ; | 203 let Pair(v, mut z_tmp) = f.differential(Pair(&μ, &z)); |
| 204 let mut τv = v * τ; | |
| 204 let μ_base = μ.clone(); | 205 let μ_base = μ.clone(); |
| 205 | 206 |
| 206 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. | 207 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. |
| 207 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( | 208 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( |
| 208 &mut μ, | 209 &mut μ, |
| 231 &μ_base, | 232 &μ_base, |
| 232 τ, | 233 τ, |
| 233 ε, | 234 ε, |
| 234 ins, | 235 ins, |
| 235 ®, | 236 ®, |
| 236 is_fb.then_some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))), | 237 Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))), |
| 237 ); | 238 ); |
| 238 } | 239 } |
| 239 | 240 |
| 240 // Prune spikes with zero weight. | 241 // Prune spikes with zero weight. |
| 241 stats.pruned += prune_with_stats(&mut μ); | 242 stats.pruned += prune_with_stats(&mut μ); |
| 242 | 243 |
| 243 // Do z variable primal update | 244 // Do z variable primal update |
| 244 let mut z_new = τz; | 245 opKz_adj.apply_add(&mut z_tmp, &y); |
| 245 opKz_adj.gemv(&mut z_new, -σ_p, &y, -σ_p / τ); | 246 z_tmp.axpy(1.0, &z, -σ_p); |
| 246 z_new = fnR.prox(σ_p, z_new + &z); | 247 let z_new = fnR.prox(σ_p, z_tmp); |
| 248 //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z); | |
| 247 // Do dual update | 249 // Do dual update |
| 248 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] | 250 z.axpy(1.0 + ω, &z_new, -ω); |
| 249 opKz.gemv(&mut y, σ_d * (1.0 + ω), &z_new, 1.0); | 251 //z = (z - &z_new) * (-ω) + &z_new; |
| 250 // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b | 252 opKz.gemv(&mut y, σ_d, z, 1.0); |
| 251 opKz.gemv(&mut y, -σ_d * ω, z, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b | |
| 252 y = starH.prox(σ_d, y); | 253 y = starH.prox(σ_d, y); |
| 253 z = z_new; | 254 z = z_new; |
| 254 | 255 |
| 255 // Update step length parameters | 256 // Update step length parameters |
| 256 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); | 257 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); |
