Mercurial > repos > pointsource_algs / changeset

--- a/src/fb.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/fb.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -79,17 +79,19 @@
 
 use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::measures::merging::SpikeMerging;
-use crate::measures::DiscreteMeasure;
+use crate::measures::{DiscreteMeasure, Radon};
 use crate::plot::Plotter;
 pub use crate::prox_penalty::{InsertionConfig, ProxPenalty, StepLengthBound};
 use crate::regularisation::RegTerm;
 use crate::types::*;
+use alg_tools::bounds::MinMaxMapping;
 use alg_tools::error::DynResult;
 use alg_tools::instance::Instance;
 use alg_tools::iterate::AlgIteratorFactory;
 use alg_tools::linops::Space;
 use alg_tools::mapping::DifferentiableMapping;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::norms::Norm;
 use colored::Colorize;
 use measures::DeltaMeasure;
 use numeric_literals::replace_float_literals;
@@ -168,9 +170,9 @@
     U: Space + Clone,
     F: Float + ToNalgebraRealField,
     I: AlgIteratorFactory<IterInfo<F>>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Dat: DifferentiableMapping<DiscreteMeasure<U, F>, Codomain = F>,
-    Dat::DerivativeDomain: ClosedMul<F>,
+    Dat::DerivativeDomain: ClosedMul<F> + MinMaxMapping<U, F>,
     Reg: RegTerm<U, F>,
     P: ProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
     Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
@@ -187,23 +189,26 @@
     let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
 
     // Statistics
-    let full_stats = |μ: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
-        value: f.apply(μ) + reg.apply(μ),
-        n_spikes: μ.len(),
-        ε,
-        //postprocessing: config.postprocessing.then(|| μ.clone()),
-        ..stats
-    };
+    let full_stats =
+        |μ: &DiscreteMeasure<U, F>, τv: &mut Dat::DerivativeDomain, ε, stats| IterInfo {
+            value: f.apply(μ) + reg.apply(μ),
+            n_spikes: μ.len(),
+            μ_norm: μ.norm(Radon),
+            ε,
+            discrepancy: Some(reg.discrepancy(τv, τ, ε, μ, config)),
+            //postprocessing: config.postprocessing.then(|| μ.clone()),
+            ..stats
+        };
     let mut stats = IterInfo::new();
 
+    // Calculate smooth part of surrogate model.
+    let mut τv = f.differential(&μ) * τ;
+
     // Run the algorithm
