pointsource_algs: comparison src/sliding

-:c5d8bd1a7728
+:6316d68b58af
 use std::iter::Iterator;
 use alg_tools::iterate::AlgIteratorFactory;
 use alg_tools::euclidean::Euclidean;
 use alg_tools::mapping::{Mapping, DifferentiableRealMapping, Instance};
-use alg_tools::norms::Norm;
+use alg_tools::norms::{Norm, Dist};
 use alg_tools::direct_product::Pair;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::linops::{
 BoundedLinear, AXPY, GEMV, Adjointable, IdOp,
 };
 TransportConfig,
 TransportStepLength,
 initial_transport,
 aposteriori_transport,
 };
-use crate::dataterm::{calculate_residual, calculate_residual2};
+use crate::dataterm::{
+calculate_residual2,
+calculate_residual,
+};
 /// Settings for [`pointsource_sliding_pdps_pair`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
 pub struct SlidingPDPSConfig<F : Float> {
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float> Default for SlidingPDPSConfig<F> {
 fn default() -> Self {
-let τ0 = 0.99;
 SlidingPDPSConfig {
-τ0,
+τ0 : 0.99,
-σd0 : 0.1,
+σd0 : 0.05,
 σp0 : 0.99,
-transport : Default::default(),
+transport : TransportConfig { θ0 : 0.1, ..Default::default()},
 insertion : Default::default()
 }
 }
 }
 + GEMV<F, Z>
 + Adjointable<Z, Y, AdjointCodomain = Z>,
 for<'b> KOpZ::Adjoint<'b> : GEMV<F, Y>,
 Y : AXPY<F> + Euclidean<F, Output=Y> + Clone + ClosedAdd,
 for<'b> &'b Y : Instance<Y>,
-Z : AXPY<F, Owned=Z> + Euclidean<F, Output=Z> + Clone + Norm<F, L2>,
+Z : AXPY<F, Owned=Z> + Euclidean<F, Output=Z> + Clone + Norm<F, L2> + Dist<F, L2>,
 for<'b> &'b Z : Instance<Z>,
 R : Prox<Z, Codomain=F>,
 H : Conjugable<Y, F, Codomain=F>,
 for<'b> H::Conjugate<'b> : Prox<Y>,
 {
 l: ℓ_v0 * b.norm2(),
 max_transport : 0.0,
 g : calculate_θ
 },
 None => TransportStepLength::FullyAdaptive{
-l : 0.0,
+l : F::EPSILON,
 max_transport : 0.0,
 g : calculate_θ
 },
 };
 // Acceleration is not currently supported
 // Run the algorithm
 for state in iterator.iter_init(|| full_stats(&residual, &μ, &z, ε, stats.clone())) {
 // Calculate initial transport
 let Pair(v, _) = opA.preadjoint().apply(&residual);
 //opKμ.preadjoint().apply_add(&mut v, y);
-let z_base = z.clone();
 // We want to proceed as in Example 4.12 but with v and v̆ as in §5.
 // With A(ν, z) = A_μ ν + A_z z, following Example 5.1, we have
 // P_ℳ[F'(ν, z) + Ξ(ν, z, y)]= A_ν^*[A_ν ν + A_z z] + K_μ ν = A_ν^*A(ν, z) + K_μ ν,
 // where A_ν^* becomes a multiplier.
 // This is much easier with K_μ = 0, which is the only reason why are enforcing it.
 v, &config.transport,
 );
 // Solve finite-dimensional subproblem several times until the dual variable for the
 // regularisation term conforms to the assumptions made for the transport above.
-let (maybe_d, _within_tolerances, Pair(τv̆, τz̆)) = 'adapt_transport: loop {
+let (maybe_d, _within_tolerances, mut τv̆, z_new) = 'adapt_transport: loop {
 // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b)
 let residual_μ̆ = calculate_residual2(Pair(&γ1, &z),
 Pair(&μ_base_minus_γ0, &zero_z),
 opA, b);
-let mut τv̆z = opA.preadjoint().apply(residual_μ̆ * τ);
+let Pair(mut τv̆, τz̆) = opA.preadjoint().apply(residual_μ̆ * τ);
 // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0);
 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
 let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh(
-&mut μ, &mut τv̆z.0, &γ1, Some(&μ_base_minus_γ0),
+&mut μ, &mut τv̆, &γ1, Some(&μ_base_minus_γ0),
 τ, ε, &config.insertion,
 &reg, &state, &mut stats,
 );
