diff -r afb3ce40cff6 -r 66b8ac0ebe00 src/run.rs
--- a/src/run.rs	Thu Jul 30 09:27:44 2026 -0500
+++ b/src/run.rs	Fri Jul 31 10:41:47 2026 -0500
@@ -24,10 +24,10 @@
 use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, Regularisation, SlidingRegTerm};
 use crate::seminorms::*;
 use crate::sliding_fb::{
-    pointsource_sliding_fb_reg, SlidingFBConfig, TransportConfig, TransportProxPenalty,
+    pointsource_sliding_fb_dom, SlidingFBConfig, TransportConfig, TransportProxPenalty,
 };
 use crate::sliding_pdps::{
-    pointsource_sliding_fb_pair, pointsource_sliding_pdps_pair, SlidingPDPSConfig,
+    pointsource_sliding_fb_pair_dom, pointsource_sliding_pdps_pair_dom, SlidingPDPSConfig,
 };
 use crate::subproblem::{InnerMethod, InnerSettings};
 use crate::tolerance::Tolerance;
@@ -989,6 +989,38 @@
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
+    init: (
+        &AlgorithmConfig<F>,
+        I,
+        Plot,
+        Option<RNDM<N, F>>,
+        String,
+    ),
+    cont: impl FnOnce(
+        (&AlgorithmConfig<F>, I, Plot, Option<RNDM<N, F>>, String),
+    ) -> DynResult<RNDM<N, F>>,
+) -> DynResult<RNDM<N, F>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
+    Dat::DerivativeDomain:
+        DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, 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>>,
+{
+    run_fb_dom(f, reg, prox_penalty, init, None, cont)
+}
+
+/// Runs FB-style algorithms if `alg` so requests and `prox_penalty` matches.
+///
+/// `cont` gives a continuation attempt to find the algorithm matching the description `alg`.
+pub fn run_fb_dom<F, Dat, Reg, P, I, Plot, const N: usize>(
+    f: &Dat,
+    reg: &Reg,
+    prox_penalty: &P,
     (alg, iterator, plotter, μ0, running): (
         &AlgorithmConfig<F>,
         I,
@@ -996,6 +1028,7 @@
         Option<RNDM<N, F>>,
         String,
     ),
+    domain: Option<&P::Dom>,
     cont: impl FnOnce(
         (&AlgorithmConfig<F>, I, Plot, Option<RNDM<N, F>>, String),
     ) -> DynResult<RNDM<N, F>>,
@@ -1024,7 +1057,7 @@
         }
         &AlgorithmConfig::SlidingFB(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
-            pointsource_sliding_fb_reg(f, reg, prox_penalty, algconfig, iterator, plotter, μ0)
+            pointsource_sliding_fb_dom(f, reg, prox_penalty, algconfig, iterator, plotter, μ0, domain)
         }
         _ => cont((alg, iterator, plotter, μ0, running)),
     }
@@ -1261,8 +1294,8 @@
                     save_extra,
                     (μ0, z, y),
                     |p| {
-                        run_pdps_pair(&f, &reg, &RadonSquared, &opKz, &fnR, &fnH, p, |q| {
-                            run_pdps_pair(&f, &reg, &op𝒟, &opKz, &fnR, &fnH, q, |_| {
+                        run_pdps_pair_dom(&f, &reg, &RadonSquared, &opKz, &fnR, &fnH, p, None, |q| {
+                            run_pdps_pair_dom(&f, &reg, &op𝒟, &opKz, &fnR, &fnH, q, None, |_| {
                                 Err(NotImplemented.into())
                             })
                         })
@@ -1281,8 +1314,8 @@
                     save_extra,
                     (μ0, z, y),
                     |p| {
-                        run_pdps_pair(&f, &reg, &RadonSquared, &opKz, &fnR, &fnH, p, |q| {
-                            run_pdps_pair(&f, &reg, &op𝒟, &opKz, &fnR, &fnH, q, |_| {
+                        run_pdps_pair_dom(&f, &reg, &RadonSquared, &opKz, &fnR, &fnH, p, None, |q| {
+                            run_pdps_pair_dom(&f, &reg, &op𝒟, &opKz, &fnR, &fnH, q, None, |_| {
                                 Err(NotImplemented.into())
                             })
                         })
@@ -1297,6 +1330,7 @@
 
 type MeasureZ<F, Z, const N: usize> = Pair<RNDM<N, F>, Z>;
 
+#[deprecated]
 pub fn run_pdps_pair<F, S, Dat, Reg, Z, R, Y, KOpZ, H, P, I, Plot, const N: usize>(
     f: &Dat,
     reg: &Reg,
@@ -1304,6 +1338,57 @@
     opKz: &KOpZ,
     fnR: &R,
     fnH: &H,
+    init: (
+        &AlgorithmConfig<F>,
+        I,
+        Plot,
+        (Option<RNDM<N, F>>, Z, Y),
+        String,
+    ),
+    cont: impl FnOnce(
+        (
+            &AlgorithmConfig<F>,
+            I,
+            Plot,
+            (Option<RNDM<N, F>>, Z, Y),
+            String,
+        ),
+    ) -> DynResult<Pair<RNDM<N, F>, Z>>,
+) -> DynResult<Pair<RNDM<N, F>, Z>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>
+        + BoundedCurvature<F>,
+    S: DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
+    for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
+    //Pair<S, Z>: ClosedMul<F>,
+    RNDM<N, F>: SpikeMerging<F>,
+    Reg: SlidingRegTerm<Loc<N, F>, F>,
+    P: TransportProxPenalty<Loc<N, F>, S, Reg, F>,
+    KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
+        + GEMV<F, Z>
+        + SimplyAdjointable<Z, Y, AdjointCodomain = Z>,
+    KOpZ::SimpleAdjoint: GEMV<F, Y>,
+    Y: ClosedEuclidean<F> + Clone,
+    for<'b> &'b Y: Instance<Y>,
+    Z: ClosedEuclidean<F> + Clone + ClosedMul<F> + AXPY,
+    for<'b> &'b Z: Instance<Z>,
+    R: Prox<Z, Codomain = F>,
+    H: Conjugable<Y, F, Codomain = F>,
+    for<'b> H::Conjugate<'b>: Prox<Y>,
+    Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
+{
+    run_pdps_pair_dom(f, reg, prox_penalty, opKz, fnR, fnH, init, None, cont)
+}
+
+pub fn run_pdps_pair_dom<F, S, Dat, Reg, Z, R, Y, KOpZ, H, P, I, Plot, const N: usize>(
+    f: &Dat,
+    reg: &Reg,
+    prox_penalty: &P,
+    opKz: &KOpZ,
+    fnR: &R,
+    fnH: &H,
     (alg, iterator, plotter, μ0zy, running): (
         &AlgorithmConfig<F>,
         I,
@@ -1311,6 +1396,7 @@
         (Option<RNDM<N, F>>, Z, Y),
         String,
     ),
+    domain: Option<&P::Dom>,
     cont: impl FnOnce(
         (
             &AlgorithmConfig<F>,
@@ -1365,7 +1451,7 @@
         }
         &AlgorithmConfig::SlidingPDPS(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
-            pointsource_sliding_pdps_pair(
+            pointsource_sliding_pdps_pair_dom(
                 f,
                 reg,
                 prox_penalty,
@@ -1376,6 +1462,7 @@
                 opKz,
                 fnR,
                 fnH,
+                domain,
             )
         }
         _ => cont((alg, iterator, plotter, μ0zy, running)),
@@ -1387,6 +1474,48 @@
     reg: &Reg,
     prox_penalty: &P,
     fnR: &R,
+    init: (
+        &AlgorithmConfig<F>,
+        I,
+        Plot,
+        (Option<RNDM<N, F>>, Z),
+        String,
+    ),
+    cont: impl FnOnce(
+        (
+            &AlgorithmConfig<F>,
+            I,
+            Plot,
+            (Option<RNDM<N, F>>, Z),
+            String,
+        ),
+    ) -> DynResult<Pair<RNDM<N, F>, Z>>,
+) -> DynResult<Pair<RNDM<N, F>, Z>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>
+        + BoundedCurvature<F>,
+    S: DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
+    RNDM<N, F>: SpikeMerging<F>,
+    Reg: SlidingRegTerm<Loc<N, F>, F>,
+    P: TransportProxPenalty<Loc<N, F>, S, Reg, F>,
+    for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
+    Z: ClosedEuclidean<F> + AXPY + Clone,
+    for<'b> &'b Z: Instance<Z>,
+    R: Prox<Z, Codomain = F>,
+    Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
+    // We should not need to explicitly require this:
+    for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>,
+{
+    run_fb_pair_dom(f, reg, prox_penalty, fnR, init, None, cont)
+}
+
+pub fn run_fb_pair_dom<F, S, Dat, Reg, Z, R, P, I, Plot, const N: usize>(
+    f: &Dat,
+    reg: &Reg,
+    prox_penalty: &P,
+    fnR: &R,
     (alg, iterator, plotter, μ0z, running): (
         &AlgorithmConfig<F>,
         I,
@@ -1394,6 +1523,7 @@
         (Option<RNDM<N, F>>, Z),
         String,
     ),
+    domain: Option<&P::Dom>,
     cont: impl FnOnce(
         (
             &AlgorithmConfig<F>,
@@ -1430,7 +1560,7 @@
         }
         &AlgorithmConfig::SlidingFB(ref algconfig, prox_type) if prox_type == pt => {
             print!("{running}");
-            pointsource_sliding_fb_pair(
+            pointsource_sliding_fb_pair_dom(
                 f,
                 reg,
                 prox_penalty,
@@ -1439,6 +1569,7 @@
                 plotter,
                 μ0z,
                 fnR,
+                domain
             )
         }
         _ => cont((alg, iterator, plotter, μ0z, running)),