-    for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
+    for state in iterator.iter_init(|| full_stats(&μ, &mut τv, ε, stats.clone())) {
         let maybe_μ_base = config.if_merge_now_check(&state, || μ.clone());
         let μ_base_len = μ.len();
 
-        // Calculate smooth part of surrogate model.
-        let mut τv = f.differential(&μ) * τ;
-
         // Do spike insertiona nd finite-dimensional weight optimisation
         prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
 
@@ -256,10 +261,18 @@
         let iter = state.iteration();
         stats.this_iters += 1;
 
+        // Calculate smooth part of surrogate model.
+        τv = f.differential(&μ) * τ;
+
         // Give statistics if needed
         state.if_verbose(|| {
             plotter.plot_spikes(iter, None, Some(&τv), &μ);
-            full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
+            full_stats(
+                &μ,
+                &mut τv,
+                ε,
+                std::mem::replace(&mut stats, IterInfo::new()),
+            )
         });
 
         // Update main tolerance for next iteration
@@ -284,9 +297,9 @@
     U: Space + Clone,
     F: Float + ToNalgebraRealField,
     I: AlgIteratorFactory<IterInfo<F>>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Dat: DifferentiableMapping<DiscreteMeasure<U, F>, Codomain = F>,
-    Dat::DerivativeDomain: ClosedMul<F>,
+    Dat::DerivativeDomain: ClosedMul<F> + MinMaxMapping<U, F>,
     Reg: RegTerm<U, F>,
     P: ProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
     Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
@@ -330,7 +343,7 @@
     I: AlgIteratorFactory<IterInfo<F>>,
     DiscreteMeasure<U, F>: SpikeMerging<F>,
     Dat: DifferentiableMapping<DiscreteMeasure<U, F>, Codomain = F>,
-    Dat::DerivativeDomain: ClosedMul<F>,
+    Dat::DerivativeDomain: ClosedMul<F> + MinMaxMapping<U, F>,
     Reg: RegTerm<U, F>,
     P: ProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
     Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
@@ -350,17 +363,21 @@
     let mut warned_merging = false;
 
     // Statistics
-    let full_stats = |ν: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
+    let full_stats = |ν: &DiscreteMeasure<U, F>,
+                      maybe_τv: Option<&mut Dat::DerivativeDomain>,
+                      ε,
+                      stats| IterInfo {
         value: f.apply(ν) + reg.apply(ν),
         n_spikes: ν.len(),
         ε,
+        discrepancy: maybe_τv.map(|τv| reg.discrepancy(τv, τ, ε, ν, config)),
         // postprocessing: config.postprocessing.then(|| ν.clone()),
         ..stats
     };
     let mut stats = IterInfo::new();
 
     // Run the algorithm
-    for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
+    for state in iterator.iter_init(|| full_stats(&μ, None, ε, stats.clone())) {
         let μ_base_len = μ.len();
         // Calculate smooth part of surrogate model.
         let mut τv = f.differential(&μ) * τ;
@@ -410,7 +427,12 @@
         // Give statistics if needed
         state.if_verbose(|| {
             plotter.plot_spikes(iter, None, Some(&τv), &μ_prev);
-            full_stats(&μ_prev, ε, std::mem::replace(&mut stats, IterInfo::new()))
+            full_stats(
+                &μ_prev,
+                Some(&mut τv),
+                ε,
+                std::mem::replace(&mut stats, IterInfo::new()),
+            )
         });
 
         // Update main tolerance for next iteration

--- a/src/forward_pdps.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/forward_pdps.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -11,6 +11,7 @@
 use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair};
 use crate::regularisation::RegTerm;
 use crate::types::*;
+use alg_tools::bounds::MinMaxMapping;
 use alg_tools::convex::{Conjugable, Prox, Zero};
 use alg_tools::direct_product::Pair;
 use alg_tools::error::DynResult;
@@ -22,11 +23,12 @@
 };
 use alg_tools::mapping::{DifferentiableMapping, Instance};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::L2;
+use alg_tools::norms::{Norm, L2};
 use anyhow::ensure;
-use measures::DeltaMeasure;
+use measures::{DeltaMeasure, Radon};
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
+//use std::ops::Mul;
 
 /// Settings for [`pointsource_forward_pdps_pair`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
@@ -123,9 +125,9 @@
     U: Space + Clone,
     I: AlgIteratorFactory<IterInfo<F>>,
     Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
-    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: RegTerm<U, F>,
     P: ProxPenalty<U, S, Reg, F>,
     KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
@@ -134,8 +136,8 @@
     KOpZ::SimpleAdjoint: GEMV<F, Y>,
     Y: ClosedEuclidean<F>,
     for<'b> &'b Y: Instance<Y>,
-    Z: ClosedEuclidean<F> + AXPY,
-    for<'b> &'b Z: Instance<Z>,
+    Z: ClosedEuclidean<F> + AXPY + Clone,
+    for<'b> &'b Z: Instance<Z>, // compiler overflow  + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     H: Conjugable<Y, F, Codomain = F>,
     for<'b> H::Conjugate<'b>: Prox<Y>,