+// Do z variable primal update here to able to estimate B_{v̆^k-v^{k+1}}
+let mut z_new = τz̆;
+opKz.adjoint().gemv(&mut z_new, -σ_p, &y, -σ_p/τ);
+z_new = fnR.prox(σ_p, z_new + &z);
 // A posteriori transport adaptation.
-// TODO: this does not properly treat v^{k+1} - v̆^k that depends on z^{k+1}!
 if aposteriori_transport(
 &mut γ1, &mut μ, &mut μ_base_minus_γ0, &μ_base_masses,
+Some(z_new.dist(&z, L2)),
 ε, &config.transport
 ) {
-break 'adapt_transport (maybe_d, within_tolerances, τv̆z)
+break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new)
 }
 };
 stats.untransported_fraction = Some({
 assert_eq!(μ_base_masses.len(), γ1.len());
 assert_eq!(μ_base_masses.len(), γ1.len());
 let (a, b) = stats.transport_error.unwrap_or((0.0, 0.0));
 (a + μ.dist_matching(&γ1), b + γ1.norm(Radon))
 });
-// // Merge spikes.
+// Merge spikes.
-// // This expects the prune below to prune γ.
+// This crucially expects the merge routine to be stable with respect to spike locations,
-// // TODO: This may not work correctly in all cases.
+// and not to performing any pruning. That is be to done below simultaneously for γ.
-// let ins = &config.insertion;
+let ins = &config.insertion;
-// if ins.merge_now(&state) {
+if ins.merge_now(&state) {
-//     if let SpikeMergingMethod::None = ins.merging {
+stats.merged += prox_penalty.merge_spikes_no_fitness(
-//     } else {
+&mut μ, &mut τv̆, &γ1, Some(&μ_base_minus_γ0), τ, ε, ins, &reg,
-//         stats.merged += μ.merge_spikes(ins.merging, |μ_candidate| {
+//Some(|μ̃ : &RNDM<F, N>| calculate_residual(Pair(μ̃, &z), opA, b).norm2_squared_div2()),
-//             let ν = μ_candidate.sub_matching(&γ1)-&μ_base_minus_γ0;
+);
-//             let mut d = &τv̆ + op𝒟.preapply(ν);
+}
-//             reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, ins)
-//         });
-//     }
-// }
 // Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the
 // latter needs to be pruned when μ is.
 // TODO: This could do with a two-vector Vec::retain to avoid copies.
 let μ_new = DiscreteMeasure::from_iter(μ.iter_spikes().filter(|δ| δ.α != F::ZERO).cloned());
 γ1.prune_by(|_| μ_iter.next().unwrap().α != F::ZERO);
 stats.pruned += μ.len() - μ_new.len();
 μ = μ_new;
 }
-// Do z variable primal update
-z.axpy(-σ_p/τ, τz̆, 1.0); // TODO: simplify nasty factors
-opKz.adjoint().gemv(&mut z, -σ_p, &y, 1.0);
-z = fnR.prox(σ_p, z);
 // Do dual update
 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0);    // y = y + σ_d K[(1+ω)(μ,z)^{k+1}]
-opKz.gemv(&mut y, σ_d*(1.0 + ω), &z, 1.0);
+opKz.gemv(&mut y, σ_d*(1.0 + ω), &z_new, 1.0);
 // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b
-opKz.gemv(&mut y, -σ_d*ω, z_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b
+opKz.gemv(&mut y, -σ_d*ω, z, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b
 y = starH.prox(σ_d, y);
+z = z_new;
 // Update residual
 residual = calculate_residual(Pair(&μ, &z), opA, b);
 // Update step length parameters
 (opA.apply(Pair(μ̃, &z))-b).norm2_squared_div2()
 //+ fnR.apply(z) + reg.apply(μ)
 + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z))
 };
-μ.merge_spikes_fitness(config.insertion.merging, fit, |&v| v);
+μ.merge_spikes_fitness(config.insertion.final_merging_method(), fit, |&v| v);
 μ.prune();
 Pair(μ, z)
 }

comparison: src/sliding_pdps.rs

src/sliding_pdps.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_pdps.rs

src/sliding_pdps.rs