pointsource_algs: comparison src/sliding

-:51dbb2cc25d9
+:8d1c959432c8
 //use nalgebra::{DVector, DMatrix};
 use itertools::izip;
 use std::iter::Iterator;
 use std::ops::MulAssign;
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::fb::*;
 use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess};
 use crate::measures::merging::SpikeMerging;
 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
 use crate::plot::Plotter;
 use alg_tools::instance::Space;
 use alg_tools::iterate::AlgIteratorFactory;
 use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Mapping, RealMapping};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::norms::Norm;
-use alg_tools::sets::{Set, Cube};
 use anyhow::ensure;
 use std::ops::ControlFlow;
 /// Transport settings for [`pointsource_sliding_fb_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 /// transport is generally preferred. If this is set to `true`, different
 /// transport adaptation heuristics will be used.
 const MINIMAL_PARTIAL_TRANSPORT: bool = true;
 const NEW_APPROACH: bool = true;
-pub trait TransportProxPenalty<Domain, PreadjointCodomain, Reg, F = f64>:
+pub trait TransportProxPenalty<Point, PreadjointCodomain, Reg, F = f64>:
-ProxPenalty<Domain, PreadjointCodomain, Reg, F>
+ProxPenalty<Point, PreadjointCodomain, Reg, F>
 where
 F: Float + ToNalgebraRealField,
