pointsource_algs: comparison src/sliding

-:e2953ffd4e0b
+:e9a460a0e638
 use crate::fb::*;
 use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess};
 use crate::measures::merging::SpikeMerging;
 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
 use crate::plot::Plotter;
-use crate::prox_penalty::{ProxPenalty, StepLengthBound};
+use crate::prox_penalty::{ProxPenalty, RadonSquared, StepLengthBound};
 use crate::regularisation::SlidingRegTerm;
+use crate::seminorms::DiscreteMeasureOp;
 use crate::types::*;
+use alg_tools::bounds::{Bounds, MinMaxMapping};
 use alg_tools::error::DynResult;
 use alg_tools::euclidean::Euclidean;
+use alg_tools::instance::Space;
 use alg_tools::iterate::AlgIteratorFactory;
-use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping};
+use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Mapping, RealMapping};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::norms::Norm;
 use anyhow::ensure;
 use std::ops::ControlFlow;
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
 pub struct TransportConfig<F: Float> {
 /// Transport step length $θ$ normalised to $(0, 1)$.
 pub θ0: F,
+/// Unnormalised transport step length $θ$. Overrides θ0.
+pub τθ: Option<F>,
 /// Factor in $(0, 1)$ for decreasing transport to adapt to tolerance.
 pub adaptation: F,
-/// A posteriori transport tolerance multiplier (C_pos)
+/// A posteriori transport tolerance multiplier
-pub tolerance_mult_con: F,
+pub tolerance_mult: F,
+/// Multiplier for rough estimate of ℓ_{∇v}. Should be ≥ 1.
+/// If explicit τθ is given, ℓ_{∇v} is one divided by this number times τθ.
+/// Otherwise, if  τθ relative to an estimate of ℓ_{∇v}, this number multiplies that estimate.
+pub ℓ_gradv_mult: F,
 /// maximum number of adaptation iterations, until cancelling transport.
 pub max_attempts: usize,
 /// Maximum number of failed transportations for a single source point
 pub max_fail: usize,
+/// Allow points to be transported partially.
+pub allow_partial_transport: bool,
+/// Use an alternative remainder control rule.
+pub alt_remainder_control: bool,
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F: Float> TransportConfig<F> {
 /// Check that the parameters are ok. Panics if not.
 pub fn check(&self) -> DynResult<()> {
 ensure!(self.θ0 > 0.0);
 ensure!(0.0 < self.adaptation && self.adaptation < 1.0);
-ensure!(self.tolerance_mult_con > 0.0);
+ensure!(self.tolerance_mult > 0.0);
 Ok(())
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F: Float> Default for TransportConfig<F> {
 fn default() -> Self {
 TransportConfig {
-θ0: 0.9,
+θ0: 0.99,
+τθ: None,
 adaptation: 0.9,
-tolerance_mult_con: 100.0,
+allow_partial_transport: true,
+alt_remainder_control: false,
+tolerance_mult: 1e1,
+ℓ_gradv_mult: 3.0,
 max_attempts: 2,
 max_fail: usize::MAX,
 }
 }
 }
 }
 }
 }
 /// Internal type of adaptive transport step length calculation
