Mercurial > repos > pointsource_algs / file comparison
comparison: src/sliding_pdps.rs
src/sliding_pdps.rs
- changeset 75
- 677a5fd1b014
- parent 72
- e9a460a0e638
- child 78
- 2a122736e91c
equal
deleted
inserted
replaced
| 74:df92e78cc3f4 | 75:677a5fd1b014 |
|---|---|
| 240 | 240 |
| 241 let mut attempts = 0; | 241 let mut attempts = 0; |
| 242 | 242 |
| 243 // Solve finite-dimensional subproblem several times until the dual variable for the | 243 // Solve finite-dimensional subproblem several times until the dual variable for the |
| 244 // regularisation term conforms to the assumptions made for the transport above. | 244 // regularisation term conforms to the assumptions made for the transport above. |
| 245 let (maybe_d, _within_tolerances, mut τv̆, z_new, μ̆) = 'adapt_transport: loop { | 245 let (mut τv̆, z_new, μ̆) = 'adapt_transport: loop { |
| 246 // Set initial guess for μ=μ^{k+1}. | 246 // Set initial guess for μ=μ^{k+1}. |
| 247 γ.μ̆_into(&mut μ); | 247 γ.μ̆_into(&mut μ); |
| 248 let μ̆ = μ.clone(); | 248 let μ̆ = μ.clone(); |
| 249 | 249 |
| 250 // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) | 250 // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) |
| 251 let Pair(v̆, mut z̆) = f.differential(Pair(&μ̆, &z)); | 251 let Pair(v̆, mut z̆) = f.differential(Pair(&μ̆, &z)); |
| 252 let mut τv̆ = v̆ * τ; | 252 let mut τv̆ = v̆ * τ; |
| 253 // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0); | 253 // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0); |
| 254 | 254 |
| 255 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. | 255 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. |
| 256 let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh( | 256 prox_penalty.insert_and_reweigh( |
| 257 &mut μ, | 257 &mut μ, |
| 258 &mut τv̆, | 258 &mut τv̆, |
| 259 τ, | 259 τ, |
| 260 ε, | 260 ε, |
| 261 &config.insertion, | 261 &config.insertion, |
| 284 reg, | 284 reg, |
| 285 &config.transport, | 285 &config.transport, |
| 286 &config.insertion.refinement, | 286 &config.insertion.refinement, |
| 287 &mut attempts, | 287 &mut attempts, |
| 288 ) { | 288 ) { |
| 289 break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new, μ̆); | 289 break 'adapt_transport (τv̆, z_new, μ̆); |
| 290 } | 290 } |
| 291 }; | 291 }; |
| 292 | 292 |
| 293 γ.get_transport_stats(&mut stats, &μ); | 293 γ.get_transport_stats(&mut stats, &μ); |
| 294 | |
| 295 // Do dual update | |
| 296 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] | |
| 297 z.axpy(1.0 + ω, &z_new, -ω); | |
| 298 //z = (z - &z_new) * (-ω) + &z_new; | |
| 299 opKz.gemv(&mut y, σ_d, z, 1.0); | |
| 300 y = starH.prox(σ_d, y); | |
| 301 z = z_new; | |
| 294 | 302 |
| 295 // Merge spikes. | 303 // Merge spikes. |
| 296 // This crucially expects the merge routine to be stable with respect to spike locations, | 304 // This crucially expects the merge routine to be stable with respect to spike locations, |
| 297 // and not to performing any pruning. That is be to done below simultaneously for γ. | 305 // and not to performing any pruning. That is be to done below simultaneously for γ. |
| 298 if config.insertion.merge_now(&state) { | 306 if config.insertion.merge_now(&state) { |
| 311 } | 319 } |
| 312 } | 320 } |
| 313 | 321 |
| 314 γ.prune_compat(&mut μ, &mut stats); | 322 γ.prune_compat(&mut μ, &mut stats); |
| 315 | 323 |
| 316 // Do dual update | 324 for _ in 0..config.insertion.extra_weight_optimisation_steps { |
| 317 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] | 325 // Do extra weight optimisation step heuristic |
| 318 z.axpy(1.0 + ω, &z_new, -ω); | 326 let Pair(v, mut _z_tmp) = f.differential(Pair(&μ, &z)); |
| 319 //z = (z - &z_new) * (-ω) + &z_new; | 327 τv̆ = v * τ; |
| 320 opKz.gemv(&mut y, σ_d, z, 1.0); | 328 prox_penalty.reweigh( |
| 321 y = starH.prox(σ_d, y); | 329 &mut μ, |
| 322 z = z_new; | 330 &mut τv̆, |
| 331 τ, | |
| 332 ε, | |
| 333 &config.insertion, | |
| 334 ®, | |
| 335 &state, | |
| 336 &mut stats, | |
| 337 )?; | |
| 338 | |
| 339 // opKz_adj.apply_add(&mut z_tmp, &y); | |
| 340 // z_tmp.axpy(1.0, &z, -σ_p); | |
| 341 // let z_new = fnR.prox(σ_p, z_tmp); | |
| 342 // z.axpy(1.0 + ω, &z_new, -ω); | |
| 343 // //z = (z - &z_new) * (-ω) + &z_new; | |
| 344 // opKz.gemv(&mut y, σ_d, z, 1.0); | |
| 345 // y = starH.prox(σ_d, y); | |
| 346 // z = z_new; | |
| 347 } | |
| 348 | |
| 349 if config.insertion.extra_weight_optimisation_steps > 0 { | |
| 350 γ.prune_compat(&mut μ, &mut stats); | |
| 351 } | |
| 323 | 352 |
| 324 // Update step length parameters | 353 // Update step length parameters |
| 325 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); | 354 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); |
| 326 | 355 |
| 327 // Give statistics if requested | 356 // Give statistics if requested |
| 328 let iter = state.iteration(); | 357 let iter = state.iteration(); |
| 329 stats.this_iters += 1; | 358 stats.this_iters += 1; |
| 330 | 359 |
| 331 state.if_verbose(|| { | 360 state.if_verbose(|| { |
| 332 plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv̆), &μ); | 361 plotter.plot_spikes(iter, None, Some(&τv̆), &μ); |
| 333 full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) | 362 full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) |
| 334 }); | 363 }); |
| 335 | 364 |
| 336 // Update main tolerance for next iteration | 365 // Update main tolerance for next iteration |
| 337 ε = tolerance.update(ε, iter); | 366 ε = tolerance.update(ε, iter); |
