pointsource_algs: comparison src/regularisation.rs

-:677a5fd1b014
+:b921ed0ab99b
 use crate::sliding_fb::pointsource_sliding_fb_reg;
 use crate::types::*;
 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};
 use crate::subproblem::{
 l1squared_nonneg::l1squared_nonneg, l1squared_unconstrained::l1squared_unconstrained,
 nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained,
 };
 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};
 /// The regularisation term $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$ for [`pointsource_fb_reg`] and other
 /// algorithms.
 &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>;
+M: MinMaxMapping<Domain, F>;
 /// Verify that `d` is in bounds `ε` for a merge candidate `μ`
 ///
 /// This version is s used for Radon-norm squared proximal term in
 /// [`crate::prox_penalty::radon_squared`].
 &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>,
+M: MinMaxMapping<Loc<N, F>, F>,
 {
-let τα = τ * self.α();
+let skip_by_rough_check = config.conservative && !config.merging.enabled;
-let mut g: Vec<_> = μ
+let τα = τ * self.α();
+let mut g = μ
 .iter_locations()
 .map(|ζ| F::to_nalgebra_mixed(-τv.apply(ζ)))
-.collect();
+.collect::<Vec<_>>();
-let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation {
+// Form finite-dimensional subproblem. The subproblem references to the original μ^k
-// Don't insert if existing spikes are almost as good
+// from the beginning of the iteration are all contained in the immutable c and g.
-if g.iter().all(|minus_τv| {
+// TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
--F::from_nalgebra_mixed(*minus_τv) > v_ξ + ε * config.refinement.tolerance_mult
+// problems have not yet been updated to sign change.
-}) {
+let mut y = μ.masses_vec();
-// Weight is found out by running the finite-dimensional optimisation algorithm
+let mut x = y.clone();
-// above
+// Solve finite-dimensional subproblem.
-// NOTE: cannot set α here before y is extracted
+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_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_ξ - τα)
+y.push(F::to_nalgebra_mixed(0.0));
-} else {
+x.push(F::to_nalgebra_mixed(-v_ξ - τα));
-None
+true
 }
-} else {
-None
 };
-// Optimise weights
+if opt {
-if μ.len() > 0 {
+let l = x.len();
-// Form finite-dimensional subproblem. The subproblem references to the original μ^k
+let mut x_na = DVectorViewMut::from_slice(x.as_mut_slice(), l);
-// from the beginning of the iteration are all contained in the immutable c and g.
+let y_na = DVectorView::from_slice(y.as_slice(), y.len());
-// TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
+let g_na = DVectorView::from_slice(g.as_slice(), g.len());
-// problems have not yet been updated to sign change.
+stats.inner_iters +=
-let y = μ.masses_dvector();
+l1squared_nonneg(&y_na, &g_na, τα, &mut x_na, &config.inner, inner_it);
-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);
-// 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>(
 &self,
 d: &mut 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>,
+M: MinMaxMapping<Loc<N, F>, F>,
 {
-let τα = τ * self.α();
+let skip_by_rough_check = config.conservative && !config.merging.enabled;
-let mut g: Vec<_> = μ
+let τα = τ * self.α();
+let mut g = μ
 .iter_locations()
-.map(|ζ| F::to_nalgebra_mixed(τv.apply(-ζ)))
+.map(|ζ| F::to_nalgebra_mixed(-τv.apply(ζ)))
-.collect();
+.collect::<Vec<_>>();
-let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation {
+// Form finite-dimensional subproblem. The subproblem references to the original μ^k
-// Don't insert if existing spikes are almost as good
+// from the beginning of the iteration are all contained in the immutable c and g.
-let n = v_ξ.abs();
+// TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
-if g.iter().all(|minus_τv| {
+// problems have not yet been updated to sign change.
-F::from_nalgebra_mixed(*minus_τv).abs() < n - ε * config.refinement.tolerance_mult
+let mut y = μ.masses_vec();
-}) {
+let mut x = y.clone();
-// Weight is found out by running the finite-dimensional optimisation algorithm
+// Solve finite-dimensional subproblem.
-// above
+let inner_tolerance = ε * config.inner.tolerance_mult;
-// NOTE: cannot initialise α before y is extracted.
+let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
-*μ += DeltaMeasure { x: ξ, α: 0.0 /*-(n + τα) * v_ξ.signum()*/ };
+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 {
+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(-(n + τα) * v_ξ.signum())
+y.push(F::to_nalgebra_mixed(0.0));
-} else {
+x.push(F::to_nalgebra_mixed(-v_ξ - τα));
-None
+true
 }
-} else {
-None
 };
-// Optimise weights
+if opt {
-if μ.len() > 0 {
+let l = x.len();
-// Form finite-dimensional subproblem. The subproblem references to the original μ^k
+let mut x_na = DVectorViewMut::from_slice(x.as_mut_slice(), l);
-// from the beginning of the iteration are all contained in the immutable c and g.
+let y_na = DVectorView::from_slice(y.as_slice(), y.len());
-// TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
+let g_na = DVectorView::from_slice(g.as_slice(), g.len());
-// 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, τα, &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>(
 &self,
 d: &mut M,

Mercurial > repos > pointsource_algs / file comparison

comparison: src/regularisation.rs

src/regularisation.rs