# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1785870738 18000
# Node ID 8d1c959432c8a707af192288123ab4ef8e0e756e
# Parent  51dbb2cc25d9e64171175e205266c67bc19ba508
Add AdmissibleDomain

diff -r 51dbb2cc25d9 -r 8d1c959432c8 Cargo.lock
--- a/Cargo.lock	Tue Aug 04 13:46:29 2026 -0500
+++ b/Cargo.lock	Tue Aug 04 14:12:18 2026 -0500
@@ -878,9 +878,9 @@
 
 [[package]]
 name = "rand"
-version = "0.9.2"
+version = "0.9.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "6db2770f06117d490610c7488547d543617b21bfa07796d7a12f6f1bd53850d1"
+checksum = "b9ef1d0d795eb7d84685bca4f72f3649f064e6641543d3a8c415898726a57b41"
 dependencies = [
  "rand_chacha",
  "rand_core",
diff -r 51dbb2cc25d9 -r 8d1c959432c8 Cargo.toml
--- a/Cargo.toml	Tue Aug 04 13:46:29 2026 -0500
+++ b/Cargo.toml	Tue Aug 04 14:12:18 2026 -0500
@@ -34,7 +34,7 @@
 [dependencies]
 serde = { version = "1.0", features = ["derive"] }
 num-traits = { version = "~0.2.14", features = ["std"] }
-rand = "~0.9.2"
+rand = "~0.9.5"
 colored = "~2.1.0"
 rand_distr = "~0.5.1"
 nalgebra = { version = "~0.34.0", features = ["rand-no-std"] }
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/domain.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/domain.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -0,0 +1,46 @@
+/*!
+Definition of an admissible domain.
+*/
+
+use alg_tools::instance::{Instance, Space};
+use alg_tools::loc::Loc;
+use alg_tools::sets::{Cube, Set};
+use alg_tools::types::Num;
+use rand::distr::Uniform;
+use rand::prelude::{Rng, StdRng};
+
+/// Admissible domain
+pub trait AdmissibleDomain<Point: Space>: Set<Point> {
+    // Generate a random point in the domain
+    fn random_point(&mut self) -> Option<Point> {
+        None
+    }
+}
+
+impl<F: Num + PartialOrd + Sized, const N: usize> AdmissibleDomain<Loc<N, F>> for Cube<N, F> {}
+
+///
+pub struct FullSpace;
+
+impl<A: Space> Set<A> for FullSpace {
+    fn contains<I: Instance<A>>(&self, _: I) -> bool {
+        true
+    }
+}
+
+impl<F: Num + PartialOrd + Sized, const N: usize> AdmissibleDomain<Loc<N, F>> for FullSpace {}
+
+pub struct UniformRandom<S, R: Rng = StdRng>(pub S, pub R);
+
+impl<A: Space, S: Set<A>, R: Rng> Set<A> for UniformRandom<S, R> {
+    fn contains<I: Instance<A>>(&self, item: I) -> bool {
+        self.0.contains(item)
+    }
+}
+
+impl<const N: usize, R: Rng> AdmissibleDomain<Loc<N, f64>> for UniformRandom<Cube<N, f64>, R> {
+    fn random_point(&mut self) -> Option<Loc<N, f64>> {
+        let UniformRandom(ref cube, ref mut rng) = *self;
+        Some(Loc(cube.map(|a, b| rng.sample(Uniform::new(a, b).unwrap()))))
+    }
+}
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;
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,
+    )
+}
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/lib.rs
--- a/src/lib.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/lib.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -21,6 +21,7 @@
 //use pyo3::pyclass;
 
 pub mod dataterm;
+pub mod domain;
 pub mod fb;
 pub mod forward_model;
 pub mod forward_pdps;
@@ -242,9 +243,14 @@
     pub max_scaling_iter: Option<usize>,
 
     #[arg(long)]
-    /// Extra finite dimensional steps to take
+    /// Extra finite dimensional steps to take after merge and random insertion.
+    /// This all methods support this.
     pub extra_weight_optimisation_steps: Option<usize>,
 
