diff -r f91f382ba4ae -r b9818eb24b47 src/sliding_pdps.rs
--- 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>,