@@ -158,9 +160,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)
+        (config.τ0 / l_prime, config.σp0 / l_z_prime, 0.0)
     } else {
         ensure!(
             //config.τ0 < 1.0 &&
@@ -168,10 +182,10 @@
             "Invalid step length parameters"
         );
         let σ_d = config.σd0 / nKz;
-        let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
+        let σ_p = config.σp0 / (l_z_prime + config.σd0 * nKz);
         let φ = 1.0 - config.σp0;
-        let a = 1.0 - σ_p * l_z;
-        let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l);
+        let a = 1.0 - σ_p * l_z_prime;
+        let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l_prime);
         (τ, σ_p, σ_d)
     };
 
@@ -187,24 +201,37 @@
     let starH = fnH.conjugate();
 
     // Statistics
-    let full_stats = |μ: &DiscreteMeasure<U, F>, z: &Z, ε, stats| IterInfo {
-        value: f.apply(Pair(μ, z))
+    let full_stats = |μ: &DiscreteMeasure<U, F>, v: &mut S, z: &Z, dz: &Z, ε, stats| {
+        let 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
+            + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z));
+        let discrepancy = reg.discrepancy(v, τ, ε, μ, &config.insertion);
+        let aux_discrepancy = is_fb.then(|| z.dist2(fnR.prox(σ_p, Z::clone(dz) * (-σ_p) + z)));
+        let joint_discrepancy = aux_discrepancy
+            .map(|d| (discrepancy.powi(2) + τ * (2.0 / σ_p - l_z) * d.powi(2)).sqrt());
+        IterInfo {
+            value,
+            n_spikes: μ.len(),
+            discrepancy: Some(discrepancy),
+            μ_norm: μ.norm(Radon),
+            // compiler overflow with the necesssary trait bounds:
+            //aux_discrepancy: is_fb.then(|| z.dist2(fnR.prox(σ_p, dz * (-σ_p) + z))),
+            aux_discrepancy,
+            joint_discrepancy,
+            ε,
+            // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
+            ..stats
+        }
     };
     let mut stats = IterInfo::new();
 
-    let zy_step = |mut z, mut y, mut z_tmp| {
+    let zy_step = |mut z, mut y, mut dz| {
         // Do z variable primal update
-        opKz_adj.apply_add(&mut z_tmp, &y);
-        z_tmp.axpy(1.0, &z, -σ_p);
-        let z_new = fnR.prox(σ_p, z_tmp);
-        //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z);
+        opKz_adj.apply_add(&mut dz, &y);
+        dz.axpy(1.0, &z, -σ_p);
+        let z_new = fnR.prox(σ_p, dz);
+        //let z_new = fnR.prox(σ_p, dz * (-σ_p) + &z);
         // Do dual update
         z.axpy(1.0 + ω, &z_new, -ω);
         //z = (z - &z_new) * (-ω) + &z_new;
@@ -213,11 +240,13 @@
         (z_new, y)
     };
 
+    // Calculate differential of smooth part
+    let Pair(mut v, mut dz) = f.differential(Pair(&μ, &z));
+    let mut τv = v * τ;
+
     // 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())) {
         // Calculate initial transport
-        let Pair(v, z_tmp) = f.differential(Pair(&μ, &z));
-        let mut τv = v * τ;
         let μ_base = μ.clone();
         let μ_base_len = μ_base.len();
 
@@ -233,7 +262,7 @@
             &mut stats,
         )?;
 
-        (z, y) = zy_step(z, y, z_tmp);
+        (z, y) = zy_step(z, y, dz);
 
         stats.inserted += μ.len() - μ_base.len();
 
