diff -r 4f468d35fa29 -r 7a8a55fd41c0 src/regularisation.rs --- 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, ) -> usize; - /// Approximately solve the problem - ///
$$ - /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x) - /// $$
- /// for $G$ depending on the trait implementation. - /// - /// Returns the number of iterations taken. - fn solve_findim_l1squared( - &self, - y: &DVector, - g: &DVector, - τ: F, - x: &mut DVector, - ε: F, - config: &InsertionConfig, - ) -> 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, ) -> Option<(Domain, F, bool)> where - M: MinMaxMapping, - { - 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( - &self, - d: &mut M, - τ: F, - ε: F, - skip_by_rough_check: bool, - config: &InsertionConfig, - slack: F, - ) -> Option<(Domain, F, bool)> - where M: MinMaxMapping; /// Verify that `d` is in bounds `ε` for a merge candidate `μ` @@ -206,6 +162,36 @@ where M: MinMaxMapping; + /// TODO: document this + fn target_bounds(&self, τ: F, ε: F) -> Option>; + + /// 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: RegTerm +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( + &self, + μ: &mut DiscreteMeasure, + τv: &mut M, + τ: F, + ε: F, + tolerance_violation: Option<(Domain, F, bool)>, + config: &InsertionConfig, + stats: &mut IterInfo, + ) where + M: Mapping; + /// 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; - - /// TODO: document this - fn target_bounds(&self, τ: F, ε: F) -> Option>; - - /// 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, - g: &DVector, - τ: F, - x: &mut DVector, - ε: F, - config: &InsertionConfig, - ) -> 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( + fn find_tolerance_violation( &self, d: &mut M, τ: F, ε: F, skip_by_rough_check: bool, config: &InsertionConfig, - slack: F, ) -> Option<(Loc, F, bool)> where M: MinMaxMapping, 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> { + let τα = τ * self.α(); + Some(Bounds(τα - ε, τα + ε)) + } + + fn tolerance_scaling(&self) -> F { + self.α() + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl RadonSquaredRegTerm, F> + for NonnegRadonRegTerm +{ + fn solve_oc_radonsq( + &self, + μ: &mut DiscreteMeasure, F>, + τv: &mut M, + τ: F, + ε: F, + tolerance_violation: Option<(Loc, F, bool)>, + config: &InsertionConfig, + stats: &mut IterInfo, + ) where + M: Mapping, 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( &self, d: &mut M, @@ -407,15 +439,6 @@ ) }; } - - fn target_bounds(&self, τ: F, ε: F) -> Option> { - 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, - g: &DVector, - τ: F, - x: &mut DVector, - ε: F, - config: &InsertionConfig, - ) -> 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( + fn find_tolerance_violation( &self, d: &mut M, τ: F, ε: F, skip_by_rough_check: bool, config: &InsertionConfig, - slack: F, ) -> Option<(Loc, F, bool)> where M: MinMaxMapping, 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> { + let τα = τ * self.α(); + Some(Bounds(-τα - ε, τα + ε)) + } + + fn tolerance_scaling(&self) -> F { + self.α() + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl RadonSquaredRegTerm, F> + for RadonRegTerm +{ + fn solve_oc_radonsq( + &self, + μ: &mut DiscreteMeasure, F>, + τv: &mut M, + τ: F, + ε: F, + tolerance_violation: Option<(Loc, F, bool)>, + config: &InsertionConfig, + stats: &mut IterInfo, + ) where + M: Mapping, 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( &self, d: &mut M, @@ -619,15 +702,6 @@ ) }; } - - fn target_bounds(&self, τ: F, ε: F) -> Option> { - let τα = τ * self.α(); - Some(Bounds(-τα - ε, τα + ε)) - } - - fn tolerance_scaling(&self) -> F { - self.α() - } } #[replace_float_literals(F::cast_from(literal))]