diff -r 677a5fd1b014 -r b921ed0ab99b src/regularisation.rs --- 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, stats: &mut IterInfo, ) where - M: Mapping; + M: MinMaxMapping; /// Verify that `d` is in bounds `ε` for a merge candidate `μ` /// @@ -343,57 +342,69 @@ τv: &mut M, τ: F, ε: F, - tolerance_violation: Option<(Loc, F, bool)>, config: &InsertionConfig, stats: &mut IterInfo, ) where - M: Mapping, Codomain = F>, + M: MinMaxMapping, 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::>(); + + // 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( @@ -602,59 +613,69 @@ τv: &mut M, τ: F, ε: F, - tolerance_violation: Option<(Loc, F, bool)>, config: &InsertionConfig, stats: &mut IterInfo, ) where - M: Mapping, Codomain = F>, + M: MinMaxMapping, 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::>(); - 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(