Mercurial > repos > pointsource_algs / changeset
changeset
Subproblem solver and sliding adjustments/improvements
Thu, 26 Feb 2026 11:36:22 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Thu, 26 Feb 2026 11:36:22 -0500
- branch
- dev
- changeset 63
- 7a8a55fd41c0
- parent 61
- 4f468d35fa29
- child 64
- d3be4f7ffd60
- child 65
- ec5f160c413b
Subproblem solver and sliding adjustments/improvements
--- a/src/fb.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/fb.rs Thu Feb 26 11:36:22 2026 -0500 @@ -107,11 +107,7 @@ #[replace_float_literals(F::cast_from(literal))] impl<F: Float> Default for FBConfig<F> { fn default() -> Self { - FBConfig { - τ0: 0.99, - σp0: 0.99, - insertion: Default::default(), - } + FBConfig { τ0: 0.99, σp0: 0.99, insertion: Default::default() } } } @@ -157,7 +153,7 @@ fbconfig: &FBConfig<F>, iterator: I, mut plotter: Plot, - μ0 : Option<RNDM<N, F>>, + μ0: Option<RNDM<N, F>>, ) -> DynResult<RNDM<N, F>> where F: Float + ToNalgebraRealField, @@ -196,21 +192,22 @@ // TODO: optimise τ to be applied to residual. let mut τv = f.differential(&μ) * τ; - // Save current base point - let μ_base = μ.clone(); + // Save current base point for merge + let μ_base_len = μ.len(); + let maybe_μ_base = config.merge_now(&state).then(|| μ.clone()); // Insert and reweigh - let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( - &mut μ, &mut τv, &μ_base, None, τ, ε, config, ®, &state, &mut stats, - )?; + let (maybe_d, _within_tolerances) = prox_penalty + .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; + + stats.inserted += μ.len() - μ_base_len; // Prune and possibly merge spikes - if config.merge_now(&state) { + if let Some(μ_base) = maybe_μ_base { stats.merged += prox_penalty.merge_spikes( &mut μ, &mut τv, &μ_base, - None, τ, ε, config, @@ -257,7 +254,7 @@ fbconfig: &FBConfig<F>, iterator: I, mut plotter: Plot, - μ0: Option<RNDM<N, F>> + μ0: Option<RNDM<N, F>>, ) -> DynResult<RNDM<N, F>> where F: Float + ToNalgebraRealField, @@ -298,13 +295,13 @@ // Calculate smooth part of surrogate model. let mut τv = f.differential(&μ) * τ; - // Save current base point - let μ_base = μ.clone(); + let μ_base_len = μ.len(); // Insert new spikes and reweigh - let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( - &mut μ, &mut τv, &μ_base, None, τ, ε, config, ®, &state, &mut stats, - )?; + let (maybe_d, _within_tolerances) = prox_penalty + .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; + + stats.inserted += μ.len() - μ_base_len; // (Do not) merge spikes. if config.merge_now(&state) && !warned_merging {
--- a/src/forward_model.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/forward_model.rs Thu Feb 26 11:36:22 2026 -0500 @@ -48,6 +48,7 @@ /// Based on Lemma 5.11 and Example 5.12, the helper bound for (5.15a) and (5.16a) is (3.8). /// Thus, subject to `guess` being correct, returns factor $(ℓ_F, Θ²)$ such that /// $B_{F'(μ)} dγ ≤ ℓ_F c_2$ and $⟨F'(μ+Δ)-F'(μ)|Δ⟩ ≤ Θ²|γ|(c_2)$, where $Δ=(π_♯^1-π_♯^0)γ$. +/// The latter is the firm transport Lipshitz property of (3.8b) and Lemma 5.11. /// /// This trait is supposed to be implemented by the data term $F$, in the basic case a /// [`Mapping`] from [`RNDM<N, F>`] to a [`Float`] `F`.
--- a/src/forward_pdps.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/forward_pdps.rs Thu Feb 26 11:36:22 2026 -0500 @@ -136,15 +136,15 @@ Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>, { // Check parameters - // ensure!( - // config.τ0 > 0.0 - // && config.τ0 < 1.0 - // && config.σp0 > 0.0 - // && config.σp0 < 1.0 - // && config.σd0 >= 0.0 - // && config.σp0 * config.σd0 <= 1.0, - // "Invalid step length parameters" - // ); + ensure!( + config.τ0 > 0.0 + && config.τ0 < 1.0 + && config.σp0 > 0.0 + && config.σp0 < 1.0 + && config.σd0 >= 0.0 + && config.σp0 * config.σd0 <= 1.0, + "Invalid step length parameters" + ); // Initialise iterates let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new()); @@ -206,8 +206,6 @@ let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( &mut μ, &mut τv, - &μ_base, - None, τ, ε, &config.insertion, @@ -216,6 +214,8 @@ &mut stats, )?; + stats.inserted += μ.len() - μ_base.len(); + // Merge spikes. // This crucially expects the merge routine to be stable with respect to spike locations, // and not to performing any pruning. That is be to done below simultaneously for γ. @@ -224,11 +224,8 @@ // and not to performing any pruning. That is be to done below simultaneously for γ. let ins = &config.insertion; if ins.merge_now(&state) { - stats.merged += prox_penalty.merge_spikes_no_fitness( - &mut μ, &mut τv, &μ_base, None, τ, ε, ins, - ®, - //Some(|μ̃ : &RNDM<N, F>| calculate_residual(Pair(μ̃, &z), opA, b).norm2_squared_div2()), - ); + stats.merged += + prox_penalty.merge_spikes_no_fitness(&mut μ, &mut τv, &μ_base, τ, ε, ins, ®); } // Prune spikes with zero weight.
--- a/src/frank_wolfe.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/frank_wolfe.rs Thu Feb 26 11:36:22 2026 -0500 @@ -408,7 +408,7 @@ config: &FWConfig<F>, iterator: I, mut plotter: Plot, - μ0 : Option<RNDM<N, F>>, + μ0: Option<RNDM<N, F>>, ) -> DynResult<RNDM<N, F>> where F: Float + ToNalgebraRealField,
--- a/src/lib.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/lib.rs Thu Feb 26 11:36:22 2026 -0500 @@ -50,6 +50,7 @@ } use run::{AlgorithmConfig, DefaultAlgorithm, Named, PlotLevel, RunnableExperiment}; +use subproblem::InnerMethod; use types::{ClapFloat, Float}; use DefaultAlgorithm::*; @@ -223,6 +224,26 @@ #[arg(long, value_names = &["ε", "θ", "p"])] /// Set the tolerance to ε_k = ε/(1+θk)^p tolerance: Option<Vec<F>>, + + #[arg(long)] + /// Method for solving inner optimisation problems + inner_method: Option<InnerMethod>, + + #[arg(long)] + /// Step length parameter for inner problem + inner_τ0: Option<F>, + + #[arg(long, value_names = &["τ0", "σ0"])] + /// Dual step length parameter for inner problem + inner_pdps_τσ0: Option<Vec<F>>, + + #[arg(long, value_names = &["τ", "growth"])] + /// Inner proximal point method step length and its growth + inner_pp_τ: Option<Vec<F>>, + + #[arg(long)] + /// Inner tolerance multiplier + inner_tol: Option<F>, } /// A generic entry point for binaries based on this library
--- a/src/pdps.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/pdps.rs Thu Feb 26 11:36:22 2026 -0500 @@ -147,7 +147,7 @@ pdpsconfig: &PDPSConfig<F>, iterator: I, mut plotter: Plot, - μ0 : Option<RNDM<N, F>>, + μ0: Option<RNDM<N, F>>, ) -> DynResult<RNDM<N, F>> where F: Float + ToNalgebraRealField, @@ -206,15 +206,15 @@ let μ_base = μ.clone(); // Insert and reweigh - let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( - &mut μ, &mut τv, &μ_base, None, τ, ε, config, ®, &state, &mut stats, - )?; + let (maybe_d, _within_tolerances) = prox_penalty + .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; // Prune and possibly merge spikes if config.merge_now(&state) { - stats.merged += prox_penalty - .merge_spikes_no_fitness(&mut μ, &mut τv, &μ_base, None, τ, ε, config, ®); + stats.merged += + prox_penalty.merge_spikes_no_fitness(&mut μ, &mut τv, &μ_base, τ, ε, config, ®); } + stats.inserted += μ.len() - μ_base.len(); stats.pruned += prune_with_stats(&mut μ); // Update step length parameters
--- a/src/prox_penalty.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/prox_penalty.rs Thu Feb 26 11:36:22 2026 -0500 @@ -129,17 +129,15 @@ /// May return `τv + w` for `w` a subdifferential of the regularisation term `reg`, /// as well as an indication of whether the tolerance bounds `ε` are satisfied. /// - /// `μ_base + ν_delta` is the base point, where `μ` and `μ_base` are assumed to have the same - /// spike locations, while `ν_delta` may have different locations. - /// /// `τv` is mutable to allow [`alg_tools::bounds::MinMaxMapping`] optimisation to /// refine data. Actual values of `τv` are not supposed to be mutated. + /// + /// `stats.inserted` should not be updated by implementations of this routine, + /// only `stats.inner_iters`. fn insert_and_reweigh<I>( &self, μ: &mut DiscreteMeasure<Domain, F>, τv: &mut PreadjointCodomain, - μ_base: &DiscreteMeasure<Domain, F>, - ν_delta: Option<&DiscreteMeasure<Domain, F>>, τ: F, ε: F, config: &InsertionConfig<F>, @@ -160,7 +158,6 @@ μ: &mut DiscreteMeasure<Domain, F>, τv: &mut PreadjointCodomain, μ_base: &DiscreteMeasure<Domain, F>, - ν_delta: Option<&DiscreteMeasure<Domain, F>>, τ: F, ε: F, config: &InsertionConfig<F>, @@ -177,7 +174,6 @@ μ: &mut DiscreteMeasure<Domain, F>, τv: &mut PreadjointCodomain, μ_base: &DiscreteMeasure<Domain, F>, - ν_delta: Option<&DiscreteMeasure<Domain, F>>, τ: F, ε: F, config: &InsertionConfig<F>, @@ -193,7 +189,6 @@ μ, τv, μ_base, - ν_delta, τ, ε, config,
--- a/src/prox_penalty/radon_squared.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/prox_penalty/radon_squared.rs Thu Feb 26 11:36:22 2026 -0500 @@ -8,8 +8,8 @@ use crate::dataterm::QuadraticDataTerm; use crate::forward_model::ForwardModel; use crate::measures::merging::SpikeMerging; -use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon}; -use crate::regularisation::RegTerm; +use crate::measures::{DiscreteMeasure, Radon}; +use crate::regularisation::RadonSquaredRegTerm; use crate::types::*; use alg_tools::bounds::MinMaxMapping; use alg_tools::error::DynResult; @@ -18,8 +18,6 @@ use alg_tools::linops::BoundedLinear; use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::norms::{Norm, L2}; -use anyhow::ensure; -use nalgebra::DVector; use numeric_literals::replace_float_literals; use serde::{Deserialize, Serialize}; @@ -35,7 +33,7 @@ for<'a> &'a Domain: Instance<Domain>, F: Float + ToNalgebraRealField, M: MinMaxMapping<Domain, F>, - Reg: RegTerm<Domain, F>, + Reg: RadonSquaredRegTerm<Domain, F>, DiscreteMeasure<Domain, F>: SpikeMerging<F>, { type ReturnMapping = M; @@ -48,8 +46,6 @@ &self, μ: &mut DiscreteMeasure<Domain, F>, τv: &mut M, - μ_base: &DiscreteMeasure<Domain, F>, - ν_delta: Option<&DiscreteMeasure<Domain, F>>, τ: F, ε: F, config: &InsertionConfig<F>, @@ -60,58 +56,8 @@ where I: AlgIterator, { - let mut y = μ_base.masses_dvector(); - - ensure!(μ_base.len() <= μ.len()); - - 'i_and_w: for i in 0..=1 { - // Optimise weights - if μ.len() > 0 { - // Form finite-dimensional subproblem. The subproblem references to the original μ^k - // from the beginning of the iteration are all contained in the immutable c and g. - // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional - // problems have not yet been updated to sign change. - let g̃ = DVector::from_iterator( - μ.len(), - μ.iter_locations() - .map(|ζ| -F::to_nalgebra_mixed(τv.apply(ζ))), - ); - let mut x = μ.masses_dvector(); - y.extend(std::iter::repeat(0.0.to_nalgebra_mixed()).take(0.max(x.len() - y.len()))); - assert_eq!(y.len(), x.len()); - // Solve finite-dimensional subproblem. - // TODO: This assumes that ν_delta has no common locations with μ-μ_base, to - // ignore it. - stats.inner_iters += reg.solve_findim_l1squared(&y, &g̃, τ, &mut x, ε, config); - - // Update masses of μ based on solution of finite-dimensional subproblem. - μ.set_masses_dvector(&x); - } - - if i > 0 { - // Simple debugging test to see if more inserts would be needed. Doesn't seem so. - //let n = μ.dist_matching(μ_base); - //println!("{:?}", reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n)); - break 'i_and_w; - } - - // Calculate ‖μ - μ_base‖_ℳ - // TODO: This assumes that ν_delta has no common locations with μ-μ_base. - let n = μ.dist_matching(μ_base) + ν_delta.map_or(0.0, |ν| ν.norm(Radon)); - - // Find a spike to insert, if needed. - // This only check the overall tolerances, not tolerances on support of μ-μ_base or μ, - // which are supposed to have been guaranteed by the finite-dimensional weight optimisation. - match reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n) { - None => break 'i_and_w, - Some((ξ, _v_ξ, _in_bounds)) => { - // Weight is found out by running the finite-dimensional optimisation algorithm - // above - *μ += DeltaMeasure { x: ξ, α: 0.0 }; - stats.inserted += 1; - } - }; - } + let violation = reg.find_tolerance_violation(τv, τ, ε, true, config); + reg.solve_oc_radonsq(μ, τv, τ, ε, violation, config, stats); Ok((None, true)) } @@ -121,7 +67,6 @@ μ: &mut DiscreteMeasure<Domain, F>, τv: &mut M, μ_base: &DiscreteMeasure<Domain, F>, - ν_delta: Option<&DiscreteMeasure<Domain, F>>, τ: F, ε: F, config: &InsertionConfig<F>, @@ -141,10 +86,7 @@ // after μ_candidate's extra points. // TODO: doesn't seem to work, maybe need to merge μ_base as well? // Although that doesn't seem to make sense. - let μ_radon = match ν_delta { - None => μ_candidate.sub_matching(μ_base), - Some(ν) => μ_candidate.sub_matching(μ_base) - ν, - }; + let μ_radon = μ_candidate.sub_matching(μ_base); reg.verify_merge_candidate_radonsq(τv, μ_candidate, τ, ε, &config, &μ_radon) //let n = μ_candidate.dist_matching(μ_base); //reg.find_tolerance_violation_slack(τv, τ, ε, false, config, n).is_none()
--- a/src/prox_penalty/wave.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/prox_penalty/wave.rs Thu Feb 26 11:36:22 2026 -0500 @@ -46,8 +46,6 @@ &self, μ: &mut DiscreteMeasure<Domain, F>, τv: &mut M, - μ_base: &DiscreteMeasure<Domain, F>, - ν_delta: Option<&DiscreteMeasure<Domain, F>>, τ: F, ε: F, config: &InsertionConfig<F>, @@ -67,10 +65,8 @@ _ => (config.max_insertions, !state.is_quiet()), }; - let ω0 = match ν_delta { - None => self.apply(μ_base), - Some(ν) => self.apply(μ_base + ν), - }; + let μ_base = μ.clone(); + let ω0 = self.apply(&μ_base); // Add points to support until within error tolerance or maximum insertion count reached. let mut count = 0; @@ -105,11 +101,7 @@ // Form d = τv + 𝒟μ - ω0 = τv + 𝒟(μ - μ^k) for checking the proximate optimality // conditions in the predual space, and finding new points for insertion, if necessary. - let mut d = &*τv - + match ν_delta { - None => self.preapply(μ.sub_matching(μ_base)), - Some(ν) => self.preapply(μ.sub_matching(μ_base) - ν), - }; + let mut d = &*τv + self.preapply(μ.sub_matching(&μ_base)); // If no merging heuristic is used, let's be more conservative about spike insertion, // and skip it after first round. If merging is done, being more greedy about spike @@ -135,7 +127,6 @@ // No point in optimising the weight here; the finite-dimensional algorithm is fast. *μ += DeltaMeasure { x: ξ, α: 0.0 }; count += 1; - stats.inserted += 1; }; if !within_tolerances && warn_insertions { @@ -156,7 +147,6 @@ μ: &mut DiscreteMeasure<Domain, F>, τv: &mut M, μ_base: &DiscreteMeasure<Domain, F>, - ν_delta: Option<&DiscreteMeasure<Domain, F>>, τ: F, ε: F, config: &InsertionConfig<F>, @@ -169,11 +159,7 @@ } } μ.merge_spikes(config.merging, |μ_candidate| { - let mut d = &*τv - + self.preapply(match ν_delta { - None => μ_candidate.sub_matching(μ_base), - Some(ν) => μ_candidate.sub_matching(μ_base) - ν, - }); + let mut d = &*τv + self.preapply(μ_candidate.sub_matching(μ_base)); reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, config) }) }
--- a/src/regularisation.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/regularisation.rs Thu Feb 26 11:36:22 2026 -0500 @@ -128,23 +128,6 @@ config: &InsertionConfig<F>, ) -> usize; - /// Approximately solve the problem - /// <div>$$ - /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x) - /// $$</div> - /// for $G$ depending on the trait implementation. - /// - /// Returns the number of iterations taken. - fn solve_findim_l1squared( - &self, - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, - τ: F, - x: &mut DVector<F::MixedType>, - ε: F, - config: &InsertionConfig<F>, - ) -> usize; - /// Find a point where `d` may violate the tolerance `ε`. /// /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we @@ -163,33 +146,6 @@ config: &InsertionConfig<F>, ) -> Option<(Domain, F, bool)> where - M: MinMaxMapping<Domain, F>, - { - self.find_tolerance_violation_slack(d, τ, ε, skip_by_rough_check, config, F::ZERO) - } - - /// Find a point where `d` may violate the tolerance `ε`. - /// - /// This version includes a `slack` parameter to expand the tolerances. - /// It is used for Radon-norm squared proximal term in [`crate::prox_penalty::radon_squared`]. - /// - /// If `skip_by_rough_check` is set, do not find the point if a rough check indicates that we - /// are in bounds. `ε` is the current main tolerance and `τ` a scaling factor for the - /// regulariser. - /// - /// Returns `None` if `d` is in bounds either based on the rough check, or a more precise check - /// terminating early. Otherwise returns a possibly violating point, the value of `d` there, - /// and a boolean indicating whether the found point is in bounds. - fn find_tolerance_violation_slack<M>( - &self, - d: &mut M, - τ: F, - ε: F, - skip_by_rough_check: bool, - config: &InsertionConfig<F>, - slack: F, - ) -> Option<(Domain, F, bool)> - where M: MinMaxMapping<Domain, F>; /// Verify that `d` is in bounds `ε` for a merge candidate `μ` @@ -206,6 +162,36 @@ where M: MinMaxMapping<Domain, F>; + /// TODO: document this + fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>; + + /// Returns a scaling factor for the tolerance sequence. + /// + /// Typically this is the regularisation parameter. + fn tolerance_scaling(&self) -> F; +} + +/// Regularisation term that can be used with [`crate::prox_penalty::radon_squared::RadonSquared`] +/// proximal penalty. +pub trait RadonSquaredRegTerm<Domain, F = f64>: RegTerm<Domain, F> +where + Domain: Space + Clone, + F: Float + ToNalgebraRealField, +{ + /// Adapt weights of μ, possibly insertion a new point at tolerance_violation (which should + /// be returned by [`RegTerm::find_tolerance_violation`]. + fn solve_oc_radonsq<M>( + &self, + μ: &mut DiscreteMeasure<Domain, F>, + τv: &mut M, + τ: F, + ε: F, + tolerance_violation: Option<(Domain, F, bool)>, + config: &InsertionConfig<F>, + stats: &mut IterInfo<F>, + ) where + M: Mapping<Domain, Codomain = F>; + /// Verify that `d` is in bounds `ε` for a merge candidate `μ` /// /// This version is s used for Radon-norm squared proximal term in @@ -225,14 +211,6 @@ ) -> bool where M: MinMaxMapping<Domain, F>; - - /// TODO: document this - fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>>; - - /// Returns a scaling factor for the tolerance sequence. - /// - /// Typically this is the regularisation parameter. - fn tolerance_scaling(&self) -> F; } /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`]. @@ -279,43 +257,23 @@ quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) } - fn solve_findim_l1squared( - &self, - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, - τ: F, - x: &mut DVector<F::MixedType>, - ε: F, - config: &InsertionConfig<F>, - ) -> usize { - let inner_tolerance = ε * config.inner.tolerance_mult; - let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); - l1squared_nonneg(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it) - } - #[inline] - fn find_tolerance_violation_slack<M>( + fn find_tolerance_violation<M>( &self, d: &mut M, τ: F, ε: F, skip_by_rough_check: bool, config: &InsertionConfig<F>, - slack: F, ) -> Option<(Loc<N, F>, F, bool)> where M: MinMaxMapping<Loc<N, F>, F>, { let τα = τ * self.α(); - let keep_above = -τα - slack - ε; - let minimise_below = -τα - slack - ε * config.insertion_cutoff_factor; + let keep_above = -τα - ε; + let minimise_below = -τα - ε * config.insertion_cutoff_factor; let refinement_tolerance = ε * config.refinement.tolerance_mult; - // println!( - // "keep_above: {keep_above}, rough lower bound: {}, tolerance: {ε}, slack: {slack}, τα: {τα}", - // d.bounds().lower() - // ); - // If preliminary check indicates that we are in bounds, and if it otherwise matches // the insertion strategy, skip insertion. if skip_by_rough_check && d.bounds().lower() >= keep_above { @@ -362,6 +320,80 @@ )); } + fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> { + let τα = τ * self.α(); + Some(Bounds(τα - ε, τα + ε)) + } + + fn tolerance_scaling(&self) -> F { + self.α() + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F: Float + ToNalgebraRealField, const N: usize> RadonSquaredRegTerm<Loc<N, F>, F> + for NonnegRadonRegTerm<F> +{ + fn solve_oc_radonsq<M>( + &self, + μ: &mut DiscreteMeasure<Loc<N, F>, F>, + τv: &mut M, + τ: F, + ε: F, + tolerance_violation: Option<(Loc<N, F>, F, bool)>, + config: &InsertionConfig<F>, + stats: &mut IterInfo<F>, + ) where + M: Mapping<Loc<N, F>, Codomain = F>, + { + let τα = τ * self.α(); + let mut g: Vec<_> = μ + .iter_locations() + .map(|ζ| F::to_nalgebra_mixed(-τv.apply(ζ))) + .collect(); + + let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation { + // Don't insert if existing spikes are almost as good + if g.iter().all(|minus_τv| { + -F::from_nalgebra_mixed(*minus_τv) > v_ξ + ε * config.refinement.tolerance_mult + }) { + // Weight is found out by running the finite-dimensional optimisation algorithm + // above + // NOTE: cannot set α here before y is extracted + *μ += DeltaMeasure { x: ξ, α: 0.0 /*-v_ξ - τα*/ }; + g.push(F::to_nalgebra_mixed(-v_ξ)); + Some(-v_ξ - τα) + } else { + None + } + } else { + None + }; + + // Optimise weights + if μ.len() > 0 { + // Form finite-dimensional subproblem. The subproblem references to the original μ^k + // from the beginning of the iteration are all contained in the immutable c and g. + // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional + // problems have not yet been updated to sign change. + let y = μ.masses_dvector(); + let mut x = y.clone(); + let g_na = DVector::from_vec(g); + if let (Some(β), Some(dest)) = (new_spike_initial_weight, x.as_mut_slice().last_mut()) + { + *dest = F::to_nalgebra_mixed(β); + } + // Solve finite-dimensional subproblem. + let inner_tolerance = ε * config.inner.tolerance_mult; + let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); + stats.inner_iters += + l1squared_nonneg(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it); + + // Update masses of μ based on solution of finite-dimensional subproblem. + μ.set_masses_dvector(&x); + } + } + fn verify_merge_candidate_radonsq<M>( &self, d: &mut M, @@ -407,15 +439,6 @@ ) }; } - - fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> { - let τα = τ * self.α(); - Some(Bounds(τα - ε, τα + ε)) - } - - fn tolerance_scaling(&self) -> F { - self.α() - } } #[replace_float_literals(F::cast_from(literal))] @@ -461,37 +484,22 @@ quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) } - fn solve_findim_l1squared( - &self, - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, - τ: F, - x: &mut DVector<F::MixedType>, - ε: F, - config: &InsertionConfig<F>, - ) -> usize { - let inner_tolerance = ε * config.inner.tolerance_mult; - let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); - l1squared_unconstrained(y, g, τ * self.α(), 1.0, x, &config.inner, inner_it) - } - - fn find_tolerance_violation_slack<M>( + fn find_tolerance_violation<M>( &self, d: &mut M, τ: F, ε: F, skip_by_rough_check: bool, config: &InsertionConfig<F>, - slack: F, ) -> Option<(Loc<N, F>, F, bool)> where M: MinMaxMapping<Loc<N, F>, F>, { let τα = τ * self.α(); - let keep_below = τα + slack + ε; - let keep_above = -(τα + slack) - ε; - let maximise_above = τα + slack + ε * config.insertion_cutoff_factor; - let minimise_below = -(τα + slack) - ε * config.insertion_cutoff_factor; + let keep_below = τα + ε; + let keep_above = -τα - ε; + let maximise_above = τα + ε * config.insertion_cutoff_factor; + let minimise_below = -τα - ε * config.insertion_cutoff_factor; let refinement_tolerance = ε * config.refinement.tolerance_mult; // If preliminary check indicates that we are in bonds, and if it otherwise matches @@ -568,6 +576,81 @@ )); } + fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> { + let τα = τ * self.α(); + Some(Bounds(-τα - ε, τα + ε)) + } + + fn tolerance_scaling(&self) -> F { + self.α() + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F: Float + ToNalgebraRealField, const N: usize> RadonSquaredRegTerm<Loc<N, F>, F> + for RadonRegTerm<F> +{ + fn solve_oc_radonsq<M>( + &self, + μ: &mut DiscreteMeasure<Loc<N, F>, F>, + τv: &mut M, + τ: F, + ε: F, + tolerance_violation: Option<(Loc<N, F>, F, bool)>, + config: &InsertionConfig<F>, + stats: &mut IterInfo<F>, + ) where + M: Mapping<Loc<N, F>, Codomain = F>, + { + let τα = τ * self.α(); + let mut g: Vec<_> = μ + .iter_locations() + .map(|ζ| F::to_nalgebra_mixed(τv.apply(-ζ))) + .collect(); + + let new_spike_initial_weight = if let Some((ξ, v_ξ, _in_bounds)) = tolerance_violation { + // Don't insert if existing spikes are almost as good + let n = v_ξ.abs(); + if g.iter().all(|minus_τv| { + F::from_nalgebra_mixed(*minus_τv).abs() < n - ε * config.refinement.tolerance_mult + }) { + // Weight is found out by running the finite-dimensional optimisation algorithm + // above + // NOTE: cannot initialise α before y is extracted. + *μ += DeltaMeasure { x: ξ, α: 0.0 /*-(n + τα) * v_ξ.signum()*/ }; + g.push(F::to_nalgebra_mixed(-v_ξ)); + Some(-(n + τα) * v_ξ.signum()) + } else { + None + } + } else { + None + }; + + // Optimise weights + if μ.len() > 0 { + // Form finite-dimensional subproblem. The subproblem references to the original μ^k + // from the beginning of the iteration are all contained in the immutable c and g. + // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional + // problems have not yet been updated to sign change. + let y = μ.masses_dvector(); + let mut x = y.clone(); + if let (Some(β), Some(dest)) = (new_spike_initial_weight, x.as_mut_slice().last_mut()) + { + *dest = F::to_nalgebra_mixed(β); + } + let g_na = DVector::from_vec(g); + // Solve finite-dimensional subproblem. + let inner_tolerance = ε * config.inner.tolerance_mult; + let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); + stats.inner_iters += + l1squared_unconstrained(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it); + + // Update masses of μ based on solution of finite-dimensional subproblem. + μ.set_masses_dvector(&x); + } + } + fn verify_merge_candidate_radonsq<M>( &self, d: &mut M, @@ -619,15 +702,6 @@ ) }; } - - fn target_bounds(&self, τ: F, ε: F) -> Option<Bounds<F>> { - let τα = τ * self.α(); - Some(Bounds(-τα - ε, τα + ε)) - } - - fn tolerance_scaling(&self) -> F { - self.α() - } } #[replace_float_literals(F::cast_from(literal))]
--- a/src/run.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/run.rs Thu Feb 26 11:36:22 2026 -0500 @@ -99,6 +99,17 @@ radius: cli.merge_radius.unwrap_or(g.radius), interp: cli.merge_interp.unwrap_or(g.interp), }; + let override_inner = |g: InnerSettings<F>| InnerSettings { + method: cli.inner_method.unwrap_or(g.method), + fb_τ0: cli.inner_τ0.unwrap_or(g.fb_τ0), + pdps_τσ0: cli + .inner_pdps_τσ0 + .as_ref() + .map_or(g.pdps_τσ0, |τσ0| (τσ0[0], τσ0[1])), + pp_τ: cli.inner_pp_τ.as_ref().map_or(g.pp_τ, |τ| (τ[0], τ[1])), + tolerance_mult: cli.inner_tol.unwrap_or(g.tolerance_mult), + ..g + }; let override_fb_generic = |g: InsertionConfig<F>| InsertionConfig { bootstrap_insertions: cli .bootstrap_insertions @@ -108,6 +119,7 @@ merging: override_merging(g.merging), final_merging: cli.final_merging.unwrap_or(g.final_merging), fitness_merging: cli.fitness_merging.unwrap_or(g.fitness_merging), + inner: override_inner(g.inner), tolerance: cli .tolerance .as_ref()
--- a/src/sliding_fb.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/sliding_fb.rs Thu Feb 26 11:36:22 2026 -0500 @@ -9,11 +9,12 @@ //use nalgebra::{DVector, DMatrix}; use itertools::izip; use std::iter::Iterator; +use std::ops::MulAssign; use crate::fb::*; use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess}; use crate::measures::merging::SpikeMerging; -use crate::measures::{DiscreteMeasure, Radon, RNDM}; +use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM}; use crate::plot::Plotter; use crate::prox_penalty::{ProxPenalty, StepLengthBound}; use crate::regularisation::SlidingRegTerm; @@ -25,6 +26,7 @@ use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::norms::Norm; use anyhow::ensure; +use std::ops::ControlFlow; /// Transport settings for [`pointsource_sliding_fb_reg`]. #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] @@ -36,6 +38,8 @@ 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, } #[replace_float_literals(F::cast_from(literal))] @@ -52,7 +56,7 @@ #[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 } + TransportConfig { θ0: 0.9, adaptation: 0.9, tolerance_mult_con: 100.0, max_attempts: 2 } } } @@ -98,158 +102,417 @@ FullyAdaptive { l: F, max_transport: F, g: G }, } -/// 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, const N: usize>( - γ1: &mut RNDM<N, F>, - μ: &mut RNDM<N, F>, - τ: F, - θ_or_adaptive: &mut TransportStepLength<F, G>, - v: D, -) -> (Vec<F>, RNDM<N, F>) -where - F: Float + ToNalgebraRealField, - G: Fn(F, F) -> F, - D: DifferentiableRealMapping<N, F>, -{ - use TransportStepLength::*; +#[derive(Clone, Debug, Serialize)] +pub struct SingleTransport<const N: usize, F: Float> { + /// Source point + x: Loc<N, F>, + /// Target point + y: Loc<N, F>, + /// Original mass + α_μ_orig: F, + /// Transported mass + α_γ: F, + /// Helper for pruning + prune: bool, +} + +#[derive(Clone, Debug, Serialize)] +pub struct Transport<const N: usize, F: Float> { + vec: Vec<SingleTransport<N, F>>, +} - // Save current base point and shift μ to new positions. Idea is that - // μ_base(_masses) = μ^k (vector of masses) - // μ_base_minus_γ0 = μ^k - π_♯^0γ^{k+1} - // γ1 = π_♯^1γ^{k+1} - // μ = μ^{k+1} - let μ_base_masses: Vec<F> = μ.iter_masses().collect(); - let mut μ_base_minus_γ0 = μ.clone(); // Weights will be set in the loop below. - // Construct μ^{k+1} and π_♯^1γ^{k+1} initial candidates - //let mut sum_norm_dv = 0.0; - let γ_prev_len = γ1.len(); - assert!(μ.len() >= γ_prev_len); - γ1.extend(μ[γ_prev_len..].iter().cloned()); +/// Whether partiall 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 ALLOW_PARTIAL_TRANSPORT: bool = true; +const MINIMAL_PARTIAL_TRANSPORT: bool = true; + +impl<const N: usize, F: Float> Transport<N, F> { + pub(crate) fn new() -> Self { + Transport { vec: Vec::new() } + } - // Calculate initial transport and step length. - // First calculate initial transported weights - for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) { - // If old transport has opposing sign, the new transport will be none. - ρ.α = if (ρ.α > 0.0 && δ.α < 0.0) || (ρ.α < 0.0 && δ.α > 0.0) { - 0.0 - } else { - δ.α - }; + pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<N, F>> { + self.vec.iter() + } + + pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &'_ mut SingleTransport<N, F>> { + self.vec.iter_mut() + } + + pub(crate) fn extend<I>(&mut self, it: I) + where + I: IntoIterator<Item = SingleTransport<N, F>>, + { + self.vec.extend(it) + } + + pub(crate) fn len(&self) -> usize { + self.vec.len() } - // Calculate transport rays. - match *θ_or_adaptive { - Fixed(θ) => { - let θτ = τ * θ; - for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) { - ρ.x = δ.x - v.differential(&δ.x) * (ρ.α.signum() * θτ); + // pub(crate) fn dist_matching(&self, μ: &RNDM<N, F>) -> F { + // self.iter() + // .zip(μ.iter_spikes()) + // .map(|(ρ, δ)| (ρ.α_γ - δ.α).abs()) + // .sum() + // } + + /// Construct `μ̆`, replacing the contents of `μ`. + #[replace_float_literals(F::cast_from(literal))] + 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; } } - AdaptiveMax { l: ℓ_F, ref mut max_transport, g: ref calculate_θ } => { - *max_transport = max_transport.max(γ1.norm(Radon)); - let θτ = τ * calculate_θ(ℓ_F, *max_transport); - for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) { - ρ.x = δ.x - v.differential(&δ.x) * (ρ.α.signum() * θτ); + + // Then source points with partial transport + let mut i = self.len(); + if ALLOW_PARTIAL_TRANSPORT { + // This can cause the number of points to explode, so cannot have partial transport. + for ρ in self.iter() { + let α = ρ.α_μ_orig - ρ.α_γ; + if ρ.α_γ.abs() > F::EPSILON && α != 0.0 { + let δ = DeltaMeasure { α, x: ρ.x }; + if i < μ.len() { + μ[i] = δ; + } else { + μ.push(δ) + } + i += 1; + } } } - FullyAdaptive { l: ref mut adaptive_ℓ_F, ref mut max_transport, g: ref calculate_θ } => { - *max_transport = max_transport.max(γ1.norm(Radon)); - let mut θ = calculate_θ(*adaptive_ℓ_F, *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 izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) { - let dv_x = v.differential(&δ.x); - let g = &dv_x * (ρ.α.signum() * θ * τ); - ρ.x = δ.x - g; - let n = g.norm2(); - if n >= F::EPSILON { - // Estimate Lipschitz factor of ∇v - let this_ℓ_F = (dv_x - v.differential(&ρ.x)).norm2() / n; - *adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F); - θ = calculate_θ(*adaptive_ℓ_F, *max_transport); - changes = true + μ.truncate(i); + } + + /// 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<G, D>( + &mut self, + μ: &RNDM<N, F>, + _τ: F, + τθ_or_adaptive: &mut TransportStepLength<F, G>, + v: D, + ) 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 (ρ.α_γ > 0.0 && δ.α < 0.0) || (ρ.α_γ < 0.0 && δ.α > 0.0) { + 0.0 + } else { + δ.α + } + } + + let γ_prev_len = self.len(); + assert!(μ.len() >= γ_prev_len); + self.extend(μ[γ_prev_len..].iter().map(|δ| SingleTransport { + x: δ.x, + y: δ.x, // Just something, will be filled properly in the next phase + α_μ_orig: δ.α, + α_γ: δ.α, + prune: false, + })); + + // Calculate transport rays. + match *τθ_or_adaptive { + Fixed(θ) => { + for ρ in self.iter_mut() { + ρ.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() * θτ); + } + } + FullyAdaptive { + l: ref mut adaptive_ℓ_F, + ref mut max_transport, + g: ref calculate_θτ, + } => { + *max_transport = max_transport.max(self.norm(Radon)); + let mut θτ = calculate_θτ(*adaptive_ℓ_F, *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); + let g = &dv_x * (ρ.α_γ.signum() * θτ); + ρ.y = ρ.x - g; + let n = g.norm2(); + if n >= F::EPSILON { + // Estimate Lipschitz factor of ∇v + let this_ℓ_F = (dv_x - v.differential(&ρ.y)).norm2() / n; + *adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F); + θτ = calculate_θτ(*adaptive_ℓ_F, *max_transport); + changes = true + } } - } - if !changes { - break; + if !changes { + break; + } } } } } - // 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 - } else { - δ.α = β; + /// A posteriori transport adaptation. + #[replace_float_literals(F::cast_from(literal))] + pub(crate) fn aposteriori_transport<D>( + &mut self, + μ: &RNDM<N, F>, + μ̆: &RNDM<N, F>, + _v: &mut D, + extra: Option<F>, + ε: F, + tconfig: &TransportConfig<F>, + attempts: &mut usize, + ) -> bool + where + D: DifferentiableRealMapping<N, F>, + { + *attempts += 1; + + // 1. If π_♯^1γ^{k+1} = γ1 has non-zero mass at some point y, but μ = μ^{k+1} does not, + // then the ansatz ∇w̃_x(y) = w^{k+1}(y) may not be satisfied. So set the mass of γ1 + // at that point to zero, and retry. + let mut all_ok = true; + for (δ, ρ) in izip!(μ.iter_spikes(), self.iter_mut()) { + if δ.α == 0.0 && ρ.α_γ != 0.0 { + all_ok = false; + ρ.α_γ = 0.0; + } } + + // 2. Through bounding ∫ B_ω(y, z) dλ(x, y, z). + // through the estimate ≤ C ‖Δ‖‖γ^{k+1}‖ for Δ := μ^{k+1}-μ̆^k + // which holds for some some C if the convolution kernel in 𝒟 has Lipschitz gradient. + let nγ = self.norm(Radon); + let nΔ = μ.dist_matching(&μ̆) + extra.unwrap_or(0.0); + let t = ε * tconfig.tolerance_mult_con; + if nγ * nΔ > t && *attempts >= tconfig.max_attempts { + all_ok = false; + } else if nγ * nΔ > t { + // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1, + // this will guarantee that eventually ‖γ‖ decreases sufficiently that we + // will not enter here. + //*self *= tconfig.adaptation * t / (nγ * nΔ); + + // We want a consistent behaviour that has the potential to set many weights to zero. + // Therefore, we find the smallest uniform reduction `chg_one`, subtracted + // from all weights, that achieves total `adapt` adaptation. + let adapt_to = tconfig.adaptation * t / nΔ; + let reduction_target = nγ - adapt_to; + assert!(reduction_target > 0.0); + if ALLOW_PARTIAL_TRANSPORT { + if MINIMAL_PARTIAL_TRANSPORT { + // This reduces weights of transport, starting from … until `adapt` is + // exhausted. It will, therefore, only ever cause one extrap point insertion + // at the sources, unlike “full” partial transport. + //let refs = self.vec.iter_mut().collect::<Vec<_>>(); + //refs.sort_by(|ρ1, ρ2| ρ1.α_γ.abs().partial_cmp(&ρ2.α_γ.abs()).unwrap()); + // let mut it = refs.into_iter(); + // + // Maybe sort by differential norm + // let mut refs = self + // .vec + // .iter_mut() + // .map(|ρ| { + // let val = v.differential(&ρ.x).norm2_squared(); + // (ρ, val) + // }) + // .collect::<Vec<_>>(); + // refs.sort_by(|(_, v1), (_, v2)| v2.partial_cmp(&v1).unwrap()); + // let mut it = refs.into_iter().map(|(ρ, _)| ρ); + let mut it = self.vec.iter_mut().rev(); + let _unused = it.try_fold(reduction_target, |left, ρ| { + let w = ρ.α_γ.abs(); + if left <= w { + ρ.α_γ = ρ.α_γ.signum() * (w - left); + ControlFlow::Break(()) + } else { + ρ.α_γ = 0.0; + ControlFlow::Continue(left - w) + } + }); + } else { + // This version equally reduces all weights. It causes partial transport, which + // has the problem that that we need to then adapt weights in both start and + // end points, in insert_and_reweigh, somtimes causing the number of spikes μ + // to explode. + let mut abs_weights = self + .vec + .iter() + .map(|ρ| ρ.α_γ.abs()) + .filter(|t| *t > F::EPSILON) + .collect::<Vec<F>>(); + abs_weights.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let n = abs_weights.len(); + // Cannot have partial transport; can cause spike count explosion + let chg = abs_weights.into_iter().zip((1..=n).rev()).try_fold( + 0.0, + |smaller_total, (w, m)| { + let mf = F::cast_from(m); + let reduction = w * mf + smaller_total; + if reduction >= reduction_target { + ControlFlow::Break((reduction_target - smaller_total) / mf) + } else { + ControlFlow::Continue(smaller_total + w) + } + }, + ); + match chg { + ControlFlow::Continue(_) => self.vec.iter_mut().for_each(|δ| δ.α_γ = 0.0), + ControlFlow::Break(chg_one) => self.vec.iter_mut().for_each(|ρ| { + let t = ρ.α_γ.abs(); + if t > 0.0 { + if ALLOW_PARTIAL_TRANSPORT { + let new = (t - chg_one).max(0.0); + ρ.α_γ = ρ.α_γ.signum() * new; + } + } + }), + } + } + } else { + // This version zeroes smallest weights, avoiding partial transport. + let mut abs_weights_idx = self + .vec + .iter() + .map(|ρ| ρ.α_γ.abs()) + .zip(0..) + .filter(|(w, _)| *w >= 0.0) + .collect::<Vec<(F, usize)>>(); + abs_weights_idx.sort_by(|(a, _), (b, _)| a.partial_cmp(b).unwrap()); + + let mut left = reduction_target; + + for (w, i) in abs_weights_idx { + left -= w; + let ρ = &mut self.vec[i]; + ρ.α_γ = 0.0; + if left < 0.0 { + break; + } + } + } + + all_ok = false + } + + if !all_ok && *attempts >= tconfig.max_attempts { + for ρ in self.iter_mut() { + ρ.α_γ = 0.0; + } + } + + all_ok } - // Calculate μ^k-π_♯^0γ^{k+1} and v̆ = A_*(A[μ_transported + μ_transported_base]-b) - μ_base_minus_γ0.set_masses( - μ_base_masses + + /// Returns $‖μ\^k - π\_♯\^0γ\^{k+1}‖$ + pub(crate) fn μ0_minus_γ0_radon(&self) -> F { + self.vec.iter().map(|ρ| (ρ.α_μ_orig - ρ.α_γ).abs()).sum() + } + + /// Returns $∫ c_2 d|γ|$ + #[replace_float_literals(F::cast_from(literal))] + pub(crate) fn c2integral(&self) -> F { + self.vec .iter() - .zip(γ1.iter_masses()) - .map(|(&a, b)| a - b), - ); - (μ_base_masses, μ_base_minus_γ0) + .map(|ρ| ρ.y.dist2_squared(&ρ.x) / 2.0 * ρ.α_γ.abs()) + .sum() + } + + #[replace_float_literals(F::cast_from(literal))] + pub(crate) fn get_transport_stats(&self, stats: &mut IterInfo<F>, μ: &RNDM<N, 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(); + + let transp = stats.get_transport_mut(); + + transp.dist = { + let (a, b) = transp.dist; + (a + self.c2integral(), b + self.norm(Radon)) + }; + transp.untransported_fraction = { + let (a, b) = transp.untransported_fraction; + let source = self.iter().map(|ρ| ρ.α_μ_orig.abs()).sum(); + (a + self.μ0_minus_γ0_radon(), b + source) + }; + transp.transport_error = { + let (a, b) = transp.transport_error; + //(a + self.dist_matching(&μ), b + self.norm(Radon)) + + // This ignores points that have been not transported at all, to only calculate + // destnation error; untransported_fraction accounts for not being able to transport + // at all. + self.iter() + .zip(μ.iter_spikes()) + .fold((a, b), |(a, b), (ρ, δ)| { + let transported = ρ.α_γ.abs(); + if transported > F::EPSILON { + (a + (ρ.α_γ - δ.α).abs(), b + transported) + } else { + (a, b) + } + }) + }; + } + + /// 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 RNDM<N, F>, stats: &mut IterInfo<F>) { + assert!(self.vec.len() <= μ.len()); + let old_len = μ.len(); + for (ρ, δ) in self.vec.iter_mut().zip(μ.iter_spikes()) { + ρ.prune = !(δ.α.abs() > F::EPSILON); + } + μ.prune_by(|δ| δ.α.abs() > F::EPSILON); + stats.pruned += old_len - μ.len(); + self.vec.retain(|ρ| !ρ.prune); + assert!(self.vec.len() <= μ.len()); + } } -/// A posteriori transport adaptation. -#[replace_float_literals(F::cast_from(literal))] -pub(crate) fn aposteriori_transport<F, const N: usize>( - γ1: &mut RNDM<N, F>, - μ: &mut RNDM<N, F>, - μ_base_minus_γ0: &mut RNDM<N, F>, - μ_base_masses: &Vec<F>, - extra: Option<F>, - ε: F, - tconfig: &TransportConfig<F>, -) -> bool -where - F: Float + ToNalgebraRealField, -{ - // 1. If π_♯^1γ^{k+1} = γ1 has non-zero mass at some point y, but μ = μ^{k+1} does not, - // then the ansatz ∇w̃_x(y) = w^{k+1}(y) may not be satisfied. So set the mass of γ1 - // at that point to zero, and retry. - let mut all_ok = true; - for (α_μ, α_γ1) in izip!(μ.iter_masses(), γ1.iter_masses_mut()) { - if α_μ == 0.0 && *α_γ1 != 0.0 { - all_ok = false; - *α_γ1 = 0.0; +impl<const N: usize, F: Float> Norm<Radon, F> for Transport<N, F> { + fn norm(&self, _: Radon) -> F { + self.iter().map(|ρ| ρ.α_γ.abs()).sum() + } +} + +impl<const N: usize, F: Float> MulAssign<F> for Transport<N, F> { + fn mul_assign(&mut self, factor: F) { + for ρ in self.iter_mut() { + ρ.α_γ *= factor; } } - - // 2. Through bounding ∫ B_ω(y, z) dλ(x, y, z). - // through the estimate ≤ C ‖Δ‖‖γ^{k+1}‖ for Δ := μ^{k+1}-μ^k-(π_♯^1-π_♯^0)γ^{k+1}, - // which holds for some some C if the convolution kernel in 𝒟 has Lipschitz gradient. - let nγ = γ1.norm(Radon); - let nΔ = μ_base_minus_γ0.norm(Radon) + μ.dist_matching(&γ1) + extra.unwrap_or(0.0); - let t = ε * tconfig.tolerance_mult_con; - if nγ * nΔ > t { - // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1, - // this will guarantee that eventually ‖γ‖ decreases sufficiently that we - // will not enter here. - *γ1 *= tconfig.adaptation * t / (nγ * nΔ); - all_ok = false - } - - if !all_ok { - // Update weights for μ_base_minus_γ0 = μ^k - π_♯^0γ^{k+1} - μ_base_minus_γ0.set_masses( - μ_base_masses - .iter() - .zip(γ1.iter_masses()) - .map(|(&a, b)| a - b), - ); - } - - all_ok } /// Iteratively solve the pointsource localisation problem using sliding forward-backward @@ -284,7 +547,7 @@ // Initialise iterates let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new()); - let mut γ1 = DiscreteMeasure::new(); + let mut γ = Transport::new(); // Set up parameters // let opAnorm = opA.opnorm_bound(Radon, L2); @@ -293,24 +556,27 @@ //let ℓ = opA.transport.lipschitz_factor(L2Squared) * max_transport; let ℓ = 0.0; let τ = config.τ0 / prox_penalty.step_length_bound(&f)?; - let (maybe_ℓ_F, maybe_transport_lip) = f.curvature_bound_components(config.guess); - let transport_lip = maybe_transport_lip?; - let calculate_θ = |ℓ_F, max_transport| { - let ℓ_r = transport_lip * max_transport; - config.transport.θ0 / (τ * (ℓ + ℓ_F + ℓ_r)) - }; - let mut θ_or_adaptive = match maybe_ℓ_F { - //Some(ℓ_F0) => TransportStepLength::Fixed(calculate_θ(ℓ_F0 * b.norm2(), 0.0)), - Ok(ℓ_F) => TransportStepLength::AdaptiveMax { - l: ℓ_F, // TODO: could estimate computing the real reesidual - max_transport: 0.0, - g: calculate_θ, - }, - Err(_) => TransportStepLength::FullyAdaptive { - l: 10.0 * F::EPSILON, // Start with something very small to estimate differentials - max_transport: 0.0, - g: calculate_θ, - }, + + let mut θ_or_adaptive = match f.curvature_bound_components(config.guess) { + (_, Err(_)) => TransportStepLength::Fixed(config.transport.θ0), + (maybe_ℓ_F, Ok(transport_lip)) => { + let calculate_θτ = move |ℓ_F, max_transport| { + let ℓ_r = transport_lip * max_transport; + config.transport.θ0 / (ℓ + ℓ_F + ℓ_r) + }; + match maybe_ℓ_F { + Ok(ℓ_F) => TransportStepLength::AdaptiveMax { + l: ℓ_F, // TODO: could estimate computing the real reesidual + max_transport: 0.0, + g: calculate_θτ, + }, + Err(_) => TransportStepLength::FullyAdaptive { + l: 10.0 * F::EPSILON, // Start with something very small to estimate differentials + max_transport: 0.0, + g: calculate_θτ, + }, + } + } }; // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled // by τ compared to the conditional gradient approach. @@ -331,25 +597,28 @@ for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) { // Calculate initial transport let v = f.differential(&μ); - let (μ_base_masses, mut μ_base_minus_γ0) = - initial_transport(&mut γ1, &mut μ, τ, &mut θ_or_adaptive, v); + γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v); + + 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. - let (maybe_d, _within_tolerances, mut τv̆) = 'adapt_transport: loop { + let (maybe_d, _within_tolerances, mut τv̆, μ̆) = 'adapt_transport: loop { + // Set initial guess for μ=μ^{k+1}. + γ.μ̆_into(&mut μ); + let μ̆ = μ.clone(); + // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) - //let residual_μ̆ = calculate_residual2(&γ1, &μ_base_minus_γ0, opA, b); + //let residual_μ̆ = calculate_residual2(&γ1, &μ0_minus_γ0, opA, b); // TODO: this could be optimised by doing the differential like the // old residual2. - let μ̆ = &γ1 + &μ_base_minus_γ0; - let mut τv̆ = f.differential(μ̆) * τ; + // NOTE: This assumes that μ = γ1 + let mut τv̆ = f.differential(&μ̆) * τ; // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh( &mut μ, &mut τv̆, - &γ1, - Some(&μ_base_minus_γ0), τ, ε, &config.insertion, @@ -359,59 +628,33 @@ )?; // A posteriori transport adaptation. - if aposteriori_transport( - &mut γ1, - &mut μ, - &mut μ_base_minus_γ0, - &μ_base_masses, - None, - ε, - &config.transport, - ) { - break 'adapt_transport (maybe_d, within_tolerances, τv̆); + if γ.aposteriori_transport(&μ, &μ̆, &mut τv̆, None, ε, &config.transport, &mut attempts) + { + break 'adapt_transport (maybe_d, within_tolerances, τv̆, μ̆); } + + stats.get_transport_mut().readjustment_iters += 1; }; - stats.untransported_fraction = Some({ - assert_eq!(μ_base_masses.len(), γ1.len()); - let (a, b) = stats.untransported_fraction.unwrap_or((0.0, 0.0)); - let source = μ_base_masses.iter().map(|v| v.abs()).sum(); - (a + μ_base_minus_γ0.norm(Radon), b + source) - }); - stats.transport_error = Some({ - assert_eq!(μ_base_masses.len(), γ1.len()); - let (a, b) = stats.transport_error.unwrap_or((0.0, 0.0)); - (a + μ.dist_matching(&γ1), b + γ1.norm(Radon)) - }); + γ.get_transport_stats(&mut stats, &μ); // Merge spikes. // This crucially expects the merge routine to be stable with respect to spike locations, // and not to performing any pruning. That is be to done below simultaneously for γ. - let ins = &config.insertion; - if ins.merge_now(&state) { + if config.insertion.merge_now(&state) { stats.merged += prox_penalty.merge_spikes( &mut μ, &mut τv̆, - &γ1, - Some(&μ_base_minus_γ0), + &μ̆, τ, ε, - ins, + &config.insertion, ®, Some(|μ̃: &RNDM<N, F>| f.apply(μ̃)), ); } - // Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the - // latter needs to be pruned when μ is. - // TODO: This could do with a two-vector Vec::retain to avoid copies. - let μ_new = DiscreteMeasure::from_iter(μ.iter_spikes().filter(|δ| δ.α != F::ZERO).cloned()); - if μ_new.len() != μ.len() { - let mut μ_iter = μ.iter_spikes(); - γ1.prune_by(|_| μ_iter.next().unwrap().α != F::ZERO); - stats.pruned += μ.len() - μ_new.len(); - μ = μ_new; - } + γ.prune_compat(&mut μ, &mut stats); let iter = state.iteration(); stats.this_iters += 1;
--- a/src/sliding_pdps.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/sliding_pdps.rs Thu Feb 26 11:36:22 2026 -0500 @@ -6,13 +6,11 @@ use crate::fb::*; use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess}; use crate::measures::merging::SpikeMerging; -use crate::measures::{DiscreteMeasure, Radon, RNDM}; +use crate::measures::{DiscreteMeasure, RNDM}; use crate::plot::Plotter; use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair}; use crate::regularisation::SlidingRegTerm; -use crate::sliding_fb::{ - aposteriori_transport, initial_transport, SlidingFBConfig, TransportConfig, TransportStepLength, -}; +use crate::sliding_fb::{SlidingFBConfig, Transport, TransportConfig, TransportStepLength}; use crate::types::*; use alg_tools::convex::{Conjugable, Prox, Zero}; use alg_tools::direct_product::Pair; @@ -24,13 +22,10 @@ }; use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Instance}; use alg_tools::nalgebra_support::ToNalgebraRealField; -use alg_tools::norms::{Norm, L2}; +use alg_tools::norms::L2; use anyhow::ensure; use numeric_literals::replace_float_literals; use serde::{Deserialize, Serialize}; -//use colored::Colorize; -//use nalgebra::{DVector, DMatrix}; -use std::iter::Iterator; /// Settings for [`pointsource_sliding_pdps_pair`]. #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] @@ -148,7 +143,7 @@ // Initialise iterates let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new()); - let mut γ1 = DiscreteMeasure::new(); + let mut γ = Transport::new(); //let zero_z = z.similar_origin(); // Set up parameters @@ -186,22 +181,25 @@ let κ = τ * σ_d * ψ / ((1.0 - β) * ψ - τ * σ_d * bigM); // The factor two in the manuscript disappears due to the definition of 𝚹 being // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2. - let (maybe_ℓ_F, maybe_transport_lip) = f.curvature_bound_components(config.guess); - let transport_lip = maybe_transport_lip?; - let calculate_θ = |ℓ_F, max_transport| { - let ℓ_r = transport_lip * max_transport; - config.transport.θ0 / (τ * (ℓ + ℓ_F + ℓ_r) + κ * bigθ * max_transport) - }; - let mut θ_or_adaptive = match maybe_ℓ_F { - // We assume that the residual is decreasing. - Ok(ℓ_F) => TransportStepLength::AdaptiveMax { - l: ℓ_F, // TODO: could estimate computing the real reesidual - max_transport: 0.0, - g: calculate_θ, - }, - Err(_) => { - TransportStepLength::FullyAdaptive { - l: F::EPSILON, max_transport: 0.0, g: calculate_θ + + let mut θ_or_adaptive = match f.curvature_bound_components(config.guess) { + (_, Err(_)) => TransportStepLength::Fixed(config.transport.θ0), + (maybe_ℓ_F, Ok(transport_lip)) => { + let calculate_θτ = move |ℓ_F, max_transport| { + let ℓ_r = transport_lip * max_transport; + config.transport.θ0 / ((ℓ + ℓ_F + ℓ_r) + κ * bigθ * max_transport / τ) + }; + match maybe_ℓ_F { + Ok(ℓ_F) => TransportStepLength::AdaptiveMax { + l: ℓ_F, // TODO: could estimate computing the real reesidual + max_transport: 0.0, + g: calculate_θτ, + }, + Err(_) => TransportStepLength::FullyAdaptive { + l: F::EPSILON, // Start with something very small to estimate differentials + max_transport: 0.0, + g: calculate_θτ, + }, } } }; @@ -243,26 +241,25 @@ //dbg!(&μ); - let (μ_base_masses, mut μ_base_minus_γ0) = - initial_transport(&mut γ1, &mut μ, τ, &mut θ_or_adaptive, v); + γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v); + + 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. - let (maybe_d, _within_tolerances, mut τv̆, z_new) = 'adapt_transport: loop { + let (maybe_d, _within_tolerances, mut τv̆, z_new, μ̆) = 'adapt_transport: loop { + // Set initial guess for μ=μ^{k+1}. + γ.μ̆_into(&mut μ); + let μ̆ = μ.clone(); + // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) - // let residual_μ̆ = - // calculate_residual2(Pair(&γ1, &z), Pair(&μ_base_minus_γ0, &zero_z), opA, b); - // let Pair(mut τv̆, τz̆) = opA.preadjoint().apply(residual_μ̆ * τ); - // TODO: might be able to optimise the measure sum working as calculate_residual2 above. - let Pair(mut τv̆, τz̆) = f.differential(Pair(&γ1 + &μ_base_minus_γ0, &z)) * τ; + let Pair(mut τv̆, τz̆) = f.differential(Pair(&μ̆, &z)) * τ; // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0); // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh( &mut μ, &mut τv̆, - &γ1, - Some(&μ_base_minus_γ0), τ, ε, &config.insertion, @@ -277,59 +274,37 @@ z_new = fnR.prox(σ_p, z_new + &z); // A posteriori transport adaptation. - if aposteriori_transport( - &mut γ1, - &mut μ, - &mut μ_base_minus_γ0, - &μ_base_masses, + if γ.aposteriori_transport( + &μ, + &μ̆, + &mut τv̆, Some(z_new.dist2(&z)), ε, &config.transport, + &mut attempts, ) { - break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new); + break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new, μ̆); } }; - stats.untransported_fraction = Some({ - assert_eq!(μ_base_masses.len(), γ1.len()); - let (a, b) = stats.untransported_fraction.unwrap_or((0.0, 0.0)); - let source = μ_base_masses.iter().map(|v| v.abs()).sum(); - (a + μ_base_minus_γ0.norm(Radon), b + source) - }); - stats.transport_error = Some({ - assert_eq!(μ_base_masses.len(), γ1.len()); - let (a, b) = stats.transport_error.unwrap_or((0.0, 0.0)); - (a + μ.dist_matching(&γ1), b + γ1.norm(Radon)) - }); + γ.get_transport_stats(&mut stats, &μ); // Merge spikes. // This crucially expects the merge routine to be stable with respect to spike locations, // and not to performing any pruning. That is be to done below simultaneously for γ. - let ins = &config.insertion; - if ins.merge_now(&state) { + if config.insertion.merge_now(&state) { stats.merged += prox_penalty.merge_spikes_no_fitness( &mut μ, &mut τv̆, - &γ1, - Some(&μ_base_minus_γ0), + &μ̆, τ, ε, - ins, + &config.insertion, ®, - //Some(|μ̃ : &RNDM<N, F>| calculate_residual(Pair(μ̃, &z), opA, b).norm2_squared_div2()), ); } - // Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the - // latter needs to be pruned when μ is. - // TODO: This could do with a two-vector Vec::retain to avoid copies. - let μ_new = DiscreteMeasure::from_iter(μ.iter_spikes().filter(|δ| δ.α != F::ZERO).cloned()); - if μ_new.len() != μ.len() { - let mut μ_iter = μ.iter_spikes(); - γ1.prune_by(|_| μ_iter.next().unwrap().α != F::ZERO); - stats.pruned += μ.len() - μ_new.len(); - μ = μ_new; - } + γ.prune_compat(&mut μ, &mut stats); // Do dual update // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}]
--- a/src/subproblem.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/subproblem.rs Thu Feb 26 11:36:22 2026 -0500 @@ -2,25 +2,21 @@ Iterative algorithms for solving the finite-dimensional subproblem. */ -use serde::{Serialize, Deserialize}; +use alg_tools::iterate::{AlgIteratorOptions, Verbose}; +use clap::ValueEnum; use numeric_literals::replace_float_literals; -use alg_tools::iterate::{ - AlgIteratorOptions, - Verbose, - -}; +use serde::{Deserialize, Serialize}; use crate::types::*; +pub mod l1squared_nonneg; +pub mod l1squared_unconstrained; pub mod nonneg; pub mod unconstrained; -pub mod l1squared_unconstrained; -pub mod l1squared_nonneg; - /// Method for solving finite-dimensional subproblems. /// Not all outer methods necessarily support all options. -#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] +#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug, ValueEnum)] #[allow(dead_code)] pub enum InnerMethod { /// Forward-backward @@ -29,43 +25,48 @@ SSN, /// PDPS PDPS, + /// Proximal point method + PP, } /// Settings for the solution of finite-dimensional subproblems #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] -pub struct InnerSettings<F : Float> { +pub struct InnerSettings<F: Float> { /// Method - pub method : InnerMethod, + pub method: InnerMethod, /// Proportional step length ∈ [0, 1) for `InnerMethod::FB`. - pub fb_τ0 : F, + pub fb_τ0: F, /// Proportional primal and dual step lengths for `InnerMethod::PDPS`. - pub pdps_τσ0 : (F, F), + pub pdps_τσ0: (F, F), + /// Proximal point step length parameter and its growth + pub pp_τ: (F, F), /// Fraction of `tolerance` given to inner algorithm - pub tolerance_mult : F, + pub tolerance_mult: F, /// Iterator options #[serde(flatten)] - pub iterator_options : AlgIteratorOptions, + pub iterator_options: AlgIteratorOptions, } #[replace_float_literals(F::cast_from(literal))] -impl<F : Float> Default for InnerSettings<F> { +impl<F: Float> Default for InnerSettings<F> { fn default() -> Self { InnerSettings { - fb_τ0 : 0.99, - pdps_τσ0 : (1.98, 0.5), - iterator_options : AlgIteratorOptions { + fb_τ0: 0.99, + pdps_τσ0: (1.98, 0.5), + pp_τ: (100.0, 100.0), + iterator_options: AlgIteratorOptions { // max_iter cannot be very small, as initially FB needs many iterations, although // on later invocations even one or two tends to be enough - max_iter : 2000, + max_iter: 2000, // verbose_iter affects testing of sufficient convergence, so we set it to // a small value… - verbose_iter : Verbose::Every(1), + verbose_iter: Verbose::Every(1), // … but don't print out anything - quiet : true, - .. Default::default() + quiet: true, + ..Default::default() }, - method : InnerMethod::SSN, - tolerance_mult : 0.01, + method: InnerMethod::SSN, + tolerance_mult: 0.01, } } }
--- a/src/subproblem/l1squared_nonneg.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/subproblem/l1squared_nonneg.rs Thu Feb 26 11:36:22 2026 -0500 @@ -2,221 +2,214 @@ Iterative algorithms for solving the finite-dimensional subproblem with constraint. */ +use itertools::izip; use nalgebra::DVector; use numeric_literals::replace_float_literals; -use itertools::izip; //use std::iter::zip; use std::cmp::Ordering::*; -use alg_tools::iterate::{ - AlgIteratorFactory, - AlgIteratorState, -}; +use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState}; use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::norms::{Dist, L1}; -use alg_tools::nanleast::NaNLeast; +use super::l1squared_unconstrained::l1squared_prox; +use super::nonneg::nonneg_soft_thresholding; +use super::{InnerMethod, InnerSettings}; use crate::types::*; -use super::{ - InnerMethod, - InnerSettings -}; -use super::nonneg::nonneg_soft_thresholding; -use super::l1squared_unconstrained::l1squared_prox; /// Return maximum of `dist` and distnce of inteval `[lb, ub]` to zero. #[replace_float_literals(F::cast_from(literal))] -pub(super) fn max_interval_dist_to_zero<F : Float>(dist : F, lb : F, ub : F) -> F { +pub(super) fn max_interval_dist_to_zero<F: Float>(dist: F, lb: F, ub: F) -> F { if lb < 0.0 { if ub > 0.0 { dist } else { dist.max(-ub) } - } else /* ub ≥ 0.0*/ { + } else + /* ub ≥ 0.0*/ + { dist.max(lb) } } -/// Returns the ∞-norm minimal subdifferential of $x ↦ (β/2)|x-y|_1^2 - g^⊤ x + λ\|x\|₁ +δ_{≥}(x)$ at $x$. +/// Returns the ∞-norm minimal subdifferential of $x ↦ (β/2)|x-y|_1^2 - g^⊤ x + λ\|x\|₁ +δ_{≥0}(x)$ at $x$. /// /// `v` will be modified and cannot be trusted to contain useful values afterwards. #[replace_float_literals(F::cast_from(literal))] -fn min_subdifferential<F : Float + nalgebra::RealField>( - y : &DVector<F>, - x : &DVector<F>, - g : &DVector<F>, - λ : F, - β : F +fn min_subdifferential<F: Float + nalgebra::RealField>( + y: &DVector<F>, + x: &DVector<F>, + g: &DVector<F>, + λ: F, + β: F, ) -> F { let mut val = 0.0; - let tmp = β*y.dist(x, L1); + let tmp = β * y.dist(x, L1); for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) { let (mut lb, mut ub) = (-g_i, -g_i); match x_i.partial_cmp(y_i) { - Some(Greater) => { lb += tmp; ub += tmp }, - Some(Less) => { lb -= tmp; ub -= tmp }, - Some(Equal) => { lb -= tmp; ub += tmp }, - None => {}, + Some(Greater) => { + lb += tmp; + ub += tmp + } + Some(Less) => { + lb -= tmp; + ub -= tmp + } + Some(Equal) => { + lb -= tmp; + ub += tmp + } + None => {} } match x_i.partial_cmp(&0.0) { - Some(Greater) => { lb += λ; ub += λ }, + Some(Greater) => { + lb += λ; + ub += λ + } // Less should not happen - Some(Less|Equal) => { lb = F::NEG_INFINITY; ub += λ }, - None => {}, + Some(Less | Equal) => { + lb = F::NEG_INFINITY; + ub += λ + } + None => {} }; val = max_interval_dist_to_zero(val, lb, ub); } val } +// #[replace_float_literals(F::cast_from(literal))] +// fn lbd_soft_thresholding<F: Float>(v: F, λ: F, b: F) -> F { +// match (b >= 0.0, v >= b) { +// (true, false) => b, +// (true, true) => b.max(v - λ), // soft-to-b from above +// (false, true) => super::unconstrained::soft_thresholding(v, λ), +// (false, false) => 0.0.min(b.max(v + λ)), // soft-to-0 with lower bound +// } +// } + +/// Calculates $\prox\_f(x)$ for $f(x)=λ\abs{x-y} + δ\_{≥0}(x)$. +/// This is the first output. The second output is the first output $-y$, i..e, the prox +/// of $f\_0$ for the shift function $f\_0(z) = f(z+y) = λ\abs{z} + δ\_{≥-y}(z)$ +/// satisfying +/// $$ +/// \prox_f(x) = \prox\_{f\_0}(x - y) + y, +/// $$ +/// which is also used internally. +/// The third output indicates whether the output is “locked” to either $0$ (resp. $-y$ in +/// the reformulation), or $y$ ($0$ in the reformulation). +/// The fourth output indicates the next highest $λ$ where such locking would happen. #[replace_float_literals(F::cast_from(literal))] -fn lbd_soft_thresholding<F : Float>(v : F, λ : F, b : F) -> F -{ - match (b >= 0.0, v >= b) { - (true, false) => b, - (true, true) => b.max(v - λ), // soft-to-b from above - (false, true) => super::unconstrained::soft_thresholding(v, λ), - (false, false) => 0.0.min(b.max(v - λ)), // soft-to-0 with lower bound +#[inline] +fn shifted_nonneg_soft_thresholding<F: Float>(x: F, y: F, λ: F) -> (F, F, bool, Option<F>) { + let z = x - y; // The shift to f_0 + if -y >= 0.0 { + if x > λ { + (x - λ, z - λ, false, Some(x)) + } else { + (0.0, -y, true, None) + } + } else if z < 0.0 { + if λ < -x { + (0.0, -y, true, Some(-x)) + } else if z + λ < 0.0 { + (x + λ, z + λ, false, Some(-z)) + } else { + (y, 0.0, true, None) + } + } else if z > λ { + (x - λ, z - λ, false, Some(z)) + } else { + (y, 0.0, true, None) } } -/// Calculate $prox_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}_1^2 + δ_{≥0}(x)$. +/// Calculate $\prox\_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}\_1\^2 + δ\_{≥0}(x)$. /// /// To derive an algorithm for this, we can use -/// $prox_f(x) = prox_{f_0}(x - y) - y$ for -/// $f_0(z)=\frac{β}{2}\norm{z}_1^2 + δ_{≥-y}(z)$. -/// Now, the optimality conditions for $w = prox_{f_0}(x)$ are +/// $$ +/// \prox_f(x) = \prox\_{f\_0}(x - y) + y +/// \quad\text{for}\quad +/// f\_0(z)=\frac{β}{2}\norm{z}\_1\^2 + δ\_{≥-y}(z). +/// $$ +/// Now, the optimality conditions for $w = \prox\_{f\_0}(x)$ are /// $$\tag{*} -/// x ∈ w + β\norm{w}_1\sign w + N_{≥ -y}(w). +/// x ∈ w + β\norm{w}\_1\sign w + N\_{≥ -y}(w). /// $$ -/// If we know $\norm{w}_1$, then this is easily solved by lower-bounded soft-thresholding. +/// If we know $\norm{w}\_1$, then this is easily solved by lower-bounded soft-thresholding. /// We find this by sorting the elements by the distance to the 'locked' lower-bounded /// soft-thresholding target ($0$ or $-y_i$). -/// Then we loop over this sorted vector, increasing our estimate of $\norm{w}_1$ as we decide -/// that the soft-thresholding parameter `β\norm{w}_1` has to be such that the passed elements +/// Then we loop over this sorted vector, increasing our estimate of $\norm{w}\_1$ as we decide +/// that the soft-thresholding parameter $β\norm{w}\_1$ has to be such that the passed elements /// will reach their locked value (after which they cannot change anymore, for a larger /// soft-thresholding parameter. This has to be slightly more fine-grained for account -/// for the case that $-y_i<0$ and $x_i < -y_i$. +/// for the case that $-y\_i<0$ and $x\_i < -y\_i$. /// -/// Indeed, we denote by x' and w' the subset of elements such that w_i ≠ 0 and w_i > -y_i, -/// we can calculate by applying $⟨\cdot, \sign w'⟩$ to the corresponding lines of (*) that +/// Indeed, denoting by $x'$ and $w'$ the subset of elements such that $w\_i ≠ 0$ and +/// $w\_i > -y\_i$, we can calculate by applying $⟨\cdot, \sign w'⟩$ to the corresponding +/// lines of (*) that /// $$ -/// \norm{x'} = \norm{w'} + β \norm{w}_1 m. -/// $$ -/// Having a value for $t = \norm{w}-\norm{w'}$, we can then calculate +/// \norm{x'} = \norm{w'} + β \norm{w}\_1 m, /// $$ -/// \norm{x'} - t = (1+β m)\norm{w}_1, +/// where $m$ is the number of unlocked components. +/// We can now calculate the mass $t = \norm{w}-\norm{w'}$ of locked components, and so obtain /// $$ -/// from where we can calculate the soft-thresholding parameter $λ=β\norm{w}_1$. +/// \norm{x'} + t = (1+β m)\norm{w}\_1, +/// $$ +/// from where we can calculate the soft-thresholding parameter $λ=β\norm{w}\_1$. /// Since we do not actually know the unlocked elements, but just loop over all the possibilities /// for them, we have to check that $λ$ is above the current lower bound for this parameter /// (`shift` in the code), and below a value that would cause changes in the locked set /// (`max_shift` in the code). #[replace_float_literals(F::cast_from(literal))] -pub(super) fn l1squared_nonneg_prox<F :Float + nalgebra::RealField>( - sorted : &mut Vec<(F, F, F, Option<(F, F)>)>, - x : &mut DVector<F>, - y : &DVector<F>, - β : F +pub fn l1squared_nonneg_prox<F: Float + nalgebra::RealField>( + x: &mut DVector<F>, + y: &DVector<F>, + β: F, ) { // nalgebra double-definition bullshit workaround - //let max = alg_tools::NumTraitsFloat::max; let abs = alg_tools::NumTraitsFloat::abs; - - *x -= y; - - for (az_x_i, &x_i, &y_i) in izip!(sorted.iter_mut(), x.iter(), y.iter()) { - // The first component of each az_x_i contains the distance of x_i to the - // soft-thresholding limit. If it is negative, it is always reached. - // The second component contains the absolute value of the result for that component - // w_i of the solution, if the soft-thresholding limit is reached. - // This is stored here due to the sorting, although could otherwise be computed directly. - // Likewise the third component contains the absolute value of x_i. - // The fourth component contains an initial lower bound. - let a_i = abs(x_i); - let b = -y_i; - *az_x_i = match (b >= 0.0, x_i >= b) { - (true, false) => (x_i-b, b, a_i, None), // w_i=b, so sorting element negative! - (true, true) => (x_i-b, b, a_i, None), // soft-to-b from above - (false, true) => (a_i, 0.0, a_i, None), // soft-to-0 - (false, false) => (a_i, 0.0, a_i, Some((b, b-x_i))), // soft-to-0 with initial limit - }; - } - sorted.as_mut_slice() - .sort_unstable_by(|(a, _, _, _), (b, _, _, _)| NaNLeast(*a).cmp(&NaNLeast(*b))); + let min = alg_tools::NumTraitsFloat::min; - let mut nwlow = 0.0; - let mut shift = 0.0; - // This main loop is over different combinations of elements of the solution locked - // to the soft-thresholding lower bound (`0` or `-y_i`), in the sorted order of locking. - for (i, az_x_i) in izip!(0.., sorted.iter()) { - // This 'attempt is over different combinations of elements locked to the - // lower bound (`-y_i ≤ 0`). It calculates `max_shift` as the maximum shift that - // can be done until the locking would change (or become non-strictly-complementary). - // If the main rule (*) gives an estimate of `λ` that stays below `max_shift`, it is - // accepted. Otherwise `shift` is updated to `max_shift`, and we attempt again, - // with the non-locking set that participates in the calculation of `λ` then including - // the elements that are no longer locked to the lower bound. - 'attempt: loop { - let mut nwthis = 0.0; // contribution to ‖w‖ from elements with locking - //soft-thresholding parameter = `shift` - let mut nxmore = 0.0; // ‖x'‖ for those elements through to not be locked to - // neither the soft-thresholding limit, nor the lower bound - let mut nwlbd = 0.0; // contribution to ‖w‖ from those elements locked to their - // lower bound - let mut m = 0; - let mut max_shift = F::INFINITY; // maximal shift for which our estimate of the set of - // unlocked elements is valid. - let mut max_shift_from_lbd = false; // Whether max_shift comes from the next element - // or from a lower bound being reached. - for az_x_j in sorted[i as usize..].iter() { - if az_x_j.0 <= shift { - nwthis += az_x_j.1; - } else { - match az_x_j.3 { - Some((l, s)) if shift < s => { - if max_shift > s { - max_shift_from_lbd = true; - max_shift = s; - } - nwlbd += -l; - }, - _ => { - nxmore += az_x_j.2; - if m == 0 && max_shift > az_x_j.0 { - max_shift = az_x_j.0; - max_shift_from_lbd = false; - } - m += 1; - } - } - } + let mut λ = 0.0; + loop { + let mut w_locked = 0.0; + let mut n_unlocked = 0; // m + let mut x_prime = 0.0; + let mut max_shift = F::INFINITY; + for (&x_i, &y_i) in izip!(x.iter(), y.iter()) { + let (_, a_shifted, locked, next_lock) = shifted_nonneg_soft_thresholding(x_i, y_i, λ); + + if let Some(t) = next_lock { + max_shift = min(max_shift, t); + assert!(max_shift > λ); } - - // We need ‖x'‖ = ‖w'‖ + β m ‖w‖, i.e. ‖x'‖ - (‖w‖-‖w'‖)= (1 + β m)‖w‖. - let tmp = β*(nxmore - (nwlow + nwthis + nwlbd))/(1.0 + β*F::cast_from(m)); - if tmp > max_shift { - if max_shift_from_lbd { - shift = max_shift; - continue 'attempt; - } else { - break 'attempt - } - } else if tmp < shift { - // TODO: this should probably be an assert!(false) - break 'attempt; + if locked { + w_locked += abs(a_shifted); } else { - // success - x.zip_apply(y, |x_i, y_i| *x_i = y_i + lbd_soft_thresholding(*x_i, tmp, -y_i)); - return + n_unlocked += 1; + x_prime += abs(x_i - y_i); } } - shift = az_x_i.0; - nwlow += az_x_i.1; + + // We need ‖x'‖ = ‖w'‖ + β m ‖w‖, i.e. ‖x'‖ + (‖w‖-‖w'‖)= (1 + β m)‖w‖. + let λ_new = (x_prime + w_locked) / (1.0 / β + F::cast_from(n_unlocked)); + + if λ_new > max_shift { + λ = max_shift; + } else { + assert!(λ_new >= λ); + // success + x.zip_apply(y, |x_i, y_i| { + let (a, _, _, _) = shifted_nonneg_soft_thresholding(*x_i, y_i, λ_new); + //*x_i = y_i + lbd_soft_thresholding(*x_i, λ_new, -y_i) + *x_i = a; + }); + return; + } } - // TODO: this fallback has not been verified to be correct - x.zip_apply(y, |x_i, y_i| *x_i = y_i + lbd_soft_thresholding(*x_i, shift, -y_i)); } /// Proximal point method implementation of [`l1squared_nonneg`]. @@ -226,31 +219,31 @@ /// for potential performance improvements. #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] pub fn l1squared_nonneg_pp<F, I>( - y : &DVector<F::MixedType>, - g : &DVector<F::MixedType>, - λ_ : F, - β_ : F, - x : &mut DVector<F::MixedType>, - τ_ : F, - θ_ : F, - iterator : I + y: &DVector<F::MixedType>, + g: &DVector<F::MixedType>, + λ_: F, + β_: F, + x: &mut DVector<F::MixedType>, + τ_: F, + θ_: F, + iterator: I, ) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory<F> +where + F: Float + ToNalgebraRealField, + I: AlgIteratorFactory<F>, { let λ = λ_.to_nalgebra_mixed(); let β = β_.to_nalgebra_mixed(); let mut τ = τ_.to_nalgebra_mixed(); let θ = θ_.to_nalgebra_mixed(); - let mut tmp = std::iter::repeat((0.0, 0.0, 0.0, None)).take(x.len()).collect(); let mut iters = 0; iterator.iterate(|state| { // Primal step: x^{k+1} = prox_{(τβ/2)|.-y|_1^2+δ_{≥0}+}(x^k - τ(λ𝟙^⊤-g)) - x.apply(|x_i| *x_i -= τ*λ); + x.apply(|x_i| *x_i -= τ * λ); x.axpy(τ, g, 1.0); - l1squared_nonneg_prox(&mut tmp, x, y, τ*β); - + l1squared_nonneg_prox(x, y, τ * β); + iters += 1; // This gives O(1/N^2) rates due to monotonicity of function values. // Higher acceleration does not seem to be numerically stable. @@ -260,9 +253,7 @@ // Higher acceleration does not seem to be numerically stable. //τ + = F::cast_from(iters).to_nalgebra_mixed()*θ; - state.if_verbose(|| { - F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)) - }) + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) }); iters @@ -276,18 +267,19 @@ /// The parameter `θ` is used to multiply the rescale the operator (identity) of the PDPS model. #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] pub fn l1squared_nonneg_pdps<F, I>( - y : &DVector<F::MixedType>, - g : &DVector<F::MixedType>, - λ_ : F, - β_ : F, - x : &mut DVector<F::MixedType>, - τ_ : F, - σ_ : F, - θ_ : F, - iterator : I + y: &DVector<F::MixedType>, + g: &DVector<F::MixedType>, + λ_: F, + β_: F, + x: &mut DVector<F::MixedType>, + τ_: F, + σ_: F, + θ_: F, + iterator: I, ) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory<F> +where + F: Float + ToNalgebraRealField, + I: AlgIteratorFactory<F>, { let λ = λ_.to_nalgebra_mixed(); let β = β_.to_nalgebra_mixed(); @@ -301,21 +293,19 @@ iterator.iterate(|state| { // Primal step: x^{k+1} = prox_{(τβ/2)|.-y|_1^2}(x^k - τ (w^k - g)) - x.axpy(-τ*θ, &w, 1.0); + x.axpy(-τ * θ, &w, 1.0); x.axpy(τ, g, 1.0); - l1squared_prox(&mut tmp, x, y, τ*β); - + l1squared_prox(&mut tmp, x, y, τ * β); + // Dual step: w^{k+1} = proj_{[-∞,λ]}(w^k + σ(2x^{k+1}-x^k)) - w.axpy(2.0*σ*θ, x, 1.0); - w.axpy(-σ*θ, &xprev, 1.0); + w.axpy(2.0 * σ * θ, x, 1.0); + w.axpy(-σ * θ, &xprev, 1.0); w.apply(|w_i| *w_i = w_i.min(λ)); xprev.copy_from(x); - - iters +=1; - state.if_verbose(|| { - F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)) - }) + iters += 1; + + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) }); iters @@ -340,26 +330,27 @@ /// $$</div> #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] pub fn l1squared_nonneg_pdps_alt<F, I>( - y : &DVector<F::MixedType>, - g : &DVector<F::MixedType>, - λ_ : F, - β_ : F, - x : &mut DVector<F::MixedType>, - τ_ : F, - σ_ : F, - θ_ : F, - iterator : I + y: &DVector<F::MixedType>, + g: &DVector<F::MixedType>, + λ_: F, + β_: F, + x: &mut DVector<F::MixedType>, + τ_: F, + σ_: F, + θ_: F, + iterator: I, ) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory<F> +where + F: Float + ToNalgebraRealField, + I: AlgIteratorFactory<F>, { let λ = λ_.to_nalgebra_mixed(); let τ = τ_.to_nalgebra_mixed(); let σ = σ_.to_nalgebra_mixed(); let θ = θ_.to_nalgebra_mixed(); let β = β_.to_nalgebra_mixed(); - let σθ = σ*θ; - let τθ = τ*θ; + let σθ = σ * θ; + let τθ = τ * θ; let mut w = DVector::zeros(x.len()); let mut tmp = DVector::zeros(x.len()); let mut xprev = x.clone(); @@ -369,8 +360,8 @@ // Primal step: x^{k+1} = nonnegsoft_τλ(x^k - τ(θ w^k -g)) x.axpy(-τθ, &w, 1.0); x.axpy(τ, g, 1.0); - x.apply(|x_i| *x_i = nonneg_soft_thresholding(*x_i, τ*λ)); - + x.apply(|x_i| *x_i = nonneg_soft_thresholding(*x_i, τ * λ)); + // Dual step: with g(x) = (β/(2θ))‖x-y‖₁² and q = w^k + σ(2x^{k+1}-x^k), // we compute w^{k+1} = prox_{σg^*}(q) for // = q - σ prox_{g/σ}(q/σ) @@ -381,22 +372,19 @@ w.axpy(2.0, x, 1.0); w.axpy(-1.0, &xprev, 1.0); xprev.copy_from(&w); // use xprev as temporary variable - l1squared_prox(&mut tmp, &mut xprev, y, β/σθ); + l1squared_prox(&mut tmp, &mut xprev, y, β / σθ); w -= &xprev; w *= σ; xprev.copy_from(x); - + iters += 1; - state.if_verbose(|| { - F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)) - }) + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) }); iters } - /// This function applies an iterative method for the solution of the problem /// <div>$$ /// \min_{x ∈ ℝ^n} \frac{β}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x). @@ -405,16 +393,17 @@ /// This function returns the number of iterations taken. #[replace_float_literals(F::cast_from(literal))] pub fn l1squared_nonneg<F, I>( - y : &DVector<F::MixedType>, - g : &DVector<F::MixedType>, - λ : F, - β : F, - x : &mut DVector<F::MixedType>, - inner : &InnerSettings<F>, - iterator : I + y: &DVector<F::MixedType>, + g: &DVector<F::MixedType>, + λ: F, + β: F, + x: &mut DVector<F::MixedType>, + inner: &InnerSettings<F>, + iterator: I, ) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory<F> +where + F: Float + ToNalgebraRealField, + I: AlgIteratorFactory<F>, { match inner.method { InnerMethod::PDPS => { @@ -423,15 +412,12 @@ let normest = inner_θ; let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest); l1squared_nonneg_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator) - }, - InnerMethod::FB => { - // The Lipschitz factor of ∇[x ↦ g^⊤ x + λ∑x]=g - λ𝟙 is FB is just a proximal point - // method with on constraints on τ. We “accelerate” it by adding to τ the constant θ - // on each iteration. Exponential growth does not seem stable. - let inner_τ = inner.fb_τ0; - let inner_θ = inner_τ; + } + InnerMethod::PP | InnerMethod::FB => { + let inner_τ = inner.pp_τ.0; + let inner_θ = inner.pp_τ.1; l1squared_nonneg_pp(y, g, λ, β, x, inner_τ, inner_θ, iterator) - }, + } other => unimplemented!("${other:?} is unimplemented"), } }
--- a/src/types.rs Thu Feb 26 11:38:43 2026 -0500 +++ b/src/types.rs Thu Feb 26 11:36:22 2026 -0500 @@ -21,6 +21,32 @@ impl ClapFloat for f32 {} impl ClapFloat for f64 {} +/// Structure for storing transport statistics +#[derive(Debug, Clone, Serialize)] +pub struct TransportInfo<F: Float = f64> { + /// Tuple of (untransported mass, source mass) + pub untransported_fraction: (F, F), + /// Tuple of (|destination mass - transported_mass|, transported mass) + pub transport_error: (F, F), + /// Number of readjustment iterations for transport + pub readjustment_iters: usize, + /// ($∫ c_2 dγ , ∫ dγ$) + pub dist: (F, F), +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F: Float> TransportInfo<F> { + /// Initialise transport statistics + pub fn new() -> Self { + TransportInfo { + untransported_fraction: (0.0, 0.0), + transport_error: (0.0, 0.0), + readjustment_iters: 0, + dist: (0.0, 0.0), + } + } +} + /// Structure for storing iteration statistics #[derive(Debug, Clone, Serialize)] pub struct IterInfo<F: Float = f64> { @@ -38,10 +64,8 @@ pub pruned: usize, /// Number of inner iterations since last IterInfo statistic pub inner_iters: usize, - /// Tuple of (transported mass, source mass) - pub untransported_fraction: Option<(F, F)>, - /// Tuple of (|destination mass - untransported_mass|, transported mass) - pub transport_error: Option<(F, F)>, + /// Transport statistis + pub transport: Option<TransportInfo<F>>, /// Current tolerance pub ε: F, // /// Solve fin.dim problem for this measure to get the optimal `value`. @@ -61,10 +85,17 @@ inner_iters: 0, ε: F::NAN, // postprocessing : None, - untransported_fraction: None, - transport_error: None, + transport: None, } } + + /// Get mutable reference to transport statistics, creating it if it is `None`. + pub fn get_transport_mut(&mut self) -> &mut TransportInfo<F> { + if self.transport.is_none() { + self.transport = Some(TransportInfo::new()); + } + self.transport.as_mut().unwrap() + } } #[replace_float_literals(F::cast_from(literal))] @@ -74,7 +105,7 @@ { fn logrepr(&self) -> ColoredString { format!( - "{}\t| N = {}, ε = {:.8}, 𝔼inner_it = {}, 𝔼ins/mer/pru = {}/{}/{}{}{}", + "{}\t| N = {}, ε = {:.2e}, 𝔼inner_it = {}, 𝔼ins/mer/pru = {}/{}/{}{}", self.value.logrepr(), self.n_spikes, self.ε, @@ -82,23 +113,20 @@ self.inserted as float / self.this_iters.max(1) as float, self.merged as float / self.this_iters.max(1) as float, self.pruned as float / self.this_iters.max(1) as float, - match self.untransported_fraction { + match &self.transport { None => format!(""), - Some((a, b)) => - if b > 0.0 { - format!(", untransported {:.2}%", 100.0 * a / b) - } else { - format!("") - }, - }, - match self.transport_error { - None => format!(""), - Some((a, b)) => - if b > 0.0 { - format!(", transport error {:.2}%", 100.0 * a / b) - } else { - format!("") - }, + Some(t) => { + let (a1, b1) = t.untransported_fraction; + let (a2, b2) = t.transport_error; + let (a3, b3) = t.dist; + format!( + ", γ-un/er/d/it = {:.2}%/{:.2}%/{:.2e}/{:.2}", + if b1 > 0.0 { 100.0 * a1 / b1 } else { F::NAN }, + if b2 > 0.0 { 100.0 * a2 / b2 } else { F::NAN }, + if b3 > 0.0 { a3 / b3 } else { F::NAN }, + t.readjustment_iters as float / self.this_iters.max(1) as float, + ) + } } ) .as_str() @@ -120,10 +148,7 @@ #[replace_float_literals(F::cast_from(literal))] impl<F: Float> Default for RefinementSettings<F> { fn default() -> Self { - RefinementSettings { - tolerance_mult: 0.1, - max_steps: 50000, - } + RefinementSettings { tolerance_mult: 0.1, max_steps: 50000 } } }
