--- a/src/regularisation.rs Sun Apr 27 15:03:51 2025 -0500
+++ b/src/regularisation.rs Thu Feb 26 11:38:43 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,7 +125,7 @@
x: &mut DVector<F::MixedType>,
mA_normest: F,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> usize;
/// Approximately solve the problem
@@ -142,7 +142,7 @@
τ: F,
x: &mut DVector<F::MixedType>,
ε: F,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
) -> usize;
/// Find a point where `d` may violate the tolerance `ε`.
@@ -154,18 +154,16 @@
/// 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>,
+ M: MinMaxMapping<Domain, F>,
{
self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO)
}
@@ -182,35 +180,31 @@
/// 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>(
+ fn find_tolerance_violation_slack<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
+ d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
slack: F,
- ) -> Option<(Loc<F, N>, F, bool)>
+ ) -> 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>;
+ 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>;
+ M: MinMaxMapping<Domain, F>;
/// Verify that `d` is in bounds `ε` for a merge candidate `μ`
///
@@ -220,19 +214,17 @@
/// 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>(
+ fn verify_merge_candidate_radonsq<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ d: &mut M,
+ μ: &DiscreteMeasure<Domain, F>,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
- radon_μ: &RNDM<F, N>,
+ config: &InsertionConfig<F>,
+ radon_μ: &DiscreteMeasure<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>;
+ M: MinMaxMapping<Domain, F>;
/// TODO: document this
fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>;
@@ -244,32 +236,33 @@
}
/// 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,7 +272,7 @@
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);
@@ -293,7 +286,7 @@
τ: F,
x: &mut DVector<F::MixedType>,
ε: 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);
@@ -301,52 +294,54 @@
}
#[inline]
- fn find_tolerance_violation_slack<G, BT>(
+ fn find_tolerance_violation_slack<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
+ d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
slack: F,
- ) -> Option<(Loc<F, N>, F, bool)>
+ ) -> 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 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 {
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 +362,17 @@
));
}
- fn verify_merge_candidate_radonsq<G, BT>(
+ fn verify_merge_candidate_radonsq<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ d: &mut M,
+ μ: &RNDM<N, F>,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
- radon_μ: &RNDM<F, N>,
+ 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;
@@ -426,24 +419,21 @@
}
#[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 +445,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,7 +454,7 @@
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);
@@ -481,26 +468,24 @@
τ: F,
x: &mut DVector<F::MixedType>,
ε: 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);
l1squared_unconstrained(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it)
}
- fn find_tolerance_violation_slack<G, BT>(
+ fn find_tolerance_violation_slack<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
+ d: &mut M,
τ: F,
ε: F,
skip_by_rough_check: bool,
- config: &FBGenericConfig<F>,
+ config: &InsertionConfig<F>,
slack: F,
- ) -> Option<(Loc<F, N>, F, bool)>
+ ) -> 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 + ε;
@@ -541,18 +526,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 +568,17 @@
));
}
- fn verify_merge_candidate_radonsq<G, BT>(
+ fn verify_merge_candidate_radonsq<M>(
&self,
- d: &mut BTFN<F, G, BT, N>,
- μ: &RNDM<F, N>,
+ d: &mut M,
+ μ: &RNDM<N, F>,
τ: F,
ε: F,
- config: &FBGenericConfig<F>,
- radon_μ: &RNDM<F, N>,
+ 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;
@@ -650,24 +631,21 @@
}
#[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| {