pointsource_algs: comparison src/regularisation.rs

-:4f468d35fa29
+:7a8a55fd41c0
 mA_normest: F,
 ε: F,
 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
 /// are in bounds. `ε` is the current main tolerance and `τ` a scaling factor for the
 /// regulariser.
 ε: F,
 skip_by_rough_check: bool,
 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 `μ`
 ///
 /// `ε` is the current main tolerance and `τ` a scaling factor for the regulariser.
 ε: F,
 config: &InsertionConfig<F>,
 ) -> 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;
+}
+/// 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
 /// [`crate::prox_penalty::radon_squared`].
 config: &InsertionConfig<F>,
 radon_μ: &DiscreteMeasure<Domain, F>,
 ) -> 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`].
 pub trait SlidingRegTerm<Domain, F = f64>: RegTerm<Domain, F>
 where
 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(
-&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 keep_above = -τα - ε;
-let minimise_below = -τα - slack - ε * config.insertion_cutoff_factor;
+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 {
 None
 || d.has_lower_bound(
 keep_above,
 refinement_tolerance,
 config.refinement.max_steps,
 ));
+}
+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,
 refinement_tolerance,
 config.refinement.max_steps,
 )
 };
 }
-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> SlidingRegTerm<Loc<N, F>, F>
 for NonnegRadonRegTerm<F>
 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<M>(
-&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>(
 &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_below = τα + ε;
-let keep_above = -(τα + slack) - ε;
+let keep_above = -τα - ε;
-let maximise_above = τα + slack + ε * config.insertion_cutoff_factor;
+let maximise_above = τα + ε * config.insertion_cutoff_factor;
-let minimise_below = -(τα + slack) - ε * 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
 // the insertion strategy, skip insertion.
 if skip_by_rough_check && Bounds(keep_above, keep_below).superset(&d.bounds()) {
 || d.has_lower_bound(
 keep_above,
 refinement_tolerance,
 config.refinement.max_steps,
 ));
+}
+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,
 refinement_tolerance,
 config.refinement.max_steps,
 )
 };
 }
-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> SlidingRegTerm<Loc<N, F>, F>
 for RadonRegTerm<F>

Mercurial > repos > pointsource_algs / file comparison

comparison: src/regularisation.rs

src/regularisation.rs