Mercurial > repos > pointsource_algs / file diff

--- a/src/sliding_pdps.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/sliding_pdps.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -242,7 +242,7 @@
 
         // Solve finite-dimensional subproblem several times until the dual variable for the
         // regularisation term conforms to the assumptions made for the transport above.
-        let (maybe_d, _within_tolerances, mut τv̆, z_new, μ̆) = 'adapt_transport: loop {
+        let (mut τv̆, z_new, μ̆) = 'adapt_transport: loop {
             // Set initial guess for μ=μ^{k+1}.
             γ.μ̆_into(&mut μ);
             let μ̆ = μ.clone();
@@ -253,7 +253,7 @@
             // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0);
 
             // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
-            let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh(
+            prox_penalty.insert_and_reweigh(
                 &mut μ,
                 &mut τv̆,
                 τ,
@@ -286,12 +286,20 @@
                 &config.insertion.refinement,
                 &mut attempts,
             ) {
-                break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new, μ̆);
+                break 'adapt_transport (τv̆, z_new, μ̆);
             }
         };
 
         γ.get_transport_stats(&mut stats, &μ);
 
+        // Do dual update
+        // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0);    // y = y + σ_d K[(1+ω)(μ,z)^{k+1}]
+        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;
+
         // 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 γ.
@@ -313,13 +321,34 @@
 
         γ.prune_compat(&mut μ, &mut stats);
 
-        // Do dual update
-        // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0);    // y = y + σ_d K[(1+ω)(μ,z)^{k+1}]
-        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;
+        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;
+        }
+
+        if config.insertion.extra_weight_optimisation_steps > 0 {
+            γ.prune_compat(&mut μ, &mut stats);
+        }
 
         // Update step length parameters
         // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ);
@@ -329,7 +358,7 @@
         stats.this_iters += 1;
 
         state.if_verbose(|| {
-            plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv̆), &μ);
+            plotter.plot_spikes(iter, None, Some(&τv̆), &μ);
             full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new()))
         });