diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/fb.rs
--- a/src/fb.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/fb.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -77,6 +77,7 @@
 [`crate::subproblem::InnerSettings`] in [`InsertionConfig::inner`].
 */
 
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::measures::merging::SpikeMerging;
 use crate::measures::{DiscreteMeasure, RNDM};
 use crate::plot::Plotter;
@@ -89,6 +90,7 @@
 use alg_tools::mapping::DifferentiableMapping;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use colored::Colorize;
+use measures::DeltaMeasure;
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
 
@@ -146,7 +148,7 @@
 ///
 /// Returns the final iterate.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fb_reg<F, I, Dat, Reg, P, Plot, const N: usize>(
+pub fn pointsource_fb_dom<F, I, Dat, Reg, P, Plot, const N: usize>(
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
@@ -154,6 +156,7 @@
     iterator: I,
     mut plotter: Plot,
     μ0: Option<RNDM<N, F>>,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<RNDM<N, F>>
 where
     F: Float + ToNalgebraRealField,
@@ -188,7 +191,7 @@
 
     // Run the algorithm
     for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
-        let maybe_μ_base = config.if_merge_now(&state, |_| μ.clone());
+        let maybe_μ_base = config.if_merge_now_check(&state, || μ.clone());
         let μ_base_len = μ.len();
 
         // Calculate smooth part of surrogate model.
@@ -200,7 +203,10 @@
         stats.inserted += μ.len() - μ_base_len;
 
         // Prune and possibly merge spikes
-        if let Some(μ_base) = maybe_μ_base {
+        let did_merge = if let Some(μ_base) = maybe_μ_base {
+            // 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.if_merge_now(&state, |ins| {
                 stats.merged += prox_penalty.merge_spikes(
                     &mut μ,
@@ -211,19 +217,29 @@
                     ins,
                     &reg,
                     Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃)),
-                )
-            });
-        }
+                );
+                stats.pruned += prune_with_stats(&mut μ);
 
-        stats.pruned += prune_with_stats(&mut μ);
+                if config.random_insertions {
+                    if let Some(x) = domain.random_point() {
+                        μ += DeltaMeasure { x, α: 0.0 };
+                        stats.inserted += 1;
+                    }
+                }
 
-        // Do extra weight optimisation step heuristic
-        for _ in 0..config.extra_weight_optimisation_steps {
-            τv = f.differential(&μ) * τ;
-            prox_penalty.reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
-        }
+                // Do extra weight optimisation step heuristic
+                for _ in 1..config.extra_weight_optimisation_steps {
+                    τv = f.differential(&μ) * τ;
+                    prox_penalty
+                        .reweigh(&mut μ, &mut τv, τ, ε, &config, &reg, &state, &mut stats)?;
+                }
+                Ok(())
+            })?
+        } else {
+            None
+        };
 
-        if config.extra_weight_optimisation_steps > 0 {
+        if did_merge.is_none() {
             stats.pruned += prune_with_stats(&mut μ);
         }
 
@@ -244,6 +260,38 @@
     postprocess(μ, config, |μ̃| f.apply(μ̃) + reg.apply(μ̃))
 }
 
+#[deprecated]
+pub fn pointsource_fb_reg<F, I, Dat, Reg, P, Plot, const N: usize>(
+    f: &Dat,
+    reg: &Reg,
+    prox_penalty: &P,
+    fbconfig: &FBConfig<F>,
+    iterator: I,
+    plotter: Plot,
+    μ0: Option<RNDM<N, F>>,
+) -> DynResult<RNDM<N, F>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    RNDM<N, F>: SpikeMerging<F>,
+    Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F>,
+    Dat::DerivativeDomain: ClosedMul<F>,
+    Reg: RegTerm<Loc<N, F>, F>,
+    P: ProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+    Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
+{
+    pointsource_fb_dom(
+        f,
+        reg,
+        prox_penalty,
+        fbconfig,
+        iterator,
+        plotter,
+        μ0,
+        &mut FullSpace,
+    )
+}
+
 /// Iteratively solve the pointsource localisation problem using inertial forward-backward splitting.
 ///
 /// The settings in `config` have their [respective documentation](FBConfig). `opA` is the
@@ -312,7 +360,7 @@
 
         // (Do not) merge spikes.
         if !warned_merging {
-            config.if_merge_now(&state, |_| {
+            config.if_merge_now_check(&state, || {
                 let err = format!("Merging not supported for μFISTA");
                 println!("{}", err.red());
                 warned_merging = true;