@@ -263,7 +292,7 @@
 
             for _ in 0..config.insertion.extra_weight_optimisation_steps {
                 // Do extra weight optimisation step heuristic
-                let Pair(v, mut _z_tmp) = f.differential(Pair(&μ, &z));
+                let Pair(v, mut _dz) = f.differential(Pair(&μ, &z));
                 τv = v * τ;
                 prox_penalty.reweigh(
                     &mut μ,
@@ -276,9 +305,9 @@
                     &mut stats,
                 )?;
 
-                // opKz_adj.apply_add(&mut z_tmp, &y);
-                // z_tmp.axpy(1.0, &z, -σ_p);
-                // let z_new = fnR.prox(σ_p, z_tmp);
+                // opKz_adj.apply_add(&mut dz, &y);
+                // dz.axpy(1.0, &z, -σ_p);
+                // let z_new = fnR.prox(σ_p, dz);
                 // z.axpy(1.0 + ω, &z_new, -ω);
                 // //z = (z - &z_new) * (-ω) + &z_new;
                 // opKz.gemv(&mut y, σ_d, z, 1.0);
@@ -300,9 +329,20 @@
         let iter = state.iteration();
         stats.this_iters += 1;
 
+        // Update differential of smooth part
+        Pair(v, dz) = f.differential(Pair(&μ, &z));
+        τv = v * τ;
+
         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
@@ -354,10 +394,10 @@
     //for<'a> &'a U: Instance<U>,
     I: AlgIteratorFactory<IterInfo<F>>,
     Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
-    S: ClosedMul<F>,
+    S: ClosedMul<F> + MinMaxMapping<U, F>,
     for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
     //Pair<S, Z>: ClosedMul<F>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Reg: RegTerm<U, F>,
     P: ProxPenalty<U, S, Reg, F>,
     KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
@@ -366,8 +406,8 @@
     KOpZ::SimpleAdjoint: GEMV<F, Y>,
     Y: ClosedEuclidean<F>,
     for<'b> &'b Y: Instance<Y>,
-    Z: ClosedEuclidean<F> + AXPY,
-    for<'b> &'b Z: Instance<Z>,
+    Z: ClosedEuclidean<F> + AXPY + Clone, // clone not neede if Mul worked below
+    for<'b> &'b Z: Instance<Z>,           // compiler overflow  + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     H: Conjugable<Y, F, Codomain = F>,
     for<'b> H::Conjugate<'b>: Prox<Y>,
@@ -412,13 +452,13 @@
     for<'a> &'a U: Instance<U>,
     I: AlgIteratorFactory<IterInfo<F>>,
     Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
-    S: ClosedMul<F>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    S: ClosedMul<F> + MinMaxMapping<U, F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Reg: RegTerm<U, F>,
     P: ProxPenalty<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>,
+    for<'b> &'b Z: Instance<Z>, // compiler overflow + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
 {
@@ -464,13 +504,13 @@
     for<'a> &'a U: Instance<U>,
     I: AlgIteratorFactory<IterInfo<F>>,
     Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
-    S: ClosedMul<F>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    S: ClosedMul<F> + MinMaxMapping<U, F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Reg: RegTerm<U, F>,
     P: ProxPenalty<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>,
+    for<'b> &'b Z: Instance<Z>, // compiler overflow  + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
 {

--- a/src/pdps.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/pdps.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -53,8 +53,10 @@
 use alg_tools::linops::{Mapping, Space, AXPY};
 use alg_tools::mapping::{DataTerm, Instance};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::norms::Norm;
 use anyhow::ensure;
 use clap::ValueEnum;
+use measures::Radon;
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
 
@@ -156,7 +158,7 @@
     A: ForwardModel<DiscreteMeasure<U, F>, F>,
     for<'b> &'b A::Observable: Instance<A::Observable>,
     A::Observable: AXPY<Field = F>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Reg: RegTerm<U, F>,
     Phi: Conjugable<A::Observable, F>,
     for<'b> Phi::Conjugate<'b>: Prox<A::Observable>,
@@ -192,6 +194,7 @@
     let full_stats = |μ: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
         value: f.apply(μ) + reg.apply(μ),
         n_spikes: μ.len(),
+        μ_norm: μ.norm(Radon),
         ε,
         // postprocessing: config.postprocessing.then(|| μ.clone()),
         ..stats

--- a/src/regularisation.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/regularisation.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -148,6 +148,20 @@
     where
         M: MinMaxMapping<Domain, F>;
 
+    /// Returns the minimum distance of the subdifferential of $τ[F+R]$ to zero,
+    /// where $F$ is the smooth data term (differential given in `τv/τ`), and $R$ is
+    /// this regulariser.
+    fn discrepancy<M>(
+        &self,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        μ: &DiscreteMeasure<Domain, F>,
+        config: &InsertionConfig<F>,
+    ) -> F
+    where
+        M: MinMaxMapping<Domain, F>;
+
     /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
     ///
     /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
@@ -290,6 +304,40 @@
         }
     }
 
+    fn discrepancy<M>(
+        &self,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        μ: &DiscreteMeasure<Loc<N, F>, F>,
+        config: &InsertionConfig<F>,
+    ) -> F
+    where
+        M: MinMaxMapping<Loc<N, F>, F>,
+    {
+        let refinement_tolerance = ε * config.refinement.tolerance_mult;
+        let τα = self.α()*τ;
+
+        let μ_disrepancy = μ
+            .iter_spikes()
+            .map(|δ| {
+                let v = τv.apply(&δ.x);
+                if δ.α > 0.0 {
+                    (v + τα).abs()
+                } else {
+                    (-v - τα).max(0.0)
+                }
+            })
+            .reduce(F::max)
+            .unwrap_or(0.0);
+
+        let (_, τv_min) = τv.minimise(refinement_tolerance, config.refinement.max_steps);
+        //let rest_discrepancy = (-d_min - self.α()).max(0.0);
+        //return rest_discrepancy.max(μ_disrepancy)
+        // Simplifies to
+        return (-τv_min - τα).max(μ_disrepancy);
+    }
+
     fn verify_merge_candidate<M>(
         &self,
         d: &mut M,
@@ -532,6 +580,54 @@
         }
     }
 