+    #[arg(long)]
+    /// Whether to perform random insertions after merge. Not all methods support this.
+    pub random_insertions: Option<bool>,
+
     #[arg(long, value_names = &["ε", "θ", "p"])]
     /// Set the tolerance to ε_k = ε/(1+θk)^p
     pub tolerance: Option<Vec<F>>,
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/pdps.rs
--- a/src/pdps.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/pdps.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -212,7 +212,8 @@
         config.if_merge_now(&state, |ins| {
             stats.merged +=
                 prox_penalty.merge_spikes_no_fitness(&mut μ, &mut τv, &μ_base, τ, ε, ins, &reg);
-        });
+            Ok(())
+        })?;
         stats.inserted += μ.len() - μ_base.len();
         stats.pruned += prune_with_stats(&mut μ);
 
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/prox_penalty.rs
--- a/src/prox_penalty.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/prox_penalty.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -66,9 +66,12 @@
     /// Iterations between merging heuristic tries
     pub merge_every: usize,
 
-    /// Additional weight optimisation steps
+    /// Additional weight optimisation steps after merge and random insertion
     pub extra_weight_optimisation_steps: usize,
 
+    /// Random insertion after merge
+    pub random_insertions: bool,
+
     /// Maximum iteration for scaling hack, zero if no maximum.
     pub max_scaling_iter: usize,
 }
@@ -91,16 +94,45 @@
             merge_tolerance_mult: 2.0,
             extra_weight_optimisation_steps: 0,
             max_scaling_iter: 0,
