Mercurial > repos > pointsource_algs / changeset

--- a/src/lib.rs	Sun Jul 19 07:30:06 2026 +0200
+++ b/src/lib.rs	Mon Jul 13 12:09:57 2026 -0500
@@ -265,6 +265,10 @@
     #[arg(long)]
     /// Inner tolerance multiplier
     pub inner_tol: Option<F>,
+
+    #[arg(long)]
+    /// Use a conservative insertion strategy
+    pub conservative: Option<bool>,
 }
 
 /// A generic entry point for binaries based on this library

--- a/src/prox_penalty.rs	Sun Jul 19 07:30:06 2026 +0200
+++ b/src/prox_penalty.rs	Mon Jul 13 12:09:57 2026 -0500
@@ -66,8 +66,8 @@
     /// Iterations between merging heuristic tries
     pub merge_every: usize,
 
-    /// Additional weight optimisation steps
-    pub extra_weight_optimisation_steps: usize,
+    /// Conservative insertion strategy
+    pub conservative: bool,
 }
 
 #[replace_float_literals(F::cast_from(literal))]
@@ -87,6 +87,7 @@
             merge_every: 10,
             merge_tolerance_mult: 2.0,
             extra_weight_optimisation_steps: 0,
+            conservative: true,
         }
     }
 }

--- a/src/prox_penalty/radon_squared.rs	Sun Jul 19 07:30:06 2026 +0200
+++ b/src/prox_penalty/radon_squared.rs	Mon Jul 13 12:09:57 2026 -0500
@@ -59,9 +59,7 @@
         // If no merging heuristic is used, let's be more conservative about spike insertion,
         // and skip it after first round. If merging is done, being more greedy about spike
         // insertion also seems to improve performance.
-        let skip_by_rough_check = !config.merging.enabled;
-        let violation = reg.find_tolerance_violation(τv, τ, ε, skip_by_rough_check, config);
-        reg.solve_oc_radonsq(μ, τv, τ, ε, violation, config, stats);
+        reg.solve_oc_radonsq(μ, τv, τ, ε, config, stats);
 
         Ok(true)
     }

--- a/src/prox_penalty/wave.rs	Sun Jul 19 07:30:06 2026 +0200
+++ b/src/prox_penalty/wave.rs	Mon Jul 13 12:09:57 2026 -0500
@@ -106,11 +106,7 @@
             // If no merging heuristic is used, let's be more conservative about spike insertion,
             // and skip it after first round. If merging is done, being more greedy about spike
             // insertion also seems to improve performance.
-            let skip_by_rough_check = if config.merging.enabled {
-                false
-            } else {
-                count > 0
-            };
+            let skip_by_rough_check = config.conservative && (!config.merging.enabled || count > 0);
 
             // Find a spike to insert, if needed
             let (ξ, _v_ξ, in_bounds) =

--- a/src/regularisation.rs	Sun Jul 19 07:30:06 2026 +0200
+++ b/src/regularisation.rs	Mon Jul 13 12:09:57 2026 -0500
@@ -12,7 +12,7 @@
 use alg_tools::instance::{Instance, Space};
 use alg_tools::linops::Mapping;
 use alg_tools::loc::Loc;
-use alg_tools::norms::Norm;
+use alg_tools::norms::{Dist, Norm, L1};
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
 
@@ -23,7 +23,7 @@
 use alg_tools::bounds::{Bounds, MinMaxMapping};
 use alg_tools::iterate::AlgIteratorFactory;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
-use nalgebra::{DMatrix, DVector};
+use nalgebra::{DMatrix, DVector, DVectorView, DVectorViewMut};
 
 use std::cmp::Ordering::{Equal, Greater, Less};
 
@@ -186,11 +186,10 @@
         τv: &mut M,
         τ: F,
         ε: F,
-        tolerance_violation: Option<(Domain, F, bool)>,
         config: &InsertionConfig<F>,
         stats: &mut IterInfo<F>,
     ) where
-        M: Mapping<Domain, Codomain = F>;
+        M: MinMaxMapping<Domain, F>;
 
     /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
     ///
@@ -343,57 +342,69 @@
         τv: &mut M,
         τ: F,
         ε: F,
-        tolerance_violation: Option<(Loc<N, F>, F, bool)>,
         config: &InsertionConfig<F>,
         stats: &mut IterInfo<F>,
     ) where
