pointsource_algs: comparison src/sliding

-:759634edb536
+:8f1e57821d76
 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::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;
 /// Transport settings for [`pointsource_sliding_fb_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 }
 }
 }
 #[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
 α_γ: F,
 /// Helper for pruning
 /// Contribution to remainder (temporary variable)
 excess: F,
 }
 #[derive(Clone, Debug, Serialize)]
-pub struct Transport<Domain, F: Float> {
+pub struct Transport<U, F: Float> {
-vec: Vec<SingleTransport<Domain, F>>,
+vec: Vec<SingleTransport<U, F>>,
 }
 /// 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 MINIMAL_PARTIAL_TRANSPORT: bool = true;
 const NEW_APPROACH: bool = true;
-pub trait TransportProxPenalty<Point, PreadjointCodomain, Reg, F = f64>:
+pub trait TransportProxPenalty<U, PreadjointCodomain, Reg, F = f64>:
-ProxPenalty<Point, PreadjointCodomain, Reg, F>
+ProxPenalty<U, PreadjointCodomain, Reg, F>
 where
 F: Float + ToNalgebraRealField,
-Reg: SlidingRegTerm<Point, F>,
+Reg: SlidingRegTerm<U, F>,
-Point: Space + Clone,
+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>,
+γ: &mut Transport<U, F>,
-μ: &DiscreteMeasure<Point, F>,
+μ: &DiscreteMeasure<U, F>,
 ε: F,
 τ: F,
 τθ_or_adaptive: &mut Self::TransportStepLength,
 v: &PreadjointCodomain,
 tconfig: &TransportConfig<F>,
 );
 /// A posteriori transport adaptation.
 fn aposteriori_transport(
 &self,
-γ: &mut Transport<Point, F>,
+γ: &mut Transport<U, F>,
-μ: &DiscreteMeasure<Point, F>,
+μ: &DiscreteMeasure<U, F>,
-μ̆: &DiscreteMeasure<Point, F>,
+μ̆: &DiscreteMeasure<U, F>,
 τv̆: &mut PreadjointCodomain,
 v: &mut PreadjointCodomain,
 extra: Option<F>,
 ε: F,
 τ: F,
 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 𝒟
+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>,
 {
 type TransportStepLength = TransportStepLength<F>;
 }
 all_ok
 }
 }
+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(
+pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &'_ mut SingleTransport<U, F>> {
-&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<Loc<N, F>, F>>,
+I: IntoIterator<Item = SingleTransport<U, F>>,
 {
 self.vec.extend(it)
 }
 pub(crate) fn len(&self) -> usize {
 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.
 stats.inserted += μ.len() - self.len();
 (a, b)
 }
 })
 };
 }
+}
-/// Prune spikes with zero weight. To maintain correct ordering between μ and γ, also the
-/// latter needs to be pruned when μ is.
+#[replace_float_literals(F::cast_from(literal))]
-pub(crate) fn prune_compat(&mut self, μ: &mut RNDM<N, F>, stats: &mut IterInfo<F>) {
+impl<U: Clone, F: Float> Transport<U, F> {
-assert!(self.vec.len() <= μ.len());
+// pub(crate) fn dist_matching(&self, μ: &RNDM<N, F>) -> F {
-let old_len = μ.len();
+//     self.iter()
-for (ρ, δ) in self.vec.iter_mut().zip(μ.iter_spikes()) {
+//         .zip(μ.iter_spikes())
-ρ.retain = δ.α != F::ZERO;
+//         .map(|(ρ, δ)| (ρ.α_γ - δ.α).abs())
-}
+//         .sum()
-// Do not prune new insertions on first round.
+// }
-μ.prune_until(self.vec.len());
-stats.pruned += old_len - μ.len();
+/// Construct `μ̆`, replacing the contents of `μ`.
-self.vec.retain(|ρ| ρ.retain);
+pub(crate) fn μ̆_into(&self, μ: &mut DiscreteMeasure<U, F>) {
-assert!(self.vec.len() <= μ.len());
+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();
+}
+}
+// 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,
 v: &D,
 μ: &RNDM<N, F>,
 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<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
 ///
 /// The parametrisation is as for [`pointsource_fb_reg`], except an additional domain restriction
 /// 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>>,
+μ0: Option<DiscreteMeasure<U, F>>,
-domain: &mut impl AdmissibleDomain<Loc<N, F>>,
+domain: &mut impl AdmissibleDomain<U>,
-) -> DynResult<RNDM<N, 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>,
+DiscreteMeasure<U, F>: SpikeMerging<F>,
-Reg: SlidingRegTerm<Loc<N, F>, F>,
+Reg: SlidingRegTerm<U, F>,
-P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
-Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
+Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
 {
 // Check parameters
 ensure!(config.τ0 > 0.0, "Invalid step length parameter");
 config.transport.check()?;
 // by τ compared to the conditional gradient approach.
 let tolerance = config.insertion.tolerance * τ * reg.tolerance_scaling();
 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(),
 ε,
 // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
 ..stats
 &μ̆,
 τ,
 ε,
 ins,
 &reg,
-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);
 if config.insertion.random_insertions {
 postprocess(μ, &config.insertion, |μ̃| f.apply(μ̃) + reg.apply(μ̃))
 }
 #[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>>,
+μ0: Option<DiscreteMeasure<U, F>>,
-) -> DynResult<RNDM<N, 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>,
+DiscreteMeasure<U, F>: SpikeMerging<F>,
-Reg: SlidingRegTerm<Loc<N, F>, F>,
+Reg: SlidingRegTerm<U, F>,
-P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
-Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
+Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
 {
 pointsource_sliding_fb(
 f,
 reg,
 prox_penalty,

comparison: src/sliding_fb.rs

src/sliding_fb.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_fb.rs

src/sliding_fb.rs