+            random_insertions: false,
         }
     }
 }
 
 impl<F: Float> InsertionConfig<F> {
-    /// Check if merging should be attempted this iteration
+    /// Check if merging should be attempted this iteration.
+    /// Wrapped in a possible error, returns [`Some`] of the return value of `f`, if merging should
+    /// be performed.
     pub fn if_merge_now<I: AlgIterator, A>(
         &self,
         state: &AlgIteratorIteration<I>,
-        mut f: impl FnMut(&Self) -> A,
+        mut f: impl FnMut(&Self) -> DynResult<A>,
+    ) -> DynResult<Option<A>> {
+        if (self.merging.enabled || self.merging.scaling.is_some())
+            && state.iteration() % self.merge_every == 0
+        {
+            if self.max_scaling_iter != 0 && state.iteration() > self.max_scaling_iter {
+                if !self.merging.enabled {
+                    Ok(None)
+                } else {
+                    Ok(Some(f(&InsertionConfig {
+                        merging: SpikeMergingMethod { scaling: None, ..self.merging },
+                        ..*self
+                    })?))
+                }
+            } else {
+                Ok(Some(f(self)?))
+            }
+        } else {
+            Ok(None)
+        }
+    }
+
+    // Lightweight version of [`Self::if_merge_now`].
+    pub fn if_merge_now_check<I: AlgIterator, A>(
+        &self,
+        state: &AlgIteratorIteration<I>,
+        mut f: impl FnMut() -> A,
     ) -> Option<A> {
         if (self.merging.enabled || self.merging.scaling.is_some())
             && state.iteration() % self.merge_every == 0
@@ -109,13 +141,10 @@
                 if !self.merging.enabled {
                     None
                 } else {
-                    Some(f(&InsertionConfig {
-                        merging: SpikeMergingMethod { scaling: None, ..self.merging },
-                        ..*self
-                    }))
+                    Some(f())
                 }
             } else {
-                Some(f(self))
+                Some(f())
             }
         } else {
             None
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/regularisation.rs
--- a/src/regularisation.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/regularisation.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -3,11 +3,11 @@
 */
 
 #[allow(unused_imports)] // Used by documentation.
-use crate::fb::pointsource_fb_reg;
+use crate::fb::pointsource_fb_dom;
 use crate::fb::InsertionConfig;
 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
 #[allow(unused_imports)] // Used by documentation.
-use crate::sliding_fb::pointsource_sliding_fb_dom;
+use crate::sliding_fb::pointsource_sliding_fb;
 use crate::types::*;
 use alg_tools::instance::{Instance, Space};
 use alg_tools::linops::Mapping;
@@ -27,7 +27,7 @@
 
 use std::cmp::Ordering::{Equal, Greater, Less};
 
-/// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other
+/// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_dom`] and other
 /// algorithms.
 ///
 /// The only member of the struct is the regularisation parameter α.
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/run.rs
--- a/src/run.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/run.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -2,7 +2,8 @@
 This module provides [`RunnableExperiment`] for running chosen algorithms on a chosen experiment.
 */
 
-use crate::fb::{pointsource_fb_reg, pointsource_fista_reg, FBConfig, InsertionConfig};
+use crate::domain::{AdmissibleDomain, FullSpace};
+use crate::fb::{pointsource_fb_dom, pointsource_fista_reg, FBConfig, InsertionConfig};
 use crate::forward_model::sensor_grid::{
     //SensorGridBTFN,
     Sensor,
@@ -11,7 +12,9 @@
     Spread,
 };
 use crate::forward_model::*;
-use crate::forward_pdps::{pointsource_fb_pair, pointsource_forward_pdps_pair, ForwardPDPSConfig};
+use crate::forward_pdps::{
+    pointsource_fb_pair_dom, pointsource_forward_pdps_pair_dom, ForwardPDPSConfig,
+};
 use crate::frank_wolfe::{pointsource_fw_reg, FWConfig, FWVariant, RegTermFW};
 use crate::kernels::*;
 use crate::measures::merging::{SpikeMerging, SpikeMergingMethod};
@@ -24,7 +27,7 @@
 use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, Regularisation, SlidingRegTerm};
 use crate::seminorms::*;
 use crate::sliding_fb::{
-    pointsource_sliding_fb_dom, SlidingFBConfig, TransportConfig, TransportProxPenalty,
+    pointsource_sliding_fb, SlidingFBConfig, TransportConfig, TransportProxPenalty,
 };
 use crate::sliding_pdps::{
     pointsource_sliding_fb_pair_dom, pointsource_sliding_pdps_pair_dom, SlidingPDPSConfig,
@@ -127,6 +130,7 @@
             extra_weight_optimisation_steps: cli
                 .extra_weight_optimisation_steps
                 .unwrap_or(g.extra_weight_optimisation_steps),
+            random_insertions: cli.random_insertions.unwrap_or(g.random_insertions),
             tolerance: cli
                 .tolerance
                 .as_ref()
@@ -989,13 +993,7 @@
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
-    init: (
-        &AlgorithmConfig<F>,
-        I,
-        Plot,
-        Option<RNDM<N, F>>,
-        String,
-    ),
+    init: (&AlgorithmConfig<F>, I, Plot, Option<RNDM<N, F>>, String),
     cont: impl FnOnce(
         (&AlgorithmConfig<F>, I, Plot, Option<RNDM<N, F>>, String),
     ) -> DynResult<RNDM<N, F>>,
@@ -1011,7 +1009,7 @@
     P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
     Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
 {
-    run_fb_dom(f, reg, prox_penalty, init, None, cont)
+    run_fb_dom(f, reg, prox_penalty, init, &mut FullSpace, cont)
 }
 
 /// Runs FB-style algorithms if `alg` so requests and `prox_penalty` matches.
@@ -1028,7 +1026,7 @@
         Option<RNDM<N, F>>,
         String,
     ),
-    domain: Option<&P::Dom>,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
     cont: impl FnOnce(
         (&AlgorithmConfig<F>, I, Plot, Option<RNDM<N, F>>, String),
     ) -> DynResult<RNDM<N, F>>,
@@ -1049,7 +1047,16 @@
     match alg {
         &AlgorithmConfig::FB(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
-            pointsource_fb_reg(f, reg, prox_penalty, algconfig, iterator, plotter, μ0)
+            pointsource_fb_dom(
+                f,
+                reg,
+                prox_penalty,
+                algconfig,
+                iterator,
+                plotter,
+                μ0,
+                domain,
+            )
         }
         &AlgorithmConfig::FISTA(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
@@ -1057,7 +1064,16 @@
         }
         &AlgorithmConfig::SlidingFB(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
-            pointsource_sliding_fb_dom(f, reg, prox_penalty, algconfig, iterator, plotter, μ0, domain)
+            pointsource_sliding_fb(
+                f,
+                reg,
+                prox_penalty,
+                algconfig,
+                iterator,
+                plotter,
+                μ0,
+                domain,
+            )
         }
         _ => cont((alg, iterator, plotter, μ0, running)),
     }
