diff -r 759634edb536 -r 8f1e57821d76 src/sliding_fb.rs
--- a/src/sliding_fb.rs	Tue Aug 04 20:11:10 2026 -0500
+++ b/src/sliding_fb.rs	Thu Aug 06 11:14:15 2026 -0500
@@ -24,11 +24,11 @@
 use alg_tools::bounds::{Bounds, MinMaxMapping};
 use alg_tools::error::DynResult;
 use alg_tools::euclidean::Euclidean;
-use alg_tools::instance::Space;
+use alg_tools::instance::{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::norms::{Dist, Norm, L2};
 use anyhow::ensure;
 use std::ops::ControlFlow;
 
@@ -193,11 +193,11 @@
 }
 
 #[derive(Clone, Debug, Serialize)]
-pub struct SingleTransport<Domain, F: Float> {
+pub struct SingleTransport<U, F: Float> {
     /// Source point
-    x: Domain,
+    x: U,
     /// Target point
-    y: Domain,
+    y: U,
     /// Original mass
     α_μ_orig: F,
     /// Transported mass
@@ -211,8 +211,8 @@
 }
 
 #[derive(Clone, Debug, Serialize)]
-pub struct Transport<Domain, F: Float> {
-    vec: Vec<SingleTransport<Domain, F>>,
+pub struct Transport<U, F: Float> {
+    vec: Vec<SingleTransport<U, F>>,
 }
 
 /// Whether partially transported points are allowed.
@@ -223,21 +223,21 @@
 const MINIMAL_PARTIAL_TRANSPORT: bool = true;
 const NEW_APPROACH: bool = true;
 
-pub trait TransportProxPenalty<Point, PreadjointCodomain, Reg, F = f64>:
-    ProxPenalty<Point, PreadjointCodomain, Reg, F>
+pub trait TransportProxPenalty<U, PreadjointCodomain, Reg, F = f64>:
+    ProxPenalty<U, PreadjointCodomain, Reg, F>
 where
     F: Float + ToNalgebraRealField,
-    Reg: SlidingRegTerm<Point, F>,
-    Point: Space + Clone,
+    Reg: SlidingRegTerm<U, F>,
+    U: Space + Clone,
 {
     type TransportStepLength;
 
     /// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
     /// with step lengh τ and transport step length `θ_or_adaptive`.
-    fn initial_transport<Dom: AdmissibleDomain<Point>>(
+    fn initial_transport<Dom: AdmissibleDomain<U>>(
         &self,
-        γ: &mut Transport<Point, F>,
-        μ: &DiscreteMeasure<Point, F>,
+        γ: &mut Transport<U, F>,
+        μ: &DiscreteMeasure<U, F>,
         ε: F,
         τ: F,
         τθ_or_adaptive: &mut Self::TransportStepLength,
@@ -249,9 +249,9 @@
     /// A posteriori transport adaptation.
     fn aposteriori_transport(
         &self,
-        γ: &mut Transport<Point, F>,
-        μ: &DiscreteMeasure<Point, F>,
-        μ̆: &DiscreteMeasure<Point, F>,
+        γ: &mut Transport<U, F>,
+        μ: &DiscreteMeasure<U, F>,
+        μ̆: &DiscreteMeasure<U, F>,
         τv̆: &mut PreadjointCodomain,
         v: &mut PreadjointCodomain,
         extra: Option<F>,
@@ -377,14 +377,12 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F, M, Reg, 𝒟, O, const N: usize> TransportProxPenalty<Loc<N, F>, M, Reg, F> for 𝒟
+impl<F, M, Reg, 𝒟, const N: usize> TransportProxPenalty<Loc<N, F>, M, Reg, F> for 𝒟
 where
     F: Float + ToNalgebraRealField,
     𝒟: DiscreteMeasureOp<Loc<N, F>, F>,
     𝒟::Codomain: RealMapping<N, F>,
     M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>,
-    for<'a> &'a M: std::ops::Add<𝒟::PreCodomain, Output = O>,
-    O: MinMaxMapping<Loc<N, F>, F>,
     Reg: SlidingRegTerm<Loc<N, F>, F>,
     Self: ProxPenalty<Loc<N, F>, M, Reg, F>,
     RNDM<N, F>: SpikeMerging<F>,
@@ -608,25 +606,37 @@
     }
 }
 
+impl<U, F: Float> Norm<Radon, F> for Transport<U, F> {
+    fn norm(&self, _: Radon) -> F {
+        self.iter().map(|ρ| ρ.α_γ.abs()).sum()
+    }
+}
+
+impl<U, F: Float> MulAssign<F> for Transport<U, F> {
+    fn mul_assign(&mut self, factor: F) {
+        for ρ in self.iter_mut() {
+            ρ.α_γ *= factor;
+        }
+    }
+}
+
 #[replace_float_literals(F::cast_from(literal))]
-impl<const N: usize, F: Float> Transport<Loc<N, F>, F> {
+impl<U, F: Float> Transport<U, F> {
     pub(crate) fn new() -> Self {
         Transport { vec: Vec::new() }
     }
 
-    pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<Loc<N, F>, F>> {
+    pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<U, F>> {
         self.vec.iter()
     }
 
-    pub(crate) fn iter_mut(
-        &mut self,
-    ) -> impl Iterator<Item = &'_ mut SingleTransport<Loc<N, F>, F>> {
+    pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &'_ mut SingleTransport<U, F>> {
         self.vec.iter_mut()
     }
 
     pub(crate) fn extend<I>(&mut self, it: I)
     where
-        I: IntoIterator<Item = SingleTransport<Loc<N, F>, F>>,
+        I: IntoIterator<Item = SingleTransport<U, F>>,
     {
         self.vec.extend(it)
     }
@@ -635,62 +645,44 @@
         self.vec.len()
     }
 
-    // pub(crate) fn dist_matching(&self, μ: &RNDM<N, F>) -> F {
-    //     self.iter()
-    //         .zip(μ.iter_spikes())
-    //         .map(|(ρ, δ)| (ρ.α_γ - δ.α).abs())
-    //         .sum()
-    // }
-
-    /// Construct `μ̆`, replacing the contents of `μ`.
-    pub(crate) fn μ̆_into(&self, μ: &mut RNDM<N, F>) {
-        assert!(self.len() <= μ.len());
-
-        // First transported points
-        for (δ, ρ) in izip!(μ.iter_spikes_mut(), self.iter()) {
-            if ρ.α_γ.abs() > 0.0 {
-                // Transport – transported point
-                δ.α = ρ.α_γ;
-                δ.x = ρ.y;
-            } else {
-                // No transport – original point
-                δ.α = ρ.α_μ_orig;
-                δ.x = ρ.x;
-            }
-        }
-
-        // Then source points with partial transport
-        let mut i = self.len();
-        // This can cause the number of points to explode, so cannot have partial transport.
-        for ρ in self.iter() {
-            let α = ρ.α_μ_orig - ρ.α_γ;
-            if ρ.α_γ.abs() > F::EPSILON && α.abs() > F::EPSILON {
-                let δ = DeltaMeasure { α, x: ρ.x };
-                if i < μ.len() {
-                    μ[i] = δ;
-                } else {
-                    μ.push(δ)
-                }
-                i += 1;
-            }
-        }
-        μ.truncate(i);
-    }
-
     /// Returns $‖μ\^k - π\_♯\^0γ\^{k+1}‖$
     pub(crate) fn μ0_minus_γ0_radon(&self) -> F {
         self.vec.iter().map(|ρ| (ρ.α_μ_orig - ρ.α_γ).abs()).sum()
     }
 
+    /// Prune spikes with zero weight. To maintain correct ordering between μ and γ, also the
+    /// latter needs to be pruned when μ is.
+    pub(crate) fn prune_compat(&mut self, μ: &mut DiscreteMeasure<U, F>, stats: &mut IterInfo<F>) {
+        assert!(self.vec.len() <= μ.len());
+        let old_len = μ.len();
+        for (ρ, δ) in self.vec.iter_mut().zip(μ.iter_spikes()) {
+            ρ.retain = δ.α != F::ZERO;
+        }
+        // Do not prune new insertions on first round.
+        μ.prune_until(self.vec.len());
+        stats.pruned += old_len - μ.len();
+        self.vec.retain(|ρ| ρ.retain);
+        assert!(self.vec.len() <= μ.len());
+    }
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<U, F> Transport<U, F>
+where
+    U: Dist<L2, F>,
+    F: Float,
+    for<'a> &'a U: Instance<U>,
+{
     /// Returns $∫ c_2 d|γ|$
     pub(crate) fn c2integral(&self) -> F {
+        let sq2 = |x| (x * x) / 2.0;
         self.vec
             .iter()
-            .map(|ρ| ρ.y.dist2_squared(&ρ.x) / 2.0 * ρ.α_γ.abs())
+            .map(|ρ| sq2(ρ.y.dist(&ρ.x, L2)) * ρ.α_γ.abs())
             .sum()
     }
 
-    pub(crate) fn get_transport_stats(&self, stats: &mut IterInfo<F>, μ: &RNDM<N, F>) {
+    pub(crate) fn get_transport_stats(&self, stats: &mut IterInfo<F>, μ: &DiscreteMeasure<U, F>) {
         // TODO: This doesn't take into account μ[i].α becoming zero in the latest tranport
         // attempt, for i < self.len(), when a corresponding source term also exists with index
         // j ≥ self.len(). For now, we let that be reflected in the prune count.
@@ -726,22 +718,55 @@
                 })
         };
     }
+}
 
-    /// Prune spikes with zero weight. To maintain correct ordering between μ and γ, also the
-    /// latter needs to be pruned when μ is.
-    pub(crate) fn prune_compat(&mut self, μ: &mut RNDM<N, F>, stats: &mut IterInfo<F>) {
-        assert!(self.vec.len() <= μ.len());
-        let old_len = μ.len();
-        for (ρ, δ) in self.vec.iter_mut().zip(μ.iter_spikes()) {
-            ρ.retain = δ.α != F::ZERO;
+#[replace_float_literals(F::cast_from(literal))]
+impl<U: Clone, F: Float> Transport<U, F> {
+    // pub(crate) fn dist_matching(&self, μ: &RNDM<N, F>) -> F {
+    //     self.iter()
+    //         .zip(μ.iter_spikes())
+    //         .map(|(ρ, δ)| (ρ.α_γ - δ.α).abs())
+    //         .sum()
+    // }
+
+    /// Construct `μ̆`, replacing the contents of `μ`.
+    pub(crate) fn μ̆_into(&self, μ: &mut DiscreteMeasure<U, F>) {
+        assert!(self.len() <= μ.len());
+
+        // First transported points
+        for (δ, ρ) in izip!(μ.iter_spikes_mut(), self.iter()) {
+            if ρ.α_γ.abs() > 0.0 {
+                // Transport – transported point
+                δ.α = ρ.α_γ;
+                δ.x = ρ.y.clone();
+            } else {
+                // No transport – original point
+                δ.α = ρ.α_μ_orig;
+                δ.x = ρ.x.clone();
+            }
         }
-        // Do not prune new insertions on first round.
-        μ.prune_until(self.vec.len());
-        stats.pruned += old_len - μ.len();
-        self.vec.retain(|ρ| ρ.retain);
-        assert!(self.vec.len() <= μ.len());
+
+        // Then source points with partial transport
+        let mut i = self.len();
+        // This can cause the number of points to explode, so cannot have partial transport.
+        for ρ in self.iter() {
+            let α = ρ.α_μ_orig - ρ.α_γ;
+            if ρ.α_γ.abs() > F::EPSILON && α.abs() > F::EPSILON {
+                let δ = DeltaMeasure { α, x: ρ.x.clone() };
+                if i < μ.len() {
+                    μ[i] = δ;
+                } else {
+                    μ.push(δ)
+                }
+                i += 1;
+            }
+        }
+        μ.truncate(i);
     }
+}
 
+#[replace_float_literals(F::cast_from(literal))]
+impl<const N: usize, F: Float> Transport<Loc<N, F>, F> {
     /// Helper for initial transport. Called from [`TransportProxPenalty::initial_transport`].
     pub(super) fn do_init_transport<D, Dom: AdmissibleDomain<Loc<N, F>>>(
         &mut self,
@@ -955,20 +980,6 @@
     }
 }
 
-impl<const N: usize, F: Float> Norm<Radon, F> for Transport<Loc<N, F>, F> {
-    fn norm(&self, _: Radon) -> F {
-        self.iter().map(|ρ| ρ.α_γ.abs()).sum()
-    }
-}
-
-impl<const N: usize, F: Float> MulAssign<F> for Transport<Loc<N, F>, F> {
-    fn mul_assign(&mut self, factor: F) {
-        for ρ in self.iter_mut() {
-            ρ.α_γ *= factor;
-        }
-    }
-}
-
 /// Iteratively solve the pointsource localisation problem using sliding forward-backward
 /// splitting
 ///
@@ -976,27 +987,29 @@
 /// may be given.
 /// Inertia is currently not supported.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_sliding_fb<F, I, Dat, Reg, Plot, P, const N: usize>(
+pub fn pointsource_sliding_fb<F, U, I, Dat, Reg, Plot, P>(
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
     config: &SlidingFBConfig<F>,
     iterator: I,
     mut plotter: Plot,
-    μ0: Option<RNDM<N, F>>,
-    domain: &mut impl AdmissibleDomain<Loc<N, F>>,
-) -> DynResult<RNDM<N, F>>
+    μ0: Option<DiscreteMeasure<U, F>>,
+    domain: &mut impl AdmissibleDomain<U>,
+) -> DynResult<DiscreteMeasure<U, F>>
 where
     F: Float + ToNalgebraRealField,
+    U: Space + Clone + Dist<L2, F>,
+    for<'a> &'a U: Instance<U>,
     I: AlgIteratorFactory<IterInfo<F>>,
-    Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
+    Dat: DifferentiableMapping<DiscreteMeasure<U, F>, Codomain = F> + BoundedCurvature<F>,
     Dat::DerivativeDomain:
-        DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
+        DifferentiableMapping<U, Codomain = F> + ClosedMul<F> + MinMaxMapping<U, F>,
     //for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
-    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>>,
+    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    Reg: SlidingRegTerm<U, F>,
+    P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+    Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
 {
     // Check parameters
     ensure!(config.τ0 > 0.0, "Invalid step length parameter");
@@ -1025,7 +1038,7 @@
     let mut ε = tolerance.initial();
 
     // Statistics
-    let full_stats = |μ: &RNDM<N, F>, ε, stats| IterInfo {
+    let full_stats = |μ: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
         value: f.apply(μ) + reg.apply(μ),
         n_spikes: μ.len(),
         ε,
@@ -1114,7 +1127,7 @@
                 ε,
                 ins,
                 &reg,
-                Some(|μ̃: &RNDM<N, F>| f.apply(μ̃) + reg.apply(μ̃)),
+                Some(|μ̃: &DiscreteMeasure<U, F>| f.apply(μ̃) + reg.apply(μ̃)),
             );
             // Prune. Do not prune new insertions on first round.
             γ.prune_compat(&mut μ, &mut stats);
@@ -1167,26 +1180,28 @@
 
 #[deprecated]
 /// Version of `pointsource_sliding_fb` without domain restriction
-pub fn pointsource_sliding_fb_reg<F, I, Dat, Reg, Plot, P, const N: usize>(
+pub fn pointsource_sliding_fb_reg<F, U, I, Dat, Reg, Plot, P>(
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
     config: &SlidingFBConfig<F>,
     iterator: I,
     plotter: Plot,
-    μ0: Option<RNDM<N, F>>,
-) -> DynResult<RNDM<N, F>>
+    μ0: Option<DiscreteMeasure<U, F>>,
+) -> DynResult<DiscreteMeasure<U, F>>
 where
     F: Float + ToNalgebraRealField,
+    U: Space + Clone + Dist<L2, F>,
+    for<'a> &'a U: Instance<U>,
     I: AlgIteratorFactory<IterInfo<F>>,
-    Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
+    Dat: DifferentiableMapping<DiscreteMeasure<U, F>, Codomain = F> + BoundedCurvature<F>,
     Dat::DerivativeDomain:
-        DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
+        DifferentiableMapping<U, Codomain = F> + ClosedMul<F> + MinMaxMapping<U, F>,
     //for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
-    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>>,
+    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    Reg: SlidingRegTerm<U, F>,
+    P: TransportProxPenalty<U, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+    Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, DiscreteMeasure<U, F>>,
 {
     pointsource_sliding_fb(
         f,