Mercurial > repos > pointsource_algs / changeset
changeset
Merging adjustments, parameter tuning, etc.
Thu, 23 Jan 2025 23:34:05 +0100
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Thu, 23 Jan 2025 23:34:05 +0100
- branch
- dev
- changeset 39
- 6316d68b58af
- parent 37
- c5d8bd1a7728
- child 40
- 896b42b5ac1a
Merging adjustments, parameter tuning, etc.
--- a/Cargo.lock Thu Jan 23 23:35:28 2025 +0100 +++ b/Cargo.lock Thu Jan 23 23:34:05 2025 +0100 @@ -136,6 +136,12 @@ checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" [[package]] +name = "base64" +version = "0.22.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b3254f16251a8381aa12e40e3c4d2f0199f8c6508fbecb9d91f575e0fbb8c6" + +[[package]] name = "bitflags" version = "2.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -313,12 +319,63 @@ ] [[package]] +name = "darling" +version = "0.20.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6f63b86c8a8826a49b8c21f08a2d07338eec8d900540f8630dc76284be802989" +dependencies = [ + "darling_core", + "darling_macro", +] + +[[package]] +name = "darling_core" +version = "0.20.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95133861a8032aaea082871032f5815eb9e98cef03fa916ab4500513994df9e5" +dependencies = [ + "fnv", + "ident_case", + "proc-macro2", + "quote", + "strsim", + "syn 2.0.93", +] + +[[package]] +name = "darling_macro" +version = "0.20.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d336a2a514f6ccccaa3e09b02d41d35330c07ddf03a62165fcec10bb561c7806" +dependencies = [ + "darling_core", + "quote", + "syn 2.0.93", +] + +[[package]] +name = "deranged" +version = "0.3.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b42b6fa04a440b495c8b04d0e71b707c585f83cb9cb28cf8cd0d976c315e31b4" +dependencies = [ + "powerfmt", + "serde", +] + +[[package]] name = "either" version = "1.13.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" [[package]] +name = "equivalent" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5443807d6dff69373d433ab9ef5378ad8df50ca6298caf15de6e52e24aaf54d5" + +[[package]] name = "errno" version = "0.3.5" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -338,6 +395,12 @@ ] [[package]] +name = "fnv" +version = "1.0.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1" + +[[package]] name = "getrandom" version = "0.2.10" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -349,12 +412,30 @@ ] [[package]] +name = "hashbrown" +version = "0.12.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8a9ee70c43aaf417c914396645a0fa852624801b24ebb7ae78fe8272889ac888" + +[[package]] +name = "hashbrown" +version = "0.15.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bf151400ff0baff5465007dd2f3e717f3fe502074ca563069ce3a6629d07b289" + +[[package]] name = "heck" version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" [[package]] +name = "hex" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7f24254aa9a54b5c858eaee2f5bccdb46aaf0e486a595ed5fd8f86ba55232a70" + +[[package]] name = "iana-time-zone" version = "0.1.61" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -378,6 +459,34 @@ ] [[package]] +name = "ident_case" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39" + +[[package]] +name = "indexmap" +version = "1.9.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bd070e393353796e801d209ad339e89596eb4c8d430d18ede6a1cced8fafbd99" +dependencies = [ + "autocfg", + "hashbrown 0.12.3", + "serde", +] + +[[package]] +name = "indexmap" +version = "2.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62f822373a4fe84d4bb149bf54e584a7f4abec90e072ed49cda0edea5b95471f" +dependencies = [ + "equivalent", + "hashbrown 0.15.2", + "serde", +] + +[[package]] name = "is_terminal_polyfill" version = "1.70.1" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -525,6 +634,12 @@ ] [[package]] +name = "num-conv" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "51d515d32fb182ee37cda2ccdcb92950d6a3c2893aa280e540671c2cd0f3b1d9" + +[[package]] name = "num-integer" version = "0.1.46" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -615,9 +730,16 @@ "regex", "serde", "serde_json", + "serde_with", ] [[package]] +name = "powerfmt" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "439ee305def115ba05938db6eb1644ff94165c5ab5e9420d1c1bcedbba909391" + +[[package]] name = "ppv-lite86" version = "0.2.17" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -802,6 +924,36 @@ ] [[package]] +name = "serde_with" +version = "3.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6b6f7f2fcb69f747921f79f3926bd1e203fce4fef62c268dd3abfb6d86029aa" +dependencies = [ + "base64", + "chrono", + "hex", + "indexmap 1.9.3", + "indexmap 2.7.0", + "serde", + "serde_derive", + "serde_json", + "serde_with_macros", + "time", +] + +[[package]] +name = "serde_with_macros" +version = "3.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8d00caa5193a3c8362ac2b73be6b9e768aa5a4b2f721d8f4b339600c3cb51f8e" +dependencies = [ + "darling", + "proc-macro2", + "quote", + "syn 2.0.93", +] + +[[package]] name = "shlex" version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" @@ -859,6 +1011,37 @@ ] [[package]] +name = "time" +version = "0.3.37" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "35e7868883861bd0e56d9ac6efcaaca0d6d5d82a2a7ec8209ff492c07cf37b21" +dependencies = [ + "deranged", + "itoa", + "num-conv", + "powerfmt", + "serde", + "time-core", + "time-macros", +] + +[[package]] +name = "time-core" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ef927ca75afb808a4d64dd374f00a2adf8d0fcff8e7b184af886c3c87ec4a3f3" + +[[package]] +name = "time-macros" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2834e6017e3e5e4b9834939793b282bc03b37a3336245fa820e35e233e2a85de" +dependencies = [ + "num-conv", + "time-core", +] + +[[package]] name = "typenum" version = "1.17.0" source = "registry+https://github.com/rust-lang/crates.io-index"
--- a/Cargo.toml Thu Jan 23 23:35:28 2025 +0100 +++ b/Cargo.toml Thu Jan 23 23:34:05 2025 +0100 @@ -43,6 +43,7 @@ serde_json = "~1.0.85" chrono = { version = "~0.4.23", features = ["alloc", "std", "serde"] } anyhow = "1.0.95" +serde_with = { version = "3.11.0", features = ["macros"] } [build-dependencies] regex = "~1.11.0"
--- a/src/experiments.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/experiments.rs Thu Jan 23 23:34:05 2025 +0100 @@ -13,10 +13,9 @@ use alg_tools::error::DynResult; use alg_tools::norms::Linfinity; -use crate::ExperimentOverrides; +use crate::{ExperimentOverrides, AlgorithmOverrides}; use crate::kernels::*; use crate::kernels::SupportProductFirst as Prod; -use crate::pdps::PDPSConfig; use crate::types::*; use crate::run::{ RunnableExperiment, @@ -24,8 +23,6 @@ ExperimentBiased, Named, DefaultAlgorithm, - AlgorithmConfig, - ProxTerm }; //use crate::fb::FBGenericConfig; use crate::rand_distr::{SerializableNormal, SaltAndPepper}; @@ -64,9 +61,12 @@ /// Two dimensions, “fast” spread, 1-norm data fidelity #[clap(name = "2d_l1_fast")] Experiment2D_L1_Fast, - /// One dimension, “fast” spread, 2-norm-squared data fidelity with extra TV-regularised bias + /// One dimension, “fast” spread, 2-norm-squared data fidelity with extra TV-regularised bias #[clap(name = "1d_tv_fast")] Experiment1D_TV_Fast, + /// Two dimensions, “fast” spread, 2-norm-squared data fidelity with extra TV-regularised bias + #[clap(name = "2d_tv_fast")] + Experiment2D_TV_Fast, } macro_rules! make_float_constant { @@ -146,25 +146,59 @@ make_float_constant!(HatBias = 0.05); // We use a different step length for PDPS in 2D experiments - let pdps_2d = || { - let τ0 = 3.0; - PDPSConfig { - τ0, - σ0 : 0.99 / τ0, + // let pdps_2d = (DefaultAlgorithm::PDPS, + // AlgorithmOverrides { + // tau0 : Some(3.0), + // sigma0 : Some(0.99 / 3.0), + // .. Default::default() + // } + // ); + // let radon_pdps_2d = (DefaultAlgorithm::RadonPDPS, + // AlgorithmOverrides { + // tau0 : Some(3.0), + // sigma0 : Some(0.99 / 3.0), + // .. Default::default() + // } + // ); + let sliding_fb_cut_gaussian = (DefaultAlgorithm::SlidingFB, + AlgorithmOverrides { + theta0 : Some(0.3), .. Default::default() } - }; - let defaults_2d = HashMap::from([ - (DefaultAlgorithm::PDPS, AlgorithmConfig::PDPS(pdps_2d(), ProxTerm::Wave)), - (DefaultAlgorithm::RadonPDPS, AlgorithmConfig::PDPS(pdps_2d(), ProxTerm::RadonSquared)) - ]); - + ); + // let higher_cpos = |alg| (alg, + // AlgorithmOverrides { + // transport_tolerance_pos : Some(1000.0), + // .. Default::default() + // } + // ); + let higher_cpos_merging = |alg| (alg, + AlgorithmOverrides { + transport_tolerance_pos : Some(1000.0), + merge : Some(true), + fitness_merging : Some(true), + .. Default::default() + } + ); + 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), + merge : Some(true), + fitness_merging : Some(true), + .. Default::default() + } + ); // We add a hash of the experiment name to the configured // noise seed to not use the same noise for different experiments. let mut h = DefaultHasher::new(); name.hash(&mut h); let noise_seed = cli.noise_seed.unwrap_or(BASE_SEED) + h.finish(); + let default_merge_radius = 0.01; + use DefaultExperiment::*; Ok(match self { Experiment1D => { @@ -173,7 +207,7 @@ Box::new(Named { name, data : ExperimentV2 { domain : [[0.0, 1.0]].into(), sensor_count : [N_SENSORS_1D], - regularisation : Regularisation::NonnegRadon(cli.alpha.unwrap_or(0.09)), + regularisation : Regularisation::NonnegRadon(cli.alpha.unwrap_or(0.08)), noise_distr : SerializableNormal::new(0.0, cli.variance.unwrap_or(0.2))?, dataterm : DataTerm::L2Squared, μ_hat : MU_TRUE_1D_BASIC.into(), @@ -182,7 +216,12 @@ kernel : Prod(AutoConvolution(spread_cutoff), base_spread), kernel_plot_width, noise_seed, - algorithm_defaults: HashMap::new(), + default_merge_radius, + algorithm_overrides: HashMap::from([ + sliding_fb_cut_gaussian, + higher_cpos_merging(DefaultAlgorithm::RadonFB), + higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), + ]), }}) }, Experiment1DFast => { @@ -199,7 +238,11 @@ kernel : base_spread, kernel_plot_width, noise_seed, - algorithm_defaults: HashMap::new(), + default_merge_radius, + algorithm_overrides: HashMap::from([ + higher_cpos_merging(DefaultAlgorithm::RadonFB), + higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), + ]), }}) }, Experiment2D => { @@ -217,7 +260,14 @@ kernel : Prod(AutoConvolution(spread_cutoff), base_spread), kernel_plot_width, noise_seed, - algorithm_defaults: defaults_2d, + default_merge_radius, + algorithm_overrides: HashMap::from([ + sliding_fb_cut_gaussian, + higher_cpos_merging(DefaultAlgorithm::RadonFB), + higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), + // pdps_2d, + // radon_pdps_2d + ]), }}) }, Experiment2DFast => { @@ -234,7 +284,13 @@ kernel : base_spread, kernel_plot_width, noise_seed, - algorithm_defaults: defaults_2d, + default_merge_radius, + algorithm_overrides: HashMap::from([ + // pdps_2d, + // radon_pdps_2d + higher_cpos_merging(DefaultAlgorithm::RadonFB), + higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), + ]), }}) }, Experiment1D_L1 => { @@ -255,7 +311,8 @@ kernel : Prod(AutoConvolution(spread_cutoff), base_spread), kernel_plot_width, noise_seed, - algorithm_defaults: HashMap::new(), + default_merge_radius, + algorithm_overrides: HashMap::new(), }}) }, Experiment1D_L1_Fast => { @@ -275,7 +332,8 @@ kernel : base_spread, kernel_plot_width, noise_seed, - algorithm_defaults: HashMap::new(), + default_merge_radius, + algorithm_overrides: HashMap::new(), }}) }, Experiment2D_L1 => { @@ -296,7 +354,11 @@ kernel : Prod(AutoConvolution(spread_cutoff), base_spread), kernel_plot_width, noise_seed, - algorithm_defaults: defaults_2d, + default_merge_radius, + algorithm_overrides: HashMap::from([ + // pdps_2d, + // radon_pdps_2d + ]), }}) }, Experiment2D_L1_Fast => { @@ -316,7 +378,11 @@ kernel : base_spread, kernel_plot_width, noise_seed, - algorithm_defaults: defaults_2d, + default_merge_radius, + algorithm_overrides: HashMap::from([ + // pdps_2d, + // radon_pdps_2d + ]), }}) }, Experiment1D_TV_Fast => { @@ -339,7 +405,39 @@ kernel : base_spread, kernel_plot_width, noise_seed, - algorithm_defaults: HashMap::new(), + default_merge_radius, + algorithm_overrides: HashMap::from([ + higher_cpos_merging(DefaultAlgorithm::RadonForwardPDPS), + higher_cpos_merging(DefaultAlgorithm::RadonSlidingPDPS), + ]), + }, + }}) + }, + Experiment2D_TV_Fast => { + let base_spread = HatConv { radius : Hat1 }; + Box::new(Named { name, data : ExperimentBiased { + λ : 0.005, + bias : MappingSum::new([ + Weighted::new(1.0, BallCharacteristic{ center : [0.3, 0.3].into(), radius : 0.2 }), + Weighted::new(0.5, BallCharacteristic{ center : [0.6, 0.6].into(), radius : 0.3 }), + ]), + base : ExperimentV2 { + domain : [[0.0, 1.0]; 2].into(), + sensor_count : [N_SENSORS_2D; 2], + regularisation : Regularisation::NonnegRadon(cli.alpha.unwrap_or(0.06)), + noise_distr : SerializableNormal::new(0.0, cli.variance.unwrap_or(0.15))?, //0.25 + dataterm : DataTerm::L2Squared, + μ_hat : MU_TRUE_2D_BASIC.into(), + sensor : BallIndicator { r : SensorWidth2D, exponent : Linfinity }, + spread : base_spread, + kernel : base_spread, + kernel_plot_width, + noise_seed, + default_merge_radius, + algorithm_overrides: HashMap::from([ + much_higher_cpos_merging_steptune(DefaultAlgorithm::RadonForwardPDPS), + much_higher_cpos_merging_steptune(DefaultAlgorithm::RadonSlidingPDPS), + ]), }, }}) },
--- a/src/fb.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/fb.rs Thu Jan 23 23:34:05 2025 +0100 @@ -93,10 +93,7 @@ DiscreteMeasure, RNDM, }; -use crate::measures::merging::{ - SpikeMergingMethod, - SpikeMerging, -}; +use crate::measures::merging::SpikeMerging; use crate::forward_model::{ ForwardModel, AdjointProductBoundedBy, @@ -164,7 +161,7 @@ RNDM<F, N> : SpikeMerging<F>, for<'a> &'a RNDM<F, N> : Instance<RNDM<F, N>>, { - μ.merge_spikes_fitness(config.merging, + μ.merge_spikes_fitness(config.final_merging_method(), |μ̃| dataterm.calculate_fit_op(μ̃, opA, b), |&v| v); μ.prune(); @@ -255,7 +252,10 @@ // Prune and possibly merge spikes if config.merge_now(&state) { - stats.merged += prox_penalty.merge_spikes(&mut μ, &mut τv, &μ_base, τ, ε, config, ®); + stats.merged += prox_penalty.merge_spikes( + &mut μ, &mut τv, &μ_base, None, τ, ε, config, ®, + Some(|μ̃ : &RNDM<F, N>| L2Squared.calculate_fit_op(μ̃, opA, b)), + ); } stats.pruned += prune_with_stats(&mut μ); @@ -363,15 +363,10 @@ ); // (Do not) merge spikes. - if config.merge_now(&state) { - match config.merging { - SpikeMergingMethod::None => { }, - _ => if !warned_merging { - let err = format!("Merging not supported for μFISTA"); - println!("{}", err.red()); - warned_merging = true; - } - } + if config.merge_now(&state) && !warned_merging { + let err = format!("Merging not supported for μFISTA"); + println!("{}", err.red()); + warned_merging = true; } // Update inertial prameters @@ -387,6 +382,9 @@ // stored in μ_prev. let n_before_prune = μ.len(); μ.pruning_sub(1.0 + θ, θ, &mut μ_prev); + //let μ_new = (&μ * (1.0 + θ)).sub_matching(&(&μ_prev * θ)); + // μ_prev = μ; + // μ = μ_new; debug_assert!(μ.len() <= n_before_prune); stats.pruned += n_before_prune - μ.len();
--- a/src/forward_model/sensor_grid.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/forward_model/sensor_grid.rs Thu Jan 23 23:34:05 2025 +0100 @@ -94,7 +94,8 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N>, - BT : SensorGridBT<F, S, P, N>, { + BT : SensorGridBT<F, S, P, N>, +{ domain : Cube<F, N>, sensor_count : [usize; N], sensor : S, @@ -109,7 +110,7 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, - /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>*/ { +{ /// Create a new sensor grid. /// @@ -148,7 +149,8 @@ BT : SensorGridBT<F, S, P, N>, S : Sensor<F, N>, P : Spread<F, N>, - Convolution<S, P> : Spread<F, N> { + Convolution<S, P> : Spread<F, N> +{ /// Return the grid of sensor locations. pub fn grid(&self) -> LinGrid<F, N> { @@ -190,7 +192,6 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N>, - //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, { type Codomain = DVector<F>; @@ -211,7 +212,6 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N>, - //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> { } @@ -222,7 +222,6 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N>, - //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> { fn gemv<I : Instance<RNDM<F, N>>>( @@ -268,8 +267,8 @@ BT : SensorGridBT<F, S, P, N, Agg=Bounds<F>>, S : Sensor<F, N>, P : Spread<F, N>, - Convolution<S, P> : Spread<F, N>, - ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> { + Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N> +{ /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped /// with the Radon norm, and $ℝ^n$ with the Euclidean norm. @@ -278,11 +277,14 @@ // |Aμ|_2 = sup_{|z|_2 ≤ 1} ⟨z,Αμ⟩ = sup_{|z|_2 ≤ 1} ⟨A^*z|μ⟩ // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ // = sup_{|z|_2 ≤ 1} |∑ φ(· - x_i)z_i|_∞ |μ|_ℳ - // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ ∑ |z_i| |μ|_ℳ - // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ √n |z|_2 |μ|_ℳ - // = |φ|_∞ √n |μ|_ℳ. + // ≤ sup_{|z|_2 ≤ 1} |φ(y)| ∑_{i:th sensor active at y}|z_i| |μ|_ℳ + // where the supremum of |∑ φ(· - x_i)z_i|_∞ is reached at y + // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ √N_ψ |z|_2 |μ|_ℳ + // where N_ψ is the maximum number of sensors that overlap, and + // |z|_2 is restricted to the active sensors. + // = |φ|_∞ √N_ψ |μ|_ℳ. // Hence - let n = F::cast_from(self.n_sensors()); + let n = self.max_overlapping(); self.base_sensor.bounds().uniform() * n.sqrt() } } @@ -297,9 +299,8 @@ BT : SensorGridBT<F, S, P, N>, S : Sensor<F, N>, P : Spread<F, N>, - Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, - /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, - Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ { + Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N> +{ type PreadjointCodomain = BTFN<F, SensorGridSupportGenerator<F, S, P, N>, BT, N>; type Preadjoint<'a> = SensorGridPreadjoint<'a, Self, F, N> where Self : 'a; @@ -317,8 +318,7 @@ P : Spread<F, N>, Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType=F> + DifferentiableMapping<Loc<F,N>> + LocalAnalysis<F, BT::Agg, N>, for<'b> <Convolution<S, P> as DifferentiableMapping<Loc<F,N>>>::Differential<'b> : Lipschitz<L2, FloatType=F>, - /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, - Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ { +{ type FloatType = F; @@ -345,7 +345,8 @@ pub struct SensorGridSupportGenerator<F, S, P, const N : usize> where F : Float, S : Sensor<F, N>, - P : Spread<F, N> { + P : Spread<F, N> +{ base_sensor : Convolution<S, P>, grid : LinGrid<F, N>, weights : DVector<F> @@ -355,7 +356,8 @@ where F : Float, S : Sensor<F, N>, P : Spread<F, N>, - Convolution<S, P> : Spread<F, N> { + Convolution<S, P> : Spread<F, N> +{ #[inline] fn construct_sensor(&self, id : usize, w : F) -> Weighted<ShiftedSensor<F, S, P, N>, F> { @@ -375,7 +377,8 @@ where F : Float, S : Sensor<F, N>, P : Spread<F, N>, - Convolution<S, P> : Spread<F, N> { + Convolution<S, P> : Spread<F, N> +{ type Id = usize; type SupportType = Weighted<ShiftedSensor<F, S, P, N>, F>; type AllDataIter<'a> = MapX<'a, Zip<RangeFrom<usize>, @@ -407,8 +410,7 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, - /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, - Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ { +{ type Observable = DVector<F>; fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError { @@ -428,9 +430,8 @@ BT : SensorGridBT<F, S, P, N>, S : Sensor<F, N>, P : Spread<F, N>, - Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, - /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, - Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ { + Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N> +{ fn findim_quadratic_model( &self, @@ -465,7 +466,8 @@ P : Spread<F, N>, Convolution<S, P> : Spread<F, N>, K : SimpleConvolutionKernel<F, N>, - AutoConvolution<P> : BoundedBy<F, K> { + AutoConvolution<P> : BoundedBy<F, K> +{ type FloatType = F; @@ -496,7 +498,8 @@ BT : SensorGridBT<F, S, P, N>, S : Sensor<F, N>, P : Spread<F, N>, - Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType = F> { + Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType = F> +{ type FloatType = F; fn transport_lipschitz_factor(&self, L2Squared : L2Squared) -> Self::FloatType { @@ -604,8 +607,7 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>, - //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, - /*Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ { +{ type Codomain = SensorGridBTFN<F, S, P, BT, N>; @@ -627,8 +629,4 @@ S : Sensor<F, N>, P : Spread<F, N>, Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>, - /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, - Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ { - -} - +{ }
--- a/src/forward_pdps.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/forward_pdps.rs Thu Jan 23 23:34:05 2025 +0100 @@ -51,10 +51,9 @@ #[replace_float_literals(F::cast_from(literal))] impl<F : Float> Default for ForwardPDPSConfig<F> { fn default() -> Self { - let τ0 = 0.99; ForwardPDPSConfig { - τ0, - σd0 : 0.1, + τ0 : 0.99, + σd0 : 0.05, σp0 : 0.99, insertion : Default::default() } @@ -172,7 +171,6 @@ for state in iterator.iter_init(|| full_stats(&residual, &μ, &z, ε, stats.clone())) { // Calculate initial transport let Pair(mut τv, τz) = opA.preadjoint().apply(residual * τ); - let z_base = z.clone(); let μ_base = μ.clone(); // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. @@ -182,34 +180,34 @@ ®, &state, &mut stats, ); - // // Merge spikes. - // // This expects the prune below to prune γ. - // // TODO: This may not work correctly in all cases. - // let ins = &config.insertion; - // if ins.merge_now(&state) { - // if let SpikeMergingMethod::None = ins.merging { - // } else { - // stats.merged += μ.merge_spikes(ins.merging, |μ_candidate| { - // let ν = μ_candidate.sub_matching(&γ1)-&μ_base_minus_γ0; - // let mut d = &τv̆ + op𝒟.preapply(ν); - // reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, ins) - // }); - // } - // } + // 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 γ. + // 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) { + stats.merged += prox_penalty.merge_spikes_no_fitness( + &mut μ, &mut τv, &μ_base, None, τ, ε, ins, ®, + //Some(|μ̃ : &RNDM<F, N>| calculate_residual(Pair(μ̃, &z), opA, b).norm2_squared_div2()), + ); + } // Prune spikes with zero weight. stats.pruned += prune_with_stats(&mut μ); // Do z variable primal update - z.axpy(-σ_p/τ, τz, 1.0); // TODO: simplify nasty factors - opKz.adjoint().gemv(&mut z, -σ_p, &y, 1.0); - z = fnR.prox(σ_p, z); + let mut z_new = τz; + opKz.adjoint().gemv(&mut z_new, -σ_p, &y, -σ_p/τ); + z_new = fnR.prox(σ_p, z_new + &z); // Do dual update // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] - opKz.gemv(&mut y, σ_d*(1.0 + ω), &z, 1.0); + opKz.gemv(&mut y, σ_d*(1.0 + ω), &z_new, 1.0); // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b - opKz.gemv(&mut y, -σ_d*ω, z_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b + opKz.gemv(&mut y, -σ_d*ω, z, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b y = starH.prox(σ_d, y); + z = z_new; // Update residual residual = calculate_residual(Pair(&μ, &z), opA, b); @@ -236,7 +234,7 @@ + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z)) }; - μ.merge_spikes_fitness(config.insertion.merging, fit, |&v| v); + μ.merge_spikes_fitness(config.insertion.final_merging_method(), fit, |&v| v); μ.prune(); Pair(μ, z) }
--- a/src/frank_wolfe.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/frank_wolfe.rs Thu Jan 23 23:34:05 2025 +0100 @@ -108,7 +108,7 @@ refinement : Default::default(), inner : Default::default(), variant : FWVariant::FullyCorrective, - merging : Default::default(), + merging : SpikeMergingMethod { enabled : true, ..Default::default() }, } } }
--- a/src/main.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/main.rs Thu Jan 23 23:34:05 2025 +0100 @@ -16,6 +16,7 @@ use clap::Parser; use serde::{Serialize, Deserialize}; use serde_json; +use serde_with::skip_serializing_none; use itertools::Itertools; use std::num::NonZeroUsize; @@ -56,12 +57,12 @@ AlgorithmConfig, }; use experiments::DefaultExperiment; -use measures::merging::SpikeMergingMethod; use DefaultExperiment::*; use DefaultAlgorithm::*; /// Command line parameters -#[derive(Parser, Debug, Serialize)] +#[skip_serializing_none] +#[derive(Parser, Debug, Serialize, Default, Clone)] #[clap( about = env!("CARGO_PKG_DESCRIPTION"), author = env!("CARGO_PKG_AUTHORS"), @@ -96,7 +97,7 @@ /// Not all algorithms are available for all the experiments. /// In particular, only PDPS is available for the experiments with L¹ data term. #[arg(value_enum, value_name = "ALGORITHM", long, short = 'a', - default_values_t = [FB, FISTA, PDPS, SlidingFB, FW, FWRelax])] + default_values_t = [FB, PDPS, SlidingFB, FW, RadonFB])] algorithm : Vec<DefaultAlgorithm>, /// Saved algorithm configration(s) to use on the experiments @@ -119,6 +120,10 @@ /// Number of threads. Overrides the maximum number. num_threads : Option<usize>, + #[arg(long, default_value_t = false)] + /// Load saved value ranges (if exists) to do partial update. + load_valuerange : bool, + #[clap(flatten, next_help_heading = "Experiment overrides")] /// Experiment setup overrides experiment_overrides : ExperimentOverrides<float>, @@ -129,7 +134,8 @@ } /// Command line experiment setup overrides -#[derive(Parser, Debug, Serialize, Deserialize)] +#[skip_serializing_none] +#[derive(Parser, Debug, Serialize, Deserialize, Default, Clone)] pub struct ExperimentOverrides<F : ClapFloat> { #[arg(long)] /// Regularisation parameter override. @@ -152,7 +158,8 @@ } /// Command line algorithm parametrisation overrides -#[derive(Parser, Debug, Serialize, Deserialize)] +#[skip_serializing_none] +#[derive(Parser, Debug, Serialize, Deserialize, Default, Clone)] pub struct AlgorithmOverrides<F : ClapFloat> { #[arg(long, value_names = &["COUNT", "EACH"])] /// Override bootstrap insertion iterations for --algorithm. @@ -187,12 +194,12 @@ theta0 : Option<F>, #[arg(long)] - /// Transport toleranced wrt. ω - transport_tolerance_omega : Option<F>, + /// A priori transport tolerance multiplier (C_pri) + transport_tolerance_pri : Option<F>, #[arg(long)] - /// Transport toleranced wrt. ∇v - transport_tolerance_dv : Option<F>, + /// A posteriori transport tolerance multiplier (C_pos) + transport_tolerance_pos : Option<F>, #[arg(long)] /// Transport adaptation factor. Must be in (0, 1). @@ -218,18 +225,26 @@ /// Only affects FB, FISTA, and PDPS. merge_every : Option<usize>, - #[arg(value_enum, long)]//, value_parser = SpikeMergingMethod::<float>::value_parser())] - /// Merging strategy - /// - /// Either the string "none", or a radius value for heuristic merging. - merging : Option<SpikeMergingMethod<F>>, + #[arg(long)] + /// Enable merging (default: determined by algorithm) + merge : Option<bool>, + + #[arg(long)] + /// Merging radius (default: determined by experiment) + merge_radius : Option<F>, - #[arg(value_enum, long)]//, value_parser = SpikeMergingMethod::<float>::value_parser())] - /// Final merging strategy - /// - /// Either the string "none", or a radius value for heuristic merging. - /// Only affects FB, FISTA, and PDPS. - final_merging : Option<SpikeMergingMethod<F>>, + #[arg(long)] + /// Interpolate when merging (default : determined by algorithm) + merge_interp : Option<bool>, + + #[arg(long)] + /// Enable final merging (default: determined by algorithm) + final_merging : Option<bool>, + + #[arg(long)] + /// Enable fitness-based merging for relevant FB-type methods. + /// This has worse convergence guarantees that merging based on optimality conditions. + fitness_merging : Option<bool>, #[arg(long, value_names = &["ε", "θ", "p"])] /// Set the tolerance to ε_k = ε/(1+θk)^p @@ -266,11 +281,12 @@ let mut algs : Vec<Named<AlgorithmConfig<float>>> = cli.algorithm .iter() - .map(|alg| alg.to_named( - experiment.algorithm_defaults(*alg) - .unwrap_or_else(|| alg.default_config()) - .cli_override(&cli.algoritm_overrides) - )) + .map(|alg| { + let cfg = alg.default_config() + .cli_override(&experiment.algorithm_overrides(*alg)) + .cli_override(&cli.algoritm_overrides); + alg.to_named(cfg) + }) .collect(); for filename in cli.saved_algorithm.iter() { let f = std::fs::File::open(filename).unwrap();
--- a/src/measures/discrete.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/measures/discrete.rs Thu Jan 23 23:34:05 2025 +0100 @@ -265,44 +265,29 @@ impl<Domain : Clone, F : Float> DiscreteMeasure<Domain, F> { /// Computes `μ1 ← θ * μ1 - ζ * μ2`, pruning entries where both `μ1` (`self`) and `μ2` have - // zero weight. `μ2` will contain copy of pruned original `μ1` without arithmetic performed. - /// **This expects `self` and `μ2` to have matching coordinates in each index**. + // zero weight. `μ2` will contain a pruned copy of pruned original `μ1` without arithmetic + /// performed. **This expects `self` and `μ2` to have matching coordinates in each index**. // `μ2` can be than `self`, but not longer. pub fn pruning_sub(&mut self, θ : F, ζ : F, μ2 : &mut Self) { - let mut μ2_get = 0; - let mut μ2_insert = 0; - self.spikes.retain_mut(|&mut DeltaMeasure{ α : ref mut α_ref, ref x }| { - // Get weight of spike in μ2, zero if out of bounds. - let β = μ2.spikes.get(μ2_get).map_or(F::ZERO, DeltaMeasure::get_mass); - μ2_get += 1; - - if *α_ref == F::ZERO && β == F::ZERO { - // Prune - true + for δ in &self[μ2.len()..] { + μ2.push(DeltaMeasure{ x : δ.x.clone(), α : F::ZERO}); + } + debug_assert_eq!(self.len(), μ2.len()); + let mut dest = 0; + for i in 0..self.len() { + let α = self[i].α; + let α_new = θ * α - ζ * μ2[i].α; + if dest < i { + μ2[dest] = DeltaMeasure{ x : self[i].x.clone(), α }; + self[dest] = DeltaMeasure{ x : self[i].x.clone(), α : α_new }; } else { - // Save self weight - let α = *α_ref; - // Modify self - *α_ref = θ * α - ζ * β; - // Make copy of old self weight in μ2 - let δ = DeltaMeasure{ α, x : x.clone() }; - match μ2.spikes.get_mut(μ2_insert) { - Some(replace) => { - *replace = δ; - }, - None => { - debug_assert_eq!(μ2.len(), μ2_insert); - μ2.spikes.push(δ); - }, - } - μ2_insert += 1; - // Keep - false + μ2[i].α = α; + self[i].α = α_new; } - }); - // Truncate μ2 to same length as self. - μ2.spikes.truncate(μ2_insert); - debug_assert_eq!(μ2.len(), self.len()); + dest += 1; + } + self.spikes.truncate(dest); + μ2.spikes.truncate(dest); } } @@ -327,17 +312,17 @@ self.set_masses(x.iter().map(|&α| F::from_nalgebra_mixed(α))); } - /// Extracts the masses of the spikes as a [`Vec`]. - pub fn masses_vec(&self) -> Vec<F::MixedType> { - self.iter_masses() - .map(|α| α.to_nalgebra_mixed()) - .collect() - } + // /// Extracts the masses of the spikes as a [`Vec`]. + // pub fn masses_vec(&self) -> Vec<F::MixedType> { + // self.iter_masses() + // .map(|α| α.to_nalgebra_mixed()) + // .collect() + // } - /// Sets the masses of the spikes from the values of a [`Vec`]. - pub fn set_masses_vec(&mut self, x : &Vec<F::MixedType>) { - self.set_masses(x.iter().map(|&α| F::from_nalgebra_mixed(α))); - } + // /// Sets the masses of the spikes from the values of a [`Vec`]. + // pub fn set_masses_vec(&mut self, x : &Vec<F::MixedType>) { + // self.set_masses(x.iter().map(|&α| F::from_nalgebra_mixed(α))); + // } } // impl<Domain, F :Num> Index<usize> for DiscreteMeasure<Domain, F> {
--- a/src/measures/merging.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/measures/merging.rs Thu Jan 23 23:34:05 2025 +0100 @@ -19,63 +19,24 @@ /// Spike merging heuristic selection #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] #[allow(dead_code)] -pub enum SpikeMergingMethod<F> { - /// Try to merge spikes within a given radius of each other, averaging the location - HeuristicRadius(F), - /// Try to merge spikes within a given radius of each other, attempting original locations - HeuristicRadiusNoInterp(F), - /// No merging - None, +pub struct SpikeMergingMethod<F> { + // Merging radius + pub(crate) radius : F, + // Enabled + pub(crate) enabled : bool, + // Interpolate merged points + pub(crate) interp : bool, } -// impl<F : Float> SpikeMergingMethod<F> { -// /// This is for [`clap`] to display command line help. -// pub fn value_parser() -> PossibleValuesParser { -// PossibleValuesParser::new([ -// PossibleValue::new("none").help("No merging"), -// PossibleValue::new("<radius>").help("Heuristic merging within indicated radius") -// ]) -// } -// } - -impl<F : ClapFloat> std::fmt::Display for SpikeMergingMethod<F> { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> Result<(), std::fmt::Error> { - match self { - Self::None => write!(f, "none"), - Self::HeuristicRadius(r) => write!(f, "i:{}", r), - Self::HeuristicRadiusNoInterp(r) => write!(f, "n:{}", r), - } - } -} - -impl<F : ClapFloat> std::str::FromStr for SpikeMergingMethod<F> { - type Err = F::Err; - - fn from_str(s: &str) -> Result<Self, Self::Err> { - if s == "none" { - Ok(Self::None) - } else { - let mut subs = s.split(':'); - match subs.next() { - None => Ok(Self::HeuristicRadius(F::from_str(s)?)), - Some(t) if t == "n" => match subs.next() { - None => Err(core::num::dec2flt::pfe_invalid()), - Some(v) => Ok(Self::HeuristicRadiusNoInterp(F::from_str(v)?)) - }, - Some(t) if t == "i" => match subs.next() { - None => Err(core::num::dec2flt::pfe_invalid()), - Some(v) => Ok(Self::HeuristicRadius(F::from_str(v)?)) - }, - Some(v) => Ok(Self::HeuristicRadius(F::from_str(v)?)) - } - } - } -} #[replace_float_literals(F::cast_from(literal))] impl<F : Float> Default for SpikeMergingMethod<F> { fn default() -> Self { - SpikeMergingMethod::HeuristicRadius(0.02) + SpikeMergingMethod{ + radius : 0.01, + enabled : false, + interp : true, + } } } @@ -90,16 +51,16 @@ /// an arbitrary value. This method will return that value for the *last* accepted merge, or /// [`None`] if no merge was accepted. /// - /// This method is stable with respect to spike locations: on merge, the weight of existing - /// spikes is set to zero, and a new one inserted at the end of the spike vector. + /// This method is stable with respect to spike locations: on merge, the weights of existing + /// removed spikes is set to zero, new ones inserted at the end of the spike vector. + /// They merge may also be performed by increasing the weights of the existing spikes, + /// without inserting new spikes. fn merge_spikes<G>(&mut self, method : SpikeMergingMethod<F>, accept : G) -> usize where G : FnMut(&'_ Self) -> bool { - match method { - SpikeMergingMethod::HeuristicRadius(ρ) => - self.do_merge_spikes_radius(ρ, true, accept), - SpikeMergingMethod::HeuristicRadiusNoInterp(ρ) => - self.do_merge_spikes_radius(ρ, false, accept), - SpikeMergingMethod::None => 0, + if method.enabled { + self.do_merge_spikes_radius(method.radius, method.interp, accept) + } else { + 0 } }
--- a/src/pdps.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/pdps.rs Thu Jan 23 23:34:05 2025 +0100 @@ -80,6 +80,8 @@ L2Squared, L1 }; +use crate::measures::merging::SpikeMergingMethod; + /// Acceleration #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, ValueEnum, Debug)] @@ -135,12 +137,15 @@ #[replace_float_literals(F::cast_from(literal))] impl<F : Float> Default for PDPSConfig<F> { fn default() -> Self { - let τ0 = 0.5; + let τ0 = 5.0; PDPSConfig { τ0, σ0 : 0.99/τ0, acceleration : Acceleration::Partial, - generic : Default::default(), + generic : FBGenericConfig { + merging : SpikeMergingMethod { enabled : true, ..Default::default() }, + .. Default::default() + }, } } } @@ -285,7 +290,9 @@ // Prune and possibly merge spikes if config.merge_now(&state) { - stats.merged += prox_penalty.merge_spikes(&mut μ, &mut τv, &μ_base, τ, ε, config, ®); + stats.merged += prox_penalty.merge_spikes_no_fitness( + &mut μ, &mut τv, &μ_base, None, τ, ε, config, ®, + ); } stats.pruned += prune_with_stats(&mut μ);
--- a/src/prox_penalty.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/prox_penalty.rs Thu Jan 23 23:34:05 2025 +0100 @@ -12,17 +12,12 @@ AlgIteratorIteration, AlgIterator, }; -use crate::measures::{ - RNDM, -}; +use crate::measures::RNDM; use crate::types::{ RefinementSettings, IterInfo, }; -use crate::subproblem::{ - InnerSettings, - InnerMethod, -}; +use crate::subproblem::InnerSettings; use crate::tolerance::Tolerance; use crate::measures::merging::SpikeMergingMethod; use crate::regularisation::RegTerm; @@ -66,8 +61,11 @@ /// Tolerance multiplier for merges pub merge_tolerance_mult : F, - /// Spike merging method after the last step - pub final_merging : SpikeMergingMethod<F>, + /// Merge spikes after last step (even if merging not generally enabled) + pub final_merging : bool, + + /// Use fitness as merging criterion. Implies worse convergence guarantees. + pub fitness_merging : bool, /// Iterations between merging heuristic tries pub merge_every : usize, @@ -86,13 +84,10 @@ max_insertions : 100, //bootstrap_insertions : None, bootstrap_insertions : Some((10, 1)), - inner : InnerSettings { - method : InnerMethod::Default, - .. Default::default() - }, - merging : SpikeMergingMethod::None, - //merging : Default::default(), - final_merging : Default::default(), + inner : Default::default(), + merging : Default::default(), + final_merging : true, + fitness_merging : false, merge_every : 10, merge_tolerance_mult : 2.0, // postprocessing : false, @@ -103,7 +98,12 @@ impl<F : Float> FBGenericConfig<F> { /// Check if merging should be attempted this iteration pub fn merge_now<I : AlgIterator>(&self, state : &AlgIteratorIteration<I>) -> bool { - state.iteration() % self.merge_every == 0 + self.merging.enabled && state.iteration() % self.merge_every == 0 + } + + /// Returns the final merging method + pub fn final_merging_method(&self) -> SpikeMergingMethod<F> { + SpikeMergingMethod{ enabled : self.final_merging, ..self.merging} } } @@ -145,14 +145,45 @@ /// Merge spikes, if possible. + /// + /// Either optimality condition merging or objective value (fitness) merging + /// may be used, the latter only if `fitness` is provided and `config.fitness_merging` + /// is set. fn merge_spikes( &self, μ : &mut RNDM<F, N>, τv : &mut PreadjointCodomain, μ_base : &RNDM<F, N>, + ν_delta: Option<&RNDM<F, N>>, τ : F, ε : F, config : &FBGenericConfig<F>, reg : &Reg, + fitness : Option<impl Fn(&RNDM<F, N>) -> F>, ) -> usize; + + /// Merge spikes, if possible. + /// + /// Unlike [`merge_spikes`], this variant only supports optimality condition based merging + #[inline] + fn merge_spikes_no_fitness( + &self, + μ : &mut RNDM<F, N>, + τv : &mut PreadjointCodomain, + μ_base : &RNDM<F, N>, + ν_delta: Option<&RNDM<F, N>>, + τ : F, + ε : F, + config : &FBGenericConfig<F>, + reg : &Reg, + ) -> usize { + /// This is a hack to create a `None` of same type as a `Some` + // for the `impl Fn` parameter of `merge_spikes`. + #[inline] + fn into_none<T>(_ : Option<T>) -> Option<T>{ + None + } + self.merge_spikes(μ, τv, μ_base, ν_delta, τ, ε, config, reg, + into_none(Some(|_ : &RNDM<F, N>| F::ZERO))) + } }
--- a/src/prox_penalty/radon_squared.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/prox_penalty/radon_squared.rs Thu Jan 23 23:34:05 2025 +0100 @@ -12,7 +12,7 @@ AlgIteratorIteration, AlgIterator }; -use alg_tools::norms::L2; +use alg_tools::norms::{L2, Norm}; use alg_tools::linops::Mapping; use alg_tools::bisection_tree::{ BTFN, @@ -75,10 +75,10 @@ where I : AlgIterator { - assert!(ν_delta.is_none(), "Transport not implemented for Radon-squared prox term"); + let mut y = μ_base.masses_dvector(); - let mut y = μ_base.masses_vec(); - + assert!(μ_base.len() <= μ.len()); + 'i_and_w: for i in 0..=1 { // Optimise weights if μ.len() > 0 { @@ -90,10 +90,12 @@ μ.iter_locations() .map(|ζ| - F::to_nalgebra_mixed(τv.apply(ζ)))); let mut x = μ.masses_dvector(); - // Ugly hack because DVector::push doesn't push but copies. - let yvec = DVector::from_column_slice(y.as_slice()); + 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. - stats.inner_iters += reg.solve_findim_l1squared(&yvec, &g̃, τ, &mut x, ε, config); + // 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); @@ -107,7 +109,8 @@ } // Calculate ‖μ - μ_base‖_ℳ - let n = μ.dist_matching(μ_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 μ, @@ -118,8 +121,6 @@ // Weight is found out by running the finite-dimensional optimisation algorithm // above *μ += DeltaMeasure { x : ξ, α : 0.0 }; - //*μ_base += DeltaMeasure { x : ξ, α : 0.0 }; - y.push(0.0.to_nalgebra_mixed()); stats.inserted += 1; } }; @@ -133,12 +134,20 @@ μ : &mut RNDM<F, N>, τv : &mut BTFN<F, GA, BTA, N>, μ_base : &RNDM<F, N>, + ν_delta: Option<&RNDM<F, N>>, τ : F, ε : F, config : &FBGenericConfig<F>, reg : &Reg, + fitness : Option<impl Fn(&RNDM<F, N>) -> F>, ) -> usize { + if config.fitness_merging { + if let Some(f) = fitness { + return μ.merge_spikes_fitness(config.merging, f, |&v| v) + .1 + } + } μ.merge_spikes(config.merging, |μ_candidate| { // Important: μ_candidate's new points are afterwards, // and do not conflict with μ_base. @@ -147,7 +156,10 @@ // 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 = μ_candidate.sub_matching(μ_base); + let μ_radon = match ν_delta { + None => μ_candidate.sub_matching(μ_base), + Some(ν) => μ_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 Jan 23 23:35:28 2025 +0100 +++ b/src/prox_penalty/wave.rs Thu Jan 23 23:34:05 2025 +0100 @@ -36,7 +36,6 @@ use crate::types::{ IterInfo, }; -use crate::measures::merging::SpikeMergingMethod; use crate::regularisation::RegTerm; use super::{ProxPenalty, FBGenericConfig}; @@ -74,7 +73,6 @@ I : AlgIterator { - // TODO: is this inefficient to do in every iteration? let op𝒟norm = self.opnorm_bound(Radon, Linfinity); // Maximum insertion count and measure difference calculation depend on insertion style. @@ -98,7 +96,7 @@ // problems have not yet been updated to sign change. let à = self.findim_matrix(μ.iter_locations()); let g̃ = DVector::from_iterator(μ.len(), - μ.iter_locations() + μ.iter_locations() .map(|ζ| ω0.apply(ζ) - τv.apply(ζ)) .map(F::to_nalgebra_mixed)); let mut x = μ.masses_dvector(); @@ -127,7 +125,7 @@ // 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 // insertion also seems to improve performance. - let skip_by_rough_check = if let SpikeMergingMethod::None = config.merging { + let skip_by_rough_check = if config.merging.enabled { false } else { count > 0 @@ -168,14 +166,25 @@ μ : &mut RNDM<F, N>, τv : &mut BTFN<F, GA, BTA, N>, μ_base : &RNDM<F, N>, + ν_delta: Option<&RNDM<F, N>>, τ : F, ε : F, config : &FBGenericConfig<F>, reg : &Reg, + fitness : Option<impl Fn(&RNDM<F, N>) -> F>, ) -> usize { + if config.fitness_merging { + if let Some(f) = fitness { + return μ.merge_spikes_fitness(config.merging, f, |&v| v) + .1 + } + } μ.merge_spikes(config.merging, |μ_candidate| { - let mut d = &*τv + self.preapply(μ_candidate.sub_matching(μ_base)); + let mut d = &*τv + self.preapply(match ν_delta { + None => μ_candidate.sub_matching(μ_base), + Some(ν) => μ_candidate.sub_matching(μ_base) - ν, + }); reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, config) }) }
--- a/src/run.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/run.rs Thu Jan 23 23:34:05 2025 +0100 @@ -56,7 +56,7 @@ use crate::kernels::*; use crate::types::*; use crate::measures::*; -use crate::measures::merging::SpikeMerging; +use crate::measures::merging::{SpikeMerging,SpikeMergingMethod}; use crate::forward_model::*; use crate::forward_model::sensor_grid::{ SensorGrid, @@ -95,7 +95,7 @@ pointsource_fw_reg, //WeightOptim, }; -//use crate::subproblem::InnerSettings; +use crate::subproblem::{InnerSettings, InnerMethod}; use crate::seminorms::*; use crate::plot::*; use crate::{AlgorithmOverrides, CommandLineArgs}; @@ -146,6 +146,13 @@ impl<F : ClapFloat> AlgorithmConfig<F> { /// Override supported parameters based on the command line. pub fn cli_override(self, cli : &AlgorithmOverrides<F>) -> Self { + let override_merging = |g : SpikeMergingMethod<F>| { + SpikeMergingMethod { + enabled : cli.merge.unwrap_or(g.enabled), + radius : cli.merge_radius.unwrap_or(g.radius), + interp : cli.merge_interp.unwrap_or(g.interp), + } + }; let override_fb_generic = |g : FBGenericConfig<F>| { FBGenericConfig { bootstrap_insertions : cli.bootstrap_insertions @@ -153,8 +160,9 @@ .map_or(g.bootstrap_insertions, |n| Some((n[0], n[1]))), merge_every : cli.merge_every.unwrap_or(g.merge_every), - merging : cli.merging.clone().unwrap_or(g.merging), - final_merging : cli.final_merging.clone().unwrap_or(g.final_merging), + 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), tolerance: cli.tolerance.as_ref().map(unpack_tolerance).unwrap_or(g.tolerance), .. g } @@ -162,8 +170,8 @@ let override_transport = |g : TransportConfig<F>| { TransportConfig { θ0 : cli.theta0.unwrap_or(g.θ0), - tolerance_ω: cli.transport_tolerance_omega.unwrap_or(g.tolerance_ω), - tolerance_dv: cli.transport_tolerance_dv.unwrap_or(g.tolerance_dv), + 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 } @@ -189,7 +197,7 @@ .. pdps }, prox), FW(fw) => FW(FWConfig { - merging : cli.merging.clone().unwrap_or(fw.merging), + merging : override_merging(fw.merging), tolerance : cli.tolerance.as_ref().map(unpack_tolerance).unwrap_or(fw.tolerance), .. fw }), @@ -282,6 +290,14 @@ /// Returns the algorithm configuration corresponding to the algorithm shorthand pub fn default_config<F : Float>(&self) -> AlgorithmConfig<F> { use DefaultAlgorithm::*; + let radon_insertion = FBGenericConfig { + merging : SpikeMergingMethod{ interp : false, .. Default::default() }, + inner : InnerSettings { + method : InnerMethod::PDPS, // SSN not implemented + .. Default::default() + }, + .. Default::default() + }; match *self { FB => AlgorithmConfig::FB(Default::default(), ProxTerm::Wave), FISTA => AlgorithmConfig::FISTA(Default::default(), ProxTerm::Wave), @@ -297,12 +313,30 @@ // Radon variants - RadonFB => AlgorithmConfig::FB(Default::default(), ProxTerm::RadonSquared), - RadonFISTA => AlgorithmConfig::FISTA(Default::default(), ProxTerm::RadonSquared), - RadonPDPS => AlgorithmConfig::PDPS(Default::default(), ProxTerm::RadonSquared), - RadonSlidingFB => AlgorithmConfig::SlidingFB(Default::default(), ProxTerm::RadonSquared), - RadonSlidingPDPS => AlgorithmConfig::SlidingPDPS(Default::default(), ProxTerm::RadonSquared), - RadonForwardPDPS => AlgorithmConfig::ForwardPDPS(Default::default(), ProxTerm::RadonSquared), + RadonFB => AlgorithmConfig::FB( + FBConfig{ generic : radon_insertion, ..Default::default() }, + ProxTerm::RadonSquared + ), + RadonFISTA => AlgorithmConfig::FISTA( + FBConfig{ generic : radon_insertion, ..Default::default() }, + ProxTerm::RadonSquared + ), + RadonPDPS => AlgorithmConfig::PDPS( + PDPSConfig{ generic : radon_insertion, ..Default::default() }, + ProxTerm::RadonSquared + ), + RadonSlidingFB => AlgorithmConfig::SlidingFB( + SlidingFBConfig{ insertion : radon_insertion, ..Default::default() }, + ProxTerm::RadonSquared + ), + RadonSlidingPDPS => AlgorithmConfig::SlidingPDPS( + SlidingPDPSConfig{ insertion : radon_insertion, ..Default::default() }, + ProxTerm::RadonSquared + ), + RadonForwardPDPS => AlgorithmConfig::ForwardPDPS( + ForwardPDPSConfig{ insertion : radon_insertion, ..Default::default() }, + ProxTerm::RadonSquared + ), } } @@ -340,6 +374,12 @@ Iter, } +impl Default for PlotLevel { + fn default() -> Self { + Self::Data + } +} + type DefaultBT<F, const N : usize> = BT< DynamicDepth, F, @@ -418,7 +458,7 @@ /// Struct for experiment configurations #[derive(Debug, Clone, Serialize)] pub struct ExperimentV2<F, NoiseDistr, S, K, P, const N : usize> -where F : Float, +where F : Float + ClapFloat, [usize; N] : Serialize, NoiseDistr : Distribution<F>, S : Sensor<F, N>, @@ -448,13 +488,14 @@ /// Data term pub dataterm : DataTerm, /// A map of default configurations for algorithms - #[serde(skip)] - pub algorithm_defaults : HashMap<DefaultAlgorithm, AlgorithmConfig<F>>, + pub algorithm_overrides : HashMap<DefaultAlgorithm, AlgorithmOverrides<F>>, + /// Default merge radius + pub default_merge_radius : F, } #[derive(Debug, Clone, Serialize)] pub struct ExperimentBiased<F, NoiseDistr, S, K, P, B, const N : usize> -where F : Float, +where F : Float + ClapFloat, [usize; N] : Serialize, NoiseDistr : Distribution<F>, S : Sensor<F, N>, @@ -477,7 +518,7 @@ algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError; /// Return algorithm default config - fn algorithm_defaults(&self, alg : DefaultAlgorithm) -> Option<AlgorithmConfig<F>>; + fn algorithm_overrides(&self, alg : DefaultAlgorithm) -> AlgorithmOverrides<F>; } /// Helper function to print experiment start message and save setup. @@ -494,8 +535,8 @@ let Named { name : experiment_name, data } = experiment; println!("{}\n{}", - format!("Performing experiment {}…", experiment_name).cyan(), - format!("{:?}", data).bright_black()); + format!("Performing experiment {}…", experiment_name).cyan(), + format!("Experiment settings: {}", serde_json::to_string(&data)?).bright_black()); // Set up output directory let prefix = format!("{}/{}/", cli.outdir, experiment_name); @@ -525,7 +566,7 @@ >; /// Helper function to run all algorithms on an experiment. -fn do_runall<F : Float, Z, const N : usize>( +fn do_runall<F : Float + for<'b> Deserialize<'b>, Z, const N : usize>( experiment_name : &String, prefix : &String, cli : &CommandLineArgs, @@ -547,7 +588,7 @@ let iterator_options = AlgIteratorOptions{ max_iter : cli.max_iter, verbose_iter : cli.verbose_iter - .map_or(Verbose::Logarithmic(10), + .map_or(Verbose::LogarithmicCap{base : 10, cap : 2}, |n| Verbose::Every(n)), quiet : cli.quiet, }; @@ -565,8 +606,8 @@ let running = if !cli.quiet { format!("{}\n{}\n{}\n", format!("Running {} on experiment {}…", alg_name, experiment_name).cyan(), - format!("{:?}", iterator_options).bright_black(), - format!("{:?}", alg).bright_black()) + format!("Iteration settings: {}", serde_json::to_string(&iterator_options)?).bright_black(), + format!("Algorithm settings: {}", serde_json::to_string(&alg)?).bright_black()) } else { "".to_string() }; @@ -614,16 +655,17 @@ save_extra(mkname(""), z)?; //logger.write_csv(mkname("log.txt"))?; logs.push((mkname("log.txt"), logger)); - } + } - save_logs(logs) + save_logs(logs, format!("{prefix}valuerange.json"), cli.load_valuerange) } #[replace_float_literals(F::cast_from(literal))] impl<F, NoiseDistr, S, K, P, /*PreadjointCodomain, */ const N : usize> RunnableExperiment<F> for Named<ExperimentV2<F, NoiseDistr, S, K, P, N>> where - F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F>, + F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F> + + Default + for<'b> Deserialize<'b>, [usize; N] : Serialize, S : Sensor<F, N> + Copy + Serialize + std::fmt::Debug, P : Spread<F, N> + Copy + Serialize + std::fmt::Debug, @@ -655,8 +697,11 @@ // RadonSquared : ProxPenalty<F, PreadjointCodomain, NonnegRadonRegTerm<F>, N>, { - fn algorithm_defaults(&self, alg : DefaultAlgorithm) -> Option<AlgorithmConfig<F>> { - self.data.algorithm_defaults.get(&alg).cloned() + fn algorithm_overrides(&self, alg : DefaultAlgorithm) -> AlgorithmOverrides<F> { + AlgorithmOverrides { + merge_radius : Some(self.data.default_merge_radius), + .. self.data.algorithm_overrides.get(&alg).cloned().unwrap_or(Default::default()) + } } fn runall(&self, cli : &CommandLineArgs, @@ -887,7 +932,8 @@ impl<F, NoiseDistr, S, K, P, B, /*PreadjointCodomain,*/ const N : usize> RunnableExperiment<F> for Named<ExperimentBiased<F, NoiseDistr, S, K, P, B, N>> where - F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F>, + F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F> + + Default + for<'b> Deserialize<'b>, [usize; N] : Serialize, S : Sensor<F, N> + Copy + Serialize + std::fmt::Debug, P : Spread<F, N> + Copy + Serialize + std::fmt::Debug, @@ -920,8 +966,11 @@ // RadonSquared : ProxPenalty<F, PreadjointCodomain, NonnegRadonRegTerm<F>, N>, { - fn algorithm_defaults(&self, alg : DefaultAlgorithm) -> Option<AlgorithmConfig<F>> { - self.data.base.algorithm_defaults.get(&alg).cloned() + fn algorithm_overrides(&self, alg : DefaultAlgorithm) -> AlgorithmOverrides<F> { + AlgorithmOverrides { + merge_radius : Some(self.data.base.default_merge_radius), + .. self.data.base.algorithm_overrides.get(&alg).cloned().unwrap_or(Default::default()) + } } fn runall(&self, cli : &CommandLineArgs, @@ -1077,16 +1126,31 @@ } } +#[derive(Copy, Clone, Debug, Serialize, Deserialize)] +struct ValueRange<F : Float> { + ini : F, + min : F, +} + +impl<F : Float> ValueRange<F> { + fn expand_with(self, other : Self) -> Self { + ValueRange { + ini : self.ini.max(other.ini), + min : self.min.min(other.min), + } + } +} /// Calculative minimum and maximum values of all the `logs`, and save them into /// corresponding file names given as the first elements of the tuples in the vectors. -fn save_logs<F : Float, const N : usize>( - logs : Vec<(String, Logger<Timed<IterInfo<F, N>>>)> +fn save_logs<F : Float + for<'b> Deserialize<'b>, const N : usize>( + logs : Vec<(String, Logger<Timed<IterInfo<F, N>>>)>, + valuerange_file : String, + load_valuerange : bool, ) -> DynError { // Process logs for relative values println!("{}", "Processing logs…"); - // Find minimum value and initial value within a single log let proc_single_log = |log : &Logger<Timed<IterInfo<F, N>>>| { let d = log.data(); @@ -1094,15 +1158,22 @@ .map(|i| i.data.value) .reduce(NumTraitsFloat::min); d.first() - .map(|i| i.data.value) - .zip(mi) + .map(|i| i.data.value) + .zip(mi) + .map(|(ini, min)| ValueRange{ ini, min }) }; // Find minimum and maximum value over all logs - let (v_ini, v_min) = logs.iter() - .filter_map(|&(_, ref log)| proc_single_log(log)) - .reduce(|(i1, m1), (i2, m2)| (i1.max(i2), m1.min(m2))) - .ok_or(anyhow!("No algorithms found"))?; + let mut v = logs.iter() + .filter_map(|&(_, ref log)| proc_single_log(log)) + .reduce(|v1, v2| v1.expand_with(v2)) + .ok_or(anyhow!("No algorithms found"))?; + + // Load existing range + if load_valuerange && std::fs::metadata(&valuerange_file).is_ok() { + let data = std::fs::read_to_string(&valuerange_file)?; + v = v.expand_with(serde_json::from_str(&data)?); + } let logmap = |Timed { cpu_time, iter, data }| { let IterInfo { @@ -1128,7 +1199,7 @@ // }, // _ => value, // }; - let relative_value = (value - v_min)/(v_ini - v_min); + let relative_value = (value - v.min)/(v.ini - v.min); CSVLog { iter, value, @@ -1145,6 +1216,8 @@ println!("{}", "Saving logs …".green()); + serde_json::to_writer_pretty(std::fs::File::create(&valuerange_file)?, &v)?; + for (name, logger) in logs { logger.map(logmap).write_csv(name)?; }
--- a/src/sliding_fb.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/sliding_fb.rs Thu Jan 23 23:34:05 2025 +0100 @@ -35,7 +35,7 @@ use crate::regularisation::SlidingRegTerm; use crate::dataterm::{ L2Squared, - //DataTerm, + DataTerm, calculate_residual, calculate_residual2, }; @@ -49,10 +49,10 @@ pub θ0 : F, /// Factor in $(0, 1)$ for decreasing transport to adapt to tolerance. pub adaptation : F, - /// Transport tolerance wrt. ω - pub tolerance_ω : F, - /// Transport tolerance wrt. ∇v - pub tolerance_dv : F, + /// A priori transport tolerance multiplier (C_pri) + pub tolerance_mult_pri : F, + /// A posteriori transport tolerance multiplier (C_pos) + pub tolerance_mult_pos : F, } #[replace_float_literals(F::cast_from(literal))] @@ -61,8 +61,8 @@ pub fn check(&self) { assert!(self.θ0 > 0.0); assert!(0.0 < self.adaptation && self.adaptation < 1.0); - assert!(self.tolerance_dv > 0.0); - assert!(self.tolerance_ω > 0.0); + assert!(self.tolerance_mult_pri > 0.0); + assert!(self.tolerance_mult_pos > 0.0); } } @@ -70,10 +70,10 @@ impl<F : Float> Default for TransportConfig<F> { fn default() -> Self { TransportConfig { - θ0 : 0.01, + θ0 : 0.4, adaptation : 0.9, - tolerance_ω : 1000.0, // TODO: no idea what this should be - tolerance_dv : 1000.0, // TODO: no idea what this should be + tolerance_mult_pos : 100.0, + tolerance_mult_pri : 1000.0, } } } @@ -181,20 +181,26 @@ FullyAdaptive{ l : ref mut adaptive_ℓ_v, ref mut max_transport, g : ref calculate_θ } => { *max_transport = max_transport.max(γ1.norm(Radon)); let mut θ = calculate_θ(*adaptive_ℓ_v, *max_transport); - loop { - let θτ = τ * θ; + // 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); - ρ.x = δ.x - &dv_x * (ρ.α.signum() * θτ); - // Estimate Lipschitz factor of ∇v - let this_ℓ_v = (dv_x - v.differential(&ρ.x)).norm2(); - *adaptive_ℓ_v = adaptive_ℓ_v.max(this_ℓ_v); + let g = &dv_x * (ρ.α.signum() * θ * τ); + ρ.x = δ.x - g; + let n = g.norm2(); + if n >= F::EPSILON { + // Estimate Lipschitz factor of ∇v + let this_ℓ_v = (dv_x - v.differential(&ρ.x)).norm2() / n; + *adaptive_ℓ_v = adaptive_ℓ_v.max(this_ℓ_v); + θ = calculate_θ(*adaptive_ℓ_v, *max_transport); + changes = true + } } - let new_θ = calculate_θ(*adaptive_ℓ_v / tconfig.adaptation, *max_transport); - if new_θ <= θ { + if !changes { break } - θ = new_θ; } } } @@ -209,7 +215,7 @@ if n <= 0.0 || nr <= 0.0 { break } - let reduction_needed = nr - (ε * tconfig.tolerance_dv / n); + let reduction_needed = nr - (ε * tconfig.tolerance_mult_pri / n); if reduction_needed <= 0.0 { break } @@ -239,7 +245,7 @@ if n <= 0.0 || nr <= 0.0 { break } - let reduction_needed = nr - (ε * tconfig.tolerance_dv / n); + let reduction_needed = nr - (ε * tconfig.tolerance_mult_pri / n); if reduction_needed <= 0.0 { break } @@ -288,6 +294,7 @@ μ : &mut RNDM<F, N>, μ_base_minus_γ0 : &mut RNDM<F, N>, μ_base_masses : &Vec<F>, + extra : Option<F>, ε : F, tconfig : &TransportConfig<F> ) -> bool @@ -308,8 +315,8 @@ // 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); - let t = ε * tconfig.tolerance_ω; + let nΔ = μ_base_minus_γ0.norm(Radon) + μ.dist_matching(&γ1) + extra.unwrap_or(0.0); + let t = ε * tconfig.tolerance_mult_pos; if nγ*nΔ > t { // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1, // this will guarantee that eventually ‖γ‖ decreases sufficiently that we @@ -379,9 +386,9 @@ // We only estimate w (the uniform Lipschitz for of v), if we also estimate ℓ_v // (the uniform Lipschitz factor of ∇v). // We assume that the residual is decreasing. - Some(ℓ_v0) => TransportStepLength::Fixed(calculate_θ(ℓ_v0 * residual.norm2(), 0.0)), + Some(ℓ_v0) => TransportStepLength::Fixed(calculate_θ(ℓ_v0 * b.norm2(), 0.0)), None => TransportStepLength::FullyAdaptive { - l : 0.0, + l : 10.0 * F::EPSILON, // Start with something very small to estimate differentials max_transport : 0.0, g : calculate_θ }, @@ -415,7 +422,7 @@ // 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, τv̆) = 'adapt_transport: loop { + let (maybe_d, _within_tolerances, mut τv̆) = 'adapt_transport: loop { // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) let residual_μ̆ = calculate_residual2(&γ1, &μ_base_minus_γ0, opA, b); let mut τv̆ = opA.preadjoint().apply(residual_μ̆ * τ); @@ -430,6 +437,7 @@ // 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̆) @@ -448,20 +456,16 @@ (a + μ.dist_matching(&γ1), b + γ1.norm(Radon)) }); - // // Merge spikes. - // // This expects the prune below to prune γ. - // // TODO: This may not work correctly in all cases. - // let ins = &config.insertion; - // if ins.merge_now(&state) { - // if let SpikeMergingMethod::None = ins.merging { - // } else { - // stats.merged += μ.merge_spikes(ins.merging, |μ_candidate| { - // let ν = μ_candidate.sub_matching(&γ1)-&μ_base_minus_γ0; - // let mut d = &τv̆ + op𝒟.preapply(ν); - // reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, ins) - // }); - // } - // } + // 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) { + stats.merged += prox_penalty.merge_spikes( + &mut μ, &mut τv̆, &γ1, Some(&μ_base_minus_γ0), τ, ε, ins, ®, + Some(|μ̃ : &RNDM<F, N>| L2Squared.calculate_fit_op(μ̃, opA, b)), + ); + } // Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the // latter needs to be pruned when μ is.
--- a/src/sliding_pdps.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/sliding_pdps.rs Thu Jan 23 23:34:05 2025 +0100 @@ -12,7 +12,7 @@ use alg_tools::iterate::AlgIteratorFactory; use alg_tools::euclidean::Euclidean; use alg_tools::mapping::{Mapping, DifferentiableRealMapping, Instance}; -use alg_tools::norms::Norm; +use alg_tools::norms::{Norm, Dist}; use alg_tools::direct_product::Pair; use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::linops::{ @@ -45,7 +45,11 @@ initial_transport, aposteriori_transport, }; -use crate::dataterm::{calculate_residual, calculate_residual2}; +use crate::dataterm::{ + calculate_residual2, + calculate_residual, +}; + /// Settings for [`pointsource_sliding_pdps_pair`]. #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] @@ -66,12 +70,11 @@ #[replace_float_literals(F::cast_from(literal))] impl<F : Float> Default for SlidingPDPSConfig<F> { fn default() -> Self { - let τ0 = 0.99; SlidingPDPSConfig { - τ0, - σd0 : 0.1, + τ0 : 0.99, + σd0 : 0.05, σp0 : 0.99, - transport : Default::default(), + transport : TransportConfig { θ0 : 0.1, ..Default::default()}, insertion : Default::default() } } @@ -134,7 +137,7 @@ for<'b> KOpZ::Adjoint<'b> : GEMV<F, Y>, Y : AXPY<F> + Euclidean<F, Output=Y> + Clone + ClosedAdd, for<'b> &'b Y : Instance<Y>, - Z : AXPY<F, Owned=Z> + Euclidean<F, Output=Z> + Clone + Norm<F, L2>, + Z : AXPY<F, Owned=Z> + Euclidean<F, Output=Z> + Clone + Norm<F, L2> + Dist<F, L2>, for<'b> &'b Z : Instance<Z>, R : Prox<Z, Codomain=F>, H : Conjugable<Y, F, Codomain=F>, @@ -202,7 +205,7 @@ g : calculate_θ }, None => TransportStepLength::FullyAdaptive{ - l : 0.0, + l : F::EPSILON, max_transport : 0.0, g : calculate_θ }, @@ -234,7 +237,6 @@ // Calculate initial transport let Pair(v, _) = opA.preadjoint().apply(&residual); //opKμ.preadjoint().apply_add(&mut v, y); - let z_base = z.clone(); // We want to proceed as in Example 4.12 but with v and v̆ as in §5. // With A(ν, z) = A_μ ν + A_z z, following Example 5.1, we have // P_ℳ[F'(ν, z) + Ξ(ν, z, y)]= A_ν^*[A_ν ν + A_z z] + K_μ ν = A_ν^*A(ν, z) + K_μ ν, @@ -250,28 +252,33 @@ // 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, Pair(τv̆, τz̆)) = 'adapt_transport: loop { + let (maybe_d, _within_tolerances, mut τv̆, z_new) = 'adapt_transport: loop { // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) let residual_μ̆ = calculate_residual2(Pair(&γ1, &z), Pair(&μ_base_minus_γ0, &zero_z), opA, b); - let mut τv̆z = opA.preadjoint().apply(residual_μ̆ * τ); + let Pair(mut τv̆, τz̆) = opA.preadjoint().apply(residual_μ̆ * τ); // 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̆z.0, &γ1, Some(&μ_base_minus_γ0), + &mut μ, &mut τv̆, &γ1, Some(&μ_base_minus_γ0), τ, ε, &config.insertion, ®, &state, &mut stats, ); + // Do z variable primal update here to able to estimate B_{v̆^k-v^{k+1}} + let mut z_new = τz̆; + opKz.adjoint().gemv(&mut z_new, -σ_p, &y, -σ_p/τ); + z_new = fnR.prox(σ_p, z_new + &z); + // A posteriori transport adaptation. - // TODO: this does not properly treat v^{k+1} - v̆^k that depends on z^{k+1}! if aposteriori_transport( &mut γ1, &mut μ, &mut μ_base_minus_γ0, &μ_base_masses, + Some(z_new.dist(&z, L2)), ε, &config.transport ) { - break 'adapt_transport (maybe_d, within_tolerances, τv̆z) + break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new) } }; @@ -287,20 +294,16 @@ (a + μ.dist_matching(&γ1), b + γ1.norm(Radon)) }); - // // Merge spikes. - // // This expects the prune below to prune γ. - // // TODO: This may not work correctly in all cases. - // let ins = &config.insertion; - // if ins.merge_now(&state) { - // if let SpikeMergingMethod::None = ins.merging { - // } else { - // stats.merged += μ.merge_spikes(ins.merging, |μ_candidate| { - // let ν = μ_candidate.sub_matching(&γ1)-&μ_base_minus_γ0; - // let mut d = &τv̆ + op𝒟.preapply(ν); - // reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, ins) - // }); - // } - // } + // 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) { + stats.merged += prox_penalty.merge_spikes_no_fitness( + &mut μ, &mut τv̆, &γ1, Some(&μ_base_minus_γ0), τ, ε, ins, ®, + //Some(|μ̃ : &RNDM<F, N>| 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. @@ -313,16 +316,13 @@ μ = μ_new; } - // Do z variable primal update - z.axpy(-σ_p/τ, τz̆, 1.0); // TODO: simplify nasty factors - opKz.adjoint().gemv(&mut z, -σ_p, &y, 1.0); - z = fnR.prox(σ_p, z); // Do dual update // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] - opKz.gemv(&mut y, σ_d*(1.0 + ω), &z, 1.0); + opKz.gemv(&mut y, σ_d*(1.0 + ω), &z_new, 1.0); // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b - opKz.gemv(&mut y, -σ_d*ω, z_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b + opKz.gemv(&mut y, -σ_d*ω, z, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b y = starH.prox(σ_d, y); + z = z_new; // Update residual residual = calculate_residual(Pair(&μ, &z), opA, b); @@ -349,7 +349,7 @@ + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z)) }; - μ.merge_spikes_fitness(config.insertion.merging, fit, |&v| v); + μ.merge_spikes_fitness(config.insertion.final_merging_method(), fit, |&v| v); μ.prune(); Pair(μ, z) }
--- a/src/subproblem.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/subproblem.rs Thu Jan 23 23:34:05 2025 +0100 @@ -18,7 +18,8 @@ pub mod l1squared_nonneg; -/// Method for solving finite-dimensional subproblems +/// Method for solving finite-dimensional subproblems. +/// Not all outer methods necessarily support all options. #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] #[allow(dead_code)] pub enum InnerMethod { @@ -28,8 +29,6 @@ SSN, /// PDPS PDPS, - /// Problem-specific choice (in practise FB or PDPS, whichever is implemented) - Default, } /// Settings for the solution of finite-dimensional subproblems @@ -65,7 +64,7 @@ quiet : true, .. Default::default() }, - method : InnerMethod::Default, + method : InnerMethod::SSN, tolerance_mult : 0.01, } }
--- a/src/subproblem/l1squared_nonneg.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/subproblem/l1squared_nonneg.rs Thu Jan 23 23:34:05 2025 +0100 @@ -405,14 +405,14 @@ I : AlgIteratorFactory<F> { match inner.method { - InnerMethod::PDPS | InnerMethod::Default => { + InnerMethod::PDPS => { let inner_θ = 0.001; // Estimate of ‖K‖ for K=θ\Id. 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 /*| InnerMethod::Default*/ => { + 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. @@ -420,6 +420,6 @@ let inner_θ = inner_τ; l1squared_nonneg_pp(y, g, λ, β, x, inner_τ, inner_θ, iterator) }, - _ => unimplemented!("{:?}", inner.method), + other => unimplemented!("${other:?} is unimplemented"), } }
--- a/src/subproblem/l1squared_unconstrained.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/subproblem/l1squared_unconstrained.rs Thu Jan 23 23:34:05 2025 +0100 @@ -188,8 +188,8 @@ let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest); match inner.method { - InnerMethod::PDPS | InnerMethod::Default => + InnerMethod::PDPS => l1squared_unconstrained_pdps(y, g, λ, β, x, inner_τ, inner_σ, iterator), - _ => unimplemented!(), + other => unimplemented!("${other:?} is unimplemented"), } }
--- a/src/subproblem/nonneg.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/subproblem/nonneg.rs Thu Jan 23 23:34:05 2025 +0100 @@ -326,7 +326,7 @@ let inner_τ = inner.fb_τ0 / mA_normest; match inner.method { - InnerMethod::FB | InnerMethod::Default => + InnerMethod::FB => quadratic_nonneg_fb(mA, g, λ, x, inner_τ, iterator), InnerMethod::SSN => quadratic_nonneg_ssn(mA, g, λ, x, inner_τ, iterator).unwrap_or_else(|e| { @@ -334,6 +334,6 @@ let ins = InnerSettings::<F>::default(); quadratic_nonneg_fb(mA, g, λ, x, inner_τ, ins.iterator_options) }), - InnerMethod::PDPS => unimplemented!(), + other => unimplemented!("${other:?} is unimplemented"), } }
--- a/src/subproblem/unconstrained.rs Thu Jan 23 23:35:28 2025 +0100 +++ b/src/subproblem/unconstrained.rs Thu Jan 23 23:34:05 2025 +0100 @@ -278,7 +278,7 @@ let inner_τ = inner.fb_τ0 / mA_normest; match inner.method { - InnerMethod::FB | InnerMethod::Default => + InnerMethod::FB => quadratic_unconstrained_fb(mA, g, λ, x, inner_τ, iterator), InnerMethod::SSN => quadratic_unconstrained_ssn(mA, g, λ, x, inner_τ, iterator).unwrap_or_else(|e| { @@ -286,6 +286,6 @@ let ins = InnerSettings::<F>::default(); quadratic_unconstrained_fb(mA, g, λ, x, inner_τ, ins.iterator_options) }), - InnerMethod::PDPS => unimplemented!(), + other => unimplemented!("${other:?} is unimplemented"), } }
