Mercurial > repos > pointsource_algs / file diff

--- a/src/regularisation.rs	Fri Aug 14 10:38:27 2026 -0500
+++ b/src/regularisation.rs	Fri Oct 02 12:03:25 2026 -0500
@@ -148,6 +148,20 @@
     where
         M: MinMaxMapping<Domain, F>;
 
+    /// Returns the minimum distance of the subdifferential of $τ[F+R]$ to zero,
+    /// where $F$ is the smooth data term (differential given in `τv/τ`), and $R$ is
+    /// this regulariser.
+    fn discrepancy<M>(
+        &self,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        μ: &DiscreteMeasure<Domain, F>,
+        config: &InsertionConfig<F>,
+    ) -> F
+    where
+        M: MinMaxMapping<Domain, F>;
+
     /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
     ///
     /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
@@ -290,6 +304,40 @@
         }
     }
 
+    fn discrepancy<M>(
+        &self,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        μ: &DiscreteMeasure<Loc<N, F>, F>,
+        config: &InsertionConfig<F>,
+    ) -> F
+    where
+        M: MinMaxMapping<Loc<N, F>, F>,
+    {
+        let refinement_tolerance = ε * config.refinement.tolerance_mult;
+        let τα = self.α()*τ;
+
+        let μ_disrepancy = μ
+            .iter_spikes()
+            .map(|δ| {
+                let v = τv.apply(&δ.x);
+                if δ.α > 0.0 {
+                    (v + τα).abs()
+                } else {
+                    (-v - τα).max(0.0)
+                }
+            })
+            .reduce(F::max)
+            .unwrap_or(0.0);
+
+        let (_, τv_min) = τv.minimise(refinement_tolerance, config.refinement.max_steps);
+        //let rest_discrepancy = (-d_min - self.α()).max(0.0);
+        //return rest_discrepancy.max(μ_disrepancy)
+        // Simplifies to
+        return (-τv_min - τα).max(μ_disrepancy);
+    }
+
     fn verify_merge_candidate<M>(
         &self,
         d: &mut M,
@@ -532,6 +580,54 @@
         }
     }
 
+    fn discrepancy<M>(
+        &self,
+        τv: &mut M,
+        τ: F,
+        ε: F,
+        μ: &DiscreteMeasure<Loc<N, F>, F>,
+        config: &InsertionConfig<F>,
+    ) -> F
+    where
+        M: MinMaxMapping<Loc<N, F>, F>,
+    {
+        use std::cmp::Ordering::*;
+        let refinement_tolerance = ε * config.refinement.tolerance_mult;
+        let τα = τ * self.α();
+
+        let threshold = |v: F, α: F| {
+            if v > α {
+                v - α
+            } else if v < -α {
+                -v - α
+            } else {
+                0.0
+            }
+        };
+
+        let μ_disrepancy = μ
+            .iter_spikes()
+            .map(|δ| {
+                let v = τv.apply(&δ.x);
+                match δ.α.total_cmp(&0.0) {
+                    Greater => (v + τα).abs(),
+                    Less => (v - τα).abs(),
+                    Equal => threshold(v, τα),
+                }
+            })
+            .reduce(F::max)
+            .unwrap_or(0.0);
+
+        let (_, τv_min) = τv.minimise(refinement_tolerance, config.refinement.max_steps);
+        let (_, τv_max) = τv.maximise(refinement_tolerance, config.refinement.max_steps);
+        //let rest_discrepancy = (-d_min - α).max(0.0);
+        //return rest_discrepancy.max(μ_disrepancy)
+        // Simplifies to
+        return threshold(τv_min, τα)
+            .max(threshold(τv_max, τα))
+            .max(μ_disrepancy);
+    }
+
     fn verify_merge_candidate<M>(
         &self,
         d: &mut M,