-pub(crate) enum TransportStepLength<F: Float, G: Fn(F, F) -> F> {
+#[derive(Clone, Debug, Serialize, Deserialize)]
+pub enum TransportStepLength<F: Float> {
 /// Fixed, known step length
 #[allow(dead_code)]
-Fixed(F),
+Fixed { τθ: F, ℓ_gradv: F },
+/// Simple step lengths that do not depend on maximum transport
+Simple { ℓ_gradv: F, τθ0: F, ℓ_base: F },
 /// Adaptive step length, only wrt. maximum transport.
-/// Content of `l` depends on use case, while `g` calculates the step length from `l`.
+AdaptiveMax {
-AdaptiveMax { l: F, max_transport: F, g: G },
+ℓ_gradv: F,
+adaptive_max_transport: F,
+τθ0: F,
+ℓ_base: F,
+ℓ_base_max_transport: F,
+},
 /// Adaptive step length.
-/// Content of `l` depends on use case, while `g` calculates the step length from `l`.
+FullyAdaptive {
-FullyAdaptive { l: F, max_transport: F, g: G },
+adaptive_ℓ_gradv: F,
+adaptive_max_transport: F,
+τθ0: F,
+ℓ_base: F,
+ℓ_base_max_transport: F,
+},
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F: Float> TransportStepLength<F> {
+fn get_ℓ_gradv(&self) -> F {
+use TransportStepLength::*;
+match *self {
+Fixed { ℓ_gradv, .. } => ℓ_gradv,
+Simple { ℓ_gradv, .. } => ℓ_gradv,
+AdaptiveMax { ℓ_gradv, .. } => ℓ_gradv,
+FullyAdaptive { adaptive_ℓ_gradv, .. } => adaptive_ℓ_gradv,
+}
+}
+fn new(
+maybe_ℓ_gradv_est: DynResult<F>,
+tconfig: &TransportConfig<F>,
+ℓ_base: F,
+ℓ_base_max_transport: F,
+) -> Self {
+if let Some(τθ) = tconfig.τθ {
+TransportStepLength::Fixed { τθ, ℓ_gradv: 1.0 / (τθ * tconfig.ℓ_gradv_mult) }
+} else {
+match maybe_ℓ_gradv_est {
+Ok(ℓ_gradv_est) => {
+if ℓ_base_max_transport == 0.0 {
+TransportStepLength::Simple {
+ℓ_gradv: tconfig.ℓ_gradv_mult * ℓ_gradv_est,
+τθ0: tconfig.θ0,
+ℓ_base,
+}
+} else {
+TransportStepLength::AdaptiveMax {
+ℓ_gradv: tconfig.ℓ_gradv_mult * ℓ_gradv_est,
+adaptive_max_transport: 0.0,
+τθ0: tconfig.θ0,
+ℓ_base,
+ℓ_base_max_transport,
+}
+}
+}
+Err(_) => TransportStepLength::FullyAdaptive {
+adaptive_ℓ_gradv: 10.0 * F::EPSILON, // Start with something very small to estimate differentials
+adaptive_max_transport: 0.0,
+τθ0: tconfig.θ0,
+ℓ_base,
+ℓ_base_max_transport,
+},
+}
+}
+}
 }
 #[derive(Clone, Debug, Serialize)]
-pub struct SingleTransport<const N: usize, F: Float> {
+pub struct SingleTransport<Domain, F: Float> {
 /// Source point
-x: Loc<N, F>,
+x: Domain,
 /// Target point
-y: Loc<N, F>,
+y: Domain,
 /// Original mass
 α_μ_orig: F,
 /// Transported mass
 α_γ: F,
 /// Helper for pruning
-prune: bool,
+retain: bool,
 /// Fail count
 fail_count: usize,
+/// Contribution to remainder (temporary variable)
+excess: F,
 }
 #[derive(Clone, Debug, Serialize)]
-pub struct Transport<const N: usize, F: Float> {
+pub struct Transport<Domain, F: Float> {
-vec: Vec<SingleTransport<N, F>>,
+vec: Vec<SingleTransport<Domain, F>>,
 }
-/// Whether partiall transported points are allowed.
+/// Whether partially 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;
+const NEW_APPROACH: bool = true;
-impl<const N: usize, F: Float> Transport<N, F> {
+pub trait TransportProxPenalty<Domain, PreadjointCodomain, Reg, F = f64>:
+ProxPenalty<Domain, PreadjointCodomain, Reg, F>
+where
+F: Float + ToNalgebraRealField,
+Reg: SlidingRegTerm<Domain, F>,
+Domain: 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(
+&self,
+γ: &mut Transport<Domain, F>,
+μ: &DiscreteMeasure<Domain, F>,
+ε: F,
+τ: F,
+τθ_or_adaptive: &mut Self::TransportStepLength,
+v: &PreadjointCodomain,
+tconfig: &TransportConfig<F>,
+);
+/// A posteriori transport adaptation.
+fn aposteriori_transport(
+&self,
+γ: &mut Transport<Domain, F>,
+μ: &DiscreteMeasure<Domain, F>,
+μ̆: &DiscreteMeasure<Domain, F>,
+τv̆: &mut PreadjointCodomain,
+v: &mut PreadjointCodomain,
+extra: Option<F>,
+ε: F,
+τ: F,
+τθ_or_adaptive: &Self::TransportStepLength,
+reg: &Reg,
+tconfig: &TransportConfig<F>,
+rconfig: &RefinementSettings<F>,
+attempts: &mut usize,
+) -> bool;
+fn get_transport_steplength(
+&self,
+lips: (DynResult<F>, DynResult<F>, DynResult<F>),
+tconfig: &TransportConfig<F>,
+ℓ_base: F,
+ℓ_base_max_transport: F,
+) -> Self::TransportStepLength;
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F, M, Reg, const N: usize> TransportProxPenalty<Loc<N, F>, M, Reg, F> for RadonSquared
+where
+RadonSquared: ProxPenalty<Loc<N, F>, M, Reg, F>,
+F: Float + ToNalgebraRealField,
+M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>,
+Reg: SlidingRegTerm<Loc<N, F>, F>,
+RNDM<N, F>: SpikeMerging<F>,
+{
+type TransportStepLength = TransportStepLength<F>;
+fn get_transport_steplength(
+&self,
+(ℓ_F, maybe_ℓ_gradv, _maybe_transport_lip): (DynResult<F>, DynResult<F>, DynResult<F>),
+tconfig: &TransportConfig<F>,
+ℓ: F,
+ℓ_base_max_transport: F,
+) -> Self::TransportStepLength {
+TransportStepLength::new(
+maybe_ℓ_gradv,
+tconfig,
+ℓ + ℓ_F.unwrap_or(0.0),
+ℓ_base_max_transport,
+)
+}
+fn initial_transport(
+&self,
+γ: &mut Transport<Loc<N, F>, F>,
+μ: &RNDM<N, F>,
+_ε: F,
+_τ: F,
+τθ_or_adaptive: &mut TransportStepLength<F>,
+v: &M,
+tconfig: &TransportConfig<F>,
+) {
+γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig);
+}
+fn aposteriori_transport(
+&self,
+γ: &mut Transport<Loc<N, F>, F>,
+μ: &RNDM<N, F>,
+μ̆: &RNDM<N, F>,
+τv̆: &mut M,
+v: &mut M,
+_extra: Option<F>,
+ε: F,
+τ: F,
+τθ_or_adaptive: &TransportStepLength<F>,
+reg: &Reg,
+tconfig: &TransportConfig<F>,
+_rconfig: &RefinementSettings<F>,
+attempts: &mut usize,
+) -> bool {
+*attempts += 1;
+if *attempts > tconfig.max_attempts {
+// Previous round has set transport to zero
+return true;
+}
+let nΔ = μ.dist_matching(&μ̆);
+let all_ok = γ.do_new_aposteriori_transport(
+μ,
+τv̆,
+v,
+ε,
+τ,
+τθ_or_adaptive,
+reg,
+tconfig,
+|_, mass_zero, negated| {
+use std::cmp::Ordering::*;
+match (mass_zero, negated) {
+(Less, _) => -nΔ,
+(Greater, _) => nΔ,
+(Equal, false) => nΔ, // pessimistic estimated for ω
+(Equal, true) => -nΔ, // pessimistic estimated for -ω
+}
+},
+);
+if !all_ok {
+if *attempts >= tconfig.max_attempts {
+for ρ in γ.iter_mut() {
+ρ.α_γ = 0.0;
+}
+}
+} else {
+for ρ in γ.iter_mut() {
+if ρ.α_γ == 0.0 {
+ρ.fail_count += 1;
+} else {
+ρ.fail_count = 0;
+}
+}
+}
+all_ok
+}
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<F, M, Reg, 𝒟, O, 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>,
+{
+type TransportStepLength = TransportStepLength<F>;
+fn get_transport_steplength(
+&self,
+(ℓ_F, maybe_ℓ_gradv, _maybe_transport_lip): (DynResult<F>, DynResult<F>, DynResult<F>),
+tconfig: &TransportConfig<F>,
+ℓ_base: F,
+ℓ_base_max_transport: F,
+) -> Self::TransportStepLength {
+TransportStepLength::new(
+maybe_ℓ_gradv,
+tconfig,
+ℓ_base + ℓ_F.unwrap_or(0.0),
+ℓ_base_max_transport,
+)
+}
+fn initial_transport(
+&self,
+γ: &mut Transport<Loc<N, F>, F>,
+μ: &RNDM<N, F>,
+_ε: F,
+_τ: F,
+τθ_or_adaptive: &mut TransportStepLength<F>,
+v: &M,
+tconfig: &TransportConfig<F>,
+) {
+γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig);
+}
+fn aposteriori_transport(
+&self,
+γ: &mut Transport<Loc<N, F>, F>,
+μ: &RNDM<N, F>,
+μ̆: &RNDM<N, F>,
+τv̆: &mut M,
+v: &mut M,
+extra: Option<F>,
+ε: F,
+τ: F,
+τθ_or_adaptive: &TransportStepLength<F>,
+reg: &Reg,
+tconfig: &TransportConfig<F>,
+_rconfig: &RefinementSettings<F>,
+attempts: &mut usize,
+) -> bool {
+*attempts += 1;
+if *attempts > tconfig.max_attempts {
+// Previous round has set transport to zero
+return true;
+}
+let all_ok = if NEW_APPROACH {
+let ω = self.apply(μ.sub_matching(&μ̆));
+γ.do_new_aposteriori_transport(
+μ,
+τv̆,
+v,
+ε,
+τ,
+τθ_or_adaptive,
+reg,
+tconfig,
+|x, _, _| ω.apply(x),
+)
+} else {
+let mut all_ok0 = true;
+// 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.
+for (δ, ρ) in izip!(μ.iter_spikes(), γ.iter_mut()) {
+if δ.α == 0.0 && ρ.α_γ != 0.0 {
+all_ok0 = false;
+ρ.α_γ = 0.0;
+}
+// TODO: sign
+}
+// 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γ = γ.norm(Radon);
+let nΔ = μ.dist_matching(&μ̆) + extra.unwrap_or(0.0);
+let t = ε * tconfig.tolerance_mult;
+if nγ * nΔ > t && *attempts >= tconfig.max_attempts {
+all_ok0 = 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.
+//*γ *= 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 tconfig.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 = γ.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 = γ
+//     .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 = γ.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 = γ
+.vec
+.iter()
+.map(|ρ| ρ.α_γ.abs())
+.filter(|t| *t > F::EPSILON)
+.collect::<Vec<F>>();
+abs_weights.sort_by(|a, b| a.total_cmp(b));
+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(_) => γ.vec.iter_mut().for_each(|δ| δ.α_γ = 0.0),
+ControlFlow::Break(chg_one) => γ.vec.iter_mut().for_each(|ρ| {
+let t = ρ.α_γ.abs();
+if t > 0.0 {
+if tconfig.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 = γ
+.vec
+.iter()
+.map(|ρ| ρ.α_γ.abs())
+.zip(0..)
+.filter(|(w, _)| *w >= 0.0)
+.collect::<Vec<(F, usize)>>();
+abs_weights_idx.sort_by(|(a, _), (b, _)| a.total_cmp(b));
+let mut left = reduction_target;
+for (w, i) in abs_weights_idx {
+left -= w;
+let ρ = &mut γ.vec[i];
+ρ.α_γ = 0.0;
+if left < 0.0 {
+break;
+}
+}
+}
+all_ok0 = false
+}
+all_ok0
+};
+if !all_ok {
+if *attempts >= tconfig.max_attempts {
+for ρ in γ.iter_mut() {
+ρ.α_γ = 0.0;
+}
+}
+} else {
+for ρ in γ.iter_mut() {
+if ρ.α_γ == 0.0 {
+ρ.fail_count += 1;
+} else {
+ρ.fail_count = 0;
+}
+}
+}
+all_ok
+}
+}
+#[replace_float_literals(F::cast_from(literal))]
+impl<const N: usize, F: Float> Transport<Loc<N, F>, F> {
 pub(crate) fn new() -> Self {
 Transport { vec: Vec::new() }
 }
-pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<N, F>> {
+pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<Loc<N, F>, F>> {
 self.vec.iter()
 }
-pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &'_ mut SingleTransport<N, F>> {
+pub(crate) fn iter_mut(
+&mut self,
+) -> impl Iterator<Item = &'_ mut SingleTransport<Loc<N, F>, F>> {
 self.vec.iter_mut()
 }
 pub(crate) fn extend<I>(&mut self, it: I)
 where
-I: IntoIterator<Item = SingleTransport<N, F>>,
+I: IntoIterator<Item = SingleTransport<Loc<N, F>, F>>,
 {
 self.vec.extend(it)
 }
 pub(crate) fn len(&self) -> usize {
 //         .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()) {
 }
 }
 // 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.
-// This can cause the number of points to explode, so cannot have partial transport.
+for ρ in self.iter() {
-for ρ in self.iter() {
+let α = ρ.α_μ_orig - ρ.α_γ;
-let α = ρ.α_μ_orig - ρ.α_γ;
+if ρ.α_γ.abs() > F::EPSILON && α.abs() > F::EPSILON {
-if ρ.α_γ.abs() > F::EPSILON && α != 0.0 {
+let δ = DeltaMeasure { α, x: ρ.x };
-let δ = DeltaMeasure { α, x: ρ.x };
+if i < μ.len() {
-if i < μ.len() {
+μ[i] = δ;
-μ[i] = δ;
+} else {
-} else {
+μ.push(δ)
-μ.push(δ)
+}
-}
+i += 1;
-i += 1;
-}
 }
 }
 μ.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,
-tconfig: &TransportConfig<F>,
-) 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 ρ.fail_count > tconfig.max_fail {
-ρ.α_γ = 0.0
-} else {
-// 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,
-fail_count: 0,
-}));
-// Calculate transport rays.
-match *τθ_or_adaptive {
-Fixed(θ) => {
-for ρ in self.iter_mut() {
-if ρ.fail_count <= tconfig.max_fail {
-ρ.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() {
-if ρ.fail_count <= tconfig.max_fail {
-ρ.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() {
-if ρ.fail_count < tconfig.max_fail {
-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;
-}
-}
-}
-}
-}
-/// 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;
-}
-}
-for ρ in self.iter_mut() {
-if ρ.α_γ == 0.0 {
-ρ.fail_count += 1;
-} else if all_ok {
-ρ.fail_count = 0;
-}
-}
-all_ok
 }
 /// 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()
 .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();
 /// 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);
+ρ.retain = δ.α.abs() > F::EPSILON;
 }
 μ.prune_by(|δ| δ.α.abs() > F::EPSILON);
 stats.pruned += old_len - μ.len();
-self.vec.retain(|ρ| !ρ.prune);
+self.vec.retain(|ρ| ρ.retain);
 assert!(self.vec.len() <= μ.len());
 }
-}
+/// Helper for initial transport. Called from [`TransportProxPenalty::initial_transport`].
-impl<const N: usize, F: Float> Norm<Radon, F> for Transport<N, F> {
+pub(super) fn do_init_transport<D>(
+&mut self,
+v: &D,
+μ: &RNDM<N, F>,
+τθ_or_adaptive: &mut TransportStepLength<F>,
+tconfig: &TransportConfig<F>,
+) where
+D: DifferentiableRealMapping<N, F>,
+{
+use TransportStepLength::*;
+// Initialise transport structure weights
+for (δ, ρ) in izip!(μ.iter_spikes(), self.iter_mut()) {
+ρ.α_μ_orig = δ.α;
+ρ.x = δ.x;
+ρ.y = δ.x; // Later updated if no fails.
+ρ.α_γ = if ρ.fail_count > tconfig.max_fail {
+0.0
+} else {
+// 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: δ.α,
+α_γ: δ.α,
+retain: true,
+fail_count: 0,
+excess: 0.0,
+}));
+// Calculate transport rays.
+let simple_τθ = match *τθ_or_adaptive {
+Fixed { τθ, .. } => Some(τθ),
+Simple { ℓ_gradv, τθ0, ℓ_base } => Some(τθ0 / (ℓ_gradv + ℓ_base)),
+AdaptiveMax {
+ℓ_gradv,
+ref mut adaptive_max_transport,
+τθ0,
+ℓ_base,
+ℓ_base_max_transport,
+} => {
+*adaptive_max_transport = adaptive_max_transport.max(self.norm(Radon));
+Some(τθ0 / (ℓ_gradv + ℓ_base + ℓ_base_max_transport * *adaptive_max_transport))
+}
+FullyAdaptive {
+ref mut adaptive_ℓ_gradv,
+ref mut adaptive_max_transport,
+τθ0,
+ℓ_base,
+ℓ_base_max_transport,
+} => {
+*adaptive_max_transport = adaptive_max_transport.max(self.norm(Radon));
+let mut τθ = τθ0
+/ (*adaptive_ℓ_gradv + ℓ_base + ℓ_base_max_transport * *adaptive_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() {
+if ρ.fail_count < tconfig.max_fail {
+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_ℓ_gradv = (dv_x - v.differential(&ρ.y)).norm2() / n;
+*adaptive_ℓ_gradv = adaptive_ℓ_gradv.max(this_ℓ_gradv);
+τθ = τθ0
+/ (*adaptive_ℓ_gradv
++ ℓ_base
++ ℓ_base_max_transport * *adaptive_max_transport);
+changes = true
+}
+}
+}
+if !changes {
+break;
+}
+}
+None
+}
+};
+if let Some(τθ) = simple_τθ {
+for ρ in self.iter_mut() {
+if ρ.fail_count <= tconfig.max_fail {
+ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * τθ);
+}
+}
+}
+}
+/// Helper for a posteriori transport error control.
+/// Called from [`TransportProxPenalty::aposteriori_transport`].
+fn do_new_aposteriori_transport<Reg, M>(
+&mut self,
+μ: &RNDM<N, F>,
+τv̆: &mut M,
+v: &mut M,
+ε: F,
+τ: F,
+τθ_or_adaptive: &TransportStepLength<F>,
+reg: &Reg,
+tconfig: &TransportConfig<F>,
+ω: impl Fn(&Loc<N, F>, std::cmp::Ordering, bool) -> F,
+) -> bool
+where
+Reg: SlidingRegTerm<Loc<N, F>, F>,
+F: ToNalgebraRealField,
+M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>,
+{
+let τℓ_gradv2 = τ * τθ_or_adaptive.get_ℓ_gradv() / 2.0;
+let Bounds(α_lower, α_upper) = reg.subdiff_range();
+let (all_ok, m) = izip!(self.vec.iter_mut(), μ.iter_spikes()).fold(
+(true, 0.0),
+|(all_ok_so_far, total_excess), (ρ, δ)| {
+use std::cmp::Ordering::*;
+// NOTE: The tolerances ε are commented out, because they will in any case
+// be consumed by `t` below by suitably large choice of `tolerance_mult`.
+// Hence, we simply implicitly adapt the `tolerance_mult` be commenting out
+// the `ε` here. That way, `d` is in its entirely multiplied by τ, as is `t`,
+// making all factors independent of τ.
+let maybe_excess = match ρ.α_γ.total_cmp(&0.0) {
+Greater => (δ.α >= 0.0).then(|| {
+let gvx = v.differential(&ρ.x);
+let d = if tconfig.alt_remainder_control {
+let τv̆y = τv̆.apply(&ρ.y);
+(/*ε +*/τ * α_upper + ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), false))
++ (τv̆y + τ * gvx.dot(&ρ.x - &ρ.y))
+- τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x)
+} else {
+(/*2.0 * ε*/-ω(&ρ.y, δ.α.total_cmp(&ρ.α_γ), true)
++ ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), false))
++ τ * gvx.dot(&ρ.x - &ρ.y)
+- τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x)
+};
+d * ρ.α_γ
+}),
+Less => (δ.α <= 0.0).then(|| {
+let gvx = v.differential(&ρ.x);
+let d = if tconfig.alt_remainder_control {
+let τv̆y = τv̆.apply(&ρ.y);
+(/*ε*/-τ * α_lower - ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), true))
+- (τv̆y + τ * gvx.dot(&ρ.x - &ρ.y))
+- τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x)
+} else {
+(/*2.0 * ε +*/ω(&ρ.y, δ.α.total_cmp(&ρ.α_γ), false)
+- ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), true))
+- τ * gvx.dot(&ρ.x - &ρ.y)
+- τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x)
+};
+d * (-ρ.α_γ)
+}),
+Equal => Some(0.0),
+};
+match maybe_excess {
+None => {
+ρ.α_γ = 0.0;
+ρ.excess = 0.0;
+(false, total_excess)
+}
+Some(e) => {
+ρ.excess = e;
+(all_ok_so_far, total_excess + e)
+}
+}
+},
+);
+let t = τ * ε * tconfig.tolerance_mult;
+if m > t {
+let mut it = self.vec.iter_mut().rev().filter(|ρ| ρ.excess > 0.0);
+let reduction_target = m - tconfig.adaptation * t;
+let _unused = it.try_fold(reduction_target, |left, ρ| {
+let d = ρ.excess;
+if d >= left {
+if tconfig.allow_partial_transport {
+ρ.α_γ *= (d - left) / d;
+} else {
+ρ.α_γ = 0.0;
+}
+ControlFlow::Break(())
+} else {
+ρ.α_γ = 0.0;
+ControlFlow::Continue(left - d)
+}
+});
+false
+} else {
+all_ok
+}
+}
+}
+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<N, F> {
+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;
 }
 }
 ) -> DynResult<RNDM<N, F>>
 where
 F: Float + ToNalgebraRealField,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