+    fn discrepancy<M>(
+        &self,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        μ: &DiscreteMeasure<Loc<N, F>, F>,
+        config: &InsertionConfig<F>,
+    ) -> F
+    where
+        M: MinMaxMapping<Loc<N, F>, F>,
+    {
+        use std::cmp::Ordering::*;
+        let refinement_tolerance = ε * config.refinement.tolerance_mult;
+        let τα = τ * self.α();
+
+        let threshold = |v: F, α: F| {
+            if v > α {
+                v - α
+            } else if v < -α {
+                -v - α
+            } else {
+                0.0
+            }
+        };
+
+        let μ_disrepancy = μ
+            .iter_spikes()
+            .map(|δ| {
+                let v = τv.apply(&δ.x);
+                match δ.α.total_cmp(&0.0) {
+                    Greater => (v + τα).abs(),
+                    Less => (v - τα).abs(),
+                    Equal => threshold(v, τα),
+                }
+            })
+            .reduce(F::max)
+            .unwrap_or(0.0);
+
+        let (_, τv_min) = τv.minimise(refinement_tolerance, config.refinement.max_steps);
+        let (_, τv_max) = τv.maximise(refinement_tolerance, config.refinement.max_steps);
+        //let rest_discrepancy = (-d_min - α).max(0.0);
+        //return rest_discrepancy.max(μ_disrepancy)
+        // Simplifies to
+        return threshold(τv_min, τα)
+            .max(threshold(τv_max, τα))
+            .max(μ_disrepancy);
+    }
+
     fn verify_merge_candidate<M>(
         &self,
         d: &mut M,

--- a/src/run.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/run.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -55,7 +55,7 @@
 };
 use alg_tools::maputil::map3;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::{Dist, NormExponent, L1, L2};
+use alg_tools::norms::{Dist, Norm, NormExponent, L1, L2};
 use alg_tools::operator_arithmetic::Weighted;
 use alg_tools::sets::Cube;
 use alg_tools::sets::SetOrd;
@@ -73,6 +73,7 @@
 use serde_json;
 use std::collections::HashMap;
 use std::hash::Hash;
+//use std::ops::Mul;
 use std::time::Instant;
 use thiserror::Error;
 