-        M: Mapping<Loc<N, F>, Codomain = F>,
+        M: MinMaxMapping<Loc<N, F>, F>,
     {
+        let skip_by_rough_check = config.conservative && !config.merging.enabled;
+
         let τα = τ * self.α();
-        let mut g: Vec<_> = μ
+        let mut g = μ
             .iter_locations()
             .map(|ζ| F::to_nalgebra_mixed(-τv.apply(ζ)))
-            .collect();
+            .collect::<Vec<_>>();
+
+        // 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 mut y = μ.masses_vec();
+        let mut x = y.clone();
+        // Solve finite-dimensional subproblem.
+        let inner_tolerance = ε * config.inner.tolerance_mult;
+        let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
 
-        let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation {
-            // Don't insert if existing spikes are almost as good
-            if g.iter().all(|minus_τv| {
-                -F::from_nalgebra_mixed(*minus_τv) > v_ξ + ε * config.refinement.tolerance_mult
-            }) {
-                // Weight is found out by running the finite-dimensional optimisation algorithm
-                // above
-                // NOTE: cannot set α here before y is extracted
+        let slack = if μ.len() > 0 && config.conservative {
+            let l = x.len();
+            let mut x_na = DVectorViewMut::from_slice(x.as_mut_slice(), l);
+            let y_na = DVectorView::from_slice(y.as_slice(), y.len());
+            let g_na = DVectorView::from_slice(g.as_slice(), g.len());
+            stats.inner_iters +=
+                l1squared_nonneg(&y_na, &g_na, τα, &mut x_na, &config.inner, inner_it);
+            F::from_nalgebra_mixed(y_na.dist(&x_na, L1))
+        } else {
+            0.0
+        };
+
+        let opt = match self.find_tolerance_violation(τv, τ, ε + slack, skip_by_rough_check, config)
+        {
+            None => !config.conservative,
+            Some((ξ, v_ξ, _in_bounds)) => {
+                // Don't insert if existing spikes are almost as good
+                // if g.iter().all(|minus_τv| {
+                //    -F::from_nalgebra_mixed(*minus_τv) > v_ξ + ε * config.refinement.tolerance_mult
+                //}) {
+
                 *μ += DeltaMeasure { x: ξ, α: 0.0 /*-v_ξ - τα*/ };
                 g.push(F::to_nalgebra_mixed(-v_ξ));
-                Some(-v_ξ - τα)
-            } else {
-                None
+                y.push(F::to_nalgebra_mixed(0.0));
+                x.push(F::to_nalgebra_mixed(-v_ξ - τα));
+                true
             }
-        } else {
-            None
         };
 
-        // Optimise weights
-        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.
-            // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
-            // problems have not yet been updated to sign change.
-            let y = μ.masses_dvector();
-            let mut x = y.clone();
-            let g_na = DVector::from_vec(g);
-            if let (Some(β), Some(dest)) = (new_spike_initial_weight, x.as_mut_slice().last_mut())
-            {
-                *dest = F::to_nalgebra_mixed(β);
-            }
-            // Solve finite-dimensional subproblem.
-            let inner_tolerance = ε * config.inner.tolerance_mult;
-            let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-            stats.inner_iters += l1squared_nonneg(&y, &g_na, τα, &mut x, &config.inner, inner_it);
+        if opt {
+            let l = x.len();
+            let mut x_na = DVectorViewMut::from_slice(x.as_mut_slice(), l);
+            let y_na = DVectorView::from_slice(y.as_slice(), y.len());
+            let g_na = DVectorView::from_slice(g.as_slice(), g.len());
+            stats.inner_iters +=
+                l1squared_nonneg(&y_na, &g_na, τα, &mut x_na, &config.inner, inner_it);
+        }
 
-            // Update masses of μ based on solution of finite-dimensional subproblem.
-            μ.set_masses_dvector(&x);
-        }
+        // Update masses of μ based on solution of finite-dimensional subproblem.
+        μ.set_masses_vec(&x);
     }
 
     fn verify_merge_candidate_radonsq<M>(
@@ -602,59 +613,69 @@
         τv: &mut M,
         τ: F,
         ε: F,
-        tolerance_violation: Option<(Loc<N, F>, F, bool)>,
         config: &InsertionConfig<F>,
         stats: &mut IterInfo<F>,
     ) where
-        M: Mapping<Loc<N, F>, Codomain = F>,
+        M: MinMaxMapping<Loc<N, F>, F>,
     {
+        let skip_by_rough_check = config.conservative && !config.merging.enabled;
+
         let τα = τ * self.α();
-        let mut g: Vec<_> = μ
+        let mut g = μ
             .iter_locations()
-            .map(|ζ| F::to_nalgebra_mixed(τv.apply(-ζ)))
-            .collect();
+            .map(|ζ| F::to_nalgebra_mixed(-τv.apply(ζ)))
+            .collect::<Vec<_>>();
 
-        let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation {
-            // Don't insert if existing spikes are almost as good
-            let n = v_ξ.abs();
-            if g.iter().all(|minus_τv| {
-                F::from_nalgebra_mixed(*minus_τv).abs() < n - ε * config.refinement.tolerance_mult
-            }) {
-                // Weight is found out by running the finite-dimensional optimisation algorithm
-                // above
-                // NOTE: cannot initialise α before y is extracted.
-                *μ += DeltaMeasure { x: ξ, α: 0.0 /*-(n + τα) * v_ξ.signum()*/ };
-                g.push(F::to_nalgebra_mixed(-v_ξ));
-                Some(-(n + τα) * v_ξ.signum())
-            } else {
-                None
-            }
+        // 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 mut y = μ.masses_vec();
+        let mut x = y.clone();
+        // Solve finite-dimensional subproblem.
+        let inner_tolerance = ε * config.inner.tolerance_mult;
+        let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
+
+        let slack = if μ.len() > 0 && config.conservative {
+            let l = x.len();
+            let mut x_na = DVectorViewMut::from_slice(x.as_mut_slice(), l);
+            let y_na = DVectorView::from_slice(y.as_slice(), y.len());
+            let g_na = DVectorView::from_slice(g.as_slice(), g.len());
+            stats.inner_iters +=
+                l1squared_unconstrained(&y_na, &g_na, τα, &mut x_na, &config.inner, inner_it);
+            F::from_nalgebra_mixed(y_na.dist(&x_na, L1))
         } else {
-            None
+            0.0
         };
 
-        // Optimise weights
-        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.
-            // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
-            // problems have not yet been updated to sign change.
-            let y = μ.masses_dvector();
-            let mut x = y.clone();
-            if let (Some(β), Some(dest)) = (new_spike_initial_weight, x.as_mut_slice().last_mut())
-            {
-                *dest = F::to_nalgebra_mixed(β);
+        let opt = match self.find_tolerance_violation(τv, τ, ε + slack, skip_by_rough_check, config)
+        {
+            None => !config.conservative,
+            Some((ξ, v_ξ, _in_bounds)) => {
+                // Don't insert if existing spikes are almost as good
+                // if g.iter().all(|minus_τv| {
+                //    -F::from_nalgebra_mixed(*minus_τv) > v_ξ + ε * config.refinement.tolerance_mult
+                //}) {
+
+                *μ += DeltaMeasure { x: ξ, α: 0.0 /*-v_ξ - τα*/ };
+                g.push(F::to_nalgebra_mixed(-v_ξ));
+                y.push(F::to_nalgebra_mixed(0.0));
+                x.push(F::to_nalgebra_mixed(-v_ξ - τα));
+                true
             }
-            let g_na = DVector::from_vec(g);
-            // Solve finite-dimensional subproblem.
-            let inner_tolerance = ε * config.inner.tolerance_mult;
-            let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
+        };
+
+        if opt {
+            let l = x.len();
+            let mut x_na = DVectorViewMut::from_slice(x.as_mut_slice(), l);
+            let y_na = DVectorView::from_slice(y.as_slice(), y.len());
+            let g_na = DVectorView::from_slice(g.as_slice(), g.len());
             stats.inner_iters +=
-                l1squared_unconstrained(&y, &g_na, τα, &mut x, &config.inner, inner_it);
+                l1squared_unconstrained(&y_na, &g_na, τα, &mut x_na, &config.inner, inner_it);
+        }
 
-            // Update masses of μ based on solution of finite-dimensional subproblem.
-            μ.set_masses_dvector(&x);
-        }
+        // Update masses of μ based on solution of finite-dimensional subproblem.
+        μ.set_masses_vec(&x);
     }
 
     fn verify_merge_candidate_radonsq<M>(

--- a/src/run.rs	Sun Jul 19 07:30:06 2026 +0200
+++ b/src/run.rs	Mon Jul 13 12:09:57 2026 -0500
@@ -131,6 +131,7 @@
                 .as_ref()
                 .map(unpack_tolerance)
                 .unwrap_or(g.tolerance),
+            conservative: cli.conservative.unwrap_or(g.conservative),
             ..g
         };
         let override_transport = |g: TransportConfig<F>| TransportConfig {