--- a/src/sliding_fb.rs Thu Jan 23 23:48:52 2025 +0100
+++ b/src/sliding_fb.rs Sun Jan 26 11:29:57 2025 -0500
@@ -15,7 +15,6 @@
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};
@@ -49,8 +48,6 @@
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,
}
@@ -61,7 +58,6 @@
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);
}
}
@@ -73,7 +69,6 @@
θ0 : 0.4,
adaptation : 0.9,
tolerance_mult_pos : 100.0,
- tolerance_mult_pri : 1000.0,
}
}
}
@@ -113,25 +108,19 @@
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'(μ)`
+/// with step lengh τ and transport step length `θ_or_adaptive`.
#[replace_float_literals(F::cast_from(literal))]
-pub(crate) fn initial_transport<F, G, D, Observable, const N : usize>(
+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>,
{
@@ -161,9 +150,7 @@
};
};
- // A priori transport adaptation based on bounding 2 ‖A‖ ‖A(γ₁-γ₀)‖‖γ‖ by scaling γ.
- // 1. Calculate transport rays.
- // If the Lipschitz factor of the values v=∇F(μ) are not known, estimate it.
+ // Calculate transport rays.
match *θ_or_adaptive {
Fixed(θ) => {
let θτ = τ * θ;
@@ -205,73 +192,6 @@
}
}
- // 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 {
@@ -375,7 +295,7 @@
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;
@@ -415,9 +335,7 @@
// Calculate initial transport
let v = opA.preadjoint().apply(residual);
let (μ_base_masses, mut μ_base_minus_γ0) = initial_transport(
- &mut γ1, &mut μ, |ν| opA.apply(ν),
- ε, τ, &mut θ_or_adaptive, opAnorm,
- v, &config.transport,
+ &mut γ1, &mut μ, τ, &mut θ_or_adaptive, v
);
// Solve finite-dimensional subproblem several times until the dual variable for the