Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_pdps.rs

src/forward_pdps.rs

branch
dev
changeset 66
fe47ad484deb
parent 65
ec5f160c413b
equal deleted inserted replaced
65:ec5f160c413b 66:fe47ad484deb
150 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new()); 150 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
151 151
152 // Set up parameters 152 // Set up parameters
153 let bigM = 0.0; //opKμ.adjoint_product_bound(prox_penalty).unwrap().sqrt(); 153 let bigM = 0.0; //opKμ.adjoint_product_bound(prox_penalty).unwrap().sqrt();
154 let nKz = opKz.opnorm_bound(L2, L2)?; 154 let nKz = opKz.opnorm_bound(L2, L2)?;
155 let is_fb = nKz == 0.0;
155 let idOpZ = IdOp::new(); 156 let idOpZ = IdOp::new();
156 let opKz_adj = opKz.adjoint(); 157 let opKz_adj = opKz.adjoint();
157 let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?; 158 let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?;
158 // We need to satisfy 159 // We need to satisfy
159 // 160 //
162 // with 1 > σ_p L_z and 1 > τ L. 163 // with 1 > σ_p L_z and 1 > τ L.
163 // 164 //
164 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition 165 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition
165 // ^^^^^ < 1. then we solve τ from the rest. 166 // ^^^^^ < 1. then we solve τ from the rest.
166 // If opKZ is the zero operator, then we set σ_d = 0 for τ to be calculated correctly below. 167 // If opKZ is the zero operator, then we set σ_d = 0 for τ to be calculated correctly below.
167 let σ_d = if nKz == 0.0 { 0.0 } else { config.σd0 / nKz }; 168 let σ_d = if is_fb { 0.0 } else { config.σd0 / nKz };
168 let σ_p = config.σp0 / (l_z + config.σd0 * nKz); 169 let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
169 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0} 170 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0}
170 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L) 171 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L)
171 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0}) 172 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0})
172 let φ = 1.0 - config.σp0; 173 let φ = 1.0 - config.σp0;
222 // Merge spikes. 223 // Merge spikes.
223 // This crucially expects the merge routine to be stable with respect to spike locations, 224 // This crucially expects the merge routine to be stable with respect to spike locations,
224 // and not to performing any pruning. That is be to done below simultaneously for γ. 225 // and not to performing any pruning. That is be to done below simultaneously for γ.
225 let ins = &config.insertion; 226 let ins = &config.insertion;
226 if ins.merge_now(&state) { 227 if ins.merge_now(&state) {
227 stats.merged += 228 stats.merged += prox_penalty.merge_spikes(
228 prox_penalty.merge_spikes_no_fitness(&mut μ, &mut τv, &μ_base, τ, ε, ins, &reg); 229 &mut μ,
230 &mut τv,
231 &μ_base,
232 τ,
233 ε,
234 ins,
235 &reg,
236 is_fb.then_some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))),
237 );
229 } 238 }
230 239
231 // Prune spikes with zero weight. 240 // Prune spikes with zero weight.
232 stats.pruned += prune_with_stats(&mut μ); 241 stats.pruned += prune_with_stats(&mut μ);
233 242

Mercurial Repository: pointsource_algs

mercurial