pointsource_algs: comparison src/sliding

-:896b42b5ac1a
+:b6bdb6cb4d44
 use alg_tools::iterate::AlgIteratorFactory;
 use alg_tools::euclidean::Euclidean;
 use alg_tools::mapping::{Mapping, DifferentiableRealMapping, Instance};
 use alg_tools::norms::Norm;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::L2;
 use crate::types::*;
 use crate::measures::{DiscreteMeasure, Radon, RNDM};
 use crate::measures::merging::SpikeMerging;
 use crate::forward_model::{
 pub struct TransportConfig<F : Float> {
 /// Transport step length $θ$ normalised to $(0, 1)$.
 pub θ0 : F,
 /// Factor in $(0, 1)$ for decreasing transport to adapt to tolerance.
 pub adaptation : F,
-/// A priori transport tolerance multiplier (C_pri)
-pub tolerance_mult_pri : F,
 /// A posteriori transport tolerance multiplier (C_pos)
 pub tolerance_mult_pos : F,
 }
 #[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) {
 assert!(self.θ0 > 0.0);
 assert!(0.0 < self.adaptation && self.adaptation < 1.0);
-assert!(self.tolerance_mult_pri > 0.0);
 assert!(self.tolerance_mult_pos > 0.0);
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
 fn default() -> Self {
 TransportConfig {
 θ0 : 0.4,
 adaptation : 0.9,
 tolerance_mult_pos : 100.0,
-tolerance_mult_pri : 1000.0,
 }
 }
 }
 /// Settings for [`pointsource_sliding_fb_reg`].
 /// Adaptive step length.
 /// Content of `l` depends on use case, while `g` calculates the step length from `l`.
 FullyAdaptive{ l : F, max_transport : F, g : G },
 }
-/// Constrution and a priori transport adaptation.
+/// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
-#[replace_float_literals(F::cast_from(literal))]
+/// with step lengh τ and transport step length `θ_or_adaptive`.
-pub(crate) fn initial_transport<F, G, D, Observable, const N : usize>(
+#[replace_float_literals(F::cast_from(literal))]
+pub(crate) fn initial_transport<F, G, D, const N : usize>(
 γ1 : &mut RNDM<F, N>,
 μ : &mut RNDM<F, N>,
-opAapply : impl Fn(&RNDM<F, N>) -> Observable,
-ε : F,
 τ : F,
 θ_or_adaptive : &mut TransportStepLength<F, G>,
-opAnorm : F,
 v : D,
-tconfig : &TransportConfig<F>
 ) -> (Vec<F>, RNDM<F, N>)
 where
 F : Float + ToNalgebraRealField,
 G : Fn(F, F) -> F,
-Observable : Euclidean<F, Output=Observable>,
-for<'a> &'a Observable : Instance<Observable>,
-//for<'b> A::Preadjoint<'b> : LipschitzValues<FloatType=F>,
 D : DifferentiableRealMapping<F, N>,
 {
 use TransportStepLength::*;
 } else {
 δ.α
 };
 };
-// A priori transport adaptation based on bounding 2 ‖A‖ ‖A(γ₁-γ₀)‖‖γ‖ by scaling γ.
+// Calculate transport rays.
-// 1. Calculate transport rays.
-//    If the Lipschitz factor of the values v=∇F(μ) are not known, estimate it.
 match *θ_or_adaptive {
 Fixed(θ) => {
 let θτ = τ * θ;
 for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) {
 ρ.x = δ.x - v.differential(&δ.x) * (ρ.α.signum() * θτ);
 }
 }
 }
 }
-// 2. Adjust transport mass, if needed.
-// This tries to remove the smallest transport masses first.
-if true {
-// Alternative 1 : subtract same amount from all transport rays until reaching zero
-loop {
-let nr =γ1.norm(Radon);
-let n = τ * 2.0 * opAnorm * (opAapply(&*γ1)-opAapply(&*μ)).norm2();
-if n <= 0.0 || nr <= 0.0 {
-break
-}
-let reduction_needed = nr - (ε * tconfig.tolerance_mult_pri / n);
-if reduction_needed <= 0.0 {
-break
-}
-let (min_nonzero, n_nonzero) = γ1.iter_masses()
-.map(|α| α.abs())
-.filter(|α| *α > F::EPSILON)
-.fold((F::INFINITY, 0), |(a, n), b| (a.min(b), n+1));
-assert!(n_nonzero > 0);
-// Reduction that can be done in all nonzero spikes simultaneously
-let h = (reduction_needed / F::cast_from(n_nonzero)).min(min_nonzero);
-for (δ, ρ) in izip!(μ.iter_spikes_mut(), γ1.iter_spikes_mut()) {
-ρ.α = ρ.α.signum() * (ρ.α.abs() - h).max(0.0);
-δ.α = ρ.α;
-}
-if min_nonzero * F::cast_from(n_nonzero) >= reduction_needed {
-break
-}
-}
-} else {
-// Alternative 2: first reduce transport rays with greater effect based on differential.
-// This is a an inefficient quick-and-dirty implementation.
-loop {
-let nr = γ1.norm(Radon);
-let a = opAapply(&*γ1)-opAapply(&*μ);
-let na = a.norm2();
-let n = τ * 2.0 * opAnorm * na;
-if n <= 0.0 || nr <= 0.0 {
-break
-}
-let reduction_needed = nr - (ε * tconfig.tolerance_mult_pri / n);
-if reduction_needed <= 0.0 {
-break
-}
-let mut max_d = 0.0;
-let mut max_d_ind = 0;
-for (δ, ρ, i) in izip!(μ.iter_spikes_mut(), γ1.iter_spikes(), 0..) {
-// Calculate differential of  ‖A(γ₁-γ₀)‖‖γ‖  wrt. each spike
-let s = δ.α.signum();
-// TODO: this is very inefficient implementation due to the limitations
-// of the closure parameters.
-let δ1 = DiscreteMeasure::from([(ρ.x, s)]);
-let δ2 = DiscreteMeasure::from([(δ.x, s)]);
-let a_part = opAapply(&δ1)-opAapply(&δ2);
-let d = a.dot(&a_part)/na * nr + 2.0 * na;
-if d > max_d {
-max_d = d;
-max_d_ind = i;
-}
-}
-// Just set mass to zero for transport ray with greater differential
-assert!(max_d > 0.0);
-γ1[max_d_ind].α = 0.0;
-μ[max_d_ind].α = 0.0;
-}
-}
 // 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
 let mut μ = DiscreteMeasure::new();
 let mut γ1 = DiscreteMeasure::new();
 let mut residual = -b; // Has to equal $Aμ-b$.
 // Set up parameters
-let opAnorm = opA.opnorm_bound(Radon, L2);
+// 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 / opA.adjoint_product_bound(prox_penalty).unwrap();
 // Run the algorithm
 for state in iterator.iter_init(|| full_stats(&residual, &μ, ε, stats.clone())) {
 // Calculate initial transport
 let v = opA.preadjoint().apply(residual);
 let (μ_base_masses, mut μ_base_minus_γ0) = initial_transport(
-&mut γ1, &mut μ, |ν| opA.apply(ν),
+&mut γ1, &mut μ, τ, &mut θ_or_adaptive, v
-ε, τ, &mut θ_or_adaptive, opAnorm,
-v, &config.transport,
 );
 // 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 {

comparison: src/sliding_fb.rs

src/sliding_fb.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_fb.rs

src/sliding_fb.rs