pointsource_algs: comparison src/sliding

-:f91f382ba4ae
+:b9818eb24b47
 use alg_tools::linops::{
 BoundedLinear, IdOp, SimplyAdjointable, StaticEuclideanOriginGenerator, ZeroOp, AXPY, GEMV,
 };
 use alg_tools::mapping::DifferentiableMapping;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::{Dist, L2};
+use alg_tools::norms::{Dist, Norm, L2};
 use anyhow::ensure;
-use measures::DeltaMeasure;
+use measures::{DeltaMeasure, Radon};
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
 /// Settings for [`pointsource_sliding_pdps_pair`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 U: Space + Clone + Dist<L2, F>,
 for<'a> &'a U: Instance<U>,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>
 + BoundedCurvature<F>,
-S: ClosedMul<F>,
+S: ClosedMul<F> + MinMaxMapping<U, F>,
 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
-DiscreteMeasure<U, F>: SpikeMerging<F>,
+DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
 Reg: SlidingRegTerm<U, F>,
 P: TransportProxPenalty<U, S, Reg, F>,
 // KOpM : Linear<DiscreteMeasure<U, F>, Codomain=Y>
 //     + GEMV<F, DiscreteMeasure<U, F>>
 //     + Preadjointable<
 let nKz = opKz.opnorm_bound(L2, L2)?;
 let is_fb = nKz == 0.0;
 let idOpZ = IdOp::new();
 let opKz_adj = opKz.adjoint();
 let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?;
+let l_z_prime = if l_z == 0.0 {
+1.0
+} else {
+ensure!(config.σp0 < 1.0);
+l_z
+};
+let l_prime = if l == 0.0 {
+1.0
+} else {
+ensure!(config.τ0 < 1.0);
+l
+};
 let (τ, σ_p, σ_d, κ) = if is_fb {
-(config.τ0 / l, config.σp0 / l_z, 0.0, 0.0)
+(config.τ0 / l_prime, config.σp0 / l_z_prime, 0.0, 0.0)
 } else {
 ensure!(
 //config.τ0 < 1.0 &&
 config.σp0 < 1.0 && config.σd0 > 0.0 && config.σp0 * config.σd0 <= 1.0,
 "Invalid step length parameters"
 //
 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition
 // ^^^^^ < 1. then we solve τ from  the rest.
 // If opKZ is the zero operator, then we set σ_d = 0 for τ to be calculated correctly below.
 let σ_d = config.σd0 / nKz;
-let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
+let σ_p = config.σp0 / (l_z_prime + config.σd0 * nKz);
 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0}
 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L)
 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0})
 let φ = 1.0 - config.σp0;
-let a = 1.0 - σ_p * l_z;
+let a = 1.0 - σ_p * l_z_prime;
 let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l);
-let ψ = 1.0 - τ * l;
+let ψ = 1.0 - τ * l_prime;
 let β = σ_p * config.σd0 * nKz / a; // σ_p * σ_d * (nKz * nK_z) / a;
 ensure!(β < 1.0);
 // Now we need κ‖K_μ(π_♯^1 - π_♯^0)γ‖^2 ≤ (1/θ - τ[ℓ_F + ℓ]) ∫ c_2 dγ for κ defined as:
 //  The factor two in the manuscript disappears due to the definition of 𝚹 being
 // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2.
 let mut ε = tolerance.initial();
 let starH = fnH.conjugate();
 // Statistics
-let full_stats = |μ: &DiscreteMeasure<U, F>, z: &Z, ε, stats| IterInfo {
+let full_stats = |μ: &DiscreteMeasure<U, F>, v: &mut S, z: &Z, dz: Z, ε, stats| {
-value: f.apply(Pair(μ, z))
+let discrepancy = reg.discrepancy(v, τ, ε, μ, &config.insertion);
-+ fnR.apply(z)
+let aux_discrepancy = is_fb.then(|| z.dist2(fnR.prox(σ_p, dz * (-σ_p) + z)));
-+ reg.apply(μ)
+let joint_discrepancy = aux_discrepancy
-+ fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)),
+.map(|d| (discrepancy.powi(2) + τ * (2.0 / σ_p - l_z) * d.powi(2)).sqrt());
-n_spikes: μ.len(),
+IterInfo {
-ε,
+value: f.apply(Pair(μ, z))
-// postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
++ fnR.apply(z)
-..stats
++ reg.apply(μ)
++ fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)),
+n_spikes: μ.len(),
+discrepancy: Some(discrepancy),
+μ_norm: μ.norm(Radon),
+aux_discrepancy,
+joint_discrepancy,
+ε,
+// postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
+..stats
+}
 };
 let mut stats = IterInfo::new();
+let Pair(mut v, mut dz) = f.differential(Pair(&μ, &z));
 // Run the algorithm
-for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) {
+for state in iterator.iter_init(|| full_stats(&μ, &mut v, &z, dz, ε, stats.clone())) {
-let Pair(mut v, _) = f.differential(Pair(&μ, &z));
 //opKμ.preadjoint().apply_add(&mut v, y);
 // 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.
 // Give statistics if requested
 let iter = state.iteration();
 stats.this_iters += 1;
+Pair(v, dz) = f.differential(Pair(&μ, &z));
 state.if_verbose(|| {
 plotter.plot_spikes(iter, None, Some(&τv̆), &μ);
-full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new()))
+full_stats(
+&μ,
+&mut v,
+&z,
+dz,
+ε,
+std::mem::replace(&mut stats, IterInfo::new()),
+)
 });
 // Update main tolerance for next iteration
 ε = tolerance.update(ε, iter);
 }
 U: Space + Clone + Dist<L2, F>,
 for<'a> &'a U: Instance<U>,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>
 + BoundedCurvature<F>,
-S: ClosedMul<F>,
+S: ClosedMul<F> + MinMaxMapping<U, F>,
-DiscreteMeasure<U, F>: SpikeMerging<F>,
+DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
 Reg: SlidingRegTerm<U, F>,
 P: TransportProxPenalty<U, S, Reg, F>,
 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
 Z: ClosedEuclidean<F> + AXPY + Clone,
 for<'b> &'b Z: Instance<Z>,
 U: Space + Clone + Dist<L2, F>,
 for<'a> &'a U: Instance<U>,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>
 + BoundedCurvature<F>,
-S: ClosedMul<F>,
+S: ClosedMul<F> + MinMaxMapping<U, F>,
 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
-DiscreteMeasure<U, F>: SpikeMerging<F>,
+DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
 Reg: SlidingRegTerm<U, F>,
 P: TransportProxPenalty<U, S, Reg, F>,
 KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
 + GEMV<F, Z>
 + SimplyAdjointable<Z, Y, AdjointCodomain = Z>,
 for<'a> &'a U: Instance<U>,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>
 + BoundedCurvature<F>,
 S: DifferentiableMapping<U, Codomain = F> + ClosedMul<F> + MinMaxMapping<U, F>,
-DiscreteMeasure<U, F>: SpikeMerging<F>,
+DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
 Reg: SlidingRegTerm<U, F>,
 P: TransportProxPenalty<U, S, Reg, F>,
 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
 Z: ClosedEuclidean<F> + AXPY + Clone,
 for<'b> &'b Z: Instance<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