--- a/src/sliding_fb.rs Tue Aug 04 13:46:29 2026 -0500
+++ b/src/sliding_fb.rs Tue Aug 04 14:12:18 2026 -0500
@@ -11,6 +11,7 @@
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;
@@ -28,7 +29,6 @@
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;
@@ -223,36 +223,35 @@
const MINIMAL_PARTIAL_TRANSPORT: bool = true;
const NEW_APPROACH: bool = true;
-pub trait TransportProxPenalty<Domain, PreadjointCodomain, Reg, F = f64>:
- ProxPenalty<Domain, PreadjointCodomain, Reg, F>
+pub trait TransportProxPenalty<Point, PreadjointCodomain, Reg, F = f64>:
+ ProxPenalty<Point, PreadjointCodomain, Reg, F>
where
F: Float + ToNalgebraRealField,
- Reg: SlidingRegTerm<Domain, F>,
- Domain: Space + Clone,
+ Reg: SlidingRegTerm<Point, F>,
+ 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>,
- μ: &DiscreteMeasure<Domain, F>,
+ γ: &mut Transport<Point, 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>,
- μ: &DiscreteMeasure<Domain, F>,
- μ̆: &DiscreteMeasure<Domain, F>,
+ γ: &mut Transport<Point, F>,
+ μ: &DiscreteMeasure<Point, F>,
+ μ̆: &DiscreteMeasure<Point, F>,
τv̆: &mut PreadjointCodomain,
v: &mut PreadjointCodomain,
extra: Option<F>,
@@ -284,7 +283,6 @@
RNDM<N, F>: SpikeMerging<F>,
{
type TransportStepLength = TransportStepLength<F>;
- type Dom = Cube<N, F>;
fn get_transport_steplength(
&self,
@@ -301,7 +299,7 @@
)
}
- fn initial_transport(
+ fn initial_transport<Dom: AdmissibleDomain<Loc<N, F>>>(
&self,
γ: &mut Transport<Loc<N, F>, F>,
μ: &RNDM<N, F>,
@@ -310,7 +308,7 @@
τθ_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);
}
@@ -392,7 +390,6 @@
RNDM<N, F>: SpikeMerging<F>,
{
type TransportStepLength = TransportStepLength<F>;
- type Dom = Cube<N, F>;
fn get_transport_steplength(
&self,
@@ -409,7 +406,7 @@
)
}
- fn initial_transport(
+ fn initial_transport<Dom: AdmissibleDomain<Loc<N, F>>>(
&self,
γ: &mut Transport<Loc<N, F>, F>,
μ: &RNDM<N, F>,
@@ -418,7 +415,7 @@
τθ_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);
}
@@ -745,13 +742,13 @@
}
/// 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>,
{
@@ -819,7 +816,7 @@
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
@@ -845,7 +842,7 @@
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 };
}
}
}
@@ -978,7 +975,7 @@
/// 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,
@@ -986,7 +983,7 @@
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,
@@ -1049,7 +1046,7 @@
&mut τθ_or_adaptive,
&v,
&config.transport,
- domain
+ domain,
);
let mut attempts = 0;
@@ -1107,8 +1104,8 @@
// 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 m = prox_penalty.merge_spikes(
+ let did_merge = config.insertion.if_merge_now(&state, |ins| {
+ stats.merged += prox_penalty.merge_spikes(
&mut μ,
&mut τv̆,
&μ̆,
@@ -1116,32 +1113,35 @@
ε,
ins,
®,
- Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃) ),
+ Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃)),
);
- //if m > 0 {
- stats.merged += m;
- //v = f.differential(&μ);
- //}
- });
+ γ.prune_compat(&mut μ, &mut stats);
- γ.prune_compat(&mut μ, &mut stats);
+ if config.insertion.random_insertions {
+ if let Some(x) = domain.random_point() {
+ μ += DeltaMeasure { x, α: 0.0 };
+ stats.inserted += 1;
+ }
+ }
- // Do extra weight optimisation step heuristic
- for _ in 1..config.insertion.extra_weight_optimisation_steps {
- τv̆ = f.differential(&μ) * τ;
- prox_penalty.reweigh(
- &mut μ,
- &mut τv̆,
- τ,
- ε,
- &config.insertion,
- ®,
- &state,
- &mut stats,
- )?;
- }
+ // Do extra weight optimisation step heuristic
+ for _ in 1..config.insertion.extra_weight_optimisation_steps {
+ τv̆ = f.differential(&μ) * τ;
+ prox_penalty.reweigh(
+ &mut μ,
+ &mut τv̆,
+ τ,
+ ε,
+ &config.insertion,
+ ®,
+ &state,
+ &mut stats,
+ )?;
+ }
+ Ok(())
+ })?;
- if config.insertion.extra_weight_optimisation_steps > 0 {
+ if did_merge.is_none() {
γ.prune_compat(&mut μ, &mut stats);
}
@@ -1185,5 +1185,14 @@
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,
+ )
}