pointsource_algs: comparison src/sliding

-:95bb12bdb6ac
+:00d0881f89a6
 pub adaptation: F,
 /// A posteriori transport tolerance multiplier (C_pos)
 pub tolerance_mult_con: 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,
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F: Float> TransportConfig<F> {
 /// Check that the parameters are ok. Panics if not.
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F: Float> Default for TransportConfig<F> {
 fn default() -> Self {
-TransportConfig { θ0: 0.9, adaptation: 0.9, tolerance_mult_con: 100.0, max_attempts: 2 }
+TransportConfig {
+θ0: 0.9,
+adaptation: 0.9,
+tolerance_mult_con: 100.0,
+max_attempts: 2,
+max_fail: usize::MAX,
+}
 }
 }
 /// Settings for [`pointsource_sliding_fb_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 α_μ_orig: F,
 /// Transported mass
 α_γ: F,
 /// Helper for pruning
 prune: bool,
+/// Fail count
+fail_count: usize,
 }
 #[derive(Clone, Debug, Serialize)]
 pub struct Transport<const N: usize, F: Float> {
 vec: Vec<SingleTransport<N, F>>,
 &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 old transport has opposing sign, the new transport will be none.
+if ρ.fail_count > tconfig.max_fail {
-ρ.α_γ = if (ρ.α_γ > 0.0 && δ.α < 0.0) || (ρ.α_γ < 0.0 && δ.α > 0.0) {
+ρ.α_γ = 0.0
-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);
 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() {
-ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * θ);
+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() {
-ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * θτ);
+if ρ.fail_count <= tconfig.max_fail {
+ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * θτ);
+}
 }
 }
 FullyAdaptive {
 l: ref mut adaptive_ℓ_F,
 ref mut 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() {
-let dv_x = v.differential(&ρ.x);
+if ρ.fail_count < tconfig.max_fail {
-let g = &dv_x * (ρ.α_γ.signum() * θτ);
+let dv_x = v.differential(&ρ.x);
-ρ.y = ρ.x - g;
+let g = &dv_x * (ρ.α_γ.signum() * θτ);
-let n = g.norm2();
+ρ.y = ρ.x - g;
-if n >= F::EPSILON {
+let n = g.norm2();
-// Estimate Lipschitz factor of ∇v
+if n >= F::EPSILON {
-let this_ℓ_F = (dv_x - v.differential(&ρ.y)).norm2() / n;
+// Estimate Lipschitz factor of ∇v
-*adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F);
+let this_ℓ_F = (dv_x - v.differential(&ρ.y)).norm2() / n;
-θτ = calculate_θτ(*adaptive_ℓ_F, *max_transport);
+*adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F);
-changes = true
+θτ = calculate_θτ(*adaptive_ℓ_F, *max_transport);
+changes = true
+}
 }
 }
 if !changes {
 break;
 }
 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 {
 // Run the algorithm
 for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
 // Calculate initial transport
 let v = f.differential(&μ);
-γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v);
+γ.initial_transport(&μ, τ, &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.

comparison: src/sliding_fb.rs

src/sliding_fb.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_fb.rs

src/sliding_fb.rs