--- a/src/regularisation.rs Tue Apr 08 13:31:39 2025 -0500
+++ b/src/regularisation.rs Fri May 08 16:47:58 2026 -0500
@@ -4,12 +4,12 @@
#[allow(unused_imports)] // Used by documentation.
use crate::fb::pointsource_fb_reg;
-use crate::fb::FBGenericConfig;
-use crate::measures::{DeltaMeasure, Radon, RNDM};
+use crate::fb::InsertionConfig;
+use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
#[allow(unused_imports)] // Used by documentation.
use crate::sliding_fb::pointsource_sliding_fb_reg;
use crate::types::*;
-use alg_tools::instance::Instance;
+use alg_tools::instance::{Instance, Space};
use alg_tools::linops::Mapping;
use alg_tools::loc::Loc;
use alg_tools::norms::Norm;
@@ -20,9 +20,7 @@
l1squared_nonneg::l1squared_nonneg, l1squared_unconstrained::l1squared_unconstrained,
nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained,
};
-use alg_tools::bisection_tree::{
- BTSearch, Bounded, Bounds, LocalAnalysis, P2Minimise, SupportGenerator, BTFN,
-};
+use alg_tools::bounds::{Bounds, MinMaxMapping};
use alg_tools::iterate::AlgIteratorFactory;
use alg_tools::nalgebra_support::ToNalgebraRealField;
use nalgebra::{DMatrix, DVector};
@@ -34,7 +32,7 @@
///
/// The only member of the struct is the regularisation parameter α.
#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
-pub struct NonnegRadonRegTerm<F: Float>(pub F /* α */);
+pub struct NonnegRadonRegTerm<F: Float = f64>(pub F /* α */);
impl<'a, F: Float> NonnegRadonRegTerm<F> {
/// Returns the regularisation parameter
@@ -44,12 +42,12 @@
}
}
-impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for NonnegRadonRegTerm<F> {
+impl<'a, F: Float, const N: usize> Mapping<RNDM<N, F>> for NonnegRadonRegTerm<F> {
type Codomain = F;
fn apply<I>(&self, μ: I) -> F
where
- I: Instance<RNDM<F, N>>,
+ I: Instance<RNDM<N, F>>,
{
self.α() * μ.eval(|x| x.norm(Radon))
}
@@ -59,7 +57,7 @@
///
/// The only member of the struct is the regularisation parameter α.
#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
-pub struct RadonRegTerm<F: Float>(pub F /* α */);
+pub struct RadonRegTerm<F: Float = f64>(pub F /* α */);
impl<'a, F: Float> RadonRegTerm<F> {
/// Returns the regularisation parameter
@@ -69,12 +67,12 @@
}
}
-impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for RadonRegTerm<F> {
+impl<'a, F: Float, const N: usize> Mapping<RNDM<N, F>> for RadonRegTerm<F> {
type Codomain = F;
fn apply<I>(&self, μ: I) -> F
where
- I: Instance<RNDM<F, N>>,
+ I: Instance<RNDM<N, F>>,
{
self.α() * μ.eval(|x| x.norm(Radon))
}
@@ -82,19 +80,19 @@
/// Regularisation term configuration
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
-pub enum Regularisation<F: Float> {
+pub enum Regularisation<F: Float = f64> {
/// $α \\|μ\\|\_{ℳ(Ω)}$
Radon(F),
/// $α\\|μ\\|\_{ℳ(Ω)} + δ_{≥ 0}(μ)$
NonnegRadon(F),
}
-impl<'a, F: Float, const N: usize> Mapping<RNDM<F, N>> for Regularisation<F> {
+impl<'a, F: Float, const N: usize> Mapping<RNDM<N, F>> for Regularisation<F> {
type Codomain = F;
fn apply<I>(&self, μ: I) -> F
where
- I: Instance<RNDM<F, N>>,
+ I: Instance<RNDM<N, F>>,
{
match *self {
Self::Radon(α) => RadonRegTerm(α).apply(μ),
@@ -104,8 +102,10 @@
}
/// Abstraction of regularisation terms.
-pub trait RegTerm<F: Float + ToNalgebraRealField, const N: usize>:
- Mapping<RNDM<F, N>, Codomain = F>
+pub trait RegTerm<Domain, F = f64>: Mapping<DiscreteMeasure<Domain, F>, Codomain = F>
+where
+ Domain: Space + Clone,
+ F: Float + ToNalgebraRealField,
{
/// Approximately solve the problem
/// <div>$$
@@ -125,24 +125,7 @@
x: &mut DVector<F::MixedType>,
mA_normest: F,
ε: F,
- config: &FBGenericConfig<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: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> usize;
/// Find a point where `d` may violate the tolerance `ε`.
@@ -154,85 +137,30 @@
/// 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<G, BT>(
+ fn find_tolerance_violation<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
+ d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
- config: &FBGenericConfig<F>,
- ) -> Option<(Loc<F, N>, F, bool)>
+ config: &InsertionConfig<F>,
+ ) -> Option<(Domain, F, bool)>
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
- {
- 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<G, BT>(
- &self,
- d: &mut BTFN<F, G, BT, N>,
- τ: F,
- ε: F,
- skip_by_rough_check: bool,
- config: &FBGenericConfig<F>,
- slack: F,
- ) -> Option<(Loc<F, N>, F, bool)>
- where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
+ M: MinMaxMapping<Domain, F>;
/// Verify that `d` is in bounds `ε` for a merge candidate `μ`
///
/// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
- fn verify_merge_candidate<G, BT>(
+ fn verify_merge_candidate<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ d: &mut M,
+ μ: &DiscreteMeasure<Domain, F>,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> bool
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
-
- /// 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`].
- /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have
- /// same coordinates at same agreeing indices.
- ///
- /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
- fn verify_merge_candidate_radonsq<G, BT>(
- &self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
- τ: F,
- ε: F,
- config: &FBGenericConfig<F>,
- radon_μ: &RNDM<F, N>,
- ) -> bool
- where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
+ M: MinMaxMapping<Domain, F>;
/// TODO: document this
fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
@@ -243,33 +171,76 @@
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
+ /// [`crate::prox_penalty::radon_squared`].
+ /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have
+ /// same coordinates at same agreeing indices.
+ ///
+ /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
+ fn verify_merge_candidate_radonsq<M>(
+ &self,
+ d: &mut M,
+ μ: &DiscreteMeasure<Domain, F>,
+ τ: F,
+ ε: F,
+ config: &InsertionConfig<F>,
+ radon_μ: &DiscreteMeasure<Domain, F>,
+ ) -> bool
+ where
+ M: MinMaxMapping<Domain, F>;
+}
+
/// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`].
-pub trait SlidingRegTerm<F: Float + ToNalgebraRealField, const N: usize>: RegTerm<F, N> {
+pub trait SlidingRegTerm<Domain, F = f64>: RegTerm<Domain, F>
+where
+ Domain: Space + Clone,
+ F: Float + ToNalgebraRealField,
+{
/// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation
/// term, such that $-ε ≤ τw - d ≤ ε$.
- fn goodness<G, BT>(
+ fn goodness<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
- y: &Loc<F, N>,
- z: &Loc<F, N>,
+ d: &mut M,
+ μ: &DiscreteMeasure<Domain, F>,
+ y: &Domain,
+ z: &Domain,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> F
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>;
+ M: MinMaxMapping<Domain, F>;
/// Convert bound on the regulariser to a bond on the Radon norm
fn radon_norm_bound(&self, b: F) -> F;
}
#[replace_float_literals(F::cast_from(literal))]
-impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for NonnegRadonRegTerm<F>
-where
- Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<Loc<N, F>, F>
+ for NonnegRadonRegTerm<F>
{
fn solve_findim(
&self,
@@ -279,45 +250,28 @@
x: &mut DVector<F::MixedType>,
mA_normest: F,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> usize {
let inner_tolerance = ε * config.inner.tolerance_mult;
let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
}
- fn solve_findim_l1squared(
+ #[inline]
+ fn find_tolerance_violation<M>(
&self,
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
- τ: F,
- x: &mut DVector<F::MixedType>,
- ε: F,
- config: &FBGenericConfig<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<G, BT>(
- &self,
- d: &mut BTFN<F, G, BT, N>,
+ d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
- config: &FBGenericConfig<F>,
- slack: F,
- ) -> Option<(Loc<F, N>, F, bool)>
+ config: &InsertionConfig<F>,
+ ) -> Option<(Loc<N, F>, F, bool)>
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ 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;
// If preliminary check indicates that we are in bounds, and if it otherwise matches
@@ -326,27 +280,26 @@
None
} else {
// If the rough check didn't indicate no insertion needed, find minimising point.
- d.minimise_below(
+ let res = d.minimise_below(
minimise_below,
refinement_tolerance,
config.refinement.max_steps,
- )
- .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
+ );
+
+ res.map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
}
}
- fn verify_merge_candidate<G, BT>(
+ fn verify_merge_candidate<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ d: &mut M,
+ μ: &RNDM<N, F>,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> bool
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ M: MinMaxMapping<Loc<N, F>, F>,
{
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
@@ -367,19 +320,91 @@
));
}
- fn verify_merge_candidate_radonsq<G, BT>(
+ 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,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ μ: &mut DiscreteMeasure<Loc<N, F>, F>,
+ τv: &mut M,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
- radon_μ: &RNDM<F, N>,
+ 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,
+ μ: &RNDM<N, F>,
+ τ: F,
+ ε: F,
+ config: &InsertionConfig<F>,
+ radon_μ: &RNDM<N, F>,
) -> bool
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ M: MinMaxMapping<Loc<N, F>, F>,
{
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
@@ -414,36 +439,24 @@
)
};
}
-
- 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> SlidingRegTerm<F, N> for NonnegRadonRegTerm<F>
-where
- Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<Loc<N, F>, F>
+ for NonnegRadonRegTerm<F>
{
- fn goodness<G, BT>(
+ fn goodness<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- _μ: &RNDM<F, N>,
- y: &Loc<F, N>,
- z: &Loc<F, N>,
+ d: &mut M,
+ _μ: &RNDM<N, F>,
+ y: &Loc<N, F>,
+ z: &Loc<N, F>,
τ: F,
ε: F,
- _config: &FBGenericConfig<F>,
+ _config: &InsertionConfig<F>,
) -> F
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ M: MinMaxMapping<Loc<N, F>, F>,
{
let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α()));
w(z) - w(y)
@@ -455,10 +468,7 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<F, N> for RadonRegTerm<F>
-where
- Cube<F, N>: P2Minimise<Loc<F, N>, F>,
-{
+impl<F: Float + ToNalgebraRealField, const N: usize> RegTerm<Loc<N, F>, F> for RadonRegTerm<F> {
fn solve_findim(
&self,
mA: &DMatrix<F::MixedType>,
@@ -467,46 +477,29 @@
x: &mut DVector<F::MixedType>,
mA_normest: F,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> usize {
let inner_tolerance = ε * config.inner.tolerance_mult;
let inner_it = config.inner.iterator_options.stop_target(inner_tolerance);
quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it)
}
- fn solve_findim_l1squared(
+ fn find_tolerance_violation<M>(
&self,
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
- τ: F,
- x: &mut DVector<F::MixedType>,
- ε: F,
- config: &FBGenericConfig<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<G, BT>(
- &self,
- d: &mut BTFN<F, G, BT, N>,
+ d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
- config: &FBGenericConfig<F>,
- slack: F,
- ) -> Option<(Loc<F, N>, F, bool)>
+ config: &InsertionConfig<F>,
+ ) -> Option<(Loc<N, F>, F, bool)>
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ 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
@@ -541,18 +534,16 @@
}
}
- fn verify_merge_candidate<G, BT>(
+ fn verify_merge_candidate<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ d: &mut M,
+ μ: &RNDM<N, F>,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> bool
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ M: MinMaxMapping<Loc<N, F>, F>,
{
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
@@ -585,19 +576,92 @@
));
}
- fn verify_merge_candidate_radonsq<G, BT>(
+ 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,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ μ: &mut DiscreteMeasure<Loc<N, F>, F>,
+ τv: &mut M,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
- radon_μ: &RNDM<F, N>,
+ 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,
+ μ: &RNDM<N, F>,
+ τ: F,
+ ε: F,
+ config: &InsertionConfig<F>,
+ radon_μ: &RNDM<N, F>,
) -> bool
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ M: MinMaxMapping<Loc<N, F>, F>,
{
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
@@ -638,36 +702,24 @@
)
};
}
-
- 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> SlidingRegTerm<F, N> for RadonRegTerm<F>
-where
- Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+impl<F: Float + ToNalgebraRealField, const N: usize> SlidingRegTerm<Loc<N, F>, F>
+ for RadonRegTerm<F>
{
- fn goodness<G, BT>(
+ fn goodness<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- _μ: &RNDM<F, N>,
- y: &Loc<F, N>,
- z: &Loc<F, N>,
+ d: &mut M,
+ _μ: &RNDM<N, F>,
+ y: &Loc<N, F>,
+ z: &Loc<N, F>,
τ: F,
ε: F,
- _config: &FBGenericConfig<F>,
+ _config: &InsertionConfig<F>,
) -> F
where
- BT: BTSearch<F, N, Agg = Bounds<F>>,
- G: SupportGenerator<F, N, Id = BT::Data>,
- G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+ M: MinMaxMapping<Loc<N, F>, F>,
{
let α = self.α();
let w = |x| {