Mercurial > repos > pointsource_algs / file comparison
comparison: src/sliding_pdps.rs
src/sliding_pdps.rs
- changeset 97
- f91f382ba4ae
- parent 96
- 5ba25d740dfc
- child 98
- b9818eb24b47
equal
deleted
inserted
replaced
| 96:5ba25d740dfc | 97:f91f382ba4ae |
|---|---|
| 253 | 253 |
| 254 let mut attempts = 0; | 254 let mut attempts = 0; |
| 255 | 255 |
| 256 // Solve finite-dimensional subproblem several times until the dual variable for the | 256 // Solve finite-dimensional subproblem several times until the dual variable for the |
| 257 // regularisation term conforms to the assumptions made for the transport above. | 257 // regularisation term conforms to the assumptions made for the transport above. |
| 258 let (mut τv̆, z_new, μ̆) = 'adapt_transport: loop { | 258 let (mut τv̆, mut z̆, μ̆) = 'adapt_transport: loop { |
| 259 // Set initial guess for μ=μ^{k+1}. | 259 // Set initial guess for μ=μ^{k+1}. |
| 260 γ.μ̆_into(&mut μ); | 260 γ.μ̆_into(&mut μ); |
| 261 let μ̆ = μ.clone(); | 261 let μ̆ = μ.clone(); |
| 262 | 262 |
| 263 // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) | 263 // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) |
| 264 let Pair(v̆, mut z̆) = f.differential(Pair(&μ̆, &z)); | 264 let Pair(v̆, z̆) = f.differential(Pair(&μ̆, &z)); |
| 265 let mut τv̆ = v̆ * τ; | 265 let mut τv̆ = v̆ * τ; |
| 266 // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0); | 266 // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0); |
| 267 | 267 |
| 268 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. | 268 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. |
| 269 prox_penalty.insert_and_reweigh( | 269 prox_penalty.insert_and_reweigh( |
| 274 &config.insertion, | 274 &config.insertion, |
| 275 ®, | 275 ®, |
| 276 &state, | 276 &state, |
| 277 &mut stats, | 277 &mut stats, |
| 278 )?; | 278 )?; |
| 279 | |
| 280 // Do z variable primal update here to able to estimate B_{v̆^k-v^{k+1}} | |
| 281 opKz_adj.apply_add(&mut z̆, &y); | |
| 282 z̆.axpy(1.0, &z, -σ_p); | |
| 283 let z_new = fnR.prox(σ_p, z̆); | |
| 284 //let z_new = fnR.prox(σ_p, z̆ * (-σ_p) + &z); | |
| 285 | 279 |
| 286 // A posteriori transport adaptation. | 280 // A posteriori transport adaptation. |
| 287 if prox_penalty.aposteriori_transport( | 281 if prox_penalty.aposteriori_transport( |
| 288 &mut γ, | 282 &mut γ, |
| 289 &μ, | 283 &μ, |
| 297 reg, | 291 reg, |
| 298 &config.transport, | 292 &config.transport, |
| 299 &config.insertion.refinement, | 293 &config.insertion.refinement, |
| 300 &mut attempts, | 294 &mut attempts, |
| 301 ) { | 295 ) { |
| 302 break 'adapt_transport (τv̆, z_new, μ̆); | 296 break 'adapt_transport (τv̆, z̆, μ̆); |
| 303 } | 297 } |
| 304 }; | 298 }; |
| 305 | 299 |
| 306 γ.get_transport_stats(&mut stats, &μ); | 300 γ.get_transport_stats(&mut stats, &μ); |
| 307 | 301 |
| 302 // Do z variable primal update | |
| 303 opKz_adj.apply_add(&mut z̆, &y); | |
| 304 z̆.axpy(1.0, &z, -σ_p); | |
| 305 let z_new = fnR.prox(σ_p, z̆); | |
| 306 //let z_new = fnR.prox(σ_p, z̆ * (-σ_p) + &z); | |
| 308 // Do dual update | 307 // Do dual update |
| 309 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] | 308 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] |
| 310 z.axpy(1.0 + ω, &z_new, -ω); | 309 z.axpy(1.0 + ω, &z_new, -ω); |
| 311 //z = (z - &z_new) * (-ω) + &z_new; | 310 //z = (z - &z_new) * (-ω) + &z_new; |
| 312 opKz.gemv(&mut y, σ_d, z, 1.0); | 311 opKz.gemv(&mut y, σ_d, z, 1.0); |