@@ -402,6 +403,10 @@
     pruned: usize,
     this_iters: usize,
     epsilon: F,
+    discrepancy: F,
+    aux_discrepancy: F,
+    joint_discrepancy: F,
+    radon_norm: F,
 }
 
 /// Collected experiment statistics
@@ -974,7 +979,7 @@
     A::Observable: AXPY,
     Reg: SlidingRegTerm<Loc<N, F>, F>,
     P: ProxPenalty<Loc<N, F>, A::PreadjointCodomain, Reg, F> + StepLengthBoundPD<F, A, RNDM<N, F>>,
-    RNDM<N, F>: SpikeMerging<F>,
+    RNDM<N, F>: SpikeMerging<F> + Norm<Radon, F>,
     I: AlgIteratorFactory<IterInfo<F>>,
     Plot: Plotter<P::ReturnMapping, A::PreadjointCodomain, RNDM<N, F>>,
 {
@@ -1421,10 +1426,10 @@
     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>,
     //Pair<S, Z>: ClosedMul<F>,
-    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>
@@ -1434,7 +1439,7 @@
     Y: ClosedEuclidean<F> + Clone,
     for<'b> &'b Y: Instance<Y>,
     Z: ClosedEuclidean<F> + Clone + ClosedMul<F> + AXPY,
-    for<'b> &'b Z: Instance<Z>,
+    for<'b> &'b Z: Instance<Z>, // compiler overflow  + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     H: Conjugable<Y, F, Codomain = F>,
     for<'b> H::Conjugate<'b>: Prox<Y>,
@@ -1485,10 +1490,10 @@
     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>,
     //Pair<S, Z>: ClosedMul<F>,
-    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>
@@ -1498,7 +1503,7 @@
     Y: ClosedEuclidean<F> + Clone,
     for<'b> &'b Y: Instance<Y>,
     Z: ClosedEuclidean<F> + Clone + ClosedMul<F> + AXPY,
-    for<'b> &'b Z: Instance<Z>,
+    for<'b> &'b Z: Instance<Z>, // compiler overflow  + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     H: Conjugable<Y, F, Codomain = F>,
     for<'b> H::Conjugate<'b>: Prox<Y>,
@@ -1572,13 +1577,13 @@
     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>,
     Z: ClosedEuclidean<F> + AXPY + Clone,
-    for<'b> &'b Z: Instance<Z>,
+    for<'b> &'b Z: Instance<Z>, // compiler overflow  + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
 {
@@ -1615,13 +1620,13 @@
     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>,
     Z: ClosedEuclidean<F> + AXPY + Clone,
-    for<'b> &'b Z: Instance<Z>,
+    for<'b> &'b Z: Instance<Z>, // compiler overflow  + Mul<F, Output = Z>,
     R: Prox<Z, Codomain = F>,
     Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
 {
@@ -1712,8 +1717,12 @@
             inner_iters,
             merged,
             pruned,
+            discrepancy,
+            aux_discrepancy,
+            joint_discrepancy,
             //postprocessing,
             this_iters,
+            μ_norm,
             ε,
             ..
         } = data;
@@ -1742,7 +1751,11 @@
             merged,
             pruned,
             this_iters,
+            radon_norm: μ_norm,
             epsilon: ε,
+            discrepancy: discrepancy.unwrap_or(F::NEG_INFINITY),
+            aux_discrepancy: aux_discrepancy.unwrap_or(F::NEG_INFINITY),
+            joint_discrepancy: joint_discrepancy.unwrap_or(F::NEG_INFINITY),
         }
     };
 

--- a/src/sliding_fb.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/sliding_fb.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -1006,7 +1006,7 @@
     Dat::DerivativeDomain:
         DifferentiableMapping<U, Codomain = F> + ClosedMul<F> + MinMaxMapping<U, F>,
     //for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Reg: SlidingRegTerm<U, F>,
     P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
     Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