@@ -1294,11 +1310,29 @@
                     save_extra,
                     (μ0, z, y),
                     |p| {
-                        run_pdps_pair_dom(&f, &reg, &RadonSquared, &opKz, &fnR, &fnH, p, None, |q| {
-                            run_pdps_pair_dom(&f, &reg, &op𝒟, &opKz, &fnR, &fnH, q, None, |_| {
-                                Err(NotImplemented.into())
-                            })
-                        })
+                        run_pdps_pair_dom(
+                            &f,
+                            &reg,
+                            &RadonSquared,
+                            &opKz,
+                            &fnR,
+                            &fnH,
+                            p,
+                            &mut FullSpace,
+                            |q| {
+                                run_pdps_pair_dom(
+                                    &f,
+                                    &reg,
+                                    &op𝒟,
+                                    &opKz,
+                                    &fnR,
+                                    &fnH,
+                                    q,
+                                    &mut FullSpace,
+                                    |_| Err(NotImplemented.into()),
+                                )
+                            },
+                        )
                         .map(|Pair(μ, z)| (μ, z))
                     },
                 )
@@ -1314,11 +1348,29 @@
                     save_extra,
                     (μ0, z, y),
                     |p| {
-                        run_pdps_pair_dom(&f, &reg, &RadonSquared, &opKz, &fnR, &fnH, p, None, |q| {
-                            run_pdps_pair_dom(&f, &reg, &op𝒟, &opKz, &fnR, &fnH, q, None, |_| {
-                                Err(NotImplemented.into())
-                            })
-                        })
+                        run_pdps_pair_dom(
+                            &f,
+                            &reg,
+                            &RadonSquared,
+                            &opKz,
+                            &fnR,
+                            &fnH,
+                            p,
+                            &mut FullSpace,
+                            |q| {
+                                run_pdps_pair_dom(
+                                    &f,
+                                    &reg,
+                                    &op𝒟,
+                                    &opKz,
+                                    &fnR,
+                                    &fnH,
+                                    q,
+                                    &mut FullSpace,
+                                    |_| Err(NotImplemented.into()),
+                                )
+                            },
+                        )
                         .map(|Pair(μ, z)| (μ, z))
                     },
                 )
@@ -1379,7 +1431,17 @@
     for<'b> H::Conjugate<'b>: Prox<Y>,
     Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
 {
-    run_pdps_pair_dom(f, reg, prox_penalty, opKz, fnR, fnH, init, None, cont)
+    run_pdps_pair_dom(
+        f,
+        reg,
+        prox_penalty,
+        opKz,
+        fnR,
+        fnH,
+        init,
+        &mut FullSpace,
+        cont,
+    )
 }
 
 pub fn run_pdps_pair_dom<F, S, Dat, Reg, Z, R, Y, KOpZ, H, P, I, Plot, const N: usize>(
@@ -1396,7 +1458,7 @@
         (Option<RNDM<N, F>>, Z, Y),
         String,
     ),