-Dat::DerivativeDomain: DifferentiableRealMapping<N, F> + ClosedMul<F>,
+Dat::DerivativeDomain:
+DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
 //for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
 RNDM<N, F>: SpikeMerging<F>,
 Reg: SlidingRegTerm<Loc<N, F>, F>,
-P: ProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
 Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
 {
 // Check parameters
 ensure!(config.τ0 > 0.0, "Invalid step length parameter");
 config.transport.check()?;
 // Set up parameters
 // let opAnorm = opA.opnorm_bound(Radon, L2);
 //let max_transport = config.max_transport.scale
 //                    * reg.radon_norm_bound(b.norm2_squared() / 2.0);
 //let ℓ = opA.transport.lipschitz_factor(L2Squared) * max_transport;
-let ℓ = 0.0;
 let τ = config.τ0 / prox_penalty.step_length_bound(&f)?;
-let mut θ_or_adaptive = match f.curvature_bound_components(config.guess) {
+let mut τθ_or_adaptive = prox_penalty.get_transport_steplength(
-(_, Err(_)) => TransportStepLength::Fixed(config.transport.θ0),
+f.curvature_bound_components(config.guess),
-(maybe_ℓ_F, Ok(transport_lip)) => {
+&config.transport,
-let calculate_θτ = move |ℓ_F, max_transport| {
+0.0,
-let ℓ_r = transport_lip * max_transport;
+0.0,
-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.
 let tolerance = config.insertion.tolerance * τ * reg.tolerance_scaling();
 let mut ε = tolerance.initial();
 };
 let mut stats = IterInfo::new();
 // Run the algorithm
 for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
+let mut v = f.differential(&μ);
 // Calculate initial transport
-let v = f.differential(&μ);
+prox_penalty.initial_transport(
-γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v, &config.transport);
+&mut γ,
+&μ,
+ε,
+τ,
+&mut τθ_or_adaptive,
+&v,
+&config.transport,
+);
 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.
 &state,
 &mut stats,
 )?;
 // A posteriori transport adaptation.
-if γ.aposteriori_transport(&μ, &μ̆, &mut τv̆, None, ε, &config.transport, &mut attempts)
+if prox_penalty.aposteriori_transport(
-{
+&mut γ,
+&μ,
+&μ̆,
+&mut τv̆,
+&mut v,
+None,
+ε,
+τ,
+&τθ_or_adaptive,
+reg,
+&config.transport,
+&config.insertion.refinement,
+&mut attempts,
+) {
 break 'adapt_transport (maybe_d, within_tolerances, τv̆, μ̆);
 }
 stats.get_transport_mut().readjustment_iters += 1;
 };
 // 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 γ.
 if config.insertion.merge_now(&state) {
-stats.merged += prox_penalty.merge_spikes(
+let m = prox_penalty.merge_spikes(
 &mut μ,
 &mut τv̆,
 &μ̆,
 τ,
 ε,
 &config.insertion,
 &reg,
 Some(|μ̃: &RNDM<N, F>| f.apply(μ̃)),
 );
+if m > 0 {
+stats.merged += m;
+v = f.differential(&μ);
+}
 }
 γ.prune_compat(&mut μ, &mut stats);
 let iter = state.iteration();

comparison: src/sliding_fb.rs

src/sliding_fb.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_fb.rs

src/sliding_fb.rs