diff -r 759634edb536 -r 8f1e57821d76 src/frank_wolfe.rs
--- a/src/frank_wolfe.rs	Tue Aug 04 20:11:10 2026 -0500
+++ b/src/frank_wolfe.rs	Thu Aug 06 11:14:15 2026 -0500
@@ -20,7 +20,7 @@
 use crate::dataterm::QuadraticDataTerm;
 use crate::forward_model::ForwardModel;
 use crate::measures::merging::{SpikeMerging, SpikeMergingMethod};
-use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM};
+use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon};
 use crate::plot::Plotter;
 use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, RegTerm};
 #[allow(unused_imports)] // Used in documentation
@@ -29,18 +29,15 @@
 };
 use crate::tolerance::Tolerance;
 use crate::types::*;
-use alg_tools::bisection_tree::P2Minimise;
 use alg_tools::bounds::MinMaxMapping;
 use alg_tools::error::DynResult;
 use alg_tools::euclidean::Euclidean;
 use alg_tools::instance::Instance;
 use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorOptions, ValueIteratorFactory};
-use alg_tools::linops::Mapping;
-use alg_tools::loc::Loc;
+use alg_tools::linops::{Mapping, Space};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::norms::Norm;
 use alg_tools::norms::L2;
-use alg_tools::sets::Cube;
 
 /// Settings for [`pointsource_fw_reg`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
@@ -106,9 +103,9 @@
 /// Trait for finite dimensional weight optimisation.
 pub trait WeightOptim<
     F: Float + ToNalgebraRealField,
-    A: ForwardModel<RNDM<N, F>, F>,
+    U: Space + Clone + PartialEq,
+    A: ForwardModel<DiscreteMeasure<U, F>, F>,
     I: AlgIteratorFactory<F>,
-    const N: usize,
 >
 {
     /// Return a pre-initialisation struct for [`Self::optimise_weights`].
@@ -135,7 +132,7 @@
     /// Returns the number of iterations taken by the method configured in `inner`.
     fn optimise_weights<'a>(
         &self,
-        μ: &mut RNDM<N, F>,
+        μ: &mut DiscreteMeasure<U, F>,
         opA: &'a A,
         b: &A::Observable,
         findim_data: &FindimData<F>,
@@ -147,10 +144,10 @@
 /// Trait for regularisation terms supported by [`pointsource_fw_reg`].
 pub trait RegTermFW<
     F: Float + ToNalgebraRealField,
-    A: ForwardModel<RNDM<N, F>, F>,
+    U: Space + Clone + PartialEq,
+    A: ForwardModel<DiscreteMeasure<U, F>, F>,
     I: AlgIteratorFactory<F>,
-    const N: usize,
->: RegTerm<Loc<N, F>, F> + WeightOptim<F, A, I, N> + Mapping<RNDM<N, F>, Codomain = F>
+>: RegTerm<U, F> + WeightOptim<F, U, A, I> + Mapping<DiscreteMeasure<U, F>, Codomain = F>
 {
     /// With $g = A\_\*(Aμ-b)$, returns $(x, g(x))$ for $x$ a new point to be inserted
     /// into $μ$, as determined by the regulariser.
@@ -162,26 +159,26 @@
         g: &mut A::PreadjointCodomain,
         refinement_tolerance: F,
         max_steps: usize,
-    ) -> (Loc<N, F>, F);
+    ) -> (U, F);
 
     /// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen.
     fn relaxed_insert<'a>(
         &self,
-        μ: &mut RNDM<N, F>,
+        μ: &mut DiscreteMeasure<U, F>,
         g: &A::PreadjointCodomain,
         opA: &'a A,
-        ξ: Loc<N, F>,
+        ξ: U,
         v_ξ: F,
         findim_data: &FindimData<F>,
     );
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N>
-    for RadonRegTerm<F>
+impl<F: Float + ToNalgebraRealField, U, A, I> WeightOptim<F, U, A, I> for RadonRegTerm<F>
 where
+    U: Space + Clone + PartialEq,
     I: AlgIteratorFactory<F>,
-    A: FindimQuadraticModel<Loc<N, F>, F>,
+    A: FindimQuadraticModel<U, F>,
 {
     fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
         Ok(FindimData {
@@ -192,7 +189,7 @@
 
     fn optimise_weights<'a>(
         &self,
-        μ: &mut RNDM<N, F>,
+        μ: &mut DiscreteMeasure<U, F>,
         opA: &'a A,
         b: &A::Observable,
         findim_data: &FindimData<F>,
@@ -221,23 +218,24 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N> for RadonRegTerm<F>
+impl<F: Float + ToNalgebraRealField, U, A, I> RegTermFW<F, U, A, I> for RadonRegTerm<F>
 where
-    Cube<N, F>: P2Minimise<Loc<N, F>, F>,
+    U: Space + Clone + PartialEq,
+    RadonRegTerm<F>: RegTerm<U, F>,
     I: AlgIteratorFactory<F>,
-    A: FindimQuadraticModel<Loc<N, F>, F>,
-    A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
+    A: FindimQuadraticModel<U, F>,
+    A::PreadjointCodomain: MinMaxMapping<U, F>,
     for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
     // FIXME: the following *should not* be needed, they are already implied
-    RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
-    DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
+    DiscreteMeasure<U, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
+    DeltaMeasure<U, F>: Copy + Mapping<A::PreadjointCodomain, Codomain = F>,
 {
     fn find_insertion(
         &self,
         g: &mut A::PreadjointCodomain,
         refinement_tolerance: F,
         max_steps: usize,
-    ) -> (Loc<N, F>, F) {
+    ) -> (U, F) {
         let (ξmax, v_ξmax) = g.maximise(refinement_tolerance, max_steps);
         let (ξmin, v_ξmin) = g.minimise(refinement_tolerance, max_steps);
         if v_ξmin < 0.0 && -v_ξmin > v_ξmax {
@@ -249,10 +247,10 @@
 
     fn relaxed_insert<'a>(
         &self,
-        μ: &mut RNDM<N, F>,
+        μ: &mut DiscreteMeasure<U, F>,
         g: &A::PreadjointCodomain,
         opA: &'a A,
-        ξ: Loc<N, F>,
+        ξ: U,
         v_ξ: F,
         findim_data: &FindimData<F>,
     ) {
@@ -285,11 +283,11 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N>
-    for NonnegRadonRegTerm<F>
+impl<F: Float + ToNalgebraRealField, A, U, I> WeightOptim<F, U, A, I> for NonnegRadonRegTerm<F>
 where
+    U: Space + Clone + PartialEq,
     I: AlgIteratorFactory<F>,
-    A: FindimQuadraticModel<Loc<N, F>, F>,
+    A: FindimQuadraticModel<U, F>,
 {
     fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
         Ok(FindimData {
@@ -300,7 +298,7 @@
 
     fn optimise_weights<'a>(
         &self,
-        μ: &mut RNDM<N, F>,
+        μ: &mut DiscreteMeasure<U, F>,
         opA: &'a A,
         b: &A::Observable,
         findim_data: &FindimData<F>,
@@ -329,33 +327,33 @@
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N>
-    for NonnegRadonRegTerm<F>
+impl<F: Float + ToNalgebraRealField, U, A, I> RegTermFW<F, U, A, I> for NonnegRadonRegTerm<F>
 where
-    Cube<N, F>: P2Minimise<Loc<N, F>, F>,
+    U: Space + Clone + PartialEq,
+    NonnegRadonRegTerm<F>: RegTerm<U, F>,
     I: AlgIteratorFactory<F>,
-    A: FindimQuadraticModel<Loc<N, F>, F>,
-    A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
+    A: FindimQuadraticModel<U, F>,
+    A::PreadjointCodomain: MinMaxMapping<U, F>,
     for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
     // FIXME: the following *should not* be needed, they are already implied
-    RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
-    DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
+    DiscreteMeasure<U, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
+    DeltaMeasure<U, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
 {
     fn find_insertion(
         &self,
         g: &mut A::PreadjointCodomain,
         refinement_tolerance: F,
         max_steps: usize,
-    ) -> (Loc<N, F>, F) {
+    ) -> (U, F) {
         g.maximise(refinement_tolerance, max_steps)
     }
 
     fn relaxed_insert<'a>(
         &self,
-        μ: &mut RNDM<N, F>,
+        μ: &mut DiscreteMeasure<U, F>,
         g: &A::PreadjointCodomain,
         opA: &'a A,
-        ξ: Loc<N, F>,
+        ξ: U,
         v_ξ: F,
         findim_data: &FindimData<F>,
     ) {
@@ -401,25 +399,25 @@
 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
 /// save intermediate iteration states as images.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fw_reg<'a, F, I, A, Reg, Plot, const N: usize>(
-    f: &'a QuadraticDataTerm<F, RNDM<N, F>, A>,
+pub fn pointsource_fw_reg<'a, F, U, I, A, Reg, Plot>(
+    f: &'a QuadraticDataTerm<F, DiscreteMeasure<U, F>, A>,
     reg: &Reg,
-    //domain : Cube<N, F>,
     config: &FWConfig<F>,
     iterator: I,
     mut 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 + PartialEq,
     I: AlgIteratorFactory<IterInfo<F>>,
-    A: ForwardModel<RNDM<N, F>, F>,
-    A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>,
+    A: ForwardModel<DiscreteMeasure<U, F>, F>,
+    A::PreadjointCodomain: MinMaxMapping<U, F>,
     &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
-    Cube<N, F>: P2Minimise<Loc<N, F>, F>,
-    RNDM<N, F>: SpikeMerging<F>,
-    Reg: RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N>,
-    Plot: Plotter<A::PreadjointCodomain, A::PreadjointCodomain, RNDM<N, F>>,
+    //Cube<N, F>: P2Minimise<U, F>,
+    DiscreteMeasure<U, F>: SpikeMerging<F>,
+    Reg: RegTermFW<F, U, A, ValueIteratorFactory<F, AlgIteratorOptions>>,
+    Plot: Plotter<A::PreadjointCodomain, A::PreadjointCodomain, DiscreteMeasure<U, F>>,
 {
     let opA = f.operator();
     let b = f.data();
@@ -438,7 +436,7 @@
     let mut residual = f.residual(&μ);
 
     // Statistics
-    let full_stats = |residual: &A::Observable, ν: &RNDM<N, F>, ε, stats| IterInfo {
+    let full_stats = |residual: &A::Observable, ν: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
         value: residual.norm2_squared_div2() + reg.apply(ν),
         n_spikes: ν.len(),
         ε,