diff -r df92e78cc3f4 -r 677a5fd1b014 src/prox_penalty/wave.rs
--- 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(