@@ -1038,19 +1038,23 @@
     let mut ε = tolerance.initial();
 
     // Statistics
-    let full_stats = |μ: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
-        value: f.apply(μ) + reg.apply(μ),
-        n_spikes: μ.len(),
-        ε,
-        // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
-        ..stats
-    };
+    let full_stats =
+        |μ: &DiscreteMeasure<U, F>, v: &mut Dat::DerivativeDomain, ε, stats| IterInfo {
+            value: f.apply(μ) + reg.apply(μ),
+            n_spikes: μ.len(),
+            discrepancy: Some(τ * reg.discrepancy(v, 1.0, ε / τ, μ, &config.insertion)),
+            μ_norm: μ.norm(Radon),
+            ε,
+            // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
+            ..stats
+        };
     let mut stats = IterInfo::new();
 
+    // Construct initial smooth part of surrogate model
+    let mut v = f.differential(&μ);
+
     // Run the algorithm
-    for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
-        let mut v = f.differential(&μ);
-
+    for state in iterator.iter_init(|| full_stats(&μ, &mut v, ε, stats.clone())) {
         // Calculate initial transport
         prox_penalty.initial_transport(
             &mut γ,
@@ -1164,10 +1168,18 @@
         let iter = state.iteration();
         stats.this_iters += 1;
 
+        // Update smooth part of surrogate model
+        v = f.differential(&μ);
+
         // Give statistics if requested
         state.if_verbose(|| {
-            plotter.plot_spikes(iter, None, Some(&τv̆), &μ);
-            full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
+            plotter.plot_spikes(iter, None, Some(&v), &μ);
+            full_stats(
+                &μ,
+                &mut v,
+                ε,
+                std::mem::replace(&mut stats, IterInfo::new()),
+            )
         });
 
         // Update main tolerance for next iteration
@@ -1198,7 +1210,7 @@
     Dat::DerivativeDomain:
         DifferentiableMapping<U, Codomain = F> + ClosedMul<F> + MinMaxMapping<U, F>,
     //for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
-    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
     Reg: SlidingRegTerm<U, F>,
     P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
     Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,

--- 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>,

--- a/src/types.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/types.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -66,10 +66,18 @@
     pub inner_iters: usize,
     /// Transport statistis
     pub transport: Option<TransportInfo<F>>,
+    /// Measure component discrepancy
+    pub discrepancy: Option<F>,
+    /// Auxiliary parameter discrepancy
+    pub aux_discrepancy: Option<F>,
+    /// Discrepancy of both components
+    pub joint_discrepancy: Option<F>,
     /// Current tolerance
     pub ε: F,
     // /// Solve fin.dim problem for this measure to get the optimal `value`.
     // pub postprocessing : Option<RNDM<N, F>>,
+    /// Norm of the measure
+    pub μ_norm: F,
 }
 
 impl<F: Float> IterInfo<F> {
@@ -86,6 +94,10 @@
             ε: F::NAN,
             // postprocessing : None,
             transport: None,
+            discrepancy: None,
+            aux_discrepancy: None,
+            joint_discrepancy: None,
+            μ_norm: F::ZERO,
         }
     }
 
@@ -105,8 +117,15 @@
 {
     fn logrepr(&self) -> ColoredString {
         format!(
-            "{}\t| N = {}, ε = {:.2e}, 𝔼inner_it = {}, 𝔼ins/mer/pru = {}/{}/{}{}",
+            "\t{:.10}{}| N = {}, ε = {:.2e}, 𝔼inner_it = {}, 𝔼ins/mer/pru = {}/{}/{}{}",
             self.value.logrepr(),
+            match &self.joint_discrepancy {
+                Some(d) => format!(" ∂_joint={:.10}", d),
+                None => match &self.discrepancy {
+                    None => format!(""),
+                    Some(d) => format!(" ∂={:.10}", d),
+                },
+            },
             self.n_spikes,
             self.ε,
             self.inner_iters as float / self.this_iters.max(1) as float,