--- a/src/sliding_pdps.rs Fri Aug 14 10:38:27 2026 -0500
+++ b/src/sliding_pdps.rs Fri Oct 02 12:03:25 2026 -0500
@@ -25,9 +25,9 @@
};
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};
@@ -108,9 +108,9 @@
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>
@@ -159,9 +159,21 @@
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 &&
@@ -178,14 +190,14 @@
// ^^^^^ < 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:
@@ -214,22 +226,32 @@
let starH = fnH.conjugate();
// Statistics
- let full_stats = |μ: &DiscreteMeasure<U, F>, z: &Z, ε, stats| IterInfo {
- value: f.apply(Pair(μ, z))
- + fnR.apply(z)
- + reg.apply(μ)
- + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)),
- n_spikes: μ.len(),
- ε,
- // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
- ..stats
+ let full_stats = |μ: &DiscreteMeasure<U, F>, v: &mut S, z: &Z, dz: Z, ε, stats| {
+ let discrepancy = reg.discrepancy(v, τ, ε, μ, &config.insertion);
+ let aux_discrepancy = is_fb.then(|| z.dist2(fnR.prox(σ_p, dz * (-σ_p) + z)));
+ let joint_discrepancy = aux_discrepancy
+ .map(|d| (discrepancy.powi(2) + τ * (2.0 / σ_p - l_z) * d.powi(2)).sqrt());
+ IterInfo {
+ value: f.apply(Pair(μ, z))
+ + fnR.apply(z)
+ + 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())) {
- let Pair(mut v, _) = f.differential(Pair(&μ, &z));
-
+ for state in iterator.iter_init(|| full_stats(&μ, &mut v, &z, dz, ε, stats.clone())) {
//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
@@ -373,9 +395,18 @@
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
@@ -418,8 +449,8 @@
I: AlgIteratorFactory<IterInfo<F>>,
Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>
+ BoundedCurvature<F>,
- S: ClosedMul<F>,
- DiscreteMeasure<U, F>: SpikeMerging<F>,
+ S: ClosedMul<F> + MinMaxMapping<U, 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>,
@@ -488,9 +519,9 @@
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>
@@ -543,7 +574,7 @@
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>,