--- a/src/regularisation.rs Thu Feb 26 11:38:43 2026 -0500
+++ b/src/regularisation.rs Thu Feb 26 11:36:22 2026 -0500
@@ -128,23 +128,6 @@
config: &InsertionConfig<F>,
) -> usize;
- /// Approximately solve the problem
- /// <div>$$
- /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x)
- /// $$</div>
- /// for $G$ depending on the trait implementation.
- ///
- /// Returns the number of iterations taken.
- fn solve_findim_l1squared(
- &self,
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
- τ: F,
- x: &mut DVector<F::MixedType>,
- ε: F,
- config: &InsertionConfig<F>,
- ) -> usize;
-
/// Find a point where `d` may violate the tolerance `ε`.
///
/// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
@@ -163,33 +146,6 @@
config: &InsertionConfig<F>,
) -> Option<(Domain, F, bool)>
where
- M: MinMaxMapping<Domain, F>,
- {
- self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO)
- }
-
- /// Find a point where `d` may violate the tolerance `ε`.
- ///
- /// This version includes a `slack` parameter to expand the tolerances.
- /// It is used for Radon-norm squared proximal term in [`crate::prox_penalty::radon_squared`].
- ///
- /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we
- /// are in bounds. `ε` is the current main tolerance and `τ` a scaling factor for the
- /// regulariser.
- ///
- /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check
- /// terminating early. Otherwise returns a possibly violating point, the value of `d` there,
- /// and a boolean indicating whether the found point is in bounds.
- fn find_tolerance_violation_slack<M>(
- &self,
- d: &mut M,
- τ: F,
- ε: F,
- skip_by_rough_check: bool,
- config: &InsertionConfig<F>,
- slack: F,
- ) -> Option<(Domain, F, bool)>
- where
M: MinMaxMapping<Domain, F>;
/// Verify that `d` is in bounds `ε` for a merge candidate `μ`
@@ -206,6 +162,36 @@
where
M: MinMaxMapping<Domain, F>;
+ /// TODO: document this
+ fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
+
+ /// Returns a scaling factor for the tolerance sequence.
+ ///
+ /// Typically this is the regularisation parameter.
+ fn tolerance_scaling(&self) -> F;
+}
+
+/// Regularisation term that can be used with [`crate::prox_penalty::radon_squared::RadonSquared`]
+/// proximal penalty.
+pub trait RadonSquaredRegTerm<Domain, F = f64>: RegTerm<Domain, F>
+where
+ Domain: Space + Clone,
+ F: Float + ToNalgebraRealField,
+{
+ /// Adapt weights of μ, possibly insertion a new point at tolerance_violation (which should
+ /// be returned by [`RegTerm::find_tolerance_violation`].
+ fn solve_oc_radonsq<M>(
+ &self,
+ μ: &mut DiscreteMeasure<Domain, F>,
+ τv: &mut M,
+ τ: F,
+ ε: F,
+ tolerance_violation: Option<(Domain, F, bool)>,
+ config: &InsertionConfig<F>,
+ stats: &mut IterInfo<F>,
+ ) where
+ M: Mapping<Domain, Codomain = F>;
+
/// Verify that `d` is in bounds `ε` for a merge candidate `μ`
///
/// This version is s used for Radon-norm squared proximal term in
@@ -225,14 +211,6 @@
) -> bool
where
M: MinMaxMapping<Domain, F>;
-
- /// TODO: document this
- fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
-
- /// Returns a scaling factor for the tolerance sequence.
- ///
- /// Typically this is the regularisation parameter.
- fn tolerance_scaling(&self) -> F;
}
/// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`].
@@ -279,43 +257,23 @@
quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
}
- fn solve_findim_l1squared(
- &self,
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
- τ: F,
- x: &mut DVector<F::MixedType>,
- ε: F,
- config: &InsertionConfig<F>,
- ) -> usize {
- let inner_tolerance = ε * config.inner.tolerance_mult;
- let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
- l1squared_nonneg(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
- }
-
#[inline]
- fn find_tolerance_violation_slack<M>(
+ fn find_tolerance_violation<M>(
&self,
d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
config: &InsertionConfig<F>,
- slack: F,
) -> Option<(Loc<N, F>, F, bool)>
where
M: MinMaxMapping<Loc<N, F>, F>,
{
let τα = τ * self.α();
- let keep_above = -τα - slack - ε;
- let minimise_below = -τα - slack - ε * config.insertion_cutoff_factor;
+ let keep_above = -τα - ε;
+ let minimise_below = -τα - ε * config.insertion_cutoff_factor;
let refinement_tolerance = ε * config.refinement.tolerance_mult;
- // println!(
- // "keep_above: {keep_above}, rough lower bound: {}, tolerance: {ε}, slack: {slack}, τα: {τα}",
- // d.bounds().lower()
- // );
-
// If preliminary check indicates that we are in bounds, and if it otherwise matches
// the insertion strategy, skip insertion.
if skip_by_rough_check && d.bounds().lower() >= keep_above {
@@ -362,6 +320,80 @@
));
}
+ fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
+ let τα = τ * self.α();
+ Some(Bounds(τα - ε, τα + ε))
+ }
+
+ fn tolerance_scaling(&self) -> F {
+ self.α()
+ }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F: Float + ToNalgebraRealField, const N: usize> RadonSquaredRegTerm<Loc<N, F>, F>
+ for NonnegRadonRegTerm<F>
+{
+ fn solve_oc_radonsq<M>(
+ &self,
+ μ: &mut DiscreteMeasure<Loc<N, F>, F>,
+ τ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>,
+ {
+ let τα = τ * self.α();
+ let mut g: Vec<_> = μ
+ .iter_locations()
+ .map(|ζ| F::to_nalgebra_mixed(-τv.apply(ζ)))
+ .collect();
+
+ 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
+ *μ += DeltaMeasure { x: ξ, α: 0.0 /*-v_ξ - τα*/ };
+ g.push(F::to_nalgebra_mixed(-v_ξ));
+ Some(-v_ξ - τα)
+ } else {
+ None
+ }
+ } 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, τα, 1.0, &mut x, &config.inner, inner_it);
+
+ // Update masses of μ based on solution of finite-dimensional subproblem.
+ μ.set_masses_dvector(&x);
+ }
+ }
+
fn verify_merge_candidate_radonsq<M>(
&self,
d: &mut M,
@@ -407,15 +439,6 @@
)
};
}
-
- fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
- let τα = τ * self.α();
- Some(Bounds(τα - ε, τα + ε))
- }
-
- fn tolerance_scaling(&self) -> F {
- self.α()
- }
}
#[replace_float_literals(F::cast_from(literal))]
@@ -461,37 +484,22 @@
quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
}
- fn solve_findim_l1squared(
- &self,
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
- τ: F,
- x: &mut DVector<F::MixedType>,
- ε: F,
- config: &InsertionConfig<F>,
- ) -> usize {
- let inner_tolerance = ε * config.inner.tolerance_mult;
- let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
- l1squared_unconstrained(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
- }
-
- fn find_tolerance_violation_slack<M>(
+ fn find_tolerance_violation<M>(
&self,
d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
config: &InsertionConfig<F>,
- slack: F,
) -> Option<(Loc<N, F>, F, bool)>
where
M: MinMaxMapping<Loc<N, F>, F>,
{
let τα = τ * self.α();
- let keep_below = τα + slack + ε;
- let keep_above = -(τα + slack) - ε;
- let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
- let minimise_below = -(τα + slack) - ε * config.insertion_cutoff_factor;
+ let keep_below = τα + ε;
+ let keep_above = -τα - ε;
+ let maximise_above = τα + ε * config.insertion_cutoff_factor;
+ let minimise_below = -τα - ε * config.insertion_cutoff_factor;
let refinement_tolerance = ε * config.refinement.tolerance_mult;
// If preliminary check indicates that we are in bonds, and if it otherwise matches
@@ -568,6 +576,81 @@
));
}
+ fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
+ let τα = τ * self.α();
+ Some(Bounds(-τα - ε, τα + ε))
+ }
+
+ fn tolerance_scaling(&self) -> F {
+ self.α()
+ }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F: Float + ToNalgebraRealField, const N: usize> RadonSquaredRegTerm<Loc<N, F>, F>
+ for RadonRegTerm<F>
+{
+ fn solve_oc_radonsq<M>(
+ &self,
+ μ: &mut DiscreteMeasure<Loc<N, F>, F>,
+ τ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>,
+ {
+ let τα = τ * self.α();
+ let mut g: Vec<_> = μ
+ .iter_locations()
+ .map(|ζ| F::to_nalgebra_mixed(τv.apply(-ζ)))
+ .collect();
+
+ 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
+ }
+ } 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();
+ if let (Some(β), Some(dest)) = (new_spike_initial_weight, x.as_mut_slice().last_mut())
+ {
+ *dest = F::to_nalgebra_mixed(β);
+ }
+ 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);
+ stats.inner_iters +=
+ l1squared_unconstrained(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it);
+
+ // Update masses of μ based on solution of finite-dimensional subproblem.
+ μ.set_masses_dvector(&x);
+ }
+ }
+
fn verify_merge_candidate_radonsq<M>(
&self,
d: &mut M,
@@ -619,15 +702,6 @@
)
};
}
-
- fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> {
- let τα = τ * self.α();
- Some(Bounds(-τα - ε, τα + ε))
- }
-
- fn tolerance_scaling(&self) -> F {
- self.α()
- }
}
#[replace_float_literals(F::cast_from(literal))]