$$ @@ -125,24 +125,7 @@ x: &mut DVector, mA_normest: F, ε: F, - config: &FBGenericConfig, - ) -> usize; - - /// Approximately solve the problem - ///

$$ - /// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + τ G(x) - /// $$

- /// for $G$ depending on the trait implementation. - /// - /// Returns the number of iterations taken. - fn solve_findim_l1squared( - &self, - y: &DVector, - g: &DVector, - τ: F, - x: &mut DVector, - ε: F, - config: &FBGenericConfig, + config: &InsertionConfig, ) -> usize; /// Find a point where `d` may violate the tolerance `ε`. @@ -154,85 +137,30 @@ /// 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( + fn find_tolerance_violation( &self, - d: &mut BTFN, + d: &mut M, τ: F, ε: F, skip_by_rough_check: bool, - config: &FBGenericConfig, - ) -> Option<(Loc, F, bool)> + config: &InsertionConfig, + ) -> Option<(Domain, F, bool)> where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, - { - 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( - &self, - d: &mut BTFN, - τ: F, - ε: F, - skip_by_rough_check: bool, - config: &FBGenericConfig, - slack: F, - ) -> Option<(Loc, F, bool)> - where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>; + M: MinMaxMapping; /// Verify that `d` is in bounds `ε` for a merge candidate `μ` /// /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser. - fn verify_merge_candidate( + fn verify_merge_candidate( &self, - d: &mut BTFN, - μ: &RNDM, + d: &mut M, + μ: &DiscreteMeasure, τ: F, ε: F, - config: &FBGenericConfig, + config: &InsertionConfig, ) -> bool where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>; - - /// Verify that `d` is in bounds `ε` for a merge candidate `μ` - /// - /// This version is s used for Radon-norm squared proximal term in - /// [`crate::prox_penalty::radon_squared`]. - /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have - /// same coordinates at same agreeing indices. - /// - /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser. - fn verify_merge_candidate_radonsq( - &self, - d: &mut BTFN, - μ: &RNDM, - τ: F, - ε: F, - config: &FBGenericConfig, - radon_μ: &RNDM, - ) -> bool - where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>; + M: MinMaxMapping; /// TODO: document this fn target_bounds(&self, τ: F, ε: F) -> Option>; @@ -243,33 +171,76 @@ fn tolerance_scaling(&self) -> F; } +/// Regularisation term that can be used with [`crate::prox_penalty::radon_squared::RadonSquared`] +/// proximal penalty. +pub trait RadonSquaredRegTerm: RegTerm +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( + &self, + μ: &mut DiscreteMeasure, + τv: &mut M, + τ: F, + ε: F, + tolerance_violation: Option<(Domain, F, bool)>, + config: &InsertionConfig, + stats: &mut IterInfo, + ) where + M: Mapping; + + /// Verify that `d` is in bounds `ε` for a merge candidate `μ` + /// + /// This version is s used for Radon-norm squared proximal term in + /// [`crate::prox_penalty::radon_squared`]. + /// The [measures][crate::measures::DiscreteMeasure] `μ` and `radon_μ` are supposed to have + /// same coordinates at same agreeing indices. + /// + /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser. + fn verify_merge_candidate_radonsq( + &self, + d: &mut M, + μ: &DiscreteMeasure, + τ: F, + ε: F, + config: &InsertionConfig, + radon_μ: &DiscreteMeasure, + ) -> bool + where + M: MinMaxMapping; +} + /// Abstraction of regularisation terms for [`pointsource_sliding_fb_reg`]. -pub trait SlidingRegTerm: RegTerm { +pub trait SlidingRegTerm: RegTerm +where + Domain: Space + Clone, + F: Float + ToNalgebraRealField, +{ /// Calculate $τ[w(z) - w(y)]$ for some w in the subdifferential of the regularisation /// term, such that $-ε ≤ τw - d ≤ ε$. - fn goodness( + fn goodness( &self, - d: &mut BTFN, - μ: &RNDM, - y: &Loc, - z: &Loc, + d: &mut M, + μ: &DiscreteMeasure, + y: &Domain, + z: &Domain, τ: F, ε: F, - config: &FBGenericConfig, + config: &InsertionConfig, ) -> F where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>; + M: MinMaxMapping; /// Convert bound on the regulariser to a bond on the Radon norm fn radon_norm_bound(&self, b: F) -> F; } #[replace_float_literals(F::cast_from(literal))] -impl RegTerm for NonnegRadonRegTerm -where - Cube: P2Minimise, F>, +impl RegTerm, F> + for NonnegRadonRegTerm { fn solve_findim( &self, @@ -279,45 +250,28 @@ x: &mut DVector, mA_normest: F, ε: F, - config: &FBGenericConfig, + config: &InsertionConfig, ) -> usize { let inner_tolerance = ε * config.inner.tolerance_mult; let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); quadratic_nonneg(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) } - fn solve_findim_l1squared( + #[inline] + fn find_tolerance_violation( &self, - y: &DVector, - g: &DVector, - τ: F, - x: &mut DVector, - ε: F, - config: &FBGenericConfig, - ) -> 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( - &self, - d: &mut BTFN, + d: &mut M, τ: F, ε: F, skip_by_rough_check: bool, - config: &FBGenericConfig, - slack: F, - ) -> Option<(Loc, F, bool)> + config: &InsertionConfig, + ) -> Option<(Loc, F, bool)> where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, 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; // If preliminary check indicates that we are in bounds, and if it otherwise matches @@ -326,27 +280,26 @@ None } else { // If the rough check didn't indicate no insertion needed, find minimising point. - d.minimise_below( + let res = d.minimise_below( minimise_below, refinement_tolerance, config.refinement.max_steps, - ) - .map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above)) + ); + + res.map(|(ξ, v_ξ)| (ξ, v_ξ, v_ξ >= keep_above)) } } - fn verify_merge_candidate( + fn verify_merge_candidate( &self, - d: &mut BTFN, - μ: &RNDM, + d: &mut M, + μ: &RNDM, τ: F, ε: F, - config: &FBGenericConfig, + config: &InsertionConfig, ) -> bool where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, F>, { let τα = τ * self.α(); let refinement_tolerance = ε * config.refinement.tolerance_mult; @@ -367,19 +320,91 @@ )); } - fn verify_merge_candidate_radonsq( + fn target_bounds(&self, τ: F, ε: F) -> Option> { + let τα = τ * self.α(); + Some(Bounds(τα - ε, τα + ε)) + } + + fn tolerance_scaling(&self) -> F { + self.α() + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl RadonSquaredRegTerm, F> + for NonnegRadonRegTerm +{ + fn solve_oc_radonsq( &self, - d: &mut BTFN, - μ: &RNDM, + μ: &mut DiscreteMeasure, F>, + τv: &mut M, τ: F, ε: F, - config: &FBGenericConfig, - radon_μ: &RNDM, + tolerance_violation: Option<(Loc, F, bool)>, + config: &InsertionConfig, + stats: &mut IterInfo, + ) where + M: Mapping, 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( + &self, + d: &mut M, + μ: &RNDM, + τ: F, + ε: F, + config: &InsertionConfig, + radon_μ: &RNDM, ) -> bool where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, F>, { let τα = τ * self.α(); let refinement_tolerance = ε * config.refinement.tolerance_mult; @@ -414,36 +439,24 @@ ) }; } - - fn target_bounds(&self, τ: F, ε: F) -> Option> { - let τα = τ * self.α(); - Some(Bounds(τα - ε, τα + ε)) - } - - fn tolerance_scaling(&self) -> F { - self.α() - } } #[replace_float_literals(F::cast_from(literal))] -impl SlidingRegTerm for NonnegRadonRegTerm -where - Cube: P2Minimise, F>, +impl SlidingRegTerm, F> + for NonnegRadonRegTerm { - fn goodness( + fn goodness( &self, - d: &mut BTFN, - _μ: &RNDM, - y: &Loc, - z: &Loc, + d: &mut M, + _μ: &RNDM, + y: &Loc, + z: &Loc, τ: F, ε: F, - _config: &FBGenericConfig, + _config: &InsertionConfig, ) -> F where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, F>, { let w = |x| 1.0.min((ε + d.apply(x)) / (τ * self.α())); w(z) - w(y) @@ -455,10 +468,7 @@ } #[replace_float_literals(F::cast_from(literal))] -impl RegTerm for RadonRegTerm -where - Cube: P2Minimise, F>, -{ +impl RegTerm, F> for RadonRegTerm { fn solve_findim( &self, mA: &DMatrix, @@ -467,46 +477,29 @@ x: &mut DVector, mA_normest: F, ε: F, - config: &FBGenericConfig, + config: &InsertionConfig, ) -> usize { let inner_tolerance = ε * config.inner.tolerance_mult; let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); quadratic_unconstrained(mA, g, τ * self.α(), x, mA_normest, &config.inner, inner_it) } - fn solve_findim_l1squared( + fn find_tolerance_violation( &self, - y: &DVector, - g: &DVector, - τ: F, - x: &mut DVector, - ε: F, - config: &FBGenericConfig, - ) -> 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( - &self, - d: &mut BTFN, + d: &mut M, τ: F, ε: F, skip_by_rough_check: bool, - config: &FBGenericConfig, - slack: F, - ) -> Option<(Loc, F, bool)> + config: &InsertionConfig, + ) -> Option<(Loc, F, bool)> where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, 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 @@ -541,18 +534,16 @@ } } - fn verify_merge_candidate( + fn verify_merge_candidate( &self, - d: &mut BTFN, - μ: &RNDM, + d: &mut M, + μ: &RNDM, τ: F, ε: F, - config: &FBGenericConfig, + config: &InsertionConfig, ) -> bool where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, F>, { let τα = τ * self.α(); let refinement_tolerance = ε * config.refinement.tolerance_mult; @@ -585,19 +576,92 @@ )); } - fn verify_merge_candidate_radonsq( + fn target_bounds(&self, τ: F, ε: F) -> Option> { + let τα = τ * self.α(); + Some(Bounds(-τα - ε, τα + ε)) + } + + fn tolerance_scaling(&self) -> F { + self.α() + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl RadonSquaredRegTerm, F> + for RadonRegTerm +{ + fn solve_oc_radonsq( &self, - d: &mut BTFN, - μ: &RNDM, + μ: &mut DiscreteMeasure, F>, + τv: &mut M, τ: F, ε: F, - config: &FBGenericConfig, - radon_μ: &RNDM, + tolerance_violation: Option<(Loc, F, bool)>, + config: &InsertionConfig, + stats: &mut IterInfo, + ) where + M: Mapping, 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( + &self, + d: &mut M, + μ: &RNDM, + τ: F, + ε: F, + config: &InsertionConfig, + radon_μ: &RNDM, ) -> bool where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, F>, { let τα = τ * self.α(); let refinement_tolerance = ε * config.refinement.tolerance_mult; @@ -638,36 +702,24 @@ ) }; } - - fn target_bounds(&self, τ: F, ε: F) -> Option> { - let τα = τ * self.α(); - Some(Bounds(-τα - ε, τα + ε)) - } - - fn tolerance_scaling(&self) -> F { - self.α() - } } #[replace_float_literals(F::cast_from(literal))] -impl SlidingRegTerm for RadonRegTerm -where - Cube: P2Minimise, F>, +impl SlidingRegTerm, F> + for RadonRegTerm { - fn goodness( + fn goodness( &self, - d: &mut BTFN, - _μ: &RNDM, - y: &Loc, - z: &Loc, + d: &mut M, + _μ: &RNDM, + y: &Loc, + z: &Loc, τ: F, ε: F, - _config: &FBGenericConfig, + _config: &InsertionConfig, ) -> F where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + M: MinMaxMapping, F>, { let α = self.α(); let w = |x| {