Mercurial > repos > pointsource_algs / changeset

--- a/src/fb.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/fb.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -188,17 +188,14 @@
 
     // Run the algorithm
     for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
+        let maybe_μ_base = config.merge_now(&state).then(|| μ.clone());
+        let μ_base_len = μ.len();
+
         // Calculate smooth part of surrogate model.
-        // TODO: optimise τ to be applied to residual.
         let mut τv = f.differential(&μ) * τ;
 
-        // Save current base point for merge
-        let μ_base_len = μ.len();
-        let maybe_μ_base = config.merge_now(&state).then(|| μ.clone());
-
-        // Insert and reweigh
-        let (maybe_d, _within_tolerances) = prox_penalty
-            .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
+        // Do spike insertiona nd finite-dimensional weight optimisation
+        prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
 
         stats.inserted += μ.len() - μ_base_len;
 
@@ -218,12 +215,22 @@
 
         stats.pruned += prune_with_stats(&mut μ);
 
+        // 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)?;
+        }
+
+        if config.extra_weight_optimisation_steps > 0 {
+            stats.pruned += prune_with_stats(&mut μ);
+        }
+
         let iter = state.iteration();
         stats.this_iters += 1;
 
         // Give statistics if needed
         state.if_verbose(|| {
-            plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ);
+            plotter.plot_spikes(iter, None, Some(&τv), &μ);
             full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
         });
 
@@ -292,14 +299,12 @@
 
     // Run the algorithm
     for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
+        let μ_base_len = μ.len();
         // Calculate smooth part of surrogate model.
         let mut τv = f.differential(&μ) * τ;
 
-        let μ_base_len = μ.len();
-
-        // Insert new spikes and reweigh
-        let (maybe_d, _within_tolerances) = prox_penalty
-            .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
+        // Do spike insertiona nd finite-dimensional weight optimisation
+        prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
 
         stats.inserted += μ.len() - μ_base_len;
 
@@ -310,6 +315,12 @@
             warned_merging = true;
         }
 
+        // 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)?;
+        }
+
         // Update inertial prameters
         let λ_prev = λ;
         λ = 2.0 * λ_prev / (λ_prev + (4.0 + λ_prev * λ_prev).sqrt());
@@ -334,7 +345,7 @@
 
         // Give statistics if needed
         state.if_verbose(|| {
-            plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ_prev);
+            plotter.plot_spikes(iter, None, Some(&τv), &μ_prev);
             full_stats(&μ_prev, ε, std::mem::replace(&mut stats, IterInfo::new()))
         });
 

--- a/src/forward_pdps.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/forward_pdps.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -197,15 +197,29 @@
     };
     let mut stats = IterInfo::new();
 
+    let zy_step = |mut z, mut y, mut z_tmp| {
+        // Do z variable primal update
+        opKz_adj.apply_add(&mut z_tmp, &y);
+        z_tmp.axpy(1.0, &z, -σ_p);
+        let z_new = fnR.prox(σ_p, z_tmp);
+        //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z);
+        // Do dual update
+        z.axpy(1.0 + ω, &z_new, -ω);
+        //z = (z - &z_new) * (-ω) + &z_new;
+        opKz.gemv(&mut y, σ_d, z, 1.0);
+        let y = starH.prox(σ_d, y);
+        (z_new, y)
+    };
+
     // Run the algorithm
     for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) {
         // Calculate initial transport
-        let Pair(v, mut z_tmp) = f.differential(Pair(&μ, &z));
+        let Pair(v, z_tmp) = f.differential(Pair(&μ, &z));
         let mut τv = v * τ;
         let μ_base = μ.clone();
 
         // 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,
             τ,
@@ -216,6 +230,8 @@
             &mut stats,
         )?;
 
+        (z, y) = zy_step(z, y, z_tmp);
+
         stats.inserted += μ.len() - μ_base.len();
 
         // Merge spikes.
@@ -241,17 +257,27 @@
         // Prune spikes with zero weight.
         stats.pruned += prune_with_stats(&mut μ);
 
-        // Do z variable primal update
-        opKz_adj.apply_add(&mut z_tmp, &y);
-        z_tmp.axpy(1.0, &z, -σ_p);
-        let z_new = fnR.prox(σ_p, z_tmp);
-        //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z);
-        // Do dual update
-        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, _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);
+        }
+
+        if config.insertion.extra_weight_optimisation_steps > 0 {
+            // Prune spikes with zero weight.
+            stats.pruned += prune_with_stats(&mut μ);
+        }
 
         // Update step length parameters
         // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ);
@@ -261,7 +287,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()))
         });
 

--- a/src/lib.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/lib.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -238,6 +238,10 @@
     /// Scaling heuristic factor
     pub scaling_heuristic: Option<F>,
 
+    #[arg(long)]
+    /// Extra finite dimensional steps to take
+    pub extra_weight_optimisation_steps: Option<usize>,
+
     #[arg(long, value_names = &["ε", "θ", "p"])]
     /// Set the tolerance to ε_k = ε/(1+θk)^p
     pub tolerance: Option<Vec<F>>,

