# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1737908997 18000
# Node ID b6bdb6cb4d442a2012b4535ef052360a9229e293
# Parent  896b42b5ac1a94b6e51912d1602db9a490b28cf6
Remove initial transport adaptation—it is not needed, after all

diff -r 896b42b5ac1a -r b6bdb6cb4d44 src/experiments.rs
--- a/src/experiments.rs	Thu Jan 23 23:48:52 2025 +0100
+++ b/src/experiments.rs	Sun Jan 26 11:29:57 2025 -0500
@@ -182,7 +182,6 @@
         );
         let much_higher_cpos_merging_steptune = |alg| (alg,
             AlgorithmOverrides {
-                transport_tolerance_pri : Some(1000.0),
                 transport_tolerance_pos : Some(10000.0),
                 sigma0 : Some(0.15),
                 theta0 : Some(0.3),
diff -r 896b42b5ac1a -r b6bdb6cb4d44 src/main.rs
--- a/src/main.rs	Thu Jan 23 23:48:52 2025 +0100
+++ b/src/main.rs	Sun Jan 26 11:29:57 2025 -0500
@@ -194,10 +194,6 @@
     theta0 : Option<F>,
 
     #[arg(long)]
-    /// A priori transport tolerance multiplier (C_pri)
-    transport_tolerance_pri : Option<F>,
-
-    #[arg(long)]
     /// A posteriori transport tolerance multiplier (C_pos)
     transport_tolerance_pos : Option<F>,
 
diff -r 896b42b5ac1a -r b6bdb6cb4d44 src/run.rs
--- a/src/run.rs	Thu Jan 23 23:48:52 2025 +0100
+++ b/src/run.rs	Sun Jan 26 11:29:57 2025 -0500
@@ -171,7 +171,6 @@
             TransportConfig {
                 θ0 : cli.theta0.unwrap_or(g.θ0),
                 tolerance_mult_pos: cli.transport_tolerance_pos.unwrap_or(g.tolerance_mult_pos),
-                tolerance_mult_pri: cli.transport_tolerance_pri.unwrap_or(g.tolerance_mult_pri),
                 adaptation: cli.transport_adaptation.unwrap_or(g.adaptation),
                 .. g
             }
diff -r 896b42b5ac1a -r b6bdb6cb4d44 src/sliding_fb.rs
--- 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
diff -r 896b42b5ac1a -r b6bdb6cb4d44 src/sliding_pdps.rs
--- a/src/sliding_pdps.rs	Thu Jan 23 23:48:52 2025 +0100
+++ b/src/sliding_pdps.rs	Sun Jan 26 11:29:57 2025 -0500
@@ -162,7 +162,7 @@
 
     // Set up parameters
     // TODO: maybe this PairNorm doesn't make sense here?
-    let opAnorm = opA.opnorm_bound(PairNorm(Radon, L2, L2), L2);
+    // let opAnorm = opA.opnorm_bound(PairNorm(Radon, L2, L2), L2);
     let bigθ = 0.0; //opKμ.transport_lipschitz_factor(L2Squared);
     let bigM = 0.0; //opKμ.adjoint_product_bound(&op𝒟).unwrap().sqrt();
     let nKz = opKz.opnorm_bound(L2, L2);
@@ -245,9 +245,7 @@
         // TODO: Write a version of initial_transport that can deal with K_μ ≠ 0.
  
         let (μ_base_masses, mut μ_base_minus_γ0) = initial_transport(
-            &mut γ1, &mut μ, |ν| opA.apply(Pair(ν, &z)),
-            ε, τ, &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