diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/forward_pdps.rs
--- a/src/forward_pdps.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/forward_pdps.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -3,6 +3,7 @@
 primal-dual proximal splitting with a forward step.
 */
 
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::fb::*;
 use crate::measures::merging::SpikeMerging;
 use crate::measures::{DiscreteMeasure, RNDM};
@@ -20,6 +21,7 @@
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::norms::L2;
 use anyhow::ensure;
+use measures::DeltaMeasure;
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
 
@@ -85,7 +87,7 @@
 ///     τ = τ_0 \frac{1 - σ_{p,0}}{(σ_d M (1-σ_p L_z) + (1 - σ_{p,0} L)}.
 /// $$
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_forward_pdps_pair<
+pub fn pointsource_forward_pdps_pair_dom<
     F,
     I,
     S,
@@ -111,6 +113,7 @@
     opKz: &KOpZ,
     fnR: &R,
     fnH: &H,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<MeasureZ<F, Z, N>>
 where
     F: Float + ToNalgebraRealField,
@@ -237,10 +240,7 @@
         // Merge spikes.
         // This crucially expects the merge routine to be stable with respect to spike locations,
         // and not to performing any pruning. That is be to done below simultaneously for γ.
-        // Merge spikes.
-        // This crucially expects the merge routine to be stable with respect to spike locations,
-        // and not to performing any pruning. That is be to done below simultaneously for γ.
-        config.insertion.if_merge_now(&state, |ins| {
+        let did_merge = config.insertion.if_merge_now(&state, |ins| {
             stats.merged += prox_penalty.merge_spikes(
                 &mut μ,
                 &mut τv,
@@ -249,32 +249,45 @@
                 ε,
                 ins,
                 &reg,
-                Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))+ reg.apply(μ̃)),
+                Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z)) + reg.apply(μ̃)),
             );
-        });
+            stats.pruned += prune_with_stats(&mut μ);
 
-        // Prune spikes with zero weight.
-        stats.pruned += prune_with_stats(&mut μ);
+            if config.insertion.random_insertions {
+                if let Some(x) = domain.random_point() {
+                    μ += DeltaMeasure { x, α: 0.0 };
+                    stats.inserted += 1;
+                }
+            }
 
-        for _ in 0..config.insertion.extra_weight_optimisation_steps {
-            // Do extra weight optimisation step heuristic
-            let Pair(v, _z_tmp) = f.differential(Pair(&μ, &z));
-            τv = v * τ;
-            prox_penalty.reweigh(
-                &mut μ,
-                &mut τv,
-                τ,
-                ε,
-                &config.insertion,
-                &reg,
-                &state,
-                &mut stats,
-            )?;
-            //(z, y) = zy_step(z, y, z_tmp);
-        }
+            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));
+                τv = v * τ;
+                prox_penalty.reweigh(
+                    &mut μ,
+                    &mut τv,
+                    τ,
+                    ε,
+                    &config.insertion,
+                    &reg,
+                    &state,
+                    &mut stats,
+                )?;
 
-        if config.insertion.extra_weight_optimisation_steps > 0 {
-            // Prune spikes with zero weight.
+                // opKz_adj.apply_add(&mut z_tmp, &y);
+                // z_tmp.axpy(1.0, &z, -σ_p);
+                // let z_new = fnR.prox(σ_p, z_tmp);
+                // z.axpy(1.0 + ω, &z_new, -ω);
+                // //z = (z - &z_new) * (-ω) + &z_new;
+                // opKz.gemv(&mut y, σ_d, z, 1.0);
+                // y = starH.prox(σ_d, y);
+                // z = z_new;
+            }
+            Ok(())
+        })?;
+
+        if did_merge.is_none() {
             stats.pruned += prune_with_stats(&mut μ);
         }
 
@@ -305,13 +318,79 @@
     Ok(Pair(μ, z))
 }
 
