src/bisection_tree/btfn.rs

branch
dev
changeset 96
962c8e346ab9
parent 63
f7b87d84864d
child 97
4e80fb049dca
--- a/src/bisection_tree/btfn.rs	Sun Apr 27 14:48:13 2025 -0500
+++ b/src/bisection_tree/btfn.rs	Sun Apr 27 15:45:40 2025 -0500
@@ -1,26 +1,24 @@
-
+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::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 +28,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 +58,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 +77,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 +91,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 +108,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 +125,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 +152,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 +188,9 @@
         }
 
         BTFN {
-            bt : bt,
-            generator : Arc::new(both),
-            _phantoms : std::marker::PhantomData,
+            bt: bt,
+            generator: Arc::new(both),
+            _phantoms: std::marker::PhantomData,
         }
     }
 }
@@ -231,7 +231,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 +274,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 +368,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 +396,7 @@
                 self.new_generator(std::ops::$trait::$fn(&self.generator))
             }
         }*/
-    }
+    };
 }
 
 make_btfn_unaryop!(Neg, neg);
@@ -405,39 +405,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 +446,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 +504,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 +535,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 +590,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 +638,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 +653,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 +700,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 +721,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 +782,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 +817,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,20 +837,28 @@
 // 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> {
-
+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()
+    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()
     }
 
     /// Maximise the `BTFN` within stated value `tolerance` subject to a lower bound.
@@ -849,19 +866,38 @@
     /// 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.")
+    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.")
     }
 
     /// 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()
+    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()
     }
 
     /// Minimise the `BTFN` within stated value `tolerance` subject to a lower bound.
@@ -869,25 +905,50 @@
     /// 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.")
+    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.")
     }
 
     /// 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.")
+    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.")
     }
 
     /// 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.")
+    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.")
     }
 }

mercurial