Mercurial > repos > pointsource_algs / file comparison
comparison: src/forward_pdps.rs
src/forward_pdps.rs
- changeset 75
- 677a5fd1b014
- parent 72
- e9a460a0e638
- child 78
- 2a122736e91c
equal
deleted
inserted
replaced
| 74:df92e78cc3f4 | 75:677a5fd1b014 |
|---|---|
| 195 // postprocessing: config.insertion.postprocessing.then(|| μ.clone()), | 195 // postprocessing: config.insertion.postprocessing.then(|| μ.clone()), |
| 196 ..stats | 196 ..stats |
| 197 }; | 197 }; |
| 198 let mut stats = IterInfo::new(); | 198 let mut stats = IterInfo::new(); |
| 199 | 199 |
| 200 let zy_step = |mut z, mut y, mut z_tmp| { | |
| 201 // Do z variable primal update | |
| 202 opKz_adj.apply_add(&mut z_tmp, &y); | |
| 203 z_tmp.axpy(1.0, &z, -σ_p); | |
| 204 let z_new = fnR.prox(σ_p, z_tmp); | |
| 205 //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z); | |
| 206 // Do dual update | |
| 207 z.axpy(1.0 + ω, &z_new, -ω); | |
| 208 //z = (z - &z_new) * (-ω) + &z_new; | |
| 209 opKz.gemv(&mut y, σ_d, z, 1.0); | |
| 210 let y = starH.prox(σ_d, y); | |
| 211 (z_new, y) | |
| 212 }; | |
| 213 | |
| 200 // Run the algorithm | 214 // Run the algorithm |
| 201 for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) { | 215 for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) { |
| 202 // Calculate initial transport | 216 // Calculate initial transport |
| 203 let Pair(v, mut z_tmp) = f.differential(Pair(&μ, &z)); | 217 let Pair(v, z_tmp) = f.differential(Pair(&μ, &z)); |
| 204 let mut τv = v * τ; | 218 let mut τv = v * τ; |
| 205 let μ_base = μ.clone(); | 219 let μ_base = μ.clone(); |
| 206 | 220 |
| 207 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. | 221 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. |
| 208 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( | 222 prox_penalty.insert_and_reweigh( |
| 209 &mut μ, | 223 &mut μ, |
| 210 &mut τv, | 224 &mut τv, |
| 211 τ, | 225 τ, |
| 212 ε, | 226 ε, |
| 213 &config.insertion, | 227 &config.insertion, |
| 214 ®, | 228 ®, |
| 215 &state, | 229 &state, |
| 216 &mut stats, | 230 &mut stats, |
| 217 )?; | 231 )?; |
| 232 | |
| 233 (z, y) = zy_step(z, y, z_tmp); | |
| 218 | 234 |
| 219 stats.inserted += μ.len() - μ_base.len(); | 235 stats.inserted += μ.len() - μ_base.len(); |
| 220 | 236 |
| 221 // Merge spikes. | 237 // Merge spikes. |
| 222 // This crucially expects the merge routine to be stable with respect to spike locations, | 238 // This crucially expects the merge routine to be stable with respect to spike locations, |
| 239 } | 255 } |
| 240 | 256 |
| 241 // Prune spikes with zero weight. | 257 // Prune spikes with zero weight. |
| 242 stats.pruned += prune_with_stats(&mut μ); | 258 stats.pruned += prune_with_stats(&mut μ); |
| 243 | 259 |
| 244 // Do z variable primal update | 260 for _ in 0..config.insertion.extra_weight_optimisation_steps { |
| 245 opKz_adj.apply_add(&mut z_tmp, &y); | 261 // Do extra weight optimisation step heuristic |
| 246 z_tmp.axpy(1.0, &z, -σ_p); | 262 let Pair(v, _z_tmp) = f.differential(Pair(&μ, &z)); |
| 247 let z_new = fnR.prox(σ_p, z_tmp); | 263 τv = v * τ; |
| 248 //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z); | 264 prox_penalty.reweigh( |
| 249 // Do dual update | 265 &mut μ, |
| 250 z.axpy(1.0 + ω, &z_new, -ω); | 266 &mut τv, |
| 251 //z = (z - &z_new) * (-ω) + &z_new; | 267 τ, |
| 252 opKz.gemv(&mut y, σ_d, z, 1.0); | 268 ε, |
| 253 y = starH.prox(σ_d, y); | 269 &config.insertion, |
| 254 z = z_new; | 270 ®, |
| 271 &state, | |
| 272 &mut stats, | |
| 273 )?; | |
| 274 //(z, y) = zy_step(z, y, z_tmp); | |
| 275 } | |
| 276 | |
| 277 if config.insertion.extra_weight_optimisation_steps > 0 { | |
| 278 // Prune spikes with zero weight. | |
| 279 stats.pruned += prune_with_stats(&mut μ); | |
| 280 } | |
| 255 | 281 |
| 256 // Update step length parameters | 282 // Update step length parameters |
| 257 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); | 283 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); |
| 258 | 284 |
| 259 // Give statistics if requested | 285 // Give statistics if requested |
| 260 let iter = state.iteration(); | 286 let iter = state.iteration(); |
| 261 stats.this_iters += 1; | 287 stats.this_iters += 1; |
| 262 | 288 |
| 263 state.if_verbose(|| { | 289 state.if_verbose(|| { |
| 264 plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ); | 290 plotter.plot_spikes(iter, None, Some(&τv), &μ); |
| 265 full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) | 291 full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) |
| 266 }); | 292 }); |
| 267 | 293 |
| 268 // Update main tolerance for next iteration | 294 // Update main tolerance for next iteration |
| 269 ε = tolerance.update(ε, iter); | 295 ε = tolerance.update(ε, iter); |