+#[deprecated]
+pub fn pointsource_forward_pdps_pair<
+    F,
+    I,
+    S,
+    Dat,
+    Reg,
+    P,
+    Z,
+    R,
+    Y,
+    /*KOpM, */ KOpZ,
+    H,
+    Plot,
+    const N: usize,
+>(
+    f: &Dat,
+    reg: &Reg,
+    prox_penalty: &P,
+    config: &ForwardPDPSConfig<F>,
+    iterator: I,
+    plotter: Plot,
+    init: (Option<RNDM<N, F>>, Z, Y),
+    //opKμ : KOpM,
+    opKz: &KOpZ,
+    fnR: &R,
+    fnH: &H,
+) -> DynResult<MeasureZ<F, Z, N>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
+    //Pair<S, Z>: ClosedMul<F>, // Doesn't really need to be closed, if make this signature more complex…
+    S: DifferentiableRealMapping<N, F> + ClosedMul<F>,
+    RNDM<N, F>: SpikeMerging<F>,
+    Reg: RegTerm<Loc<N, F>, F>,
+    P: ProxPenalty<Loc<N, F>, S, Reg, F>,
+    for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
+    KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
+        + GEMV<F, Z, Y>
+        + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>,
+    KOpZ::SimpleAdjoint: GEMV<F, Y, Z>,
+    Y: ClosedEuclidean<F>,
+    for<'b> &'b Y: Instance<Y>,
+    Z: ClosedEuclidean<F> + AXPY,
+    for<'b> &'b Z: Instance<Z>,
+    R: Prox<Z, Codomain = F>,
+    H: Conjugable<Y, F, Codomain = F>,
+    for<'b> H::Conjugate<'b>: Prox<Y>,
+    Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
+{
+    pointsource_forward_pdps_pair_dom(
+        f,
+        reg,
+        prox_penalty,
+        config,
+        iterator,
+        plotter,
+        init,
+        opKz,
+        fnR,
+        fnH,
+        &mut FullSpace,
+    )
+}
+
 /// Iteratively solve the pointsource localisation with an additional variable
 /// using forward-backward splitting.
 ///
 /// The implementation uses [`pointsource_forward_pdps_pair`] with appropriate dummy
 /// variables, operators, and functions.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fb_pair<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>(
+pub fn pointsource_fb_pair_dom<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>(
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
@@ -321,6 +400,7 @@
     (μ0, z): (Option<RNDM<N, F>>, Z),
     //opKμ : KOpM,
     fnR: &R,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<MeasureZ<F, Z, N>>
 where
     F: Float + ToNalgebraRealField,
@@ -345,7 +425,7 @@
     let &FBConfig { τ0, σp0, insertion } = config;
     let pdps_config = ForwardPDPSConfig { τ0, σp0, insertion, σd0: 0.0 };
 
-    pointsource_forward_pdps_pair(
+    pointsource_forward_pdps_pair_dom(
         f,
         reg,
         prox_penalty,
@@ -356,5 +436,48 @@
         &opKz,
         fnR,
         &fnH,
+        domain,
     )
 }
+
+#[deprecated]
+#[replace_float_literals(F::cast_from(literal))]
+pub fn pointsource_fb_pair<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>(
+    f: &Dat,
+    reg: &Reg,
+    prox_penalty: &P,
+    config: &FBConfig<F>,
+    iterator: I,
+    plotter: Plot,
+    init: (Option<RNDM<N, F>>, Z),
+    //opKμ : KOpM,
+    fnR: &R,
+) -> DynResult<MeasureZ<F, Z, N>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
+    S: DifferentiableRealMapping<N, F> + ClosedMul<F>,
+    RNDM<N, F>: SpikeMerging<F>,
+    Reg: RegTerm<Loc<N, F>, F>,
+    P: ProxPenalty<Loc<N, F>, S, Reg, F>,
+    for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
+    Z: ClosedEuclidean<F> + AXPY<Field = F> + Clone,
+    for<'b> &'b Z: Instance<Z>,
+    R: Prox<Z, Codomain = F>,
+    Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
+    // We should not need to explicitly require this:
+    for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>,
+{
+    pointsource_fb_pair_dom(
+        f,
+        reg,
+        prox_penalty,
+        config,
+        iterator,
+        plotter,
+        init,
+        fnR,
+        &mut FullSpace,
+    )
+}