--- a/src/sliding_fb.rs Thu Feb 26 11:38:43 2026 -0500
+++ b/src/sliding_fb.rs Thu Feb 26 11:36:22 2026 -0500
@@ -9,11 +9,12 @@
//use nalgebra::{DVector, DMatrix};
use itertools::izip;
use std::iter::Iterator;
+use std::ops::MulAssign;
use crate::fb::*;
use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess};
use crate::measures::merging::SpikeMerging;
-use crate::measures::{DiscreteMeasure, Radon, RNDM};
+use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
use crate::plot::Plotter;
use crate::prox_penalty::{ProxPenalty, StepLengthBound};
use crate::regularisation::SlidingRegTerm;
@@ -25,6 +26,7 @@
use alg_tools::nalgebra_support::ToNalgebraRealField;
use alg_tools::norms::Norm;
use anyhow::ensure;
+use std::ops::ControlFlow;
/// Transport settings for [`pointsource_sliding_fb_reg`].
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
@@ -36,6 +38,8 @@
pub adaptation: F,
/// A posteriori transport tolerance multiplier (C_pos)
pub tolerance_mult_con: F,
+ /// maximum number of adaptation iterations, until cancelling transport.
+ pub max_attempts: usize,
}
#[replace_float_literals(F::cast_from(literal))]
@@ -52,7 +56,7 @@
#[replace_float_literals(F::cast_from(literal))]
impl<F: Float> Default for TransportConfig<F> {
fn default() -> Self {
- TransportConfig { θ0: 0.9, adaptation: 0.9, tolerance_mult_con: 100.0 }
+ TransportConfig { θ0: 0.9, adaptation: 0.9, tolerance_mult_con: 100.0, max_attempts: 2 }
}
}
@@ -98,158 +102,417 @@
FullyAdaptive { l: F, max_transport: F, g: G },
}
-/// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
-/// with step lengh τ and transport step length `θ_or_adaptive`.
-#[replace_float_literals(F::cast_from(literal))]
-pub(crate) fn initial_transport<F, G, D, const N: usize>(
- γ1: &mut RNDM<N, F>,
- μ: &mut RNDM<N, F>,
- τ: F,
- θ_or_adaptive: &mut TransportStepLength<F, G>,
- v: D,
-) -> (Vec<F>, RNDM<N, F>)
-where
- F: Float + ToNalgebraRealField,
- G: Fn(F, F) -> F,
- D: DifferentiableRealMapping<N, F>,
-{
- use TransportStepLength::*;
+#[derive(Clone, Debug, Serialize)]
+pub struct SingleTransport<const N: usize, F: Float> {
+ /// Source point
+ x: Loc<N, F>,
+ /// Target point
+ y: Loc<N, F>,
+ /// Original mass
+ α_μ_orig: F,
+ /// Transported mass
+ α_γ: F,
+ /// Helper for pruning
+ prune: bool,
+}
+
+#[derive(Clone, Debug, Serialize)]
+pub struct Transport<const N: usize, F: Float> {
+ vec: Vec<SingleTransport<N, F>>,
+}
- // Save current base point and shift μ to new positions. Idea is that
- // μ_base(_masses) = μ^k (vector of masses)
- // μ_base_minus_γ0 = μ^k - π_♯^0γ^{k+1}
- // γ1 = π_♯^1γ^{k+1}
- // μ = μ^{k+1}
- let μ_base_masses: Vec<F> = μ.iter_masses().collect();
- let mut μ_base_minus_γ0 = μ.clone(); // Weights will be set in the loop below.
- // Construct μ^{k+1} and π_♯^1γ^{k+1} initial candidates
- //let mut sum_norm_dv = 0.0;
- let γ_prev_len = γ1.len();
- assert!(μ.len() >= γ_prev_len);
- γ1.extend(μ[γ_prev_len..].iter().cloned());
+/// Whether partiall transported points are allowed.
+///
+/// Partial transport can cause spike count explosion, so full or zero
+/// transport is generally preferred. If this is set to `true`, different
+/// transport adaptation heuristics will be used.
+const ALLOW_PARTIAL_TRANSPORT: bool = true;
+const MINIMAL_PARTIAL_TRANSPORT: bool = true;
+
+impl<const N: usize, F: Float> Transport<N, F> {
+ pub(crate) fn new() -> Self {
+ Transport { vec: Vec::new() }
+ }
- // Calculate initial transport and step length.
- // First calculate initial transported weights
- for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) {
- // If old transport has opposing sign, the new transport will be none.
- ρ.α = if (ρ.α > 0.0 && δ.α < 0.0) || (ρ.α < 0.0 && δ.α > 0.0) {
- 0.0
- } else {
- δ.α
- };
+ pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<N, F>> {
+ self.vec.iter()
+ }
+
+ pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &'_ mut SingleTransport<N, F>> {
+ self.vec.iter_mut()
+ }
+
+ pub(crate) fn extend<I>(&mut self, it: I)
+ where
+ I: IntoIterator<Item = SingleTransport<N, F>>,
+ {
+ self.vec.extend(it)
+ }
+
+ pub(crate) fn len(&self) -> usize {
+ self.vec.len()
}
- // Calculate transport rays.
- match *θ_or_adaptive {
- Fixed(θ) => {
- let θτ = τ * θ;
- for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) {
- ρ.x = δ.x - v.differential(&δ.x) * (ρ.α.signum() * θτ);
+ // pub(crate) fn dist_matching(&self, μ: &RNDM<N, F>) -> F {
+ // self.iter()
+ // .zip(μ.iter_spikes())
+ // .map(|(ρ, δ)| (ρ.α_γ - δ.α).abs())
+ // .sum()
+ // }
+
+ /// Construct `μ̆`, replacing the contents of `μ`.
+ #[replace_float_literals(F::cast_from(literal))]
+ 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;
}
}
- AdaptiveMax { l: ℓ_F, ref mut max_transport, g: ref calculate_θ } => {
- *max_transport = max_transport.max(γ1.norm(Radon));
- let θτ = τ * calculate_θ(ℓ_F, *max_transport);
- for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) {
- ρ.x = δ.x - v.differential(&δ.x) * (ρ.α.signum() * θτ);
+
+ // Then source points with partial transport
+ let mut i = self.len();
+ if ALLOW_PARTIAL_TRANSPORT {
+ // This can cause the number of points to explode, so cannot have partial transport.
+ for ρ in self.iter() {
+ let α = ρ.α_μ_orig - ρ.α_γ;
+ if ρ.α_γ.abs() > F::EPSILON && α != 0.0 {
+ let δ = DeltaMeasure { α, x: ρ.x };
+ if i < μ.len() {
+ μ[i] = δ;
+ } else {
+ μ.push(δ)
+ }
+ i += 1;
+ }
}
}
- FullyAdaptive { l: ref mut adaptive_ℓ_F, ref mut max_transport, g: ref calculate_θ } => {
- *max_transport = max_transport.max(γ1.norm(Radon));
- let mut θ = calculate_θ(*adaptive_ℓ_F, *max_transport);
- // Do two runs through the spikes to update θ, breaking if first run did not cause
- // a change.
- for _i in 0..=1 {
- let mut changes = false;
- for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) {
- let dv_x = v.differential(&δ.x);
- let g = &dv_x * (ρ.α.signum() * θ * τ);
- ρ.x = δ.x - g;
- let n = g.norm2();
- if n >= F::EPSILON {
- // Estimate Lipschitz factor of ∇v
- let this_ℓ_F = (dv_x - v.differential(&ρ.x)).norm2() / n;
- *adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F);
- θ = calculate_θ(*adaptive_ℓ_F, *max_transport);
- changes = true
+ μ.truncate(i);
+ }
+
+ /// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
+ /// with step lengh τ and transport step length `θ_or_adaptive`.
+ #[replace_float_literals(F::cast_from(literal))]
+ pub(crate) fn initial_transport<G, D>(
+ &mut self,
+ μ: &RNDM<N, F>,
+ _τ: F,
+ τθ_or_adaptive: &mut TransportStepLength<F, G>,
+ v: D,
+ ) where
+ G: Fn(F, F) -> F,
+ D: DifferentiableRealMapping<N, F>,
+ {
+ use TransportStepLength::*;
+
+ // Initialise transport structure weights
+ for (δ, ρ) in izip!(μ.iter_spikes(), self.iter_mut()) {
+ ρ.α_μ_orig = δ.α;
+ ρ.x = δ.x;
+ // If old transport has opposing sign, the new transport will be none.
+ ρ.α_γ = if (ρ.α_γ > 0.0 && δ.α < 0.0) || (ρ.α_γ < 0.0 && δ.α > 0.0) {
+ 0.0
+ } else {
+ δ.α
+ }
+ }
+
+ let γ_prev_len = self.len();
+ assert!(μ.len() >= γ_prev_len);
+ self.extend(μ[γ_prev_len..].iter().map(|δ| SingleTransport {
+ x: δ.x,
+ y: δ.x, // Just something, will be filled properly in the next phase
+ α_μ_orig: δ.α,
+ α_γ: δ.α,
+ prune: false,
+ }));
+
+ // Calculate transport rays.
+ match *τθ_or_adaptive {
+ Fixed(θ) => {
+ for ρ in self.iter_mut() {
+ ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * θ);
+ }
+ }
+ AdaptiveMax { l: ℓ_F, ref mut max_transport, g: ref calculate_θτ } => {
+ *max_transport = max_transport.max(self.norm(Radon));
+ let θτ = calculate_θτ(ℓ_F, *max_transport);
+ for ρ in self.iter_mut() {
+ ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * θτ);
+ }
+ }
+ FullyAdaptive {
+ l: ref mut adaptive_ℓ_F,
+ ref mut max_transport,
+ g: ref calculate_θτ,
+ } => {
+ *max_transport = max_transport.max(self.norm(Radon));
+ let mut θτ = calculate_θτ(*adaptive_ℓ_F, *max_transport);
+ // Do two runs through the spikes to update θ, breaking if first run did not cause
+ // a change.
+ for _i in 0..=1 {
+ let mut changes = false;
+ for ρ in self.iter_mut() {
+ let dv_x = v.differential(&ρ.x);
+ let g = &dv_x * (ρ.α_γ.signum() * θτ);
+ ρ.y = ρ.x - g;
+ let n = g.norm2();
+ if n >= F::EPSILON {
+ // Estimate Lipschitz factor of ∇v
+ let this_ℓ_F = (dv_x - v.differential(&ρ.y)).norm2() / n;
+ *adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F);
+ θτ = calculate_θτ(*adaptive_ℓ_F, *max_transport);
+ changes = true
+ }
}
- }
- if !changes {
- break;
+ if !changes {
+ break;
+ }
}
}
}
}
- // Set initial guess for μ=μ^{k+1}.
- for (δ, ρ, &β) in izip!(μ.iter_spikes_mut(), γ1.iter_spikes(), μ_base_masses.iter()) {
- if ρ.α.abs() > F::EPSILON {
- δ.x = ρ.x;
- //δ.α = ρ.α; // already set above
- } else {
- δ.α = β;
+ /// A posteriori transport adaptation.
+ #[replace_float_literals(F::cast_from(literal))]
+ pub(crate) fn aposteriori_transport<D>(
+ &mut self,
+ μ: &RNDM<N, F>,
+ μ̆: &RNDM<N, F>,
+ _v: &mut D,
+ extra: Option<F>,
+ ε: F,
+ tconfig: &TransportConfig<F>,
+ attempts: &mut usize,
+ ) -> bool
+ where
+ D: DifferentiableRealMapping<N, F>,
+ {
+ *attempts += 1;
+
+ // 1. If π_♯^1γ^{k+1} = γ1 has non-zero mass at some point y, but μ = μ^{k+1} does not,
+ // then the ansatz ∇w̃_x(y) = w^{k+1}(y) may not be satisfied. So set the mass of γ1
+ // at that point to zero, and retry.
+ let mut all_ok = true;
+ for (δ, ρ) in izip!(μ.iter_spikes(), self.iter_mut()) {
+ if δ.α == 0.0 && ρ.α_γ != 0.0 {
+ all_ok = false;
+ ρ.α_γ = 0.0;
+ }
}
+
+ // 2. Through bounding ∫ B_ω(y, z) dλ(x, y, z).
+ // through the estimate ≤ C ‖Δ‖‖γ^{k+1}‖ for Δ := μ^{k+1}-μ̆^k
+ // which holds for some some C if the convolution kernel in 𝒟 has Lipschitz gradient.
+ let nγ = self.norm(Radon);
+ let nΔ = μ.dist_matching(&μ̆) + extra.unwrap_or(0.0);
+ let t = ε * tconfig.tolerance_mult_con;
+ if nγ * nΔ > t && *attempts >= tconfig.max_attempts {
+ all_ok = false;
+ } else if nγ * nΔ > t {
+ // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1,
+ // this will guarantee that eventually ‖γ‖ decreases sufficiently that we
+ // will not enter here.
+ //*self *= tconfig.adaptation * t / (nγ * nΔ);
+
+ // We want a consistent behaviour that has the potential to set many weights to zero.
+ // Therefore, we find the smallest uniform reduction `chg_one`, subtracted
+ // from all weights, that achieves total `adapt` adaptation.
+ let adapt_to = tconfig.adaptation * t / nΔ;
+ let reduction_target = nγ - adapt_to;
+ assert!(reduction_target > 0.0);
+ if ALLOW_PARTIAL_TRANSPORT {
+ if MINIMAL_PARTIAL_TRANSPORT {
+ // This reduces weights of transport, starting from … until `adapt` is
+ // exhausted. It will, therefore, only ever cause one extrap point insertion
+ // at the sources, unlike “full” partial transport.
+ //let refs = self.vec.iter_mut().collect::<Vec<_>>();
+ //refs.sort_by(|ρ1, ρ2| ρ1.α_γ.abs().partial_cmp(&ρ2.α_γ.abs()).unwrap());
+ // let mut it = refs.into_iter();
+ //
+ // Maybe sort by differential norm
+ // let mut refs = self
+ // .vec
+ // .iter_mut()
+ // .map(|ρ| {
+ // let val = v.differential(&ρ.x).norm2_squared();
+ // (ρ, val)
+ // })
+ // .collect::<Vec<_>>();
+ // refs.sort_by(|(_, v1), (_, v2)| v2.partial_cmp(&v1).unwrap());
+ // let mut it = refs.into_iter().map(|(ρ, _)| ρ);
+ let mut it = self.vec.iter_mut().rev();
+ let _unused = it.try_fold(reduction_target, |left, ρ| {
+ let w = ρ.α_γ.abs();
+ if left <= w {
+ ρ.α_γ = ρ.α_γ.signum() * (w - left);
+ ControlFlow::Break(())
+ } else {
+ ρ.α_γ = 0.0;
+ ControlFlow::Continue(left - w)
+ }
+ });
+ } else {
+ // This version equally reduces all weights. It causes partial transport, which
+ // has the problem that that we need to then adapt weights in both start and
+ // end points, in insert_and_reweigh, somtimes causing the number of spikes μ
+ // to explode.
+ let mut abs_weights = self
+ .vec
+ .iter()
+ .map(|ρ| ρ.α_γ.abs())
+ .filter(|t| *t > F::EPSILON)
+ .collect::<Vec<F>>();
+ abs_weights.sort_by(|a, b| a.partial_cmp(b).unwrap());
+ let n = abs_weights.len();
+ // Cannot have partial transport; can cause spike count explosion
+ let chg = abs_weights.into_iter().zip((1..=n).rev()).try_fold(
+ 0.0,
+ |smaller_total, (w, m)| {
+ let mf = F::cast_from(m);
+ let reduction = w * mf + smaller_total;
+ if reduction >= reduction_target {
+ ControlFlow::Break((reduction_target - smaller_total) / mf)
+ } else {
+ ControlFlow::Continue(smaller_total + w)
+ }
+ },
+ );
+ match chg {
+ ControlFlow::Continue(_) => self.vec.iter_mut().for_each(|δ| δ.α_γ = 0.0),
+ ControlFlow::Break(chg_one) => self.vec.iter_mut().for_each(|ρ| {
+ let t = ρ.α_γ.abs();
+ if t > 0.0 {
+ if ALLOW_PARTIAL_TRANSPORT {
+ let new = (t - chg_one).max(0.0);
+ ρ.α_γ = ρ.α_γ.signum() * new;
+ }
+ }
+ }),
+ }
+ }
+ } else {
+ // This version zeroes smallest weights, avoiding partial transport.
+ let mut abs_weights_idx = self
+ .vec
+ .iter()
+ .map(|ρ| ρ.α_γ.abs())
+ .zip(0..)
+ .filter(|(w, _)| *w >= 0.0)
+ .collect::<Vec<(F, usize)>>();
+ abs_weights_idx.sort_by(|(a, _), (b, _)| a.partial_cmp(b).unwrap());
+
+ let mut left = reduction_target;
+
+ for (w, i) in abs_weights_idx {
+ left -= w;
+ let ρ = &mut self.vec[i];
+ ρ.α_γ = 0.0;
+ if left < 0.0 {
+ break;
+ }
+ }
+ }
+
+ all_ok = false
+ }
+
+ if !all_ok && *attempts >= tconfig.max_attempts {
+ for ρ in self.iter_mut() {
+ ρ.α_γ = 0.0;
+ }
+ }
+
+ all_ok
}
- // Calculate μ^k-π_♯^0γ^{k+1} and v̆ = A_*(A[μ_transported + μ_transported_base]-b)
- μ_base_minus_γ0.set_masses(
- μ_base_masses
+
+ /// Returns $‖μ\^k - π\_♯\^0γ\^{k+1}‖$
+ pub(crate) fn μ0_minus_γ0_radon(&self) -> F {
+ self.vec.iter().map(|ρ| (ρ.α_μ_orig - ρ.α_γ).abs()).sum()
+ }
+
+ /// Returns $∫ c_2 d|γ|$
+ #[replace_float_literals(F::cast_from(literal))]
+ pub(crate) fn c2integral(&self) -> F {
+ self.vec
.iter()
- .zip(γ1.iter_masses())
- .map(|(&a, b)| a - b),
- );
- (μ_base_masses, μ_base_minus_γ0)
+ .map(|ρ| ρ.y.dist2_squared(&ρ.x) / 2.0 * ρ.α_γ.abs())
+ .sum()
+ }
+
+ #[replace_float_literals(F::cast_from(literal))]
+ pub(crate) fn get_transport_stats(&self, stats: &mut IterInfo<F>, μ: &RNDM<N, 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.
+ stats.inserted += μ.len() - self.len();
+
+ let transp = stats.get_transport_mut();
+
+ transp.dist = {
+ let (a, b) = transp.dist;
+ (a + self.c2integral(), b + self.norm(Radon))
+ };
+ transp.untransported_fraction = {
+ let (a, b) = transp.untransported_fraction;
+ let source = self.iter().map(|ρ| ρ.α_μ_orig.abs()).sum();
+ (a + self.μ0_minus_γ0_radon(), b + source)
+ };
+ transp.transport_error = {
+ let (a, b) = transp.transport_error;
+ //(a + self.dist_matching(&μ), b + self.norm(Radon))
+
+ // This ignores points that have been not transported at all, to only calculate
+ // destnation error; untransported_fraction accounts for not being able to transport
+ // at all.
+ self.iter()
+ .zip(μ.iter_spikes())
+ .fold((a, b), |(a, b), (ρ, δ)| {
+ let transported = ρ.α_γ.abs();
+ if transported > F::EPSILON {
+ (a + (ρ.α_γ - δ.α).abs(), b + transported)
+ } else {
+ (a, b)
+ }
+ })
+ };
+ }
+
+ /// 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()) {
+ ρ.prune = !(δ.α.abs() > F::EPSILON);
+ }
+ μ.prune_by(|δ| δ.α.abs() > F::EPSILON);
+ stats.pruned += old_len - μ.len();
+ self.vec.retain(|ρ| !ρ.prune);
+ assert!(self.vec.len() <= μ.len());
+ }
}
-/// A posteriori transport adaptation.
-#[replace_float_literals(F::cast_from(literal))]
-pub(crate) fn aposteriori_transport<F, const N: usize>(
- γ1: &mut RNDM<N, F>,
- μ: &mut RNDM<N, F>,
- μ_base_minus_γ0: &mut RNDM<N, F>,
- μ_base_masses: &Vec<F>,
- extra: Option<F>,
- ε: F,
- tconfig: &TransportConfig<F>,
-) -> bool
-where
- F: Float + ToNalgebraRealField,
-{
- // 1. If π_♯^1γ^{k+1} = γ1 has non-zero mass at some point y, but μ = μ^{k+1} does not,
- // then the ansatz ∇w̃_x(y) = w^{k+1}(y) may not be satisfied. So set the mass of γ1
- // at that point to zero, and retry.
- let mut all_ok = true;
- for (α_μ, α_γ1) in izip!(μ.iter_masses(), γ1.iter_masses_mut()) {
- if α_μ == 0.0 && *α_γ1 != 0.0 {
- all_ok = false;
- *α_γ1 = 0.0;
+impl<const N: usize, F: Float> Norm<Radon, F> for Transport<N, F> {
+ fn norm(&self, _: Radon) -> F {
+ self.iter().map(|ρ| ρ.α_γ.abs()).sum()
+ }
+}
+
+impl<const N: usize, F: Float> MulAssign<F> for Transport<N, F> {
+ fn mul_assign(&mut self, factor: F) {
+ for ρ in self.iter_mut() {
+ ρ.α_γ *= factor;
}
}
-
- // 2. Through bounding ∫ B_ω(y, z) dλ(x, y, z).
- // through the estimate ≤ C ‖Δ‖‖γ^{k+1}‖ for Δ := μ^{k+1}-μ^k-(π_♯^1-π_♯^0)γ^{k+1},
- // which holds for some some C if the convolution kernel in 𝒟 has Lipschitz gradient.
- let nγ = γ1.norm(Radon);
- let nΔ = μ_base_minus_γ0.norm(Radon) + μ.dist_matching(&γ1) + extra.unwrap_or(0.0);
- let t = ε * tconfig.tolerance_mult_con;
- if nγ * nΔ > t {
- // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1,
- // this will guarantee that eventually ‖γ‖ decreases sufficiently that we
- // will not enter here.
- *γ1 *= tconfig.adaptation * t / (nγ * nΔ);
- all_ok = false
- }
-
- if !all_ok {
- // Update weights for μ_base_minus_γ0 = μ^k - π_♯^0γ^{k+1}
- μ_base_minus_γ0.set_masses(
- μ_base_masses
- .iter()
- .zip(γ1.iter_masses())
- .map(|(&a, b)| a - b),
- );
- }
-
- all_ok
}
/// Iteratively solve the pointsource localisation problem using sliding forward-backward
@@ -284,7 +547,7 @@
// Initialise iterates
let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
- let mut γ1 = DiscreteMeasure::new();
+ let mut γ = Transport::new();
// Set up parameters
// let opAnorm = opA.opnorm_bound(Radon, L2);
@@ -293,24 +556,27 @@
//let ℓ = opA.transport.lipschitz_factor(L2Squared) * max_transport;
let ℓ = 0.0;
let τ = config.τ0 / prox_penalty.step_length_bound(&f)?;
- let (maybe_ℓ_F, maybe_transport_lip) = f.curvature_bound_components(config.guess);
- let transport_lip = maybe_transport_lip?;
- let calculate_θ = |ℓ_F, max_transport| {
- let ℓ_r = transport_lip * max_transport;
- config.transport.θ0 / (τ * (ℓ + ℓ_F + ℓ_r))
- };
- let mut θ_or_adaptive = match maybe_ℓ_F {
- //Some(ℓ_F0) => TransportStepLength::Fixed(calculate_θ(ℓ_F0 * b.norm2(), 0.0)),
- Ok(ℓ_F) => TransportStepLength::AdaptiveMax {
- l: ℓ_F, // TODO: could estimate computing the real reesidual
- max_transport: 0.0,
- g: calculate_θ,
- },
- Err(_) => TransportStepLength::FullyAdaptive {
- l: 10.0 * F::EPSILON, // Start with something very small to estimate differentials
- max_transport: 0.0,
- g: calculate_θ,
- },
+
+ let mut θ_or_adaptive = match f.curvature_bound_components(config.guess) {
+ (_, Err(_)) => TransportStepLength::Fixed(config.transport.θ0),
+ (maybe_ℓ_F, Ok(transport_lip)) => {
+ let calculate_θτ = move |ℓ_F, max_transport| {
+ let ℓ_r = transport_lip * max_transport;
+ config.transport.θ0 / (ℓ + ℓ_F + ℓ_r)
+ };
+ match maybe_ℓ_F {
+ Ok(ℓ_F) => TransportStepLength::AdaptiveMax {
+ l: ℓ_F, // TODO: could estimate computing the real reesidual
+ max_transport: 0.0,
+ g: calculate_θτ,
+ },
+ Err(_) => TransportStepLength::FullyAdaptive {
+ l: 10.0 * F::EPSILON, // Start with something very small to estimate differentials
+ max_transport: 0.0,
+ g: calculate_θτ,
+ },
+ }
+ }
};
// We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
// by τ compared to the conditional gradient approach.
@@ -331,25 +597,28 @@
for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
// Calculate initial transport
let v = f.differential(&μ);
- let (μ_base_masses, mut μ_base_minus_γ0) =
- initial_transport(&mut γ1, &mut μ, τ, &mut θ_or_adaptive, v);
+ γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v);
+
+ let mut attempts = 0;
// Solve finite-dimensional subproblem several times until the dual variable for the
// regularisation term conforms to the assumptions made for the transport above.
- let (maybe_d, _within_tolerances, mut τv̆) = 'adapt_transport: loop {
+ let (maybe_d, _within_tolerances, mut τv̆, μ̆) = 'adapt_transport: loop {
+ // Set initial guess for μ=μ^{k+1}.
+ γ.μ̆_into(&mut μ);
+ let μ̆ = μ.clone();
+
// Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b)
- //let residual_μ̆ = calculate_residual2(&γ1, &μ_base_minus_γ0, opA, b);
+ //let residual_μ̆ = calculate_residual2(&γ1, &μ0_minus_γ0, opA, b);
// TODO: this could be optimised by doing the differential like the
// old residual2.
- let μ̆ = &γ1 + &μ_base_minus_γ0;
- let mut τv̆ = f.differential(μ̆) * τ;
+ // NOTE: This assumes that μ = γ1
+ let mut τv̆ = f.differential(&μ̆) * τ;
// Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh(
&mut μ,
&mut τv̆,
- &γ1,
- Some(&μ_base_minus_γ0),
τ,
ε,
&config.insertion,
@@ -359,59 +628,33 @@
)?;
// A posteriori transport adaptation.
- if aposteriori_transport(
- &mut γ1,
- &mut μ,
- &mut μ_base_minus_γ0,
- &μ_base_masses,
- None,
- ε,
- &config.transport,
- ) {
- break 'adapt_transport (maybe_d, within_tolerances, τv̆);
+ if γ.aposteriori_transport(&μ, &μ̆, &mut τv̆, None, ε, &config.transport, &mut attempts)
+ {
+ break 'adapt_transport (maybe_d, within_tolerances, τv̆, μ̆);
}
+
+ stats.get_transport_mut().readjustment_iters += 1;
};
- stats.untransported_fraction = Some({
- assert_eq!(μ_base_masses.len(), γ1.len());
- let (a, b) = stats.untransported_fraction.unwrap_or((0.0, 0.0));
- let source = μ_base_masses.iter().map(|v| v.abs()).sum();
- (a + μ_base_minus_γ0.norm(Radon), b + source)
- });
- stats.transport_error = Some({
- assert_eq!(μ_base_masses.len(), γ1.len());
- let (a, b) = stats.transport_error.unwrap_or((0.0, 0.0));
- (a + μ.dist_matching(&γ1), b + γ1.norm(Radon))
- });
+ γ.get_transport_stats(&mut stats, &μ);
// Merge spikes.
// This crucially expects the merge routine to be stable with respect to spike locations,
// and not to performing any pruning. That is be to done below simultaneously for γ.
- let ins = &config.insertion;
- if ins.merge_now(&state) {
+ if config.insertion.merge_now(&state) {
stats.merged += prox_penalty.merge_spikes(
&mut μ,
&mut τv̆,
- &γ1,
- Some(&μ_base_minus_γ0),
+ &μ̆,
τ,
ε,
- ins,
+ &config.insertion,
®,
Some(|μ̃: &RNDM<N, F>| f.apply(μ̃)),
);
}
- // Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the
- // latter needs to be pruned when μ is.
- // TODO: This could do with a two-vector Vec::retain to avoid copies.
- let μ_new = DiscreteMeasure::from_iter(μ.iter_spikes().filter(|δ| δ.α != F::ZERO).cloned());
- if μ_new.len() != μ.len() {
- let mut μ_iter = μ.iter_spikes();
- γ1.prune_by(|_| μ_iter.next().unwrap().α != F::ZERO);
- stats.pruned += μ.len() - μ_new.len();
- μ = μ_new;
- }
+ γ.prune_compat(&mut μ, &mut stats);
let iter = state.iteration();
stats.this_iters += 1;