--- a/src/sliding_fb.rs Tue Aug 04 20:11:10 2026 -0500
+++ b/src/sliding_fb.rs Thu Aug 06 11:14:15 2026 -0500
@@ -24,11 +24,11 @@
use alg_tools::bounds::{Bounds, MinMaxMapping};
use alg_tools::error::DynResult;
use alg_tools::euclidean::Euclidean;
-use alg_tools::instance::Space;
+use alg_tools::instance::{Instance, Space};
use alg_tools::iterate::AlgIteratorFactory;
use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Mapping, RealMapping};
use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::Norm;
+use alg_tools::norms::{Dist, Norm, L2};
use anyhow::ensure;
use std::ops::ControlFlow;
@@ -193,11 +193,11 @@
}
#[derive(Clone, Debug, Serialize)]
-pub struct SingleTransport<Domain, F: Float> {
+pub struct SingleTransport<U, F: Float> {
/// Source point
- x: Domain,
+ x: U,
/// Target point
- y: Domain,
+ y: U,
/// Original mass
α_μ_orig: F,
/// Transported mass
@@ -211,8 +211,8 @@
}
#[derive(Clone, Debug, Serialize)]
-pub struct Transport<Domain, F: Float> {
- vec: Vec<SingleTransport<Domain, F>>,
+pub struct Transport<U, F: Float> {
+ vec: Vec<SingleTransport<U, F>>,
}
/// Whether partially transported points are allowed.
@@ -223,21 +223,21 @@
const MINIMAL_PARTIAL_TRANSPORT: bool = true;
const NEW_APPROACH: bool = true;
-pub trait TransportProxPenalty<Point, PreadjointCodomain, Reg, F = f64>:
- ProxPenalty<Point, PreadjointCodomain, Reg, F>
+pub trait TransportProxPenalty<U, PreadjointCodomain, Reg, F = f64>:
+ ProxPenalty<U, PreadjointCodomain, Reg, F>
where
F: Float + ToNalgebraRealField,
- Reg: SlidingRegTerm<Point, F>,
- Point: Space + Clone,
+ Reg: SlidingRegTerm<U, F>,
+ U: Space + Clone,
{
type TransportStepLength;
/// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
/// with step lengh τ and transport step length `θ_or_adaptive`.
- fn initial_transport<Dom: AdmissibleDomain<Point>>(
+ fn initial_transport<Dom: AdmissibleDomain<U>>(
&self,
- γ: &mut Transport<Point, F>,
- μ: &DiscreteMeasure<Point, F>,
+ γ: &mut Transport<U, F>,
+ μ: &DiscreteMeasure<U, F>,
ε: F,
τ: F,
τθ_or_adaptive: &mut Self::TransportStepLength,
@@ -249,9 +249,9 @@
/// A posteriori transport adaptation.
fn aposteriori_transport(
&self,
- γ: &mut Transport<Point, F>,
- μ: &DiscreteMeasure<Point, F>,
- μ̆: &DiscreteMeasure<Point, F>,
+ γ: &mut Transport<U, F>,
+ μ: &DiscreteMeasure<U, F>,
+ μ̆: &DiscreteMeasure<U, F>,
τv̆: &mut PreadjointCodomain,
v: &mut PreadjointCodomain,
extra: Option<F>,
@@ -377,14 +377,12 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<F, M, Reg, 𝒟, O, const N: usize> TransportProxPenalty<Loc<N, F>, M, Reg, F> for 𝒟
+impl<F, M, Reg, 𝒟, const N: usize> TransportProxPenalty<Loc<N, F>, M, Reg, F> for 𝒟
where
F: Float + ToNalgebraRealField,
𝒟: DiscreteMeasureOp<Loc<N, F>, F>,
𝒟::Codomain: RealMapping<N, F>,
M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>,
- for<'a> &'a M: std::ops::Add<𝒟::PreCodomain, Output = O>,
- O: MinMaxMapping<Loc<N, F>, F>,
Reg: SlidingRegTerm<Loc<N, F>, F>,
Self: ProxPenalty<Loc<N, F>, M, Reg, F>,
RNDM<N, F>: SpikeMerging<F>,
@@ -608,25 +606,37 @@
}
}
+impl<U, F: Float> Norm<Radon, F> for Transport<U, F> {
+ fn norm(&self, _: Radon) -> F {
+ self.iter().map(|ρ| ρ.α_γ.abs()).sum()
+ }
+}
+
+impl<U, F: Float> MulAssign<F> for Transport<U, F> {
+ fn mul_assign(&mut self, factor: F) {
+ for ρ in self.iter_mut() {
+ ρ.α_γ *= factor;
+ }
+ }
+}
+
#[replace_float_literals(F::cast_from(literal))]
-impl<const N: usize, F: Float> Transport<Loc<N, F>, F> {
+impl<U, F: Float> Transport<U, F> {
pub(crate) fn new() -> Self {
Transport { vec: Vec::new() }
}
- pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<Loc<N, F>, F>> {
+ pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<U, F>> {
self.vec.iter()
}
- pub(crate) fn iter_mut(
- &mut self,
- ) -> impl Iterator<Item = &'_ mut SingleTransport<Loc<N, F>, F>> {
+ pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &'_ mut SingleTransport<U, F>> {
self.vec.iter_mut()
}
pub(crate) fn extend<I>(&mut self, it: I)
where
- I: IntoIterator<Item = SingleTransport<Loc<N, F>, F>>,
+ I: IntoIterator<Item = SingleTransport<U, F>>,
{
self.vec.extend(it)
}
@@ -635,62 +645,44 @@
self.vec.len()
}
- // pub(crate) fn dist_matching(&self, μ: &RNDM<N, F>) -> F {
- // self.iter()
- // .zip(μ.iter_spikes())
- // .map(|(ρ, δ)| (ρ.α_γ - δ.α).abs())
- // .sum()
- // }
-
- /// Construct `μ̆`, replacing the contents of `μ`.
- pub(crate) fn μ̆_into(&self, μ: &mut RNDM<N, F>) {
- assert!(self.len() <= μ.len());
-
- // First transported points
- for (δ, ρ) in izip!(μ.iter_spikes_mut(), self.iter()) {
- if ρ.α_γ.abs() > 0.0 {
- // Transport – transported point
- δ.α = ρ.α_γ;
- δ.x = ρ.y;
- } else {
- // No transport – original point
- δ.α = ρ.α_μ_orig;
- δ.x = ρ.x;
- }
- }
-
- // Then source points with partial transport
- let mut i = self.len();
- // This can cause the number of points to explode, so cannot have partial transport.
- for ρ in self.iter() {
- let α = ρ.α_μ_orig - ρ.α_γ;
- if ρ.α_γ.abs() > F::EPSILON && α.abs() > F::EPSILON {
- let δ = DeltaMeasure { α, x: ρ.x };
- if i < μ.len() {
- μ[i] = δ;
- } else {
- μ.push(δ)
- }
- i += 1;
- }
- }
- μ.truncate(i);
- }
-
/// Returns $‖μ\^k - π\_♯\^0γ\^{k+1}‖$
pub(crate) fn μ0_minus_γ0_radon(&self) -> F {
self.vec.iter().map(|ρ| (ρ.α_μ_orig - ρ.α_γ).abs()).sum()
}
+ /// Prune spikes with zero weight. To maintain correct ordering between μ and γ, also the
+ /// latter needs to be pruned when μ is.
+ pub(crate) fn prune_compat(&mut self, μ: &mut DiscreteMeasure<U, F>, stats: &mut IterInfo<F>) {
+ assert!(self.vec.len() <= μ.len());
+ let old_len = μ.len();
+ for (ρ, δ) in self.vec.iter_mut().zip(μ.iter_spikes()) {
+ ρ.retain = δ.α != F::ZERO;
+ }
+ // Do not prune new insertions on first round.
+ μ.prune_until(self.vec.len());
+ stats.pruned += old_len - μ.len();
+ self.vec.retain(|ρ| ρ.retain);
+ assert!(self.vec.len() <= μ.len());
+ }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<U, F> Transport<U, F>
+where
+ U: Dist<L2, F>,
+ F: Float,
+ for<'a> &'a U: Instance<U>,
+{
/// Returns $∫ c_2 d|γ|$
pub(crate) fn c2integral(&self) -> F {
+ let sq2 = |x| (x * x) / 2.0;
self.vec
.iter()
- .map(|ρ| ρ.y.dist2_squared(&ρ.x) / 2.0 * ρ.α_γ.abs())
+ .map(|ρ| sq2(ρ.y.dist(&ρ.x, L2)) * ρ.α_γ.abs())
.sum()
}
- pub(crate) fn get_transport_stats(&self, stats: &mut IterInfo<F>, μ: &RNDM<N, F>) {
+ pub(crate) fn get_transport_stats(&self, stats: &mut IterInfo<F>, μ: &DiscreteMeasure<U, F>) {
// TODO: This doesn't take into account μ[i].α becoming zero in the latest tranport
// attempt, for i < self.len(), when a corresponding source term also exists with index
// j ≥ self.len(). For now, we let that be reflected in the prune count.
@@ -726,22 +718,55 @@
})
};
}
+}
- /// Prune spikes with zero weight. To maintain correct ordering between μ and γ, also the
- /// latter needs to be pruned when μ is.
- pub(crate) fn prune_compat(&mut self, μ: &mut RNDM<N, F>, stats: &mut IterInfo<F>) {
- assert!(self.vec.len() <= μ.len());
- let old_len = μ.len();
- for (ρ, δ) in self.vec.iter_mut().zip(μ.iter_spikes()) {
- ρ.retain = δ.α != F::ZERO;
+#[replace_float_literals(F::cast_from(literal))]
+impl<U: Clone, F: Float> Transport<U, F> {
+ // pub(crate) fn dist_matching(&self, μ: &RNDM<N, F>) -> F {
+ // self.iter()
+ // .zip(μ.iter_spikes())
+ // .map(|(ρ, δ)| (ρ.α_γ - δ.α).abs())
+ // .sum()
+ // }
+
+ /// Construct `μ̆`, replacing the contents of `μ`.
+ pub(crate) fn μ̆_into(&self, μ: &mut DiscreteMeasure<U, F>) {
+ assert!(self.len() <= μ.len());
+
+ // First transported points
+ for (δ, ρ) in izip!(μ.iter_spikes_mut(), self.iter()) {
+ if ρ.α_γ.abs() > 0.0 {
+ // Transport – transported point
+ δ.α = ρ.α_γ;
+ δ.x = ρ.y.clone();
+ } else {
+ // No transport – original point
+ δ.α = ρ.α_μ_orig;
+ δ.x = ρ.x.clone();
+ }
}
- // Do not prune new insertions on first round.
- μ.prune_until(self.vec.len());
- stats.pruned += old_len - μ.len();
- self.vec.retain(|ρ| ρ.retain);
- assert!(self.vec.len() <= μ.len());
+
+ // Then source points with partial transport
+ let mut i = self.len();
+ // This can cause the number of points to explode, so cannot have partial transport.
+ for ρ in self.iter() {
+ let α = ρ.α_μ_orig - ρ.α_γ;
+ if ρ.α_γ.abs() > F::EPSILON && α.abs() > F::EPSILON {
+ let δ = DeltaMeasure { α, x: ρ.x.clone() };
+ if i < μ.len() {
+ μ[i] = δ;
+ } else {
+ μ.push(δ)
+ }
+ i += 1;
+ }
+ }
+ μ.truncate(i);
}
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<const N: usize, F: Float> Transport<Loc<N, F>, F> {
/// Helper for initial transport. Called from [`TransportProxPenalty::initial_transport`].
pub(super) fn do_init_transport<D, Dom: AdmissibleDomain<Loc<N, F>>>(
&mut self,
@@ -955,20 +980,6 @@
}
}
-impl<const N: usize, F: Float> Norm<Radon, F> for Transport<Loc<N, F>, F> {
- fn norm(&self, _: Radon) -> F {
- self.iter().map(|ρ| ρ.α_γ.abs()).sum()
- }
-}
-
-impl<const N: usize, F: Float> MulAssign<F> for Transport<Loc<N, F>, F> {
- fn mul_assign(&mut self, factor: F) {
- for ρ in self.iter_mut() {
- ρ.α_γ *= factor;
- }
- }
-}
-
/// Iteratively solve the pointsource localisation problem using sliding forward-backward
/// splitting
///
@@ -976,27 +987,29 @@
/// may be given.
/// Inertia is currently not supported.
#[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_sliding_fb<F, I, Dat, Reg, Plot, P, const N: usize>(
+pub fn pointsource_sliding_fb<F, U, I, Dat, Reg, Plot, P>(
f: &Dat,
reg: &Reg,
prox_penalty: &P,
config: &SlidingFBConfig<F>,
iterator: I,
mut plotter: Plot,
- μ0: Option<RNDM<N, F>>,
- domain: &mut impl AdmissibleDomain<Loc<N, F>>,
-) -> DynResult<RNDM<N, F>>
+ μ0: Option<DiscreteMeasure<U, F>>,
+ domain: &mut impl AdmissibleDomain<U>,
+) -> DynResult<DiscreteMeasure<U, F>>
where
F: Float + ToNalgebraRealField,
+ U: Space + Clone + Dist<L2, F>,
+ for<'a> &'a U: Instance<U>,
I: AlgIteratorFactory<IterInfo<F>>,
- Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
+ Dat: DifferentiableMapping<DiscreteMeasure<U, F>, Codomain = F> + BoundedCurvature<F>,
Dat::DerivativeDomain:
- DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
+ DifferentiableMapping<U, Codomain = F> + ClosedMul<F> + MinMaxMapping<U, F>,
//for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
- RNDM<N, F>: SpikeMerging<F>,
- Reg: SlidingRegTerm<Loc<N, F>, F>,
- P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
- Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
+ DiscreteMeasure<U, F>: SpikeMerging<F>,
+ Reg: SlidingRegTerm<U, F>,
+ P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+ Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
{
// Check parameters
ensure!(config.τ0 > 0.0, "Invalid step length parameter");
@@ -1025,7 +1038,7 @@
let mut ε = tolerance.initial();
// Statistics
- let full_stats = |μ: &RNDM<N, F>, ε, stats| IterInfo {
+ let full_stats = |μ: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
value: f.apply(μ) + reg.apply(μ),
n_spikes: μ.len(),
ε,
@@ -1114,7 +1127,7 @@
ε,
ins,
®,
- Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃)),
+ Some(|μ̃: &DiscreteMeasure<U, F>| f.apply(μ̃) + reg.apply(μ̃)),
);
// Prune. Do not prune new insertions on first round.
γ.prune_compat(&mut μ, &mut stats);
@@ -1167,26 +1180,28 @@
#[deprecated]
/// Version of `pointsource_sliding_fb` without domain restriction
-pub fn pointsource_sliding_fb_reg<F, I, Dat, Reg, Plot, P, const N: usize>(
+pub fn pointsource_sliding_fb_reg<F, U, I, Dat, Reg, Plot, P>(
f: &Dat,
reg: &Reg,
prox_penalty: &P,
config: &SlidingFBConfig<F>,
iterator: I,
plotter: Plot,
- μ0: Option<RNDM<N, F>>,
-) -> DynResult<RNDM<N, F>>
+ μ0: Option<DiscreteMeasure<U, F>>,
+) -> DynResult<DiscreteMeasure<U, F>>
where
F: Float + ToNalgebraRealField,
+ U: Space + Clone + Dist<L2, F>,
+ for<'a> &'a U: Instance<U>,
I: AlgIteratorFactory<IterInfo<F>>,
- Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
+ Dat: DifferentiableMapping<DiscreteMeasure<U, F>, Codomain = F> + BoundedCurvature<F>,
Dat::DerivativeDomain:
- DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
+ DifferentiableMapping<U, Codomain = F> + ClosedMul<F> + MinMaxMapping<U, F>,
//for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
- RNDM<N, F>: SpikeMerging<F>,
- Reg: SlidingRegTerm<Loc<N, F>, F>,
- P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
- Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
+ DiscreteMeasure<U, F>: SpikeMerging<F>,
+ Reg: SlidingRegTerm<U, F>,
+ P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+ Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
{
pointsource_sliding_fb(
f,