# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1745804496 18000
# Node ID 2e4517b55442207ff49474d456df3e61b8fc0633
# Parent  4e80fb049dcaabb943f7bc3467a2afa6b698a8bc# Parent  1f19c6bbf07b1305db75b5fef0a129bad63ac67e
merge stable/nightly fixes from default

diff -r 1f19c6bbf07b -r 2e4517b55442 Cargo.lock
--- a/Cargo.lock	Sun Apr 27 20:29:43 2025 -0500
+++ b/Cargo.lock	Sun Apr 27 20:41:36 2025 -0500
@@ -4,7 +4,7 @@
 
 [[package]]
 name = "alg_tools"
-version = "0.3.2"
+version = "0.4.0-dev"
 dependencies = [
  "anyhow",
  "colored",
diff -r 1f19c6bbf07b -r 2e4517b55442 Cargo.toml
--- a/Cargo.toml	Sun Apr 27 20:29:43 2025 -0500
+++ b/Cargo.toml	Sun Apr 27 20:41:36 2025 -0500
@@ -1,6 +1,6 @@
 [package]
 name = "alg_tools"
-version = "0.3.2"
+version = "0.4.0-dev"
 edition = "2021"
 rust-version = "1.85"
 authors = ["Tuomo Valkonen <tuomov@iki.fi>"]
diff -r 1f19c6bbf07b -r 2e4517b55442 src/bisection_tree/aggregator.rs
--- a/src/bisection_tree/aggregator.rs	Sun Apr 27 20:29:43 2025 -0500
+++ b/src/bisection_tree/aggregator.rs	Sun Apr 27 20:41:36 2025 -0500
@@ -2,9 +2,8 @@
 Aggregation / summarisation of information in branches of bisection trees.
 */
 
+pub use crate::bounds::Bounds;
 use crate::types::*;
-use crate::sets::Set;
-use crate::instance::Instance;
 
 /// Trait for aggregating information about a branch of a [bisection tree][super::BT].
 ///
@@ -19,180 +18,60 @@
 ///    of a function on a greater domain from bounds on subdomains
 ///    (in practise [`Cube`][crate::sets::Cube]s).
 ///
-pub trait Aggregator : Clone + Sync + Send + 'static + std::fmt::Debug {
+pub trait Aggregator: Clone + Sync + Send + 'static + std::fmt::Debug {
     /// Aggregate a new data to current state.
-    fn aggregate<I>(&mut self, aggregates : I)
-    where I : Iterator<Item=Self>;
+    fn aggregate<I>(&mut self, aggregates: I)
+    where
+        I: Iterator<Item = Self>;
 
     /// Summarise several other aggregators, resetting current state.
-    fn summarise<'a, I>(&'a mut self, aggregates : I)
-    where I : Iterator<Item=&'a Self>;
+    fn summarise<'a, I>(&'a mut self, aggregates: I)
+    where
+        I: Iterator<Item = &'a Self>;
 
     /// Create a new “empty” aggregate data.
     fn new() -> Self;
 }
 
 /// An [`Aggregator`] that doesn't aggregate anything.
-#[derive(Clone,Debug)]
+#[derive(Clone, Debug)]
 pub struct NullAggregator;
 
 impl Aggregator for NullAggregator {
-    fn aggregate<I>(&mut self, _aggregates : I)
-    where I : Iterator<Item=Self> {}
-
-    fn summarise<'a, I>(&'a mut self, _aggregates : I)
-    where I : Iterator<Item=&'a Self> {}
-
-    fn new() -> Self { NullAggregator }
-}
-
-/// Upper and lower bounds on an `F`-valued function.
-#[derive(Copy,Clone,Debug)]
-pub struct Bounds<F>(
-    /// Lower bound
-    pub F,
-    /// Upper bound
-    pub F
-);
-
-impl<F : Copy> Bounds<F> {
-    /// Returns the lower bound
-    #[inline]
-    pub fn lower(&self) -> F { self.0 }
-
-    /// Returns the upper bound
-    #[inline]
-    pub fn upper(&self) -> F { self.1 }
-}
-
-impl<F : Float> Bounds<F> {
-    /// Returns a uniform bound.
-    ///
-    /// This is maximum over the absolute values of the upper and lower bound.
-    #[inline]
-    pub fn uniform(&self) -> F {
-        let &Bounds(lower, upper) = self;
-        lower.abs().max(upper.abs())
+    fn aggregate<I>(&mut self, _aggregates: I)
+    where
+        I: Iterator<Item = Self>,
+    {
     }
 
-    /// Construct a bounds, making sure `lower` bound is less than `upper`
-    #[inline]
-    pub fn corrected(lower : F, upper : F) -> Self {
-        if lower <= upper {
-            Bounds(lower, upper)
-        } else {
-            Bounds(upper, lower)
-        }
+    fn summarise<'a, I>(&'a mut self, _aggregates: I)
+    where
+        I: Iterator<Item = &'a Self>,
+    {
     }
 
-    /// Refine the lower bound
-    #[inline]
-    pub fn refine_lower(&self, lower : F) -> Self {
-        let &Bounds(l, u) = self;
-        debug_assert!(l <= u);
-        Bounds(l.max(lower), u.max(lower))
-    }
-
-    /// Refine the lower bound
-    #[inline]
-    pub fn refine_upper(&self, upper : F) -> Self {
-        let &Bounds(l, u) = self;
-        debug_assert!(l <= u);
-        Bounds(l.min(upper), u.min(upper))
+    fn new() -> Self {
+        NullAggregator
     }
 }
 
-impl<'a, F : Float> std::ops::Add<Self> for Bounds<F> {
-    type Output = Self;
-    #[inline]
-    fn add(self, Bounds(l2, u2) : Self) -> Self::Output {
-        let Bounds(l1, u1) = self;
-        debug_assert!(l1 <= u1 && l2 <= u2);
-        Bounds(l1 + l2, u1 + u2)
-    }
-}
-
-impl<'a, F : Float> std::ops::Mul<Self> for Bounds<F> {
-    type Output = Self;
-    #[inline]
-    fn mul(self, Bounds(l2, u2) : Self) -> Self::Output {
-        let Bounds(l1, u1) = self;
-        debug_assert!(l1 <= u1 && l2 <= u2);
-        let a = l1 * l2;
-        let b = u1 * u2;
-        // The order may flip when negative numbers are involved, so need min/max
-        Bounds(a.min(b), a.max(b))
-    }
-}
-
-impl<F : Float> std::iter::Product for Bounds<F> {
+impl<F: Float> Aggregator for Bounds<F> {
     #[inline]
-    fn product<I>(mut iter: I) -> Self
-    where I: Iterator<Item = Self> {
-        match iter.next() {
-            None => Bounds(F::ZERO, F::ZERO),
-            Some(init) => iter.fold(init, |a, b| a*b)
-        }
-    }
-}
-
-impl<F : Float> Set<F> for Bounds<F> {
-    fn contains<I : Instance<F>>(&self, item : I) -> bool {
-        let v = item.own();
-        let &Bounds(l, u) = self;
-        debug_assert!(l <= u);
-        l <= v && v <= u
-    }
-}
-
-impl<F : Float> Bounds<F> {
-    /// Calculate a common bound (glb, lub) for two bounds.
-    #[inline]
-    pub fn common(&self, &Bounds(l2, u2) : &Self) -> Self {
-        let &Bounds(l1, u1) = self;
-        debug_assert!(l1 <= u1 && l2 <= u2);
-        Bounds(l1.min(l2), u1.max(u2))
-    }
-
-    /// Indicates whether `Self` is a superset of the argument bound.
-    #[inline]
-    pub fn superset(&self, &Bounds(l2, u2) : &Self) -> bool {
-        let &Bounds(l1, u1) = self;
-        debug_assert!(l1 <= u1 && l2 <= u2);
-        l1 <= l2 && u2 <= u1
-    }
-
-    /// Returns the greatest bound contained by both argument bounds, if one exists.
-    #[inline]
-    pub fn glb(&self, &Bounds(l2, u2) : &Self) -> Option<Self> {
-        let &Bounds(l1, u1) = self;
-        debug_assert!(l1 <= u1 && l2 <= u2);
-        let l = l1.max(l2);
-        let u = u1.min(u2);
-        debug_assert!(l <= u);
-        if l < u {
-            Some(Bounds(l, u))
-        } else {
-            None
-        }
-    }
-}
-
-impl<F : Float> Aggregator for Bounds<F> {
-    #[inline]
-    fn aggregate<I>(&mut self, aggregates : I)
-    where I : Iterator<Item=Self> {
+    fn aggregate<I>(&mut self, aggregates: I)
+    where
+        I: Iterator<Item = Self>,
+    {
         *self = aggregates.fold(*self, |a, b| a + b);
     }
 
     #[inline]
-    fn summarise<'a, I>(&'a mut self, mut aggregates : I)
-    where I : Iterator<Item=&'a Self> {
+    fn summarise<'a, I>(&'a mut self, mut aggregates: I)
+    where
+        I: Iterator<Item = &'a Self>,
+    {
         *self = match aggregates.next() {
             None => Bounds(F::ZERO, F::ZERO), // No parts in this cube; the function is zero
-            Some(&bounds) => {
-                aggregates.fold(bounds, |a, b| a.common(b))
-            }
+            Some(&bounds) => aggregates.fold(bounds, |a, b| a.common(b)),
         }
     }
 
diff -r 1f19c6bbf07b -r 2e4517b55442 src/bisection_tree/bt.rs
--- a/src/bisection_tree/bt.rs	Sun Apr 27 20:29:43 2025 -0500
+++ b/src/bisection_tree/bt.rs	Sun Apr 27 20:41:36 2025 -0500
@@ -1,34 +1,28 @@
-
 /*!
 Bisection tree basics, [`BT`] type and the [`BTImpl`] trait.
 */
 
-use std::slice::IterMut;
-use std::iter::once;
-use std::sync::Arc;
-use serde::{Serialize, Deserialize};
+use itertools::izip;
 pub(super) use nalgebra::Const;
-use itertools::izip;
+use serde::{Deserialize, Serialize};
+use std::iter::once;
+use std::slice::IterMut;
+use std::sync::Arc;
 
-use crate::types::{Float, Num};
-use crate::parallelism::{with_task_budget, TaskBudget};
+use super::aggregator::*;
+use super::support::*;
 use crate::coefficients::pow;
-use crate::maputil::{
-    array_init,
-    map2,
-    map2_indexed,
-    collect_into_array_unchecked
-};
+use crate::loc::Loc;
+use crate::maputil::{array_init, collect_into_array_unchecked, map2, map2_indexed};
+use crate::parallelism::{with_task_budget, TaskBudget};
 use crate::sets::Cube;
-use crate::loc::Loc;
-use super::support::*;
-use super::aggregator::*;
+use crate::types::{Float, Num};
 
 /// An enum that indicates whether a [`Node`] of a [`BT`] is uninitialised, leaf, or branch.
 ///
 /// For the type and const parametere, see the [module level documentation][super].
-#[derive(Clone,Debug)]
-pub(super) enum NodeOption<F : Num, D, A : Aggregator, const N : usize, const P : usize> {
+#[derive(Clone, Debug)]
+pub(super) enum NodeOption<F: Num, D, A: Aggregator, const N: usize, const P: usize> {
     /// Indicates an uninitilised node; may become a branch or a leaf.
     // TODO: Could optimise Uninitialised away by simply treat Leaf with an empty Vec as
     // something that can be still replaced with Branches.
@@ -44,39 +38,41 @@
 ///
 /// For the type and const parameteres, see the [module level documentation][super].
 #[derive(Clone, Debug)]
-pub struct Node<F : Num, D, A : Aggregator, const N : usize, const P : usize> {
+pub struct Node<F: Num, D, A: Aggregator, const N: usize, const P: usize> {
     /// The data or branches under the node.
-    pub(super) data : NodeOption<F, D, A, N, P>,
+    pub(super) data: NodeOption<F, D, A, N, P>,
     /// Aggregator for `data`.
-    pub(super) aggregator : A,
+    pub(super) aggregator: A,
 }
 
 /// Branching information of a [`Node`] of a [`BT`] bisection tree into `P` subnodes.
 ///
 /// For the type and const parameters, see the [module level documentation][super].
 #[derive(Clone, Debug)]
-pub(super) struct Branches<F : Num, D, A : Aggregator, const N : usize, const P : usize> {
+pub(super) struct Branches<F: Num, D, A: Aggregator, const N: usize, const P: usize> {
     /// Point for subdivision of the (unstored) [`Cube`] corresponding to the node.
-    pub(super) branch_at : Loc<F, N>,
+    pub(super) branch_at: Loc<F, N>,
     /// Subnodes
-    pub(super) nodes : [Node<F, D, A, N, P>; P],
+    pub(super) nodes: [Node<F, D, A, N, P>; P],
 }
 
 /// Dirty workaround to broken Rust drop, see [https://github.com/rust-lang/rust/issues/58068]().
-impl<F : Num, D, A : Aggregator, const N : usize, const P : usize>
-Drop for Node<F, D, A, N, P> {
+impl<F: Num, D, A: Aggregator, const N: usize, const P: usize> Drop for Node<F, D, A, N, P> {
     fn drop(&mut self) {
         use NodeOption as NO;
 
-        let process = |brc : Arc<Branches<F, D, A, N, P>>,
-                       to_drop : &mut Vec<Arc<Branches<F, D, A, N, P>>>| {
+        let process = |brc: Arc<Branches<F, D, A, N, P>>,
+                       to_drop: &mut Vec<Arc<Branches<F, D, A, N, P>>>| {
             // We only drop Branches if we have the only strong reference.
             // FIXME: update the RwLocks on Nodes.
-            Arc::try_unwrap(brc).ok().map(|branches| branches.nodes.map(|mut node| {
-                if let NO::Branches(brc2) = std::mem::replace(&mut node.data, NO::Uninitialised) {
-                    to_drop.push(brc2)
-                }
-            }));
+            Arc::try_unwrap(brc).ok().map(|branches| {
+                branches.nodes.map(|mut node| {
+                    if let NO::Branches(brc2) = std::mem::replace(&mut node.data, NO::Uninitialised)
+                    {
+                        to_drop.push(brc2)
+                    }
+                })
+            });
         };
 
         // We mark Self as NodeOption::Uninitialised, extracting the real contents.
@@ -98,9 +94,9 @@
 /// Trait for the depth of a [`BT`].
 ///
 /// This will generally be either a runtime [`DynamicDepth`] or compile-time [`Const`] depth.
-pub trait Depth : 'static + Copy + Send + Sync + std::fmt::Debug {
+pub trait Depth: 'static + Copy + Send + Sync + std::fmt::Debug {
     /// Lower depth type.
-    type Lower : Depth;
+    type Lower: Depth;
 
     /// Returns a lower depth, if there still is one.
     fn lower(&self) -> Option<Self::Lower>;
@@ -113,26 +109,26 @@
 }
 
 /// Dynamic (runtime) [`Depth`] for a [`BT`].
-#[derive(Copy,Clone,Debug,Serialize,Deserialize)]
+#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
 pub struct DynamicDepth(
     /// The depth
-    pub u8
+    pub u8,
 );
 
 impl Depth for DynamicDepth {
     type Lower = Self;
     #[inline]
     fn lower(&self) -> Option<Self> {
-        if self.0>0 {
-            Some(DynamicDepth(self.0-1))
-         } else {
+        if self.0 > 0 {
+            Some(DynamicDepth(self.0 - 1))
+        } else {
             None
         }
     }
 
     #[inline]
     fn lower_or(&self) -> Self {
-        DynamicDepth(if self.0>0 { self.0 - 1 } else { 0 })
+        DynamicDepth(if self.0 > 0 { self.0 - 1 } else { 0 })
     }
 
     #[inline]
@@ -143,9 +139,15 @@
 
 impl Depth for Const<0> {
     type Lower = Self;
-    fn lower(&self) -> Option<Self::Lower> { None }
-    fn lower_or(&self) -> Self::Lower { Const }
-    fn value(&self) -> u32 { 0 }
+    fn lower(&self) -> Option<Self::Lower> {
+        None
+    }
+    fn lower_or(&self) -> Self::Lower {
+        Const
+    }
+    fn value(&self) -> u32 {
+        0
+    }
 }
 
 macro_rules! impl_constdepth {
@@ -165,7 +167,7 @@
 /// The const parameter `P` from the [module level documentation][super] is required to satisfy
 /// `Const<P> : Branchcount<N>`.
 /// This trait is implemented for `P=pow(2, N)` for small `N`.
-pub trait BranchCount<const N : usize> {}
+pub trait BranchCount<const N: usize> {}
 macro_rules! impl_branchcount {
     ($($n:literal)*) => { $(
         impl BranchCount<$n> for Const<{pow(2, $n)}>{}
@@ -173,18 +175,19 @@
 }
 impl_branchcount!(1 2 3 4 5 6 7 8);
 
-impl<F : Float, D, A, const N : usize, const P : usize> Branches<F,D,A,N,P>
-where Const<P> : BranchCount<N>,
-      A : Aggregator
+impl<F: Float, D, A, const N: usize, const P: usize> Branches<F, D, A, N, P>
+where
+    Const<P>: BranchCount<N>,
+    A: Aggregator,
 {
     /// Returns the index in {0, …, `P`-1} for the branch to which the point `x` corresponds.
     ///
     /// This only takes the branch subdivision point $d$ into account, so is always succesfull.
     /// Thus, for this point, each branch corresponds to a quadrant of $ℝ^N$ relative to $d$.
-    fn get_node_index(&self, x : &Loc<F, N>) -> usize {
-        izip!(0..P, x.iter(), self.branch_at.iter()).map(|(i, x_i, branch_i)|
-            if x_i > branch_i { 1<<i } else { 0 }
-        ).sum()
+    fn get_node_index(&self, x: &Loc<F, N>) -> usize {
+        izip!(0..P, x.iter(), self.branch_at.iter())
+            .map(|(i, x_i, branch_i)| if x_i > branch_i { 1 << i } else { 0 })
+            .sum()
     }
 
     /// Returns the node within `Self` containing the point `x`.
@@ -192,24 +195,24 @@
     /// This only takes the branch subdivision point $d$ into account, so is always succesfull.
     /// Thus, for this point, each branch corresponds to a quadrant of $ℝ^N$ relative to $d$.
     #[inline]
-    fn get_node(&self, x : &Loc<F,N>) -> &Node<F,D,A,N,P> {
-         &self.nodes[self.get_node_index(x)]
+    fn get_node(&self, x: &Loc<F, N>) -> &Node<F, D, A, N, P> {
+        &self.nodes[self.get_node_index(x)]
     }
 }
 
 /// An iterator over the $P=2^N$ subcubes of a [`Cube`] subdivided at a point `d`.
-pub(super) struct SubcubeIter<'b, F : Float, const N : usize, const P : usize> {
-    domain : &'b Cube<F, N>,
-    branch_at : Loc<F, N>,
-    index : usize,
+pub(super) struct SubcubeIter<'b, F: Float, const N: usize, const P: usize> {
+    domain: &'b Cube<F, N>,
+    branch_at: Loc<F, N>,
+    index: usize,
 }
 
 /// Returns the `i`:th subcube of `domain` subdivided at `branch_at`.
 #[inline]
-fn get_subcube<F : Float, const N : usize>(
-    branch_at : &Loc<F, N>,
-    domain : &Cube<F, N>,
-    i : usize
+fn get_subcube<F: Float, const N: usize>(
+    branch_at: &Loc<F, N>,
+    domain: &Cube<F, N>,
+    i: usize,
 ) -> Cube<F, N> {
     map2_indexed(branch_at, domain, move |j, &branch, &[start, end]| {
         if i & (1 << j) != 0 {
@@ -217,11 +220,11 @@
         } else {
             [start, branch]
         }
-    }).into()
+    })
+    .into()
 }
 
-impl<'a, 'b, F : Float, const N : usize, const P : usize> Iterator
-for SubcubeIter<'b, F, N, P> {
+impl<'a, 'b, F: Float, const N: usize, const P: usize> Iterator for SubcubeIter<'b, F, N, P> {
     type Item = Cube<F, N>;
     #[inline]
     fn next(&mut self) -> Option<Self::Item> {
@@ -235,30 +238,33 @@
     }
 }
 
-impl<F : Float, D, A, const N : usize, const P : usize>
-Branches<F,D,A,N,P>
-where Const<P> : BranchCount<N>,
-      A : Aggregator,
-      D : 'static + Copy + Send + Sync {
-
+impl<F: Float, D, A, const N: usize, const P: usize> Branches<F, D, A, N, P>
+where
+    Const<P>: BranchCount<N>,
+    A: Aggregator,
+    D: 'static + Copy + Send + Sync,
+{
     /// Creates a new node branching structure, subdividing `domain` based on the
     /// [hint][Support::support_hint] of `support`.
-    pub(super) fn new_with<S : LocalAnalysis <F, A, N>>(
-        domain : &Cube<F,N>,
-        support : &S
+    pub(super) fn new_with<S: LocalAnalysis<F, A, N> + Support<F, N>>(
+        domain: &Cube<F, N>,
+        support: &S,
     ) -> Self {
         let hint = support.bisection_hint(domain);
         let branch_at = map2(&hint, domain, |h, r| {
-            h.unwrap_or_else(|| (r[0]+r[1])/F::TWO).max(r[0]).min(r[1])
-        }).into();
-        Branches{
-            branch_at : branch_at,
-            nodes : array_init(|| Node::new()),
+            h.unwrap_or_else(|| (r[0] + r[1]) / F::TWO)
+                .max(r[0])
+                .min(r[1])
+        })
+        .into();
+        Branches {
+            branch_at: branch_at,
+            nodes: array_init(|| Node::new()),
         }
     }
 
     /// Summarises the aggregators of these branches into `agg`
-    pub(super) fn summarise_into(&self, agg : &mut A) {
+    pub(super) fn summarise_into(&self, agg: &mut A) {
         // We need to create an array of the aggregators clones due to the RwLock.
         agg.summarise(self.nodes.iter().map(Node::get_aggregator));
     }
@@ -266,12 +272,11 @@
     /// Returns an iterator over the subcubes of `domain` subdivided at the branching point
     /// of `self`.
     #[inline]
-    pub(super) fn iter_subcubes<'b>(&self, domain : &'b Cube<F, N>)
-    -> SubcubeIter<'b, F, N, P> {
+    pub(super) fn iter_subcubes<'b>(&self, domain: &'b Cube<F, N>) -> SubcubeIter<'b, F, N, P> {
         SubcubeIter {
-            domain : domain,
-            branch_at : self.branch_at,
-            index : 0,
+            domain: domain,
+            branch_at: self.branch_at,
+            index: 0,
         }
     }
 
@@ -286,8 +291,10 @@
 
     /// Mutably iterate over all nodes and corresponding subcubes of `self`.
     #[inline]
-    pub(super) fn nodes_and_cubes_mut<'a, 'b>(&'a mut self, domain : &'b Cube<F, N>)
-    -> std::iter::Zip<IterMut<'a, Node<F,D,A,N,P>>, SubcubeIter<'b, F, N, P>> {
+    pub(super) fn nodes_and_cubes_mut<'a, 'b>(
+        &'a mut self,
+        domain: &'b Cube<F, N>,
+    ) -> std::iter::Zip<IterMut<'a, Node<F, D, A, N, P>>, SubcubeIter<'b, F, N, P>> {
         let subcube_iter = self.iter_subcubes(domain);
         self.nodes.iter_mut().zip(subcube_iter)
     }
@@ -296,12 +303,16 @@
     #[inline]
     fn recurse<'scope, 'smaller, 'refs>(
         &'smaller mut self,
-        domain : &'smaller Cube<F, N>,
-        task_budget : TaskBudget<'scope, 'refs>,
-        guard : impl Fn(&Node<F,D,A,N,P>, &Cube<F, N>) -> bool + Send + 'smaller,
-        mut f : impl for<'a> FnMut(&mut Node<F,D,A,N,P>, &Cube<F, N>, TaskBudget<'smaller, 'a>)
-                + Send + Copy + 'smaller
-    ) where 'scope : 'smaller {
+        domain: &'smaller Cube<F, N>,
+        task_budget: TaskBudget<'scope, 'refs>,
+        guard: impl Fn(&Node<F, D, A, N, P>, &Cube<F, N>) -> bool + Send + 'smaller,
+        mut f: impl for<'a> FnMut(&mut Node<F, D, A, N, P>, &Cube<F, N>, TaskBudget<'smaller, 'a>)
+            + Send
+            + Copy
+            + 'smaller,
+    ) where
+        'scope: 'smaller,
+    {
         let subs = self.nodes_and_cubes_mut(domain);
         task_budget.zoom(move |s| {
             for (node, subcube) in subs {
@@ -321,19 +332,23 @@
     ///  * `support` is the [`Support`] that is used determine with which subcubes of `domain`
     ///     (at subdivision depth `new_leaf_depth`) the data `d` is to be associated with.
     ///
-    pub(super) fn insert<'refs, 'scope, M : Depth, S : LocalAnalysis<F, A, N>>(
+    pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis<F, A, N> + Support<F, N>>(
         &mut self,
-        domain : &Cube<F,N>,
-        d : D,
-        new_leaf_depth : M,
-        support : &S,
-        task_budget : TaskBudget<'scope, 'refs>,
+        domain: &Cube<F, N>,
+        d: D,
+        new_leaf_depth: M,
+        support: &S,
+        task_budget: TaskBudget<'scope, 'refs>,
     ) {
         let support_hint = support.support_hint();
-        self.recurse(domain, task_budget,
-                     |_, subcube| support_hint.intersects(&subcube),
-                     move |node, subcube, new_budget| node.insert(subcube, d, new_leaf_depth, support,
-                                                                  new_budget));
+        self.recurse(
+            domain,
+            task_budget,
+            |_, subcube| support_hint.intersects(&subcube),
+            move |node, subcube, new_budget| {
+                node.insert(subcube, d, new_leaf_depth, support, new_budget)
+            },
+        );
     }
 
     /// Construct a new instance of the branch for a different aggregator.
@@ -344,20 +359,24 @@
     /// generator's `SupportType`.
     pub(super) fn convert_aggregator<ANew, G>(
         self,
-        generator : &G,
-        domain : &Cube<F, N>
-    ) -> Branches<F,D,ANew,N,P>
-    where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=D>,
-          G::SupportType : LocalAnalysis<F, ANew, N> {
+        generator: &G,
+        domain: &Cube<F, N>,
+    ) -> Branches<F, D, ANew, N, P>
+    where
+        ANew: Aggregator,
+        G: SupportGenerator<F, N, Id = D>,
+        G::SupportType: LocalAnalysis<F, ANew, N>,
+    {
         let branch_at = self.branch_at;
         let subcube_iter = self.iter_subcubes(domain);
-        let new_nodes = self.nodes.into_iter().zip(subcube_iter).map(|(node, subcube)| {
-            Node::convert_aggregator(node, generator, &subcube)
-        });
+        let new_nodes = self
+            .nodes
+            .into_iter()
+            .zip(subcube_iter)
+            .map(|(node, subcube)| Node::convert_aggregator(node, generator, &subcube));
         Branches {
-            branch_at : branch_at,
-            nodes : collect_into_array_unchecked(new_nodes),
+            branch_at: branch_at,
+            nodes: collect_into_array_unchecked(new_nodes),
         }
     }
 
@@ -367,30 +386,36 @@
     /// [`Support`]s. The `domain` is the cube corresponding to `self`.
     pub(super) fn refresh_aggregator<'refs, 'scope, G>(
         &mut self,
-        generator : &G,
-        domain : &Cube<F, N>,
-        task_budget : TaskBudget<'scope, 'refs>,
-    ) where G : SupportGenerator<F, N, Id=D>,
-            G::SupportType : LocalAnalysis<F, A, N> {
-        self.recurse(domain, task_budget,
-                     |_, _| true,
-                     move |node, subcube, new_budget| node.refresh_aggregator(generator, subcube,
-                                                                              new_budget));
+        generator: &G,
+        domain: &Cube<F, N>,
+        task_budget: TaskBudget<'scope, 'refs>,
+    ) where
+        G: SupportGenerator<F, N, Id = D>,
+        G::SupportType: LocalAnalysis<F, A, N>,
+    {
+        self.recurse(
+            domain,
+            task_budget,
+            |_, _| true,
+            move |node, subcube, new_budget| {
+                node.refresh_aggregator(generator, subcube, new_budget)
+            },
+        );
     }
 }
 
-impl<F : Float, D, A, const N : usize, const P : usize>
-Node<F,D,A,N,P>
-where Const<P> : BranchCount<N>,
-      A : Aggregator,
-      D : 'static + Copy + Send + Sync {
-
+impl<F: Float, D, A, const N: usize, const P: usize> Node<F, D, A, N, P>
+where
+    Const<P>: BranchCount<N>,
+    A: Aggregator,
+    D: 'static + Copy + Send + Sync,
+{
     /// Create a new node
     #[inline]
     pub(super) fn new() -> Self {
         Node {
-            data : NodeOption::Uninitialised,
-            aggregator : A::new(),
+            data: NodeOption::Uninitialised,
+            aggregator: A::new(),
         }
     }
 
@@ -407,7 +432,7 @@
 
     /// Get leaf data iterator
     #[inline]
-    pub(super) fn get_leaf_data_iter(&self, x : &Loc<F, N>) -> Option<std::slice::Iter<'_, D>> {
+    pub(super) fn get_leaf_data_iter(&self, x: &Loc<F, N>) -> Option<std::slice::Iter<'_, D>> {
         match self.data {
             NodeOption::Uninitialised => None,
             NodeOption::Leaf(ref data) => Some(data.iter()),
@@ -434,13 +459,13 @@
     /// If `self` is a [`NodeOption::Branches`], the data is passed to branches whose subcubes
     /// `support` intersects. If an [`NodeOption::Uninitialised`] node is encountered, a new leaf is
     /// created at a minimum depth of `new_leaf_depth`.
-    pub(super) fn insert<'refs, 'scope, M : Depth, S : LocalAnalysis <F, A, N>>(
+    pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis<F, A, N> + Support<F, N>>(
         &mut self,
-        domain : &Cube<F,N>,
-        d : D,
-        new_leaf_depth : M,
-        support : &S,
-        task_budget : TaskBudget<'scope, 'refs>,
+        domain: &Cube<F, N>,
+        d: D,
+        new_leaf_depth: M,
+        support: &S,
+        task_budget: TaskBudget<'scope, 'refs>,
     ) {
         match &mut self.data {
             NodeOption::Uninitialised => {
@@ -451,10 +476,10 @@
                         self.aggregator.aggregate(once(a));
                         // TODO: this is currently a dirty hard-coded heuristic;
                         // should add capacity as a parameter
-                        let mut vec = Vec::with_capacity(2*P+1);
+                        let mut vec = Vec::with_capacity(2 * P + 1);
                         vec.push(d);
                         NodeOption::Leaf(vec)
-                    },
+                    }
                     Some(lower) => {
                         let b = Arc::new({
                             let mut b0 = Branches::new_with(domain, support);
@@ -465,19 +490,19 @@
                         NodeOption::Branches(b)
                     }
                 }
-            },
+            }
             NodeOption::Leaf(leaf) => {
                 leaf.push(d);
                 let a = support.local_analysis(&domain);
                 self.aggregator.aggregate(once(a));
-            },
+            }
             NodeOption::Branches(b) => {
                 // FIXME: recursion that may cause stack overflow if the tree becomes
                 // very deep, e.g. due to [`BTSearch::search_and_refine`].
                 let bm = Arc::make_mut(b);
                 bm.insert(domain, d, new_leaf_depth.lower_or(), support, task_budget);
                 bm.summarise_into(&mut self.aggregator);
-            },
+            }
         }
     }
 
@@ -489,18 +514,19 @@
     /// generator's `SupportType`.
     pub(super) fn convert_aggregator<ANew, G>(
         mut self,
-        generator : &G,
-        domain : &Cube<F, N>
-    ) -> Node<F,D,ANew,N,P>
-    where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=D>,
-          G::SupportType : LocalAnalysis<F, ANew, N> {
-
+        generator: &G,
+        domain: &Cube<F, N>,
+    ) -> Node<F, D, ANew, N, P>
+    where
+        ANew: Aggregator,
+        G: SupportGenerator<F, N, Id = D>,
+        G::SupportType: LocalAnalysis<F, ANew, N>,
+    {
         // The mem::replace is needed due to the [`Drop`] implementation to extract self.data.
         match std::mem::replace(&mut self.data, NodeOption::Uninitialised) {
             NodeOption::Uninitialised => Node {
-                data : NodeOption::Uninitialised,
-                aggregator : ANew::new(),
+                data: NodeOption::Uninitialised,
+                aggregator: ANew::new(),
             },
             NodeOption::Leaf(v) => {
                 let mut anew = ANew::new();
@@ -510,10 +536,10 @@
                 }));
 
                 Node {
-                    data : NodeOption::Leaf(v),
-                    aggregator : anew,
+                    data: NodeOption::Leaf(v),
+                    aggregator: anew,
                 }
-            },
+            }
             NodeOption::Branches(b) => {
                 // FIXME: recursion that may cause stack overflow if the tree becomes
                 // very deep, e.g. due to [`BTSearch::search_and_refine`].
@@ -521,8 +547,8 @@
                 let mut anew = ANew::new();
                 bnew.summarise_into(&mut anew);
                 Node {
-                    data : NodeOption::Branches(Arc::new(bnew)),
-                    aggregator : anew,
+                    data: NodeOption::Branches(Arc::new(bnew)),
+                    aggregator: anew,
                 }
             }
         }
@@ -534,20 +560,22 @@
     /// [`Support`]s. The `domain` is the cube corresponding to `self`.
     pub(super) fn refresh_aggregator<'refs, 'scope, G>(
         &mut self,
-        generator : &G,
-        domain : &Cube<F, N>,
-        task_budget : TaskBudget<'scope, 'refs>,
-    ) where G : SupportGenerator<F, N, Id=D>,
-            G::SupportType : LocalAnalysis<F, A, N> {
+        generator: &G,
+        domain: &Cube<F, N>,
+        task_budget: TaskBudget<'scope, 'refs>,
+    ) where
+        G: SupportGenerator<F, N, Id = D>,
+        G::SupportType: LocalAnalysis<F, A, N>,
+    {
         match &mut self.data {
-            NodeOption::Uninitialised => { },
+            NodeOption::Uninitialised => {}
             NodeOption::Leaf(v) => {
                 self.aggregator = A::new();
-                self.aggregator.aggregate(v.iter().map(|d| {
-                    generator.support_for(*d)
-                             .local_analysis(&domain)
-                }));
-            },
+                self.aggregator.aggregate(
+                    v.iter()
+                        .map(|d| generator.support_for(*d).local_analysis(&domain)),
+                );
+            }
             NodeOption::Branches(ref mut b) => {
                 // FIXME: recursion that may cause stack overflow if the tree becomes
                 // very deep, e.g. due to [`BTSearch::search_and_refine`].
@@ -563,11 +591,13 @@
 ///
 /// This can be removed and the methods implemented directly on [`BT`] once Rust's const generics
 /// are flexible enough to allow fixing `P=pow(2, N)`.
-pub trait BTNode<F, D, A, const N : usize>
-where F : Float,
-        D : 'static + Copy,
-        A : Aggregator {
-    type Node : Clone + std::fmt::Debug;
+pub trait BTNode<F, D, A, const N: usize>
+where
+    F: Float,
+    D: 'static + Copy,
+    A: Aggregator,
+{
+    type Node: Clone + std::fmt::Debug;
 }
 
 /// Helper structure for looking up a [`Node`] without the knowledge of `P`.
@@ -580,48 +610,52 @@
 /// Basic interface to a [`BT`] bisection tree.
 ///
 /// Further routines are provided by the [`BTSearch`][super::refine::BTSearch] trait.
-pub trait BTImpl<F : Float, const N : usize> : std::fmt::Debug + Clone + GlobalAnalysis<F, Self::Agg> {
+pub trait BTImpl<F: Float, const N: usize>:
+    std::fmt::Debug + Clone + GlobalAnalysis<F, Self::Agg>
+{
     /// The data type stored in the tree
-    type Data : 'static + Copy + Send + Sync;
+    type Data: 'static + Copy + Send + Sync;
     /// The depth type of the tree
-    type Depth : Depth;
+    type Depth: Depth;
     /// The type for the [aggregate information][Aggregator] about the `Data` stored in each node
     /// of the tree.
-    type Agg : Aggregator;
+    type Agg: Aggregator;
     /// The type of the tree with the aggregator converted to `ANew`.
-    type Converted<ANew> : BTImpl<F, N, Data=Self::Data, Agg=ANew> where ANew : Aggregator;
+    type Converted<ANew>: BTImpl<F, N, Data = Self::Data, Agg = ANew>
+    where
+        ANew: Aggregator;
 
     /// Insert the data `d` into the tree for `support`.
     ///
     /// Every leaf node of the tree that intersects the `support` will contain a copy of
     /// `d`.
-    fn insert<S : LocalAnalysis<F, Self::Agg, N>>(
+    fn insert<S: LocalAnalysis<F, Self::Agg, N> + Support<F, N>>(
         &mut self,
-        d : Self::Data,
-        support : &S
+        d: Self::Data,
+        support: &S,
     );
 
     /// Construct a new instance of the tree for a different aggregator
     ///
     /// The `generator` is used to convert the data of type [`Self::Data`] contained in the tree
     /// into corresponding [`Support`]s.
-    fn convert_aggregator<ANew, G>(self, generator : &G)
-    -> Self::Converted<ANew>
-    where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=Self::Data>,
-          G::SupportType : LocalAnalysis<F, ANew, N>;
-
+    fn convert_aggregator<ANew, G>(self, generator: &G) -> Self::Converted<ANew>
+    where
+        ANew: Aggregator,
+        G: SupportGenerator<F, N, Id = Self::Data>,
+        G::SupportType: LocalAnalysis<F, ANew, N>;
 
     /// Refreshes the aggregator of the three after possible changes to the support generator.
     ///
     /// The `generator` is used to convert the data of type [`Self::Data`] contained in the tree
     /// into corresponding [`Support`]s.
-    fn refresh_aggregator<G>(&mut self, generator : &G)
-    where G : SupportGenerator<F, N, Id=Self::Data>,
-          G::SupportType : LocalAnalysis<F, Self::Agg, N>;
+    fn refresh_aggregator<G>(&mut self, generator: &G)
+    where
+        G: SupportGenerator<F, N, Id = Self::Data>,
+        G::SupportType: LocalAnalysis<F, Self::Agg, N>;
 
     /// Returns an iterator over all [`Self::Data`] items at the point `x` of the domain.
-    fn iter_at(&self, x : &Loc<F,N>) -> std::slice::Iter<'_, Self::Data>;
+    fn iter_at(&self, x: &Loc<F, N>) -> std::slice::Iter<'_, Self::Data>;
 
     /*
     /// Returns all [`Self::Data`] items at the point `x` of the domain.
@@ -629,7 +663,7 @@
     */
 
     /// Create a new tree on `domain` of indicated `depth`.
-    fn new(domain : Cube<F, N>, depth : Self::Depth) -> Self;
+    fn new(domain: Cube<F, N>, depth: Self::Depth) -> Self;
 }
 
 /// The main bisection tree structure.
@@ -637,20 +671,17 @@
 /// It should be accessed via the [`BTImpl`] trait to hide the `const P : usize` parameter until
 /// const generics are flexible enough to fix `P=pow(2, N)` and thus also get rid of
 /// the `BTNodeLookup : BTNode<F, D, A, N>` trait bound.
-#[derive(Clone,Debug)]
-pub struct BT<
-    M : Depth,
-    F : Float,
-    D : 'static + Copy,
-    A : Aggregator,
-    const N : usize,
-> where BTNodeLookup : BTNode<F, D, A, N> {
+#[derive(Clone, Debug)]
+pub struct BT<M: Depth, F: Float, D: 'static + Copy, A: Aggregator, const N: usize>
+where
+    BTNodeLookup: BTNode<F, D, A, N>,
+{
     /// The depth of the tree (initial, before refinement)
-    pub(super) depth : M,
+    pub(super) depth: M,
     /// The domain of the toplevel node
-    pub(super) domain : Cube<F, N>,
+    pub(super) domain: Cube<F, N>,
     /// The toplevel node of the tree
-    pub(super) topnode : <BTNodeLookup as BTNode<F, D, A, N>>::Node,
+    pub(super) topnode: <BTNodeLookup as BTNode<F, D, A, N>>::Node,
 }
 
 macro_rules! impl_bt {
@@ -672,7 +703,7 @@
             type Agg = A;
             type Converted<ANew> = BT<M,F,D,ANew,$n> where ANew : Aggregator;
 
-            fn insert<S: LocalAnalysis<F, A, $n>>(
+            fn insert<S: LocalAnalysis<F, A, $n> + Support<F, $n>>(
                 &mut self,
                 d : D,
                 support : &S
@@ -739,4 +770,3 @@
 }
 
 impl_bt!(1 2 3 4);
-
diff -r 1f19c6bbf07b -r 2e4517b55442 src/bisection_tree/btfn.rs
--- a/src/bisection_tree/btfn.rs	Sun Apr 27 20:29:43 2025 -0500
+++ b/src/bisection_tree/btfn.rs	Sun Apr 27 20:41:36 2025 -0500
@@ -1,26 +1,25 @@
-
+use crate::mapping::{
+    BasicDecomposition, DifferentiableImpl, DifferentiableMapping, Instance, Mapping, Space,
+};
+use crate::types::Float;
 use numeric_literals::replace_float_literals;
 use std::iter::Sum;
 use std::marker::PhantomData;
 use std::sync::Arc;
-use crate::types::Float;
-use crate::mapping::{
-    Instance, Mapping, DifferentiableImpl, DifferentiableMapping, Space,
-    BasicDecomposition,
-};
 //use crate::linops::{Apply, Linear};
-use crate::sets::Set;
-use crate::sets::Cube;
-use crate::loc::Loc;
+use super::aggregator::*;
+use super::bt::*;
+use super::either::*;
+use super::refine::*;
 use super::support::*;
-use super::bt::*;
-use super::refine::*;
-use super::aggregator::*;
-use super::either::*;
+use crate::bounds::MinMaxMapping;
 use crate::fe_model::base::RealLocalModel;
 use crate::fe_model::p2_local_model::*;
+use crate::loc::Loc;
+use crate::sets::Cube;
+use crate::sets::Set;
 
-/// Presentation for (mathematical) functions constructed as a sum of components functions with 
+/// Presentation for (mathematical) functions constructed as a sum of components functions with
 /// typically small support.
 ///
 /// The domain of the function is [`Loc`]`<F, N>`, where `F` is the type of floating point numbers,
@@ -30,36 +29,29 @@
 /// Identifiers of the components ([`SupportGenerator::Id`], usually `usize`) are stored stored
 /// in a [bisection tree][BTImpl], when one is provided as `bt`. However `bt` may also be `()`
 /// for a [`PreBTFN`] that is only useful for vector space operations with a full [`BTFN`].
-#[derive(Clone,Debug)]
-pub struct BTFN<
-    F : Float,
-    G : SupportGenerator<F, N>,
-    BT /*: BTImpl<F, N>*/,
-    const N : usize
-> /*where G::SupportType : LocalAnalysis<F, A, N>*/ {
-    bt : BT,
-    generator : Arc<G>,
-    _phantoms : PhantomData<F>,
+#[derive(Clone, Debug)]
+pub struct BTFN<F: Float, G: SupportGenerator<F, N>, BT /*: BTImpl<F, N>*/, const N: usize> /*where G::SupportType : LocalAnalysis<F, A, N>*/
+{
+    bt: BT,
+    generator: Arc<G>,
+    _phantoms: PhantomData<F>,
 }
 
-impl<F : Float, G, BT, const N : usize>
-Space for BTFN<F, G, BT, N>
+impl<F: Float, G, BT, const N: usize> Space for BTFN<F, G, BT, N>
 where
-    G : SupportGenerator<F, N, Id=BT::Data>,
-    G::SupportType : LocalAnalysis<F, BT::Agg, N>,
-    BT : BTImpl<F, N>
+    G: SupportGenerator<F, N, Id = BT::Data>,
+    G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+    BT: BTImpl<F, N>,
 {
     type Decomp = BasicDecomposition;
 }
 
-impl<F : Float, G, BT, const N : usize>
-BTFN<F, G, BT, N>
+impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N>
 where
-    G : SupportGenerator<F, N, Id=BT::Data>,
-    G::SupportType : LocalAnalysis<F, BT::Agg, N>,
-    BT : BTImpl<F, N>
+    G: SupportGenerator<F, N, Id = BT::Data>,
+    G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+    BT: BTImpl<F, N>,
 {
-
     /// Create a new BTFN from a support generator and a pre-initialised bisection tree.
     ///
     /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`.
@@ -67,15 +59,15 @@
     /// when the aggregators of the tree may need updates.
     ///
     /// See the documentation for [`BTFN`] on the role of the `generator`.
-    pub fn new(bt : BT, generator : G) -> Self {
+    pub fn new(bt: BT, generator: G) -> Self {
         Self::new_arc(bt, Arc::new(generator))
     }
 
-    fn new_arc(bt : BT, generator : Arc<G>) -> Self {
+    fn new_arc(bt: BT, generator: Arc<G>) -> Self {
         BTFN {
-            bt : bt,
-            generator : generator,
-            _phantoms : std::marker::PhantomData,
+            bt: bt,
+            generator: generator,
+            _phantoms: std::marker::PhantomData,
         }
     }
 
@@ -86,7 +78,7 @@
     /// the aggregator may be out of date.
     ///
     /// See the documentation for [`BTFN`] on the role of the `generator`.
-    pub fn new_refresh(bt : &BT, generator : G) -> Self {
+    pub fn new_refresh(bt: &BT, generator: G) -> Self {
         // clone().refresh_aggregator(…) as opposed to convert_aggregator
         // ensures that type is maintained. Due to Rc-pointer copy-on-write,
         // the effort is not significantly different.
@@ -100,11 +92,11 @@
     /// The top node of the created [`BT`] will have the given `domain`.
     ///
     /// See the documentation for [`BTFN`] on the role of the `generator`.
-    pub fn construct(domain : Cube<F, N>, depth : BT::Depth, generator : G) -> Self {
+    pub fn construct(domain: Cube<F, N>, depth: BT::Depth, generator: G) -> Self {
         Self::construct_arc(domain, depth, Arc::new(generator))
     }
 
-    fn construct_arc(domain : Cube<F, N>, depth : BT::Depth, generator : Arc<G>) -> Self {
+    fn construct_arc(domain: Cube<F, N>, depth: BT::Depth, generator: Arc<G>) -> Self {
         let mut bt = BT::new(domain, depth);
         for (d, support) in generator.all_data() {
             bt.insert(d, &support);
@@ -117,14 +109,16 @@
     /// This will construct a [`BTFN`] with the same components and generator as the (consumed)
     /// `self`, but a new `BT` with [`Aggregator`]s of type `ANew`.
     pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N>
-    where ANew : Aggregator,
-          G : SupportGenerator<F, N, Id=BT::Data>,
-          G::SupportType : LocalAnalysis<F, ANew, N> {
+    where
+        ANew: Aggregator,
+        G: SupportGenerator<F, N, Id = BT::Data>,
+        G::SupportType: LocalAnalysis<F, ANew, N>,
+    {
         BTFN::new_arc(self.bt.convert_aggregator(&*self.generator), self.generator)
     }
 
     /// Change the generator (after, e.g., a scaling of the latter).
-    fn new_generator(&self, generator : G) -> Self {
+    fn new_generator(&self, generator: G) -> Self {
         BTFN::new_refresh(&self.bt, generator)
     }
 
@@ -132,20 +126,24 @@
     fn refresh_aggregator(&mut self) {
         self.bt.refresh_aggregator(&*self.generator);
     }
-
 }
 
-impl<F : Float, G, BT, const N : usize>
-BTFN<F, G, BT, N>
-where G : SupportGenerator<F, N> {
+impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N>
+where
+    G: SupportGenerator<F, N>,
+{
     /// Change the [bisection tree][BTImpl] of the [`BTFN`] to a different one.
     ///
     /// This can be used to convert a [`PreBTFN`] to a full [`BTFN`], or the change
     /// the aggreagator; see also [`self.convert_aggregator`].
-    pub fn instantiate<
-        BTNew : BTImpl<F, N, Data=G::Id>,
-    > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<F, G, BTNew, N>
-    where G::SupportType : LocalAnalysis<F, BTNew::Agg, N>  {
+    pub fn instantiate<BTNew: BTImpl<F, N, Data = G::Id>>(
+        self,
+        domain: Cube<F, N>,
+        depth: BTNew::Depth,
+    ) -> BTFN<F, G, BTNew, N>
+    where
+        G::SupportType: LocalAnalysis<F, BTNew::Agg, N>,
+    {
         BTFN::construct_arc(domain, depth, self.generator)
     }
 }
@@ -155,31 +153,34 @@
 /// Most BTFN methods are not available, but if a BTFN is going to be summed with another
 /// before other use, it will be more efficient to not construct an unnecessary bisection tree
 /// that would be shortly dropped.
-pub type PreBTFN<F, G, const N : usize> = BTFN<F, G, (), N>;
+pub type PreBTFN<F, G, const N: usize> = BTFN<F, G, (), N>;
 
-impl<F : Float, G, const N : usize> PreBTFN<F, G, N> where G : SupportGenerator<F, N> {
-
+impl<F: Float, G, const N: usize> PreBTFN<F, G, N>
+where
+    G: SupportGenerator<F, N>,
+{
     /// Create a new [`PreBTFN`] with no bisection tree.
-    pub fn new_pre(generator : G) -> Self {
+    pub fn new_pre(generator: G) -> Self {
         BTFN {
-            bt : (),
-            generator : Arc::new(generator),
-            _phantoms : std::marker::PhantomData,
+            bt: (),
+            generator: Arc::new(generator),
+            _phantoms: std::marker::PhantomData,
         }
     }
 }
 
-impl<F : Float, G, BT, const N : usize>
-BTFN<F, G, BT, N>
-where G : SupportGenerator<F, N, Id=usize>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N>,
-      BT : BTImpl<F, N, Data=usize> {
-
+impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N>
+where
+    G: SupportGenerator<F, N, Id = usize>,
+    G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+    BT: BTImpl<F, N, Data = usize>,
+{
     /// Helper function for implementing [`std::ops::Add`].
-    fn add_another<G2>(&self, g2 : Arc<G2>) -> BTFN<F, BothGenerators<G, G2>, BT, N>
-    where G2 : SupportGenerator<F, N, Id=usize>,
-          G2::SupportType : LocalAnalysis<F, BT::Agg, N> {
-
+    fn add_another<G2>(&self, g2: Arc<G2>) -> BTFN<F, BothGenerators<G, G2>, BT, N>
+    where
+        G2: SupportGenerator<F, N, Id = usize>,
+        G2::SupportType: LocalAnalysis<F, BT::Agg, N>,
+    {
         let mut bt = self.bt.clone();
         let both = BothGenerators(Arc::clone(&self.generator), g2);
 
@@ -188,9 +189,9 @@
         }
 
         BTFN {
-            bt : bt,
-            generator : Arc::new(both),
-            _phantoms : std::marker::PhantomData,
+            bt: bt,
+            generator: Arc::new(both),
+            _phantoms: std::marker::PhantomData,
         }
     }
 }
@@ -231,7 +232,7 @@
 }
 
 make_btfn_add!(BTFN<F, G1, BT1, N>, std::convert::identity, );
-make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone, );
+make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone,);
 
 macro_rules! make_btfn_sub {
     ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => {
@@ -274,52 +275,52 @@
 }
 
 make_btfn_sub!(BTFN<F, G1, BT1, N>, std::convert::identity, );
-make_btfn_sub!(&'a BTFN<F, G1, BT1, N>, std::convert::identity, );
+make_btfn_sub!(&'a BTFN<F, G1, BT1, N>, std::convert::identity,);
 
 macro_rules! make_btfn_scalarop_rhs {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
-        impl<F : Float, G, BT, const N : usize>
-        std::ops::$trait_assign<F>
-        for BTFN<F, G, BT, N>
-        where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N> {
+        impl<F: Float, G, BT, const N: usize> std::ops::$trait_assign<F> for BTFN<F, G, BT, N>
+        where
+            BT: BTImpl<F, N>,
+            G: SupportGenerator<F, N, Id = BT::Data>,
+            G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+        {
             #[inline]
-            fn $fn_assign(&mut self, t : F) {
+            fn $fn_assign(&mut self, t: F) {
                 Arc::make_mut(&mut self.generator).$fn_assign(t);
                 self.refresh_aggregator();
             }
         }
 
-        impl<F : Float, G, BT, const N : usize>
-        std::ops::$trait<F>
-        for BTFN<F, G, BT, N>
-        where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N> {
+        impl<F: Float, G, BT, const N: usize> std::ops::$trait<F> for BTFN<F, G, BT, N>
+        where
+            BT: BTImpl<F, N>,
+            G: SupportGenerator<F, N, Id = BT::Data>,
+            G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+        {
             type Output = Self;
             #[inline]
-            fn $fn(mut self, t : F) -> Self::Output {
+            fn $fn(mut self, t: F) -> Self::Output {
                 Arc::make_mut(&mut self.generator).$fn_assign(t);
                 self.refresh_aggregator();
                 self
             }
         }
 
-        impl<'a, F : Float, G, BT, const N : usize>
-        std::ops::$trait<F>
-        for &'a BTFN<F, G, BT, N>
-        where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N>,
-              &'a G : std::ops::$trait<F,Output=G> {
+        impl<'a, F: Float, G, BT, const N: usize> std::ops::$trait<F> for &'a BTFN<F, G, BT, N>
+        where
+            BT: BTImpl<F, N>,
+            G: SupportGenerator<F, N, Id = BT::Data>,
+            G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+            &'a G: std::ops::$trait<F, Output = G>,
+        {
             type Output = BTFN<F, G, BT, N>;
             #[inline]
-            fn $fn(self, t : F) -> Self::Output {
+            fn $fn(self, t: F) -> Self::Output {
                 self.new_generator(self.generator.$fn(t))
             }
         }
-    }
+    };
 }
 
 make_btfn_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
@@ -368,12 +369,12 @@
 
 macro_rules! make_btfn_unaryop {
     ($trait:ident, $fn:ident) => {
-        impl<F : Float, G, BT, const N : usize>
-        std::ops::$trait
-        for BTFN<F, G, BT, N>
-        where BT : BTImpl<F, N>,
-              G : SupportGenerator<F, N, Id=BT::Data>,
-              G::SupportType : LocalAnalysis<F, BT::Agg, N> {
+        impl<F: Float, G, BT, const N: usize> std::ops::$trait for BTFN<F, G, BT, N>
+        where
+            BT: BTImpl<F, N>,
+            G: SupportGenerator<F, N, Id = BT::Data>,
+            G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+        {
             type Output = Self;
             #[inline]
             fn $fn(mut self) -> Self::Output {
@@ -396,7 +397,7 @@
                 self.new_generator(std::ops::$trait::$fn(&self.generator))
             }
         }*/
-    }
+    };
 }
 
 make_btfn_unaryop!(Neg, neg);
@@ -405,39 +406,38 @@
 // Apply, Mapping, Differentiate
 //
 
-impl<F : Float, G, BT, V, const N : usize> Mapping<Loc<F, N>>
-for BTFN<F, G, BT, N>
+impl<F: Float, G, BT, V, const N: usize> Mapping<Loc<F, N>> for BTFN<F, G, BT, N>
 where
-    BT : BTImpl<F, N>,
-    G : SupportGenerator<F, N, Id=BT::Data>,
-    G::SupportType : LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F, N>, Codomain = V>,
-    V : Sum + Space,
+    BT: BTImpl<F, N>,
+    G: SupportGenerator<F, N, Id = BT::Data>,
+    G::SupportType: LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F, N>, Codomain = V>,
+    V: Sum + Space,
 {
-
     type Codomain = V;
 
-    fn apply<I : Instance<Loc<F,N>>>(&self, x : I) -> Self::Codomain {
+    fn apply<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Codomain {
         let xc = x.cow();
-        self.bt.iter_at(&*xc)
-            .map(|&d| self.generator.support_for(d).apply(&*xc)).sum()
+        self.bt
+            .iter_at(&*xc)
+            .map(|&d| self.generator.support_for(d).apply(&*xc))
+            .sum()
     }
 }
 
-impl<F : Float, G, BT, V, const N : usize> DifferentiableImpl<Loc<F, N>>
-for BTFN<F, G, BT, N>
+impl<F: Float, G, BT, V, const N: usize> DifferentiableImpl<Loc<F, N>> for BTFN<F, G, BT, N>
 where
-    BT : BTImpl<F, N>,
-    G : SupportGenerator<F, N, Id=BT::Data>,
-    G::SupportType : LocalAnalysis<F, BT::Agg, N>
-        + DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>,
-    V : Sum + Space,
+    BT: BTImpl<F, N>,
+    G: SupportGenerator<F, N, Id = BT::Data>,
+    G::SupportType:
+        LocalAnalysis<F, BT::Agg, N> + DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>,
+    V: Sum + Space,
 {
-
     type Derivative = V;
 
-    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x :I) -> Self::Derivative {
+    fn differential_impl<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Derivative {
         let xc = x.cow();
-        self.bt.iter_at(&*xc)
+        self.bt
+            .iter_at(&*xc)
             .map(|&d| self.generator.support_for(d).differential(&*xc))
             .sum()
     }
@@ -447,12 +447,12 @@
 // GlobalAnalysis
 //
 
-impl<F : Float, G, BT, const N : usize> GlobalAnalysis<F, BT::Agg>
-for BTFN<F, G, BT, N>
-where BT : BTImpl<F, N>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : LocalAnalysis<F, BT::Agg, N> {
-
+impl<F: Float, G, BT, const N: usize> GlobalAnalysis<F, BT::Agg> for BTFN<F, G, BT, N>
+where
+    BT: BTImpl<F, N>,
+    G: SupportGenerator<F, N, Id = BT::Data>,
+    G::SupportType: LocalAnalysis<F, BT::Agg, N>,
+{
     #[inline]
     fn global_analysis(&self) -> BT::Agg {
         self.bt.global_analysis()
@@ -505,24 +505,23 @@
 ///
 /// `U` is the domain, generally [`Loc`]`<F, N>`, and `F` the type of floating point numbers.
 /// `Self` is generally a set of `U`, for example, [`Cube`]`<F, N>`.
-pub trait P2Minimise<U : Space, F : Float> : Set<U> {
+pub trait P2Minimise<U: Space, F: Float>: Set<U> {
     /// Minimise `g` over the set presented by `Self`.
     ///
     /// The function returns `(x, v)` where `x` is the minimiser `v` an approximation of `g(x)`.
-    fn p2_minimise<G : Fn(&U) -> F>(&self, g : G) -> (U, F);
-
+    fn p2_minimise<G: Fn(&U) -> F>(&self, g: G) -> (U, F);
 }
 
-impl<F : Float> P2Minimise<Loc<F, 1>, F> for Cube<F, 1> {
-    fn p2_minimise<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> (Loc<F, 1>, F) {
+impl<F: Float> P2Minimise<Loc<F, 1>, F> for Cube<F, 1> {
+    fn p2_minimise<G: Fn(&Loc<F, 1>) -> F>(&self, g: G) -> (Loc<F, 1>, F) {
         let interval = Simplex(self.corners());
         interval.p2_model(&g).minimise(&interval)
     }
 }
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Float> P2Minimise<Loc<F, 2>, F> for Cube<F, 2> {
-    fn p2_minimise<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> (Loc<F, 2>, F) {
+impl<F: Float> P2Minimise<Loc<F, 2>, F> for Cube<F, 2> {
+    fn p2_minimise<G: Fn(&Loc<F, 2>) -> F>(&self, g: G) -> (Loc<F, 2>, F) {
         if false {
             // Split into two triangle (simplex) with separate P2 model in each.
             // The six nodes of each triangle are the corners and the edges.
@@ -537,49 +536,46 @@
             let [vab, vbc, vca, vcd, vda] = [g(&ab), g(&bc), g(&ca), g(&cd), g(&da)];
 
             let s1 = Simplex([a, b, c]);
-            let m1 = P2LocalModel::<F, 2, 3>::new(
-                &[a, b, c, ab, bc, ca],
-                &[va, vb, vc, vab, vbc, vca]
-            );
+            let m1 =
+                P2LocalModel::<F, 2, 3>::new(&[a, b, c, ab, bc, ca], &[va, vb, vc, vab, vbc, vca]);
 
-            let r1@(_, v1) = m1.minimise(&s1);
+            let r1 @ (_, v1) = m1.minimise(&s1);
 
             let s2 = Simplex([c, d, a]);
-            let m2 = P2LocalModel::<F, 2, 3>::new(
-                &[c, d, a, cd, da, ca],
-                &[vc, vd, va, vcd, vda, vca]
-            );
+            let m2 =
+                P2LocalModel::<F, 2, 3>::new(&[c, d, a, cd, da, ca], &[vc, vd, va, vcd, vda, vca]);
+
+            let r2 @ (_, v2) = m2.minimise(&s2);
 
-            let r2@(_, v2) = m2.minimise(&s2);
-
-            if v1 < v2 { r1 } else { r2 }
+            if v1 < v2 {
+                r1
+            } else {
+                r2
+            }
         } else {
             // Single P2 model for the entire cube.
             let [a, b, c, d] = self.corners();
             let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)];
             let [e, f] = match 'r' {
-                 'm' => [(&a + &b + &c) / 3.0,    (&c + &d + &a) / 3.0],
-                 'c' => [midpoint(&a, &b),        midpoint(&a, &d)],
-                 'w' => [(&a + &b * 2.0) / 3.0,   (&a + &d * 2.0) / 3.0],
-                 'r' => {
+                'm' => [(&a + &b + &c) / 3.0, (&c + &d + &a) / 3.0],
+                'c' => [midpoint(&a, &b), midpoint(&a, &d)],
+                'w' => [(&a + &b * 2.0) / 3.0, (&a + &d * 2.0) / 3.0],
+                'r' => {
                     // Pseudo-randomise edge midpoints
                     let Loc([x, y]) = a;
-                    let tmp : f64 = (x+y).as_();
+                    let tmp: f64 = (x + y).as_();
                     match tmp.to_bits() % 4 {
-                        0 => [midpoint(&a, &b),        midpoint(&a, &d)],
-                        1 => [midpoint(&c, &d),        midpoint(&a, &d)],
-                        2 => [midpoint(&a, &b),        midpoint(&b, &c)],
-                        _ => [midpoint(&c, &d),        midpoint(&b, &c)],
+                        0 => [midpoint(&a, &b), midpoint(&a, &d)],
+                        1 => [midpoint(&c, &d), midpoint(&a, &d)],
+                        2 => [midpoint(&a, &b), midpoint(&b, &c)],
+                        _ => [midpoint(&c, &d), midpoint(&b, &c)],
                     }
-                 },
-                 _ => [self.center(),           (&a + &b) / 2.0],
+                }
+                _ => [self.center(), (&a + &b) / 2.0],
             };
             let [ve, vf] = [g(&e), g(&f)];
 
-            let m1 = P2LocalModel::<F, 2, 3>::new(
-                &[a, b, c, d, e, f],
-                &[va, vb, vc, vd, ve, vf],
-            );
+            let m1 = P2LocalModel::<F, 2, 3>::new(&[a, b, c, d, e, f], &[va, vb, vc, vd, ve, vf]);
 
             m1.minimise(self)
         }
@@ -595,44 +591,46 @@
 /// A bisection tree [`Refiner`] for maximising or minimising a [`BTFN`].
 ///
 /// The type parameter `T` should be either [`RefineMax`] or [`RefineMin`].
-struct P2Refiner<F : Float, T> {
+struct P2Refiner<F: Float, T> {
     /// The maximum / minimum should be above / below this threshold.
     /// If the threshold cannot be satisfied, the refiner will return `None`.
-    bound : Option<F>,
+    bound: Option<F>,
     /// Tolerance for function value estimation.
-    tolerance : F,
+    tolerance: F,
     /// Maximum number of steps to execute the refiner for
-    max_steps : usize,
+    max_steps: usize,
     /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes.
     #[allow(dead_code)] // `how` is just for type system purposes.
-    how : T,
+    how: T,
 }
 
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
-for P2Refiner<F, RefineMax>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
-      G : SupportGenerator<F, N>,
-      G::SupportType : Mapping<Loc<F, N>, Codomain=F>
-                       + LocalAnalysis<F, Bounds<F>, N> {
+impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for P2Refiner<F, RefineMax>
+where
+    Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+    G: SupportGenerator<F, N>,
+    G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+{
     type Result = Option<(Loc<F, N>, F)>;
     type Sorting = UpperBoundSorting<F>;
 
     fn refine(
         &self,
-        aggregator : &Bounds<F>,
-        cube : &Cube<F, N>,
-        data : &[G::Id],
-        generator : &G,
-        step : usize
+        aggregator: &Bounds<F>,
+        cube: &Cube<F, N>,
+        data: &[G::Id],
+        generator: &G,
+        step: usize,
     ) -> RefinerResult<Bounds<F>, Self::Result> {
-
-        if self.bound.map_or(false, |b| aggregator.upper() <= b + self.tolerance) {
+        if self
+            .bound
+            .map_or(false, |b| aggregator.upper() <= b + self.tolerance)
+        {
             // The upper bound is below the maximisation threshold. Don't bother with this cube.
-            return RefinerResult::Uncertain(*aggregator, None)
+            return RefinerResult::Uncertain(*aggregator, None);
         }
 
         // g gives the negative of the value of the function presented by `data` and `generator`.
-        let g = move |x : &Loc<F, N>| {
+        let g = move |x: &Loc<F, N>| {
             let f = move |&d| generator.support_for(d).apply(x);
             -data.iter().map(f).sum::<F>()
         };
@@ -641,8 +639,9 @@
         //let v = -neg_v;
         let v = -g(&x);
 
-        if step < self.max_steps && (aggregator.upper() > v + self.tolerance
-                                     /*|| aggregator.lower() > v - self.tolerance*/) {
+        if step < self.max_steps
+            && (aggregator.upper() > v + self.tolerance/*|| aggregator.lower() > v - self.tolerance*/)
+        {
             // The function isn't refined enough in `cube`, so return None
             // to indicate that further subdivision is required.
             RefinerResult::NeedRefinement
@@ -655,41 +654,46 @@
         }
     }
 
-    fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) {
+    fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
         match (*r1, r2) {
-            (Some((_, v1)), Some((_, v2))) => if v1 < v2 { *r1 = r2 }
+            (Some((_, v1)), Some((_, v2))) => {
+                if v1 < v2 {
+                    *r1 = r2
+                }
+            }
             (None, Some(_)) => *r1 = r2,
-            (_, _) => {},
+            (_, _) => {}
         }
     }
 }
 
-
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
-for P2Refiner<F, RefineMin>
-where Cube<F, N> : P2Minimise<Loc<F, N>, F>,
-      G : SupportGenerator<F, N>,
-      G::SupportType : Mapping<Loc<F, N>, Codomain=F>
-                       + LocalAnalysis<F, Bounds<F>, N> {
+impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for P2Refiner<F, RefineMin>
+where
+    Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+    G: SupportGenerator<F, N>,
+    G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+{
     type Result = Option<(Loc<F, N>, F)>;
     type Sorting = LowerBoundSorting<F>;
 
     fn refine(
         &self,
-        aggregator : &Bounds<F>,
-        cube : &Cube<F, N>,
-        data : &[G::Id],
-        generator : &G,
-        step : usize
+        aggregator: &Bounds<F>,
+        cube: &Cube<F, N>,
+        data: &[G::Id],
+        generator: &G,
+        step: usize,
     ) -> RefinerResult<Bounds<F>, Self::Result> {
-    
-        if self.bound.map_or(false, |b| aggregator.lower() >= b - self.tolerance) {
+        if self
+            .bound
+            .map_or(false, |b| aggregator.lower() >= b - self.tolerance)
+        {
             // The lower bound is above the minimisation threshold. Don't bother with this cube.
-            return RefinerResult::Uncertain(*aggregator, None)
+            return RefinerResult::Uncertain(*aggregator, None);
         }
 
         // g gives the value of the function presented by `data` and `generator`.
-        let g = move |x : &Loc<F, N>| {
+        let g = move |x: &Loc<F, N>| {
             let f = move |&d| generator.support_for(d).apply(x);
             data.iter().map(f).sum::<F>()
         };
@@ -697,8 +701,9 @@
         let (x, _v) = cube.p2_minimise(g);
         let v = g(&x);
 
-         if step < self.max_steps && (aggregator.lower() < v - self.tolerance
-                                      /*|| aggregator.upper() < v + self.tolerance*/) {
+        if step < self.max_steps
+            && (aggregator.lower() < v - self.tolerance/*|| aggregator.upper() < v + self.tolerance*/)
+        {
             // The function isn't refined enough in `cube`, so return None
             // to indicate that further subdivision is required.
             RefinerResult::NeedRefinement
@@ -717,47 +722,51 @@
         }
     }
 
-    fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) {
+    fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
         match (*r1, r2) {
-            (Some((_, v1)), Some((_, v2))) => if v1 > v2 { *r1 = r2 }
+            (Some((_, v1)), Some((_, v2))) => {
+                if v1 > v2 {
+                    *r1 = r2
+                }
+            }
             (_, Some(_)) => *r1 = r2,
-            (_, _) => {},
+            (_, _) => {}
         }
     }
 }
 
-
 /// A bisection tree [`Refiner`] for checking that a [`BTFN`] is within a stated
 //// upper or lower bound.
 ///
 /// The type parameter `T` should be either [`RefineMax`] for upper bound or [`RefineMin`]
 /// for lower bound.
 
-struct BoundRefiner<F : Float, T> {
+struct BoundRefiner<F: Float, T> {
     /// The upper/lower bound to check for
-    bound : F,
+    bound: F,
     /// Tolerance for function value estimation.
-    tolerance : F,
+    tolerance: F,
     /// Maximum number of steps to execute the refiner for
-    max_steps : usize,
+    max_steps: usize,
     #[allow(dead_code)] // `how` is just for type system purposes.
     /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes.
-    how : T,
+    how: T,
 }
 
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
-for BoundRefiner<F, RefineMax>
-where G : SupportGenerator<F, N> {
+impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for BoundRefiner<F, RefineMax>
+where
+    G: SupportGenerator<F, N>,
+{
     type Result = bool;
     type Sorting = UpperBoundSorting<F>;
 
     fn refine(
         &self,
-        aggregator : &Bounds<F>,
-        _cube : &Cube<F, N>,
-        _data : &[G::Id],
-        _generator : &G,
-        step : usize
+        aggregator: &Bounds<F>,
+        _cube: &Cube<F, N>,
+        _data: &[G::Id],
+        _generator: &G,
+        step: usize,
     ) -> RefinerResult<Bounds<F>, Self::Result> {
         if aggregator.upper() <= self.bound + self.tolerance {
             // Below upper bound within tolerances. Indicate uncertain success.
@@ -774,24 +783,25 @@
         }
     }
 
-    fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) {
+    fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
         *r1 = *r1 && r2;
     }
 }
 
-impl<F : Float, G, const N : usize> Refiner<F, Bounds<F>, G, N>
-for BoundRefiner<F, RefineMin>
-where G : SupportGenerator<F, N> {
+impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for BoundRefiner<F, RefineMin>
+where
+    G: SupportGenerator<F, N>,
+{
     type Result = bool;
     type Sorting = UpperBoundSorting<F>;
 
     fn refine(
         &self,
-        aggregator : &Bounds<F>,
-        _cube : &Cube<F, N>,
-        _data : &[G::Id],
-        _generator : &G,
-        step : usize
+        aggregator: &Bounds<F>,
+        _cube: &Cube<F, N>,
+        _data: &[G::Id],
+        _generator: &G,
+        step: usize,
     ) -> RefinerResult<Bounds<F>, Self::Result> {
         if aggregator.lower() >= self.bound - self.tolerance {
             // Above lower bound within tolerances. Indicate uncertain success.
@@ -808,7 +818,7 @@
         }
     }
 
-    fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) {
+    fn fuse_results(r1: &mut Self::Result, r2: Self::Result) {
         *r1 = *r1 && r2;
     }
 }
@@ -828,66 +838,94 @@
 // there should be a result, or new nodes above the `glb` inserted into the queue. Then the waiting
 // threads can also continue processing. If, however, numerical inaccuracy destroyes the `glb`,
 // the queue may run out, and we get “Refiner failure”.
-impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N>
-where BT : BTSearch<F, N, Agg=Bounds<F>>,
-      G : SupportGenerator<F, N, Id=BT::Data>,
-      G::SupportType : Mapping<Loc<F, N>, Codomain=F>
-                       + LocalAnalysis<F, Bounds<F>, N>,
-      Cube<F, N> : P2Minimise<Loc<F, N>, F> {
-
-    /// Maximise the `BTFN` within stated value `tolerance`.
-    ///
-    /// At most `max_steps` refinement steps are taken.
-    /// Returns the approximate maximiser and the corresponding function value.
-    pub fn maximise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) {
-        let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : None };
-        self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.").unwrap()
+impl<F: Float, G, BT, const N: usize> MinMaxMapping<F, N> for BTFN<F, G, BT, N>
+where
+    BT: BTSearch<F, N, Agg = Bounds<F>>,
+    G: SupportGenerator<F, N, Id = BT::Data>,
+    G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>,
+    Cube<F, N>: P2Minimise<Loc<F, N>, F>,
+{
+    fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F) {
+        let refiner = P2Refiner {
+            tolerance,
+            max_steps,
+            how: RefineMax,
+            bound: None,
+        };
+        self.bt
+            .search_and_refine(refiner, &self.generator)
+            .expect("Refiner failure.")
+            .unwrap()
     }
 
-    /// Maximise the `BTFN` within stated value `tolerance` subject to a lower bound.
-    ///
-    /// At most `max_steps` refinement steps are taken.
-    /// Returns the approximate maximiser and the corresponding function value when one is found
-    /// above the `bound` threshold, otherwise `None`.
-    pub fn maximise_above(&mut self, bound : F, tolerance : F, max_steps : usize)
-    -> Option<(Loc<F, N>, F)> {
-        let refiner = P2Refiner{ tolerance, max_steps, how : RefineMax, bound : Some(bound) };
-        self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.")
+    fn maximise_above(
+        &mut self,
+        bound: F,
+        tolerance: F,
+        max_steps: usize,
+    ) -> Option<(Loc<F, N>, F)> {
+        let refiner = P2Refiner {
+            tolerance,
+            max_steps,
+            how: RefineMax,
+            bound: Some(bound),
+        };
+        self.bt
+            .search_and_refine(refiner, &self.generator)
+            .expect("Refiner failure.")
     }
 
-    /// Minimise the `BTFN` within stated value `tolerance`.
-    ///
-    /// At most `max_steps` refinement steps are taken.
-    /// Returns the approximate minimiser and the corresponding function value.
-    pub fn minimise(&mut self, tolerance : F, max_steps : usize) -> (Loc<F, N>, F) {
-        let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : None };
-        self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.").unwrap()
+    fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F) {
+        let refiner = P2Refiner {
+            tolerance,
+            max_steps,
+            how: RefineMin,
+            bound: None,
+        };
+        self.bt
+            .search_and_refine(refiner, &self.generator)
+            .expect("Refiner failure.")
+            .unwrap()
     }
 
-    /// Minimise the `BTFN` within stated value `tolerance` subject to a lower bound.
-    ///
-    /// At most `max_steps` refinement steps are taken.
-    /// Returns the approximate minimiser and the corresponding function value when one is found
-    /// above the `bound` threshold, otherwise `None`.
-    pub fn minimise_below(&mut self, bound : F, tolerance : F, max_steps : usize)
-    -> Option<(Loc<F, N>, F)> {
-        let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : Some(bound) };
-        self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.")
+    fn minimise_below(
+        &mut self,
+        bound: F,
+        tolerance: F,
+        max_steps: usize,
+    ) -> Option<(Loc<F, N>, F)> {
+        let refiner = P2Refiner {
+            tolerance,
+            max_steps,
+            how: RefineMin,
+            bound: Some(bound),
+        };
+        self.bt
+            .search_and_refine(refiner, &self.generator)
+            .expect("Refiner failure.")
     }
 
-    /// Verify that the `BTFN` has a given upper `bound` within indicated `tolerance`.
-    ///
-    /// At most `max_steps` refinement steps are taken.
-    pub fn has_upper_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool {
-        let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMax };
-        self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.")
+    fn has_upper_bound(&mut self, bound: F, tolerance: F, max_steps: usize) -> bool {
+        let refiner = BoundRefiner {
+            bound,
+            tolerance,
+            max_steps,
+            how: RefineMax,
+        };
+        self.bt
+            .search_and_refine(refiner, &self.generator)
+            .expect("Refiner failure.")
     }
 
-    /// Verify that the `BTFN` has a given lower `bound` within indicated `tolerance`.
-    ///
-    /// At most `max_steps` refinement steps are taken.
-    pub fn has_lower_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool {
-        let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMin };
-        self.bt.search_and_refine(refiner, &self.generator).expect("Refiner failure.")
+    fn has_lower_bound(&mut self, bound: F, tolerance: F, max_steps: usize) -> bool {
+        let refiner = BoundRefiner {
+            bound,
+            tolerance,
+            max_steps,
+            how: RefineMin,
+        };
+        self.bt
+            .search_and_refine(refiner, &self.generator)
+            .expect("Refiner failure.")
     }
 }
diff -r 1f19c6bbf07b -r 2e4517b55442 src/bisection_tree/support.rs
--- a/src/bisection_tree/support.rs	Sun Apr 27 20:29:43 2025 -0500
+++ b/src/bisection_tree/support.rs	Sun Apr 27 20:41:36 2025 -0500
@@ -1,31 +1,29 @@
-
 /*!
 Traits for representing the support of a [`Mapping`], and analysing the mapping on a [`Cube`].
 */
-use serde::Serialize;
-use std::ops::{MulAssign,DivAssign,Neg};
-use crate::types::{Float, Num};
+use super::aggregator::Bounds;
+pub use crate::bounds::{GlobalAnalysis, LocalAnalysis};
+use crate::loc::Loc;
+use crate::mapping::{DifferentiableImpl, DifferentiableMapping, Instance, Mapping, Space};
 use crate::maputil::map2;
-use crate::mapping::{
-    Instance, Mapping, DifferentiableImpl, DifferentiableMapping, Space
-};
+use crate::norms::{Linfinity, Norm, L1, L2};
+pub use crate::operator_arithmetic::{Constant, Weighted};
 use crate::sets::Cube;
-use crate::loc::Loc;
-use super::aggregator::Bounds;
-use crate::norms::{Norm, L1, L2, Linfinity};
-pub use crate::operator_arithmetic::{Weighted, Constant};
+use crate::types::{Float, Num};
+use serde::Serialize;
+use std::ops::{DivAssign, MulAssign, Neg};
 
 /// A trait for working with the supports of [`Mapping`]s.
 ///
 /// `Mapping` is not a super-trait to allow more general use.
-pub trait Support<F : Num, const N : usize> : Sized + Sync + Send + 'static {
+pub trait Support<F: Num, const N: usize>: Sized + Sync + Send + 'static {
     /// Return a cube containing the support of the function represented by `self`.
     ///
     /// The hint may be larger than the actual support, but must contain it.
-    fn support_hint(&self) -> Cube<F,N>;
+    fn support_hint(&self) -> Cube<F, N>;
 
     /// Indicate whether `x` is in the support of the function represented by `self`.
-    fn in_support(&self, x : &Loc<F,N>) -> bool;
+    fn in_support(&self, x: &Loc<F, N>) -> bool;
 
     // Indicate whether `cube` is fully in the support of the function represented by `self`.
     //fn fully_in_support(&self, cube : &Cube<F,N>) -> bool;
@@ -41,131 +39,93 @@
     /// The default implementation returns `[None; N]`.
     #[inline]
     #[allow(unused_variables)]
-    fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] {
         [None; N]
     }
 
     /// Translate `self` by `x`.
     #[inline]
-    fn shift(self, x : Loc<F, N>) -> Shift<Self, F, N> {
-        Shift { shift : x, base_fn : self }
+    fn shift(self, x: Loc<F, N>) -> Shift<Self, F, N> {
+        Shift {
+            shift: x,
+            base_fn: self,
+        }
     }
 }
 
-/// Trait for globally analysing a property `A` of a [`Mapping`].
-///
-/// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as
-/// [`Bounds`][super::aggregator::Bounds].
-pub trait GlobalAnalysis<F : Num, A> {
-    /// Perform global analysis of the property `A` of `Self`.
-    ///
-    /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds],
-    /// this function will return global upper and lower bounds for the mapping
-    /// represented by `self`.
-    fn global_analysis(&self) -> A;
+/// Shift of [`Support`] and [`Mapping`]; output of [`Support::shift`].
+#[derive(Copy, Clone, Debug, Serialize)] // Serialize! but not implemented by Loc.
+pub struct Shift<T, F, const N: usize> {
+    shift: Loc<F, N>,
+    base_fn: T,
 }
 
-// default impl<F, A, N, L> GlobalAnalysis<F, A, N> for L
-// where L : LocalAnalysis<F, A, N> {
-//     #[inline]
-//     fn global_analysis(&self) -> Bounds<F> {
-//         self.local_analysis(&self.support_hint())
-//     }
-// }
-
-/// Trait for locally analysing a property `A` of a [`Mapping`] (implementing [`Support`])
-/// within a [`Cube`].
-///
-/// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as
-/// [`Bounds`][super::aggregator::Bounds].
-pub trait LocalAnalysis<F : Num, A, const N : usize> : GlobalAnalysis<F, A> + Support<F, N> {
-    /// Perform local analysis of the property `A` of `Self`.
-    ///
-    /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds],
-    /// this function will return upper and lower bounds within `cube` for the mapping
-    /// represented by `self`.
-    fn local_analysis(&self, cube : &Cube<F, N>) -> A;
-}
-
-/// Trait for determining the upper and lower bounds of an float-valued [`Mapping`].
-///
-/// This is a blanket-implemented alias for [`GlobalAnalysis`]`<F, Bounds<F>>`
-/// [`Mapping`] is not a supertrait to allow flexibility in the implementation of either
-/// reference or non-reference arguments.
-pub trait Bounded<F : Float> : GlobalAnalysis<F, Bounds<F>> {
-    /// Return lower and upper bounds for the values of of `self`.
-    #[inline]
-    fn bounds(&self) -> Bounds<F> {
-        self.global_analysis()
-    }
-}
-
-impl<F : Float, T : GlobalAnalysis<F, Bounds<F>>> Bounded<F> for T { }
-
-/// Shift of [`Support`] and [`Mapping`]; output of [`Support::shift`].
-#[derive(Copy,Clone,Debug,Serialize)] // Serialize! but not implemented by Loc.
-pub struct Shift<T, F, const N : usize> {
-    shift : Loc<F, N>,
-    base_fn : T,
-}
-
-impl<'a, T, V : Space, F : Float, const N : usize> Mapping<Loc<F, N>> for Shift<T,F,N>
-where T : Mapping<Loc<F, N>, Codomain=V> {
+impl<'a, T, V: Space, F: Float, const N: usize> Mapping<Loc<F, N>> for Shift<T, F, N>
+where
+    T: Mapping<Loc<F, N>, Codomain = V>,
+{
     type Codomain = V;
 
     #[inline]
-    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+    fn apply<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Codomain {
         self.base_fn.apply(x.own() - &self.shift)
     }
 }
 
-impl<'a, T, V : Space, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>> for Shift<T,F,N>
-where T : DifferentiableMapping<Loc<F, N>, DerivativeDomain=V> {
+impl<'a, T, V: Space, F: Float, const N: usize> DifferentiableImpl<Loc<F, N>> for Shift<T, F, N>
+where
+    T: DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>,
+{
     type Derivative = V;
 
     #[inline]
-    fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
+    fn differential_impl<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Derivative {
         self.base_fn.differential(x.own() - &self.shift)
     }
 }
 
-impl<'a, T, F : Float, const N : usize> Support<F,N> for Shift<T,F,N>
-where T : Support<F, N> {
+impl<'a, T, F: Float, const N: usize> Support<F, N> for Shift<T, F, N>
+where
+    T: Support<F, N>,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F,N> {
+    fn support_hint(&self) -> Cube<F, N> {
         self.base_fn.support_hint().shift(&self.shift)
     }
 
     #[inline]
-    fn in_support(&self, x : &Loc<F,N>) -> bool {
+    fn in_support(&self, x: &Loc<F, N>) -> bool {
         self.base_fn.in_support(&(x - &self.shift))
     }
-    
+
     // fn fully_in_support(&self, _cube : &Cube<F,N>) -> bool {
     //     //self.base_fn.fully_in_support(cube.shift(&vectorneg(self.shift)))
     //     todo!("Not implemented, but not used at the moment")
     // }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F,N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] {
         let base_hint = self.base_fn.bisection_hint(cube);
         map2(base_hint, &self.shift, |h, s| h.map(|z| z + *s))
     }
-
 }
 
-impl<'a, T, F : Float, const N : usize> GlobalAnalysis<F, Bounds<F>> for Shift<T,F,N>
-where T : LocalAnalysis<F, Bounds<F>, N> {
+impl<'a, T, F: Float, const N: usize> GlobalAnalysis<F, Bounds<F>> for Shift<T, F, N>
+where
+    T: LocalAnalysis<F, Bounds<F>, N>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         self.base_fn.global_analysis()
     }
 }
 
-impl<'a, T, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Shift<T,F,N>
-where T : LocalAnalysis<F, Bounds<F>, N> {
+impl<'a, T, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Shift<T, F, N>
+where
+    T: LocalAnalysis<F, Bounds<F>, N>,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+    fn local_analysis(&self, cube: &Cube<F, N>) -> Bounds<F> {
         self.base_fn.local_analysis(&cube.shift(&(-self.shift)))
     }
 }
@@ -184,33 +144,36 @@
 
 impl_shift_norm!(L1 L2 Linfinity);
 
-impl<'a, T, F : Float, C, const N : usize> Support<F,N> for Weighted<T, C>
-where T : Support<F, N>,
-      C : Constant<Type=F> {
-
+impl<'a, T, F: Float, C, const N: usize> Support<F, N> for Weighted<T, C>
+where
+    T: Support<F, N>,
+    C: Constant<Type = F>,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F,N> {
+    fn support_hint(&self) -> Cube<F, N> {
         self.base_fn.support_hint()
     }
 
     #[inline]
-    fn in_support(&self, x : &Loc<F,N>) -> bool {
+    fn in_support(&self, x: &Loc<F, N>) -> bool {
         self.base_fn.in_support(x)
     }
-    
+
     // fn fully_in_support(&self, cube : &Cube<F,N>) -> bool {
     //     self.base_fn.fully_in_support(cube)
     // }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F,N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] {
         self.base_fn.bisection_hint(cube)
     }
 }
 
-impl<'a, T, F : Float, C> GlobalAnalysis<F, Bounds<F>> for Weighted<T, C>
-where T : GlobalAnalysis<F, Bounds<F>>,
-      C : Constant<Type=F> {
+impl<'a, T, F: Float, C> GlobalAnalysis<F, Bounds<F>> for Weighted<T, C>
+where
+    T: GlobalAnalysis<F, Bounds<F>>,
+    C: Constant<Type = F>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         let Bounds(lower, upper) = self.base_fn.global_analysis();
@@ -222,11 +185,13 @@
     }
 }
 
-impl<'a, T, F : Float, C, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Weighted<T, C>
-where T : LocalAnalysis<F, Bounds<F>, N>,
-      C : Constant<Type=F> {
+impl<'a, T, F: Float, C, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Weighted<T, C>
+where
+    T: LocalAnalysis<F, Bounds<F>, N>,
+    C: Constant<Type = F>,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+    fn local_analysis(&self, cube: &Cube<F, N>) -> Bounds<F> {
         let Bounds(lower, upper) = self.base_fn.local_analysis(cube);
         debug_assert!(lower <= upper);
         match self.weight.value() {
@@ -238,31 +203,36 @@
 
 macro_rules! make_weighted_scalarop_rhs {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
-        impl<F : Float, T> std::ops::$trait_assign<F> for Weighted<T, F> {
+        impl<F: Float, T> std::ops::$trait_assign<F> for Weighted<T, F> {
             #[inline]
-            fn $fn_assign(&mut self, t : F) {
+            fn $fn_assign(&mut self, t: F) {
                 self.weight.$fn_assign(t);
             }
         }
 
-        impl<'a, F : Float, T> std::ops::$trait<F> for Weighted<T, F> {
+        impl<'a, F: Float, T> std::ops::$trait<F> for Weighted<T, F> {
             type Output = Self;
             #[inline]
-            fn $fn(mut self, t : F) -> Self {
+            fn $fn(mut self, t: F) -> Self {
                 self.weight.$fn_assign(t);
                 self
             }
         }
 
-        impl<'a, F : Float, T> std::ops::$trait<F> for &'a Weighted<T, F>
-        where T : Clone {
+        impl<'a, F: Float, T> std::ops::$trait<F> for &'a Weighted<T, F>
+        where
+            T: Clone,
+        {
             type Output = Weighted<T, F>;
             #[inline]
-            fn $fn(self, t : F) -> Self::Output {
-                Weighted { weight : self.weight.$fn(t), base_fn : self.base_fn.clone() }
+            fn $fn(self, t: F) -> Self::Output {
+                Weighted {
+                    weight: self.weight.$fn(t),
+                    base_fn: self.base_fn.clone(),
+                }
             }
         }
-    }
+    };
 }
 
 make_weighted_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
@@ -282,52 +252,60 @@
 
 impl_weighted_norm!(L1 L2 Linfinity);
 
-
 /// Normalisation of [`Support`] and [`Mapping`] to L¹ norm 1.
 ///
 /// Currently only scalar-valued functions are supported.
 #[derive(Copy, Clone, Debug, Serialize, PartialEq)]
 pub struct Normalised<T>(
     /// The base [`Support`] or [`Mapping`].
-    pub T
+    pub T,
 );
 
-impl<'a, T, F : Float, const N : usize> Mapping<Loc<F, N>> for Normalised<T>
-where T : Norm<F, L1> + Mapping<Loc<F,N>, Codomain=F> {
+impl<'a, T, F: Float, const N: usize> Mapping<Loc<F, N>> for Normalised<T>
+where
+    T: Norm<F, L1> + Mapping<Loc<F, N>, Codomain = F>,
+{
     type Codomain = F;
 
     #[inline]
-    fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+    fn apply<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Codomain {
         let w = self.0.norm(L1);
-        if w == F::ZERO { F::ZERO } else { self.0.apply(x) / w }
+        if w == F::ZERO {
+            F::ZERO
+        } else {
+            self.0.apply(x) / w
+        }
     }
 }
 
-impl<'a, T, F : Float, const N : usize> Support<F,N> for Normalised<T>
-where T : Norm<F, L1> + Support<F, N> {
-
+impl<'a, T, F: Float, const N: usize> Support<F, N> for Normalised<T>
+where
+    T: Norm<F, L1> + Support<F, N>,
+{
     #[inline]
-    fn support_hint(&self) -> Cube<F,N> {
+    fn support_hint(&self) -> Cube<F, N> {
         self.0.support_hint()
     }
 
     #[inline]
-    fn in_support(&self, x : &Loc<F,N>) -> bool {
+    fn in_support(&self, x: &Loc<F, N>) -> bool {
         self.0.in_support(x)
     }
-    
+
     // fn fully_in_support(&self, cube : &Cube<F,N>) -> bool {
     //     self.0.fully_in_support(cube)
     // }
 
     #[inline]
-    fn bisection_hint(&self, cube : &Cube<F,N>) -> [Option<F>; N] {
+    fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] {
         self.0.bisection_hint(cube)
     }
 }
 
-impl<'a, T, F : Float> GlobalAnalysis<F, Bounds<F>> for Normalised<T>
-where T : Norm<F, L1> + GlobalAnalysis<F, Bounds<F>> {
+impl<'a, T, F: Float> GlobalAnalysis<F, Bounds<F>> for Normalised<T>
+where
+    T: Norm<F, L1> + GlobalAnalysis<F, Bounds<F>>,
+{
     #[inline]
     fn global_analysis(&self) -> Bounds<F> {
         let Bounds(lower, upper) = self.0.global_analysis();
@@ -338,10 +316,12 @@
     }
 }
 
-impl<'a, T, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N> for Normalised<T>
-where T : Norm<F, L1> + LocalAnalysis<F, Bounds<F>, N> {
+impl<'a, T, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Normalised<T>
+where
+    T: Norm<F, L1> + LocalAnalysis<F, Bounds<F>, N>,
+{
     #[inline]
-    fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+    fn local_analysis(&self, cube: &Cube<F, N>) -> Bounds<F> {
         let Bounds(lower, upper) = self.0.local_analysis(cube);
         debug_assert!(lower <= upper);
         let w = self.0.norm(L1);
@@ -350,12 +330,18 @@
     }
 }
 
-impl<'a, T, F : Float> Norm<F, L1> for Normalised<T>
-where T : Norm<F, L1> {
+impl<'a, T, F: Float> Norm<F, L1> for Normalised<T>
+where
+    T: Norm<F, L1>,
+{
     #[inline]
-    fn norm(&self, _ : L1) -> F {
+    fn norm(&self, _: L1) -> F {
         let w = self.0.norm(L1);
-        if w == F::ZERO { F::ZERO } else { F::ONE }
+        if w == F::ZERO {
+            F::ZERO
+        } else {
+            F::ONE
+        }
     }
 }
 
@@ -388,24 +374,26 @@
 
 /// Generator of [`Support`]-implementing component functions based on low storage requirement
 /// [ids][`Self::Id`].
-pub trait SupportGenerator<F : Float, const N : usize>
-: MulAssign<F> + DivAssign<F> + Neg<Output=Self> + Clone + Sync + Send + 'static {
+pub trait SupportGenerator<F: Float, const N: usize>:
+    MulAssign<F> + DivAssign<F> + Neg<Output = Self> + Clone + Sync + Send + 'static
+{
     /// The identification type
-    type Id : 'static + Copy;
+    type Id: 'static + Copy;
     /// The type of the [`Support`] (often also a [`Mapping`]).
-    type SupportType : 'static + Support<F, N>;
+    type SupportType: 'static + Support<F, N>;
     /// An iterator over all the [`Support`]s of the generator.
-    type AllDataIter<'a> : Iterator<Item=(Self::Id, Self::SupportType)> where Self : 'a;
+    type AllDataIter<'a>: Iterator<Item = (Self::Id, Self::SupportType)>
+    where
+        Self: 'a;
 
     /// Returns the component identified by `id`.
     ///
     /// Panics if `id` is an invalid identifier.
-    fn support_for(&self, id : Self::Id) -> Self::SupportType;
-    
+    fn support_for(&self, id: Self::Id) -> Self::SupportType;
+
     /// Returns the number of different components in this generator.
     fn support_count(&self) -> usize;
 
     /// Returns an iterator over all pairs of `(id, support)`.
     fn all_data(&self) -> Self::AllDataIter<'_>;
 }
-
diff -r 1f19c6bbf07b -r 2e4517b55442 src/bounds.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/bounds.rs	Sun Apr 27 20:41:36 2025 -0500
@@ -0,0 +1,246 @@
+/*!
+Bounded and minimizable/maximizable mappings.
+*/
+
+use crate::instance::Instance;
+use crate::loc::Loc;
+use crate::mapping::RealMapping;
+use crate::sets::{Cube, Set};
+use crate::types::{Float, Num};
+
+/// Trait for globally analysing a property `A` of a [`Mapping`].
+///
+/// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as
+/// [`Bounds`][super::aggregator::Bounds].
+pub trait GlobalAnalysis<F: Num, A> {
+    /// Perform global analysis of the property `A` of `Self`.
+    ///
+    /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds],
+    /// this function will return global upper and lower bounds for the mapping
+    /// represented by `self`.
+    fn global_analysis(&self) -> A;
+}
+
+// default impl<F, A, N, L> GlobalAnalysis<F, A, N> for L
+// where L : LocalAnalysis<F, A, N> {
+//     #[inline]
+//     fn global_analysis(&self) -> Bounds<F> {
+//         self.local_analysis(&self.support_hint())
+//     }
+// }
+
+/// Trait for locally analysing a property `A` of a [`Mapping`] (implementing [`Support`])
+/// within a [`Cube`].
+///
+/// Typically `A` is an [`Aggregator`][super::aggregator::Aggregator] such as
+/// [`Bounds`][super::aggregator::Bounds].
+pub trait LocalAnalysis<F: Num, A, const N: usize>: GlobalAnalysis<F, A> {
+    /// Perform local analysis of the property `A` of `Self`.
+    ///
+    /// As an example, in the case of `A` being [`Bounds`][super::aggregator::Bounds],
+    /// this function will return upper and lower bounds within `cube` for the mapping
+    /// represented by `self`.
+    fn local_analysis(&self, cube: &Cube<F, N>) -> A;
+}
+
+/// Trait for determining the upper and lower bounds of an float-valued [`Mapping`].
+///
+/// This is a blanket-implemented alias for [`GlobalAnalysis`]`<F, Bounds<F>>`
+/// [`Mapping`] is not a supertrait to allow flexibility in the implementation of either
+/// reference or non-reference arguments.
+pub trait Bounded<F: Float>: GlobalAnalysis<F, Bounds<F>> {
+    /// Return lower and upper bounds for the values of of `self`.
+    #[inline]
+    fn bounds(&self) -> Bounds<F> {
+        self.global_analysis()
+    }
+}
+
+impl<F: Float, T: GlobalAnalysis<F, Bounds<F>>> Bounded<F> for T {}
+
+/// A [`RealMapping`] that provides rough bounds as well as minimisation and maximisation.
+pub trait MinMaxMapping<F: Float, const N: usize>: RealMapping<F, N> + Bounded<F> {
+    /// Maximise the mapping within stated value `tolerance`.
+    ///
+    /// At most `max_steps` refinement steps are taken.
+    /// Returns the approximate maximiser and the corresponding function value.
+    fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F);
+
+    /// Maximise the mapping within stated value `tolerance` subject to a lower bound.
+    ///
+    /// At most `max_steps` refinement steps are taken.
+    /// Returns the approximate maximiser and the corresponding function value when one is found
+    /// above the `bound` threshold, otherwise `None`.
+    fn maximise_above(
+        &mut self,
+        bound: F,
+        tolerance: F,
+        max_steps: usize,
+    ) -> Option<(Loc<F, N>, F)>;
+
+    /// Minimise the mapping within stated value `tolerance`.
+    ///
+    /// At most `max_steps` refinement steps are taken.
+    /// Returns the approximate minimiser and the corresponding function value.
+    fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F);
+
+    /// Minimise the mapping within stated value `tolerance` subject to a lower bound.
+    ///
+    /// At most `max_steps` refinement steps are taken.
+    /// Returns the approximate minimiser and the corresponding function value when one is found
+    /// above the `bound` threshold, otherwise `None`.
+    fn minimise_below(
+        &mut self,
+        bound: F,
+        tolerance: F,
+        max_steps: usize,
+    ) -> Option<(Loc<F, N>, F)>;
+
+    /// Verify that the mapping has a given upper `bound` within indicated `tolerance`.
+    ///
+    /// At most `max_steps` refinement steps are taken.
+    fn has_upper_bound(&mut self, bound: F, tolerance: F, max_steps: usize) -> bool;
+
+    /// Verify that the mapping has a given lower `bound` within indicated `tolerance`.
+    ///
+    /// At most `max_steps` refinement steps are taken.
+    fn has_lower_bound(&mut self, bound: F, tolerance: F, max_steps: usize) -> bool;
+}
+
+/// Upper and lower bounds on an `F`-valued function.
+#[derive(Copy, Clone, Debug)]
+pub struct Bounds<F>(
+    /// Lower bound
+    pub F,
+    /// Upper bound
+    pub F,
+);
+
+impl<F: Copy> Bounds<F> {
+    /// Returns the lower bound
+    #[inline]
+    pub fn lower(&self) -> F {
+        self.0
+    }
+
+    /// Returns the upper bound
+    #[inline]
+    pub fn upper(&self) -> F {
+        self.1
+    }
+}
+
+impl<F: Float> Bounds<F> {
+    /// Returns a uniform bound.
+    ///
+    /// This is maximum over the absolute values of the upper and lower bound.
+    #[inline]
+    pub fn uniform(&self) -> F {
+        let &Bounds(lower, upper) = self;
+        lower.abs().max(upper.abs())
+    }
+
+    /// Construct a bounds, making sure `lower` bound is less than `upper`
+    #[inline]
+    pub fn corrected(lower: F, upper: F) -> Self {
+        if lower <= upper {
+            Bounds(lower, upper)
+        } else {
+            Bounds(upper, lower)
+        }
+    }
+
+    /// Refine the lower bound
+    #[inline]
+    pub fn refine_lower(&self, lower: F) -> Self {
+        let &Bounds(l, u) = self;
+        debug_assert!(l <= u);
+        Bounds(l.max(lower), u.max(lower))
+    }
+
+    /// Refine the lower bound
+    #[inline]
+    pub fn refine_upper(&self, upper: F) -> Self {
+        let &Bounds(l, u) = self;
+        debug_assert!(l <= u);
+        Bounds(l.min(upper), u.min(upper))
+    }
+}
+
+impl<'a, F: Float> std::ops::Add<Self> for Bounds<F> {
+    type Output = Self;
+    #[inline]
+    fn add(self, Bounds(l2, u2): Self) -> Self::Output {
+        let Bounds(l1, u1) = self;
+        debug_assert!(l1 <= u1 && l2 <= u2);
+        Bounds(l1 + l2, u1 + u2)
+    }
+}
+
+impl<'a, F: Float> std::ops::Mul<Self> for Bounds<F> {
+    type Output = Self;
+    #[inline]
+    fn mul(self, Bounds(l2, u2): Self) -> Self::Output {
+        let Bounds(l1, u1) = self;
+        debug_assert!(l1 <= u1 && l2 <= u2);
+        let a = l1 * l2;
+        let b = u1 * u2;
+        // The order may flip when negative numbers are involved, so need min/max
+        Bounds(a.min(b), a.max(b))
+    }
+}
+
+impl<F: Float> std::iter::Product for Bounds<F> {
+    #[inline]
+    fn product<I>(mut iter: I) -> Self
+    where
+        I: Iterator<Item = Self>,
+    {
+        match iter.next() {
+            None => Bounds(F::ZERO, F::ZERO),
+            Some(init) => iter.fold(init, |a, b| a * b),
+        }
+    }
+}
+
+impl<F: Float> Set<F> for Bounds<F> {
+    fn contains<I: Instance<F>>(&self, item: I) -> bool {
+        let v = item.own();
+        let &Bounds(l, u) = self;
+        debug_assert!(l <= u);
+        l <= v && v <= u
+    }
+}
+
+impl<F: Float> Bounds<F> {
+    /// Calculate a common bound (glb, lub) for two bounds.
+    #[inline]
+    pub fn common(&self, &Bounds(l2, u2): &Self) -> Self {
+        let &Bounds(l1, u1) = self;
+        debug_assert!(l1 <= u1 && l2 <= u2);
+        Bounds(l1.min(l2), u1.max(u2))
+    }
+
+    /// Indicates whether `Self` is a superset of the argument bound.
+    #[inline]
+    pub fn superset(&self, &Bounds(l2, u2): &Self) -> bool {
+        let &Bounds(l1, u1) = self;
+        debug_assert!(l1 <= u1 && l2 <= u2);
+        l1 <= l2 && u2 <= u1
+    }
+
+    /// Returns the greatest bound contained by both argument bounds, if one exists.
+    #[inline]
+    pub fn glb(&self, &Bounds(l2, u2): &Self) -> Option<Self> {
+        let &Bounds(l1, u1) = self;
+        debug_assert!(l1 <= u1 && l2 <= u2);
+        let l = l1.max(l2);
+        let u = u1.min(u2);
+        debug_assert!(l <= u);
+        if l < u {
+            Some(Bounds(l, u))
+        } else {
+            None
+        }
+    }
+}
diff -r 1f19c6bbf07b -r 2e4517b55442 src/lib.rs
--- a/src/lib.rs	Sun Apr 27 20:29:43 2025 -0500
+++ b/src/lib.rs	Sun Apr 27 20:41:36 2025 -0500
@@ -34,6 +34,7 @@
 #[macro_use]
 pub mod loc;
 pub mod bisection_tree;
+pub mod bounds;
 pub mod coefficients;
 pub mod convex;
 pub mod direct_product;