-Reg: SlidingRegTerm<Domain, F>,
+Reg: SlidingRegTerm<Point, F>,
-Domain: Space + Clone,
+Point: Space + Clone,
 {
 type TransportStepLength;
-type Dom : Set<Domain>;
 /// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
 /// with step lengh τ and transport step length `θ_or_adaptive`.
-fn initial_transport(
+fn initial_transport<Dom: AdmissibleDomain<Point>>(
 &self,
-γ: &mut Transport<Domain, F>,
+γ: &mut Transport<Point, F>,
-μ: &DiscreteMeasure<Domain, F>,
+μ: &DiscreteMeasure<Point, F>,
 ε: F,
 τ: F,
 τθ_or_adaptive: &mut Self::TransportStepLength,
 v: &PreadjointCodomain,
 tconfig: &TransportConfig<F>,
-domain: Option<&Self::Dom>,
+domain: &Dom,
 );
 /// A posteriori transport adaptation.
 fn aposteriori_transport(
 &self,
-γ: &mut Transport<Domain, F>,
+γ: &mut Transport<Point, F>,
-μ: &DiscreteMeasure<Domain, F>,
+μ: &DiscreteMeasure<Point, F>,
-μ̆: &DiscreteMeasure<Domain, F>,
+μ̆: &DiscreteMeasure<Point, F>,
 τv̆: &mut PreadjointCodomain,
 v: &mut PreadjointCodomain,
 extra: Option<F>,
 ε: F,
 τ: F,
 M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>,
 Reg: SlidingRegTerm<Loc<N, F>, F>,
 RNDM<N, F>: SpikeMerging<F>,
 {
 type TransportStepLength = TransportStepLength<F>;
-type Dom = Cube<N, F>;
 fn get_transport_steplength(
 &self,
 (ℓ_F, maybe_ℓ_gradv, _maybe_transport_lip): (DynResult<F>, DynResult<F>, DynResult<F>),
 tconfig: &TransportConfig<F>,
 ℓ + ℓ_F.unwrap_or(0.0),
 ℓ_base_max_transport,
 )
 }
-fn initial_transport(
+fn initial_transport<Dom: AdmissibleDomain<Loc<N, F>>>(
 &self,
 γ: &mut Transport<Loc<N, F>, F>,
 μ: &RNDM<N, F>,
 _ε: F,
 _τ: F,
 τθ_or_adaptive: &mut TransportStepLength<F>,
 v: &M,
 tconfig: &TransportConfig<F>,
-domain: Option<&Cube<N, F>>,
+domain: &Dom,
 ) {
 γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig, domain);
 }
 fn aposteriori_transport(
 Reg: SlidingRegTerm<Loc<N, F>, F>,
 Self: ProxPenalty<Loc<N, F>, M, Reg, F>,
 RNDM<N, F>: SpikeMerging<F>,
 {
 type TransportStepLength = TransportStepLength<F>;
-type Dom = Cube<N, F>;
 fn get_transport_steplength(
 &self,
 (ℓ_F, maybe_ℓ_gradv, _maybe_transport_lip): (DynResult<F>, DynResult<F>, DynResult<F>),
 tconfig: &TransportConfig<F>,
 ℓ_base + ℓ_F.unwrap_or(0.0),
 ℓ_base_max_transport,
 )
 }
-fn initial_transport(
+fn initial_transport<Dom: AdmissibleDomain<Loc<N, F>>>(
 &self,
 γ: &mut Transport<Loc<N, F>, F>,
 μ: &RNDM<N, F>,
 _ε: F,
 _τ: F,
 τθ_or_adaptive: &mut TransportStepLength<F>,
 v: &M,
 tconfig: &TransportConfig<F>,
-domain: Option<&Cube<N, F>>,
+domain: &Dom,
 ) {
 γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig, domain);
 }
 fn aposteriori_transport(
 self.vec.retain(|ρ| ρ.retain);
 assert!(self.vec.len() <= μ.len());
 }
 /// Helper for initial transport. Called from [`TransportProxPenalty::initial_transport`].
-pub(super) fn do_init_transport<D>(
+pub(super) fn do_init_transport<D, Dom: AdmissibleDomain<Loc<N, F>>>(
 &mut self,
 v: &D,
 μ: &RNDM<N, F>,
 τθ_or_adaptive: &mut TransportStepLength<F>,
 tconfig: &TransportConfig<F>,
-domain: Option<&Cube<N, F>>,
+domain: &Dom,
 ) where
 D: DifferentiableRealMapping<N, F>,
 {
 use TransportStepLength::*;
 for ρ in self.iter_mut() {
 if ρ.fail_count < tconfig.max_fail {
 let dv_x = v.differential(&ρ.x);
 let g = &dv_x * (ρ.α_γ.signum() * τθ);
 let tmp = ρ.x - g;
-ρ.y = if domain.is_none_or(|d| d.contains(&tmp)) { tmp } else { ρ.x };
+ρ.y = if domain.contains(&tmp) { tmp } else { ρ.x };
 let n = g.norm2();
 if n >= F::EPSILON {
 // Estimate Lipschitz factor of ∇v
 let this_ℓ_gradv = (dv_x - v.differential(&ρ.y)).norm2() / n;
 *adaptive_ℓ_gradv = adaptive_ℓ_gradv.max(this_ℓ_gradv);
 if let Some(τθ) = simple_τθ {
 for ρ in self.iter_mut() {
 if ρ.fail_count <= tconfig.max_fail {
 let tmp = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * τθ);
-ρ.y = if domain.is_none_or(|d| d.contains(&tmp)) { tmp } else { ρ.x };
+ρ.y = if domain.contains(&tmp) { tmp } else { ρ.x };
 }
 }
 }
 }
 ///
 /// The parametrisation is as for [`pointsource_fb_reg`], except an additional domain restriction
 /// may be given.
 /// Inertia is currently not supported.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_sliding_fb_dom<F, I, Dat, Reg, Plot, P, const N: usize>(
+pub fn pointsource_sliding_fb<F, I, Dat, Reg, Plot, P, const N: usize>(
 f: &Dat,
 reg: &Reg,
 prox_penalty: &P,
 config: &SlidingFBConfig<F>,
 iterator: I,
 mut plotter: Plot,
 μ0: Option<RNDM<N, F>>,
-domain: Option<&P::Dom>,
+domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<RNDM<N, F>>
 where
 F: Float + ToNalgebraRealField,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
 ε,
 τ,
 &mut τθ_or_adaptive,
 &v,
 &config.transport,
-domain
+domain,
 );
 let mut attempts = 0;
 // Solve finite-dimensional subproblem several times until the dual variable for the
 γ.get_transport_stats(&mut stats, &μ);
 // Merge spikes.
 // This crucially expects the merge routine to be stable with respect to spike locations,
 // and not to performing any pruning. That is be to done below simultaneously for γ.
-config.insertion.if_merge_now(&state, |ins| {
+let did_merge = config.insertion.if_merge_now(&state, |ins| {
-let m = prox_penalty.merge_spikes(
+stats.merged += prox_penalty.merge_spikes(
 &mut μ,
 &mut τv̆,
 &μ̆,
 τ,
 ε,
 ins,
 &reg,
-Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃) ),
+Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃)),
 );
-//if m > 0 {
+γ.prune_compat(&mut μ, &mut stats);
-stats.merged += m;
-//v = f.differential(&μ);
+if config.insertion.random_insertions {
-//}
+if let Some(x) = domain.random_point() {
-});
+μ += DeltaMeasure { x, α: 0.0 };
+stats.inserted += 1;
-γ.prune_compat(&mut μ, &mut stats);
+}
+}
-// Do extra weight optimisation step heuristic
-for _ in 1..config.insertion.extra_weight_optimisation_steps {
+// Do extra weight optimisation step heuristic
-τv̆ = f.differential(&μ) * τ;
+for _ in 1..config.insertion.extra_weight_optimisation_steps {
-prox_penalty.reweigh(
+τv̆ = f.differential(&μ) * τ;
-&mut μ,
+prox_penalty.reweigh(
-&mut τv̆,
+&mut μ,
-τ,
+&mut τv̆,
-ε,
+τ,
-&config.insertion,
+ε,
-&reg,
+&config.insertion,
-&state,
+&reg,
-&mut stats,
+&state,
-)?;
+&mut stats,
-}
+)?;
+}
-if config.insertion.extra_weight_optimisation_steps > 0 {
+Ok(())
+})?;
+if did_merge.is_none() {
 γ.prune_compat(&mut μ, &mut stats);
 }
 let iter = state.iteration();
 stats.this_iters += 1;
 RNDM<N, F>: SpikeMerging<F>,
 Reg: SlidingRegTerm<Loc<N, F>, F>,
 P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
 Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
 {
-pointsource_sliding_fb_dom(f, reg, prox_penalty, config, iterator, plotter, μ0, None)
+pointsource_sliding_fb(
-}
+f,
+reg,
+prox_penalty,
+config,
+iterator,
+plotter,
+μ0,
+&mut FullSpace,
+)
+}

comparison: src/sliding_fb.rs

src/sliding_fb.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_fb.rs

src/sliding_fb.rs