-    domain: Option<&P::Dom>,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
     cont: impl FnOnce(
         (
             &AlgorithmConfig<F>,
@@ -1436,7 +1498,7 @@
     match alg {
         &AlgorithmConfig::ForwardPDPS(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
-            pointsource_forward_pdps_pair(
+            pointsource_forward_pdps_pair_dom(
                 f,
                 reg,
                 prox_penalty,
@@ -1447,6 +1509,7 @@
                 opKz,
                 fnR,
                 fnH,
+                domain,
             )
         }
         &AlgorithmConfig::SlidingPDPS(ref algconfig, prox_type) if prox_type == pt => {
@@ -1508,7 +1571,7 @@
     // We should not need to explicitly require this:
     for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>,
 {
-    run_fb_pair_dom(f, reg, prox_penalty, fnR, init, None, cont)
+    run_fb_pair_dom(f, reg, prox_penalty, fnR, init, &mut FullSpace, cont)
 }
 
 pub fn run_fb_pair_dom<F, S, Dat, Reg, Z, R, P, I, Plot, const N: usize>(
@@ -1523,7 +1586,7 @@
         (Option<RNDM<N, F>>, Z),
         String,
     ),
-    domain: Option<&P::Dom>,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
     cont: impl FnOnce(
         (
             &AlgorithmConfig<F>,
@@ -1556,7 +1619,17 @@
     match alg {
         &AlgorithmConfig::FB(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
-            pointsource_fb_pair(f, reg, prox_penalty, algconfig, iterator, plotter, μ0z, fnR)
+            pointsource_fb_pair_dom(
+                f,
+                reg,
+                prox_penalty,
+                algconfig,
+                iterator,
+                plotter,
+                μ0z,
+                fnR,
+                domain,
+            )
         }
         &AlgorithmConfig::SlidingFB(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
@@ -1569,7 +1642,7 @@
                 plotter,
                 μ0z,
                 fnR,
-                domain
+                domain,
             )
         }
         _ => cont((alg, iterator, plotter, μ0z, running)),
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/sliding_fb.rs
--- a/src/sliding_fb.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/sliding_fb.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -11,6 +11,7 @@
 use std::iter::Iterator;
 use std::ops::MulAssign;
 
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::fb::*;
 use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess};
 use crate::measures::merging::SpikeMerging;
@@ -28,7 +29,6 @@
 use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Mapping, RealMapping};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::norms::Norm;
-use alg_tools::sets::{Set, Cube};
 use anyhow::ensure;
 use std::ops::ControlFlow;
 
@@ -223,36 +223,35 @@
 const MINIMAL_PARTIAL_TRANSPORT: bool = true;
 const NEW_APPROACH: bool = true;
 
-pub trait TransportProxPenalty<Domain, PreadjointCodomain, Reg, F = f64>:
-    ProxPenalty<Domain, PreadjointCodomain, Reg, F>
+pub trait TransportProxPenalty<Point, PreadjointCodomain, Reg, F = f64>:
+    ProxPenalty<Point, PreadjointCodomain, Reg, F>
 where
     F: Float + ToNalgebraRealField,
-    Reg: SlidingRegTerm<Domain, F>,
-    Domain: Space + Clone,
+    Reg: SlidingRegTerm<Point, F>,
+    Point: Space + Clone,
 {
     type TransportStepLength;
-    type Dom : Set<Domain>;
 
     /// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
     /// with step lengh τ and transport step length `θ_or_adaptive`.
-    fn initial_transport(
+    fn initial_transport<Dom: AdmissibleDomain<Point>>(
         &self,
-        γ: &mut Transport<Domain, F>,
-        μ: &DiscreteMeasure<Domain, F>,
+        γ: &mut Transport<Point, F>,
+        μ: &DiscreteMeasure<Point, F>,
         ε: F,
         τ: F,
         τθ_or_adaptive: &mut Self::TransportStepLength,
         v: &PreadjointCodomain,
         tconfig: &TransportConfig<F>,
-        domain: Option<&Self::Dom>,
+        domain: &Dom,
     );
 
     /// A posteriori transport adaptation.
     fn aposteriori_transport(
         &self,
-        γ: &mut Transport<Domain, F>,
-        μ: &DiscreteMeasure<Domain, F>,
-        μ̆: &DiscreteMeasure<Domain, F>,
+        γ: &mut Transport<Point, F>,
+        μ: &DiscreteMeasure<Point, F>,
+        μ̆: &DiscreteMeasure<Point, F>,
         τv̆: &mut PreadjointCodomain,
         v: &mut PreadjointCodomain,
         extra: Option<F>,
@@ -284,7 +283,6 @@
     RNDM<N, F>: SpikeMerging<F>,
 {
     type TransportStepLength = TransportStepLength<F>;
-    type Dom = Cube<N, F>;
 
     fn get_transport_steplength(
         &self,
@@ -301,7 +299,7 @@
         )
     }
 
-    fn initial_transport(
+    fn initial_transport<Dom: AdmissibleDomain<Loc<N, F>>>(
         &self,
         γ: &mut Transport<Loc<N, F>, F>,
         μ: &RNDM<N, F>,
@@ -310,7 +308,7 @@
         τθ_or_adaptive: &mut TransportStepLength<F>,
         v: &M,
         tconfig: &TransportConfig<F>,
-        domain: Option<&Cube<N, F>>,
+        domain: &Dom,
     ) {
         γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig, domain);
     }
@@ -392,7 +390,6 @@
     RNDM<N, F>: SpikeMerging<F>,
 {
     type TransportStepLength = TransportStepLength<F>;
-    type Dom = Cube<N, F>;
 
     fn get_transport_steplength(
         &self,
@@ -409,7 +406,7 @@
         )
     }
 
-    fn initial_transport(
+    fn initial_transport<Dom: AdmissibleDomain<Loc<N, F>>>(
         &self,
         γ: &mut Transport<Loc<N, F>, F>,
         μ: &RNDM<N, F>,
@@ -418,7 +415,7 @@
         τθ_or_adaptive: &mut TransportStepLength<F>,
         v: &M,
         tconfig: &TransportConfig<F>,
-        domain: Option<&Cube<N, F>>,
+        domain: &Dom,
     ) {
         γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig, domain);
     }
@@ -745,13 +742,13 @@
     }
 
     /// Helper for initial transport. Called from [`TransportProxPenalty::initial_transport`].
-    pub(super) fn do_init_transport<D>(
+    pub(super) fn do_init_transport<D, Dom: AdmissibleDomain<Loc<N, F>>>(
         &mut self,
         v: &D,
         μ: &RNDM<N, F>,
         τθ_or_adaptive: &mut TransportStepLength<F>,
         tconfig: &TransportConfig<F>,
-        domain: Option<&Cube<N, F>>,
+        domain: &Dom,
     ) where
         D: DifferentiableRealMapping<N, F>,
     {
@@ -819,7 +816,7 @@
                             let dv_x = v.differential(&ρ.x);
                             let g = &dv_x * (ρ.α_γ.signum() * τθ);
                             let tmp = ρ.x - g;
-                            ρ.y = if domain.is_none_or(|d| d.contains(&tmp)) { tmp } else { ρ.x };
+                            ρ.y = if domain.contains(&tmp) { tmp } else { ρ.x };
                             let n = g.norm2();
                             if n >= F::EPSILON {
                                 // Estimate Lipschitz factor of ∇v
@@ -845,7 +842,7 @@
             for ρ in self.iter_mut() {
                 if ρ.fail_count <= tconfig.max_fail {
                     let tmp = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * τθ);
-                    ρ.y = if domain.is_none_or(|d| d.contains(&tmp)) { tmp } else { ρ.x };
+                    ρ.y = if domain.contains(&tmp) { tmp } else { ρ.x };
                 }
             }
         }
@@ -978,7 +975,7 @@
 /// may be given.
 /// Inertia is currently not supported.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_sliding_fb_dom<F, I, Dat, Reg, Plot, P, const N: usize>(
+pub fn pointsource_sliding_fb<F, I, Dat, Reg, Plot, P, const N: usize>(
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
@@ -986,7 +983,7 @@
     iterator: I,
     mut plotter: Plot,
     μ0: Option<RNDM<N, F>>,
-    domain: Option<&P::Dom>,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<RNDM<N, F>>
 where
     F: Float + ToNalgebraRealField,
@@ -1049,7 +1046,7 @@
             &mut τθ_or_adaptive,
             &v,
             &config.transport,
-            domain
+            domain,
         );
 
         let mut attempts = 0;
@@ -1107,8 +1104,8 @@
         // 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 m = prox_penalty.merge_spikes(
+        let did_merge = config.insertion.if_merge_now(&state, |ins| {
+            stats.merged += prox_penalty.merge_spikes(
                 &mut μ,
                 &mut τv̆,
                 &μ̆,
@@ -1116,32 +1113,35 @@
                 ε,
                 ins,
                 &reg,
-                Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃) ),
+                Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃)),
             );
-            //if m > 0 {
-            stats.merged += m;
-            //v = f.differential(&μ);
-            //}
-        });
+            γ.prune_compat(&mut μ, &mut stats);
 
-        γ.prune_compat(&mut μ, &mut stats);
+            if config.insertion.random_insertions {
+                if let Some(x) = domain.random_point() {
+                    μ += DeltaMeasure { x, α: 0.0 };
+                    stats.inserted += 1;
+                }
+            }
 
-        // Do extra weight optimisation step heuristic
-        for _ in 1..config.insertion.extra_weight_optimisation_steps {
-            τv̆ = f.differential(&μ) * τ;
-            prox_penalty.reweigh(
-                &mut μ,
-                &mut τv̆,
-                τ,
-                ε,
-                &config.insertion,
-                &reg,
-                &state,
-                &mut stats,
-            )?;
-        }
+            // Do extra weight optimisation step heuristic
+            for _ in 1..config.insertion.extra_weight_optimisation_steps {
+                τv̆ = f.differential(&μ) * τ;
+                prox_penalty.reweigh(
+                    &mut μ,
+                    &mut τv̆,
+                    τ,
+                    ε,
+                    &config.insertion,
+                    &reg,
+                    &state,
+                    &mut stats,
+                )?;
+            }
+            Ok(())
+        })?;
 
-        if config.insertion.extra_weight_optimisation_steps > 0 {
+        if did_merge.is_none() {
             γ.prune_compat(&mut μ, &mut stats);
         }
 
@@ -1185,5 +1185,14 @@
     P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
     Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
 {
-    pointsource_sliding_fb_dom(f, reg, prox_penalty, config, iterator, plotter, μ0, None)
+    pointsource_sliding_fb(
+        f,
+        reg,
+        prox_penalty,
+        config,
+        iterator,
+        plotter,
+        μ0,
+        &mut FullSpace,
+    )
 }
diff -r 51dbb2cc25d9 -r 8d1c959432c8 src/sliding_pdps.rs
--- a/src/sliding_pdps.rs	Tue Aug 04 13:46:29 2026 -0500
+++ b/src/sliding_pdps.rs	Tue Aug 04 14:12:18 2026 -0500
@@ -3,6 +3,7 @@
 primal-dual proximal splitting method.
 */
 
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::fb::*;
 use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess};
 use crate::measures::merging::SpikeMerging;
@@ -14,7 +15,6 @@
 use crate::types::*;
 use alg_tools::bounds::MinMaxMapping;
 use alg_tools::convex::{Conjugable, Prox, Zero};
-use alg_tools::norms::L2;
 use alg_tools::direct_product::Pair;
 use alg_tools::error::DynResult;
 use alg_tools::euclidean::ClosedEuclidean;
@@ -24,7 +24,9 @@
 };
 use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Instance};
 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};
 
@@ -96,7 +98,7 @@
     opKz: &KOpZ,
     fnR: &R,
     fnH: &H,
-    domain: Option<&P::Dom>,
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<MeasureZ<F, Z, N>>
 where
     F: Float + ToNalgebraRealField,
@@ -306,8 +308,8 @@
         // 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 m = prox_penalty.merge_spikes(
+        let did_merge = config.insertion.if_merge_now(&state, |ins| {
+            stats.merged += prox_penalty.merge_spikes(
                 &mut μ,
                 &mut τv̆,
                 &μ̆,
@@ -317,39 +319,43 @@
                 &reg,
                 Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z)) + reg.apply(μ̃)),
             );
-            if m > 0 {
-                stats.merged += m;
+            γ.prune_compat(&mut μ, &mut stats);
+
+            if config.insertion.random_insertions {
+                if let Some(x) = domain.random_point() {
+                    μ += DeltaMeasure { x, α: 0.0 };
+                    stats.inserted += 1;
+                }
             }
-        });
-
-        γ.prune_compat(&mut μ, &mut stats);
 
-        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,
-            )?;
+            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,
+                )?;
 
-            // 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;
-        }
+                // 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 config.insertion.extra_weight_optimisation_steps > 0 {
+        if did_merge.is_none() {
             γ.prune_compat(&mut μ, &mut stats);
         }
 
@@ -396,7 +402,7 @@
     (μ0, z): (Option<RNDM<N, F>>, Z),
     //opKμ : KOpM,
     fnR: &R,
-    domain: Option<&P::Dom>
+    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<MeasureZ<F, Z, N>>
 where
     F: Float + ToNalgebraRealField,
@@ -437,11 +443,10 @@
         &opKz,
         fnR,
         &fnH,
-        domain
+        domain,
     )
 }
 
-
 #[deprecated]
 /// Version of [`pointsource_sliding_pdps_pair`] without domain specification.
 #[replace_float_literals(F::cast_from(literal))]
@@ -517,7 +522,7 @@
         opKz,
         fnR,
         fnH,
-        None,
+        &mut FullSpace,
     )
 }
 
@@ -574,6 +579,6 @@
         &opKz,
         fnR,
         &fnH,
-        None
+        &mut FullSpace,
     )
 }