Mercurial > repos > pointsource_algs / changeset

--- a/src/forward_pdps.rs	Thu Feb 26 12:55:38 2026 -0500
+++ b/src/forward_pdps.rs	Thu Feb 26 13:05:07 2026 -0500
@@ -152,6 +152,7 @@
     // Set up parameters
     let bigM = 0.0; //opKμ.adjoint_product_bound(prox_penalty).unwrap().sqrt();
     let nKz = opKz.opnorm_bound(L2, L2)?;
+    let is_fb = nKz == 0.0;
     let idOpZ = IdOp::new();
     let opKz_adj = opKz.adjoint();
     let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?;
@@ -164,7 +165,7 @@
     // To do so, we first solve σ_p and σ_d from standard PDPS step length condition
     // ^^^^^ < 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 = if nKz == 0.0 { 0.0 } else { config.σd0 / nKz };
+    let σ_d = if is_fb { 0.0 } else { config.σd0 / nKz };
     let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
     // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0}
     // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L)
@@ -224,8 +225,16 @@
         // and not to performing any pruning. That is be to done below simultaneously for γ.
         let ins = &config.insertion;
         if ins.merge_now(&state) {
-            stats.merged +=
-                prox_penalty.merge_spikes_no_fitness(&mut μ, &mut τv, &μ_base, τ, ε, ins, &reg);
+            stats.merged += prox_penalty.merge_spikes(
+                &mut μ,
+                &mut τv,
+                &μ_base,
+                τ,
+                ε,
+                ins,
+                &reg,
+                is_fb.then_some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))),
+            );
         }
 
         // Prune spikes with zero weight.

--- a/src/sliding_pdps.rs	Thu Feb 26 12:55:38 2026 -0500
+++ b/src/sliding_pdps.rs	Thu Feb 26 13:05:07 2026 -0500
@@ -152,6 +152,7 @@
     let bigθ = 0.0; //opKμ.transport_lipschitz_factor(L2Squared);
     let bigM = 0.0; //opKμ.adjoint_product_bound(&op𝒟).unwrap().sqrt();
     let nKz = opKz.opnorm_bound(L2, L2)?;
+    let is_fb = nKz == 0.0;
     let ℓ = 0.0;
     let idOpZ = IdOp::new();
     let opKz_adj = opKz.adjoint();
@@ -166,7 +167,7 @@
     // To do so, we first solve σ_p and σ_d from standard PDPS step length condition
     // ^^^^^ < 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 = if nKz == 0.0 { 0.0 } else { config.σd0 / nKz };
+    let σ_d = if is_fb { 0.0 } else { config.σd0 / nKz };
     let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
     // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0}
     // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L)
@@ -293,7 +294,7 @@
         // 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 γ.
         if config.insertion.merge_now(&state) {
-            stats.merged += prox_penalty.merge_spikes_no_fitness(
+            stats.merged += prox_penalty.merge_spikes(
                 &mut μ,
                 &mut τv̆,
                 &μ̆,
@@ -301,6 +302,7 @@
                 ε,
                 &config.insertion,
                 &reg,
+                is_fb.then_some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))),
             );
         }