--- a/src/regularisation.rs Thu Aug 29 00:00:00 2024 -0500
+++ b/src/regularisation.rs Tue Dec 31 09:25:45 2024 -0500
@@ -5,13 +5,14 @@
use numeric_literals::replace_float_literals;
use serde::{Serialize, Deserialize};
use alg_tools::norms::Norm;
-use alg_tools::linops::Apply;
+use alg_tools::linops::Mapping;
+use alg_tools::instance::Instance;
use alg_tools::loc::Loc;
use crate::types::*;
use crate::measures::{
- DiscreteMeasure,
+ RNDM,
DeltaMeasure,
- Radon
+ Radon,
};
use crate::fb::FBGenericConfig;
#[allow(unused_imports)] // Used by documentation.
@@ -21,7 +22,6 @@
use nalgebra::{DVector, DMatrix};
use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::mapping::Mapping;
use alg_tools::bisection_tree::{
BTFN,
Bounds,
@@ -56,12 +56,13 @@
}
}
-impl<'a, F : Float, const N : usize> Apply<&'a DiscreteMeasure<Loc<F, N>, F>>
+impl<'a, F : Float, const N : usize> Mapping<RNDM<F, N>>
for NonnegRadonRegTerm<F> {
- type Output = F;
+ type Codomain = F;
- fn apply(&self, μ : &'a DiscreteMeasure<Loc<F, N>, F>) -> F {
- self.α() * μ.norm(Radon)
+ fn apply<I>(&self, μ : I) -> F
+ where I : Instance<RNDM<F, N>> {
+ self.α() * μ.eval(|x| x.norm(Radon))
}
}
@@ -81,12 +82,13 @@
}
}
-impl<'a, F : Float, const N : usize> Apply<&'a DiscreteMeasure<Loc<F, N>, F>>
+impl<'a, F : Float, const N : usize> Mapping<RNDM<F, N>>
for RadonRegTerm<F> {
- type Output = F;
+ type Codomain = F;
- fn apply(&self, μ : &'a DiscreteMeasure<Loc<F, N>, F>) -> F {
- self.α() * μ.norm(Radon)
+ fn apply<I>(&self, μ : I) -> F
+ where I : Instance<RNDM<F, N>> {
+ self.α() * μ.eval(|x| x.norm(Radon))
}
}
@@ -99,11 +101,12 @@
NonnegRadon(F),
}
-impl<'a, F : Float, const N : usize> Apply<&'a DiscreteMeasure<Loc<F, N>, F>>
+impl<'a, F : Float, const N : usize> Mapping<RNDM<F, N>>
for Regularisation<F> {
- type Output = F;
-
- fn apply(&self, μ : &'a DiscreteMeasure<Loc<F, N>, F>) -> F {
+ type Codomain = F;
+
+ fn apply<I>(&self, μ : I) -> F
+ where I : Instance<RNDM<F, N>> {
match *self {
Self::Radon(α) => RadonRegTerm(α).apply(μ),
Self::NonnegRadon(α) => NonnegRadonRegTerm(α).apply(μ),
@@ -111,9 +114,9 @@
}
}
-/// Abstraction of regularisation terms for [`generic_pointsource_fb_reg`].
+/// Abstraction of regularisation terms.
pub trait RegTerm<F : Float + ToNalgebraRealField, const N : usize>
-: for<'a> Apply<&'a DiscreteMeasure<Loc<F, N>, F>, Output = F> {
+: Mapping<RNDM<F, N>, Codomain = F> {
/// Approximately solve the problem
/// <div>$$
/// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + τ G(x)
@@ -171,8 +174,7 @@
) -> 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> {
+ 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)
}
@@ -199,8 +201,7 @@
) -> 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>;
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>;
/// Verify that `d` is in bounds `ε` for a merge candidate `μ`
@@ -209,36 +210,34 @@
fn verify_merge_candidate<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
+ μ : &RNDM<F, N>,
τ : F,
ε : F,
config : &FBGenericConfig<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>;
+ 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::radon_fb`].
- /// The [`DiscreteMeasure`]s `μ` and `radon_μ` are supposed to have same coordinates at
- /// same agreeing indices.
+ /// 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>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
+ μ : &RNDM<F, N>,
τ : F,
ε : F,
config : &FBGenericConfig<F>,
- radon_μ :&DiscreteMeasure<Loc<F, N>, 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>;
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>;
/// TODO: document this
@@ -258,7 +257,7 @@
fn goodness<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
+ μ : &RNDM<F, N>,
y : &Loc<F, N>,
z : &Loc<F, N>,
τ : F,
@@ -267,8 +266,7 @@
) -> 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>;
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N>;
/// Convert bound on the regulariser to a bond on the Radon norm
fn radon_norm_bound(&self, b : F) -> F;
@@ -323,69 +321,66 @@
) -> 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> {
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
let τα = τ * self.α();
- let keep_below = τα + slack + ε;
- let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
+ let keep_above = -τα - slack - ε;
+ let minimise_below = -τα - slack - ε * 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
// the insertion strategy, skip insertion.
- if skip_by_rough_check && d.bounds().upper() <= keep_below {
+ if skip_by_rough_check && d.bounds().lower() >= keep_above {
None
} else {
- // If the rough check didn't indicate no insertion needed, find maximising point.
- d.maximise_above(maximise_above, refinement_tolerance, config.refinement.max_steps)
- .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ <= keep_below))
+ // If the rough check didn't indicate no insertion needed, find minimising point.
+ d.minimise_below(minimise_below, refinement_tolerance, config.refinement.max_steps)
+ .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above))
}
}
fn verify_merge_candidate<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
+ μ : &RNDM<F, N>,
τ : F,
ε : F,
config : &FBGenericConfig<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> {
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
let merge_tolerance = config.merge_tolerance_mult * ε;
- let keep_below = τα + merge_tolerance;
- let keep_supp_above = τα - merge_tolerance;
+ let keep_above = -τα - merge_tolerance;
+ let keep_supp_below = -τα + merge_tolerance;
let bnd = d.bounds();
return (
- bnd.lower() >= keep_supp_above
+ bnd.upper() <= keep_supp_below
||
μ.iter_spikes().all(|&DeltaMeasure{ α, ref x }| {
- (α == 0.0) || d.apply(x) >= keep_supp_above
+ (α == 0.0) || d.apply(x) <= keep_supp_below
})
) && (
- bnd.upper() <= keep_below
+ bnd.lower() >= keep_above
||
- d.has_upper_bound(keep_below, refinement_tolerance, config.refinement.max_steps)
+ d.has_lower_bound(keep_above, refinement_tolerance, config.refinement.max_steps)
)
}
fn verify_merge_candidate_radonsq<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
+ μ : &RNDM<F, N>,
τ : F,
ε : F,
config : &FBGenericConfig<F>,
- radon_μ :&DiscreteMeasure<Loc<F, N>, 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> {
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
let merge_tolerance = config.merge_tolerance_mult * ε;
@@ -395,7 +390,8 @@
return {
μ.both_matching(radon_μ)
.all(|(α, rα, x)| {
- let v = d.apply(x);
+ let v = -d.apply(x); // TODO: observe ad hoc negation here, after minus_τv
+ // switch to τv.
let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) {
Greater => (τα, τα),
_ => (F::NEG_INFINITY, τα),
@@ -410,10 +406,10 @@
(l1 + l2 - merge_tolerance <= v) && (v <= u1 + u2 + merge_tolerance)
})
} && {
- let keep_below = τα + slack + merge_tolerance;
- bnd.upper() <= keep_below
+ let keep_above = -τα - slack - merge_tolerance;
+ bnd.lower() <= keep_above
||
- d.has_upper_bound(keep_below, refinement_tolerance, config.refinement.max_steps)
+ d.has_lower_bound(keep_above, refinement_tolerance, config.refinement.max_steps)
}
}
@@ -435,7 +431,7 @@
fn goodness<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- _μ : &DiscreteMeasure<Loc<F, N>, F>,
+ _μ : &RNDM<F, N>,
y : &Loc<F, N>,
z : &Loc<F, N>,
τ : F,
@@ -444,8 +440,7 @@
) -> 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> {
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
let w = |x| 1.0.min((ε + d.apply(x))/(τ * self.α()));
w(z) - w(y)
}
@@ -500,8 +495,7 @@
) -> 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> {
+ G::SupportType : Mapping<Loc<F, N>, Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
let τα = τ * self.α();
let keep_below = τα + slack + ε;
let keep_above = -(τα + slack) - ε;
@@ -538,15 +532,14 @@
fn verify_merge_candidate<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
+ μ : &RNDM<F, N>,
τ : F,
ε : F,
config : &FBGenericConfig<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> {
+ G::SupportType : Mapping<Loc<F, N>,Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
let merge_tolerance = config.merge_tolerance_mult * ε;
@@ -582,16 +575,15 @@
fn verify_merge_candidate_radonsq<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- μ : &DiscreteMeasure<Loc<F, N>, F>,
+ μ : &RNDM<F, N>,
τ : F,
ε : F,
config : &FBGenericConfig<F>,
- radon_μ : &DiscreteMeasure<Loc<F, N>, 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> {
+ G::SupportType : Mapping<Loc<F, N>,Codomain=F> + LocalAnalysis<F, Bounds<F>, N> {
let τα = τ * self.α();
let refinement_tolerance = ε * config.refinement.tolerance_mult;
let merge_tolerance = config.merge_tolerance_mult * ε;
@@ -647,7 +639,7 @@
fn goodness<G, BT>(
&self,
d : &mut BTFN<F, G, BT, N>,
- _μ : &DiscreteMeasure<Loc<F, N>, F>,
+ _μ : &RNDM<F, N>,
y : &Loc<F, N>,
z : &Loc<F, N>,
τ : F,