--- a/src/pdps.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/pdps.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -206,8 +206,7 @@
         let μ_base = μ.clone();
 
         // Insert and reweigh
-        let (maybe_d, _within_tolerances) = prox_penalty
-            .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
+        prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, &reg, &state, &mut stats)?;
 
         // Prune and possibly merge spikes
         if config.merge_now(&state) {
@@ -234,7 +233,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(&μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
         });
 

--- a/src/prox_penalty.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/prox_penalty.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -65,8 +65,9 @@
 
     /// Iterations between merging heuristic tries
     pub merge_every: usize,
-    // /// Save $μ$ for postprocessing optimisation
-    // pub postprocessing : bool
+
+    /// Additional weight optimisation steps
+    pub extra_weight_optimisation_steps: usize,
 }
 
 #[replace_float_literals(F::cast_from(literal))]
@@ -85,7 +86,7 @@
             fitness_merging: false,
             merge_every: 10,
             merge_tolerance_mult: 2.0,
-            // postprocessing : false,
+            extra_weight_optimisation_steps: 0,
         }
     }
 }
@@ -118,6 +119,7 @@
     Reg: RegTerm<Domain, F>,
     Domain: Space + Clone,
 {
+    /// Unused, but required as a plotter parametrisation.
     type ReturnMapping: Mapping<Domain, Codomain = F>;
 
     /// Returns the type of this proximality penalty
@@ -126,8 +128,7 @@
     /// Insert new spikes into `μ` to approximately satisfy optimality conditions
     /// with the forward step term fixed to `τv`.
     ///
-    /// May return `τv + w` for `w` a subdifferential of the regularisation term `reg`,
-    /// as well as an indication of whether the tolerance bounds `ε` are satisfied.
+    /// Returns an indication of whether the tolerance bounds `ε` are satisfied.
     ///
     /// `τv` is mutable to allow [`alg_tools::bounds::MinMaxMapping`] optimisation to
     /// refine data. Actual values of `τv` are not supposed to be mutated.
@@ -144,7 +145,23 @@
         reg: &Reg,
         state: &AlgIteratorIteration<I>,
         stats: &mut IterInfo<F>,
-    ) -> DynResult<(Option<Self::ReturnMapping>, bool)>
+    ) -> DynResult<bool>
+    where
+        I: AlgIterator;
+
+    /// A variant of [`insert_and_reweigh`] that only does finite-dimensional weight optimisation,
+    /// without inserting spikes.
+    fn reweigh<I>(
+        &self,
+        μ: &mut DiscreteMeasure<Domain, F>,
+        τv: &mut PreadjointCodomain,
+        τ: F,
+        ε: F,
+        config: &InsertionConfig<F>,
+        reg: &Reg,
+        state: &AlgIteratorIteration<I>,
+        stats: &mut IterInfo<F>,
+    ) -> DynResult<()>
     where
         I: AlgIterator;
 

--- a/src/prox_penalty/radon_squared.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/prox_penalty/radon_squared.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -52,7 +52,7 @@
         reg: &Reg,
         _state: &AlgIteratorIteration<I>,
         stats: &mut IterInfo<F>,
-    ) -> DynResult<(Option<Self::ReturnMapping>, bool)>
+    ) -> DynResult<bool>
     where
         I: AlgIterator,
     {
@@ -63,7 +63,25 @@
         let violation = reg.find_tolerance_violation(τv, τ, ε, skip_by_rough_check, config);
         reg.solve_oc_radonsq(μ, τv, τ, ε, violation, config, stats);
 
-        Ok((None, true))
+        Ok(true)
+    }
+
+    fn reweigh<I>(
+        &self,
+        μ: &mut DiscreteMeasure<Domain, F>,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        config: &InsertionConfig<F>,
+        reg: &Reg,
+        _state: &AlgIteratorIteration<I>,
+        stats: &mut IterInfo<F>,
+    ) -> DynResult<()>
+    where
+        I: AlgIterator,
+    {
+        reg.solve_oc_radonsq(μ, τv, τ, ε, None, config, stats);
+        Ok(())
     }
 
     fn merge_spikes(

--- a/src/prox_penalty/wave.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/prox_penalty/wave.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -52,7 +52,7 @@
         reg: &Reg,
         state: &AlgIteratorIteration<I>,
         stats: &mut IterInfo<F>,
-    ) -> DynResult<(Option<Self::ReturnMapping>, bool)>
+    ) -> DynResult<bool>
     where
         I: AlgIterator,
     {
@@ -70,7 +70,7 @@
 
         // Add points to support until within error tolerance or maximum insertion count reached.
         let mut count = 0;
-        let (within_tolerances, d) = 'insertion: loop {
+        let within_tolerances = 'insertion: loop {
             if μ.len() > 0 {
                 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
                 // from the beginning of the iteration are all contained in the immutable c and g.
@@ -115,13 +115,13 @@
             // Find a spike to insert, if needed
             let (ξ, _v_ξ, in_bounds) =
                 match reg.find_tolerance_violation(&mut d, τ, ε, skip_by_rough_check, config) {
-                    None => break 'insertion (true, d),
+                    None => break 'insertion true,
                     Some(res) => res,
                 };
 
             // Break if maximum insertion count reached
             if count >= max_insertions {
-                break 'insertion (in_bounds, d);
+                break 'insertion in_bounds;
             }
 
             // No point in optimising the weight here; the finite-dimensional algorithm is fast.
@@ -139,7 +139,55 @@
             println!("{}", err.red());
         }
 
-        Ok((Some(d), within_tolerances))
+        Ok(within_tolerances)
+    }
+
+    fn reweigh<I>(
+        &self,
+        μ: &mut DiscreteMeasure<Domain, F>,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        config: &InsertionConfig<F>,
+        reg: &Reg,
+        _state: &AlgIteratorIteration<I>,
+        stats: &mut IterInfo<F>,
+    ) -> DynResult<()>
+    where
+        I: AlgIterator,
+    {
+        if μ.len() > 0 {
+            let op𝒟norm = self.opnorm_bound(Radon, Linfinity)?;
+            let ω0 = self.apply(&*μ);
+
+            // Form finite-dimensional subproblem. The subproblem references to the original μ^k
+            // from the beginning of the iteration are all contained in the immutable c and g.
+            // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
+            // problems have not yet been updated to sign change.
+            let à = self.findim_matrix(μ.iter_locations());
+            let g̃ = DVector::from_iterator(
+                μ.len(),
+                μ.iter_locations()
+                    .map(|ζ| ω0.apply(ζ) - τv.apply(ζ))
+                    .map(F::to_nalgebra_mixed),
+            );
+            let mut x = μ.masses_dvector();
+
+            // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃.
+            // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩
+            // ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟Cx|_∞ ≤  sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟| |Cx|_ℳ
+            // ≤ sup_{|z|_2 ≤ 1} |z|_1 |𝒟| |x|_1 ≤ sup_{|z|_2 ≤ 1} n |z|_2 |𝒟| |x|_2
+            // = n |𝒟| |x|_2, where n is the number of points. Therefore
+            let Ã_normest = op𝒟norm * F::cast_from(μ.len());
+
+            // Solve finite-dimensional subproblem.
+            stats.inner_iters += reg.solve_findim(&Ã, &g̃, τ, &mut x, Ã_normest, ε, config);
+
+            // Update masses of μ based on solution of finite-dimensional subproblem.
+            μ.set_masses_dvector(&x);
+        }
+
+        Ok(())
     }
 
     fn merge_spikes(

--- a/src/run.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/run.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -123,6 +123,9 @@
             final_merging: cli.final_merging.unwrap_or(g.final_merging),
             fitness_merging: cli.fitness_merging.unwrap_or(g.fitness_merging),
             inner: override_inner(g.inner),
+            extra_weight_optimisation_steps: cli
+                .extra_weight_optimisation_steps
+                .unwrap_or(g.extra_weight_optimisation_steps),
             tolerance: cli
                 .tolerance
                 .as_ref()
@@ -797,7 +800,6 @@
                     DefaultAlgorithm::SlidingFB,
                     DefaultAlgorithm::RadonSlidingFB,
                     DefaultAlgorithm::FW,
-                    DefaultAlgorithm::FWRelax,
                     DefaultAlgorithm::PDPS,
                 ],
                 DataTermType::L1 => &[DefaultAlgorithm::PDPS],

--- a/src/sliding_fb.rs	Thu Jul 16 19:51:51 2026 +0300
+++ b/src/sliding_fb.rs	Sun Jul 19 07:30:06 2026 +0200
@@ -1041,7 +1041,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̆, μ̆) = 'adapt_transport: loop {
+        let (mut τv̆, μ̆) = 'adapt_transport: loop {
             // Set initial guess for μ=μ^{k+1}.
             γ.μ̆_into(&mut μ);
             let μ̆ = μ.clone();
@@ -1054,7 +1054,7 @@
             let mut τv̆ = f.differential(&μ̆) * τ;
 
             // 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̆,
                 τ,
@@ -1081,7 +1081,7 @@
                 &config.insertion.refinement,
                 &mut attempts,
             ) {
-                break 'adapt_transport (maybe_d, within_tolerances, τv̆, μ̆);
+                break 'adapt_transport (τv̆, μ̆);
             }
 
             stats.get_transport_mut().readjustment_iters += 1;
@@ -1103,20 +1103,39 @@
                 &reg,
                 Some(|μ̃: &RNDM<N, F>| f.apply(μ̃)),
             );
-            if m > 0 {
-                stats.merged += m;
-                v = f.differential(&μ);
-            }
+            //if m > 0 {
+            stats.merged += m;
+            //v = f.differential(&μ);
+            //}
         }
 
         γ.prune_compat(&mut μ, &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,
+            )?;
+        }
+
+        if config.insertion.extra_weight_optimisation_steps > 0 {
+            γ.prune_compat(&mut μ, &mut stats);
+        }
+
         let iter = state.iteration();
         stats.this_iters += 1;
 
         // Give statistics if requested
         state.if_verbose(|| {
-            plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv̆), &μ);
+            plotter.plot_spikes(iter, None, Some(&τv̆), &μ);
             full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
         });
 

--- 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()))
         });