Mercurial > repos > pointsource_algs / file diff

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