src/bisection_tree/btfn.rs


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::{Apply, Mapping, Differentiable};
//use crate::linops::{Apply, Linear};
use crate::sets::Set;
use crate::sets::Cube;
use crate::loc::Loc;
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::*;

/// 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,
/// and `N` the dimension.
///
/// The `generator` lists the component functions that have to implement [`Support`].
/// 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>,
}

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> {

    /// 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`.
    /// Use [`Self::construct`] if no preinitialised tree is available. Use [`Self::new_refresh`]
    /// 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 {
        Self::new_arc(bt, Arc::new(generator))
    }

    fn new_arc(bt : BT, generator : Arc<G>) -> Self {
        BTFN {
            bt : bt,
            generator : generator,
            _phantoms : std::marker::PhantomData,
        }
    }

    /// Create a new BTFN support generator and a pre-initialised bisection tree,
    /// cloning the tree and refreshing aggregators.
    ///
    /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`, but
    /// 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 {
        // 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.
        let mut btnew = bt.clone();
        btnew.refresh_aggregator(&generator);
        BTFN::new(btnew, generator)
    }

    /// Create a new BTFN from a support generator, domain, and depth for a new [`BT`].
    ///
    /// 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 {
        Self::construct_arc(domain, depth, Arc::new(generator))
    }

    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);
        }
        Self::new_arc(bt, generator)
    }

    /// Convert the aggregator of the [`BTFN`] to a different one.
    ///
    /// 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> {
        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 {
        BTFN::new_refresh(&self.bt, generator)
    }

    /// Refresh aggregator after updates to generator
    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> {
    /// 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>  {
        BTFN::construct_arc(domain, depth, self.generator)
    }
}

/// A BTFN with no bisection tree.
///
/// 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>;

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 {
        BTFN {
            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> {

    /// 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> {

        let mut bt = self.bt.clone();
        let both = BothGenerators(Arc::clone(&self.generator), g2);

        for (d, support) in both.all_right_data() {
            bt.insert(d, &support);
        }

        BTFN {
            bt : bt,
            generator : Arc::new(both),
            _phantoms : std::marker::PhantomData,
        }
    }
}

macro_rules! make_btfn_add {
    ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => {
        impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize>
        std::ops::Add<BTFN<F, G2, BT2, N>> for
        $lhs
        where BT1 : BTImpl<F, N, Data=usize>,
              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
              G2 : SupportGenerator<F, N, Id=usize>,
              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
              G2::SupportType : LocalAnalysis<F, BT1::Agg, N> {
            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
            #[inline]
            fn add(self, other : BTFN<F, G2, BT2, N>) -> Self::Output {
                $preprocess(self).add_another(other.generator)
            }
        }

        impl<'a, 'b, F : Float, G1, G2,  BT1, BT2, const N : usize>
        std::ops::Add<&'b BTFN<F, G2, BT2, N>> for
        $lhs
        where BT1 : BTImpl<F, N, Data=usize>,
              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
              G2 : SupportGenerator<F, N, Id=usize>,
              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
              G2::SupportType : LocalAnalysis<F, BT1::Agg, N> {

            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
            #[inline]
            fn add(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output {
                $preprocess(self).add_another(other.generator.clone())
            }
        }
    }
}

make_btfn_add!(BTFN<F, G1, BT1, N>, std::convert::identity, );
make_btfn_add!(&'a BTFN<F, G1, BT1, N>, Clone::clone, );

macro_rules! make_btfn_sub {
    ($lhs:ty, $preprocess:path, $($extra_trait:ident)?) => {
        impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize>
        std::ops::Sub<BTFN<F, G2, BT2, N>> for
        $lhs
        where BT1 : BTImpl<F, N, Data=usize>,
              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
              G2 : SupportGenerator<F, N, Id=usize>,
              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
              G2::SupportType : LocalAnalysis<F, BT1::Agg, N> {
            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
            #[inline]
            fn sub(self, other : BTFN<F, G2, BT2, N>) -> Self::Output {
                $preprocess(self).add_another(Arc::new(
                    Arc::try_unwrap(other.generator)
                        .unwrap_or_else(|arc| (*arc).clone())
                        .neg()
                ))
            }
        }

        impl<'a, 'b, F : Float, G1, G2,  BT1, BT2, const N : usize>
        std::ops::Sub<&'b BTFN<F, G2, BT2, N>> for
        $lhs
        where BT1 : BTImpl<F, N, Data=usize>,
              G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?,
              G2 : SupportGenerator<F, N, Id=usize> + Clone,
              G1::SupportType : LocalAnalysis<F, BT1::Agg, N>,
              G2::SupportType : LocalAnalysis<F, BT1::Agg, N>,
              &'b G2 : std::ops::Neg<Output=G2> {

            type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>;
            #[inline]
            fn sub(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output {
                $preprocess(self).add_another(Arc::new((*other.generator).clone().neg()))
            }
        }
    }
}

make_btfn_sub!(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> {
            #[inline]
            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> {
            type Output = Self;
            #[inline]
            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> {
            type Output = BTFN<F, G, BT, N>;
            #[inline]
            fn $fn(self, t : F) -> Self::Output {
                self.new_generator(self.generator.$fn(t))
            }
        }
    }
}

make_btfn_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
make_btfn_scalarop_rhs!(Div, div, DivAssign, div_assign);

macro_rules! make_btfn_scalarop_lhs {
    ($trait:ident, $fn:ident, $fn_assign:ident, $($f:ident)+) => { $(
        impl<G, BT, const N : usize>
        std::ops::$trait<BTFN<$f, G, BT, N>>
        for $f
        where BT : BTImpl<$f, N>,
              G : SupportGenerator<$f, N, Id=BT::Data>,
              G::SupportType : LocalAnalysis<$f, BT::Agg, N> {
            type Output = BTFN<$f, G, BT, N>;
            #[inline]
            fn $fn(self, mut a : BTFN<$f, G, BT, N>) -> Self::Output {
                Arc::make_mut(&mut a.generator).$fn_assign(self);
                a.refresh_aggregator();
                a
            }
        }

        impl<'a, G, BT, const N : usize>
        std::ops::$trait<&'a BTFN<$f, G, BT, N>>
        for $f
        where BT : BTImpl<$f, N>,
              G : SupportGenerator<$f, N, Id=BT::Data> + Clone,
              G::SupportType : LocalAnalysis<$f, BT::Agg, N>,
              // FIXME: This causes compiler overflow
              /*&'a G : std::ops::$trait<$f,Output=G>*/ {
            type Output = BTFN<$f, G, BT, N>;
            #[inline]
            fn $fn(self, a : &'a BTFN<$f, G, BT, N>) -> Self::Output {
                let mut tmp = (*a.generator).clone();
                tmp.$fn_assign(self);
                a.new_generator(tmp)
                // FIXME: Prevented by the compiler overflow above.
                //a.new_generator(a.generator.$fn(a))
            }
        }
    )+ }
}

make_btfn_scalarop_lhs!(Mul, mul, mul_assign, f32 f64);
make_btfn_scalarop_lhs!(Div, div, div_assign, f32 f64);

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> {
            type Output = Self;
            #[inline]
            fn $fn(mut self) -> Self::Output {
                self.generator = Arc::new(Arc::unwrap_or_clone(self.generator).$fn());
                self.refresh_aggregator();
                self
            }
        }

        /*impl<'a, F : Float, G, BT, const N : usize>
        std::ops::$trait
        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<Output=G> {
            type Output = BTFN<F, G, BT, N>;
            #[inline]
            fn $fn(self) -> Self::Output {
                self.new_generator(std::ops::$trait::$fn(&self.generator))
            }
        }*/
    }
}

make_btfn_unaryop!(Neg, neg);

//
// Apply, Mapping, Differentiate
//

impl<'a, F : Float, G, BT, V, const N : usize> Apply<&'a 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> + Apply<&'a Loc<F, N>, Output = V>,
      V : Sum {

    type Output = V;

    fn apply(&self, x : &'a Loc<F, N>) -> Self::Output {
        self.bt.iter_at(x)
            .map(|&d| self.generator.support_for(d).apply(x)).sum()
    }
}

impl<F : Float, G, BT, V, const N : usize> Apply<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> + Apply<Loc<F, N>, Output = V>,
      V : Sum {

    type Output = V;

    fn apply(&self, x : Loc<F, N>) -> Self::Output {
        self.bt.iter_at(&x)
            .map(|&d| self.generator.support_for(d).apply(x)).sum()
    }
}

impl<'a, F : Float, G, BT, V, W, const N : usize> Differentiable<&'a 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>
                       + Differentiable<&'a Loc<F, N>, Output = V>
                       + Apply<&'a Loc<F, N>, Output = W>,
      V : Sum,
      W : Sum {

    type Output = V;

    fn differential(&self, x : &'a Loc<F, N>) -> V {
        self.bt.iter_at(x)
            .map(|&d| self.generator.support_for(d).differential(x))
            .sum()
    }

    fn linearisation_error_gen(&self, x : &'a Loc<F, N>, y : &'a Loc<F, N>, z : &'a Loc<F, N>) -> W {
        self.bt.iter_at(x)
            .map(|&d| self.generator.support_for(d).linearisation_error_gen(x, y, z)).sum()
    }

    fn linearisation_error(&self, x : &'a Loc<F, N>, y : &'a Loc<F, N>) -> W {
        self.bt.iter_at(x)
            .map(|&d| self.generator.support_for(d).linearisation_error(x, y)).sum()
    }
}

impl<F : Float, G, BT, V, W, const N : usize> Differentiable<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>
                       + Differentiable<Loc<F, N>, Output = V>
                       + Apply<Loc<F, N>, Output = W>,
      V : Sum,
      W : Sum {

    type Output = V;

    fn differential(&self, x : Loc<F, N>) -> V {
        self.bt.iter_at(&x)
            .map(|&d| self.generator.support_for(d).differential(x))
            .sum()
    }

    fn linearisation_error_gen(&self, x : Loc<F, N>, y : Loc<F, N>, z : Loc<F, N>) -> W {
        self.bt.iter_at(&x)
            .map(|&d| self.generator.support_for(d).linearisation_error_gen(x, y, z)).sum()
    }

    fn linearisation_error(&self, x : Loc<F, N>, y : Loc<F, N>) -> W {
        self.bt.iter_at(&x)
            .map(|&d| self.generator.support_for(d).linearisation_error(x, y)).sum()
    }
}

//
// 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> {

    #[inline]
    fn global_analysis(&self) -> BT::Agg {
        self.bt.global_analysis()
    }
}

//
// Blanket implementation of BTFN as a linear functional over objects
// that are linear functionals over BTFN.
//

/*
impl<'b, X, F : Float, G, BT, const N : usize> Apply<&'b X, 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>,
      X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> {

    #[inline]
    fn apply(&self, x : &'b X) -> F {
        x.apply(&self)
    }
}

impl<X, F : Float, G, BT, const N : usize> Apply<X, 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>,
      X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> {

    #[inline]
    fn apply(&self, x : X) -> F {
        x.apply(&self)
    }
}

impl<X, F : Float, G, BT, const N : usize> Linear<X>
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>,
      X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> {
    type Codomain = F;
}
*/

/// Helper trait for performing approximate minimisation using P2 elements.
///
/// `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, 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);

}

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) {
        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.
            let [a, b, c, d] = self.corners();
            let [va, vb, vc, vd] = [g(&a), g(&b), g(&c), g(&d)];

            let ab = midpoint(&a, &b);
            let bc = midpoint(&b, &c);
            let ca = midpoint(&c, &a);
            let cd = midpoint(&c, &d);
            let da = midpoint(&d, &a);
            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 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 r2@(_, v2) = m2.minimise(&s2);

            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' => {
                    // Pseudo-randomise edge midpoints
                    let Loc([x, y]) = a;
                    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)],
                    }
                 },
                 _ => [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],
            );

            m1.minimise(self)
        }
    }
}

/// Helper type to use [`P2Refiner`] for maximisation.
struct RefineMax;

/// Helper type to use [`P2Refiner`] for minimisation.
struct RefineMin;

/// 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> {
    /// The maximum / minimum should be above / below this threshold.
    /// If the threshold cannot be satisfied, the refiner will return `None`.
    bound : Option<F>,
    /// Tolerance for function value estimation.
    tolerance : F,
    /// Maximum number of steps to execute the refiner for
    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,
}

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
    ) -> RefinerResult<Bounds<F>, Self::Result> {

        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)
        }

        // g gives the negative of the value of the function presented by `data` and `generator`.
        let g = move |x : &Loc<F, N>| {
            let f = move |&d| generator.support_for(d).apply(x);
            -data.iter().map(f).sum::<F>()
        };
        // … so the negative of the minimum is the maximm we want.
        let (x, _neg_v) = cube.p2_minimise(g);
        //let v = -neg_v;
        let v = -g(&x);

        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
        } else {
            // The data is refined enough, so return new hopefully better bounds
            // and the maximiser.
            let res = (x, v);
            let bounds = Bounds(v, v);
            RefinerResult::Uncertain(bounds, Some(res))
        }
    }

    fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) {
        match (*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> {
    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
    ) -> RefinerResult<Bounds<F>, Self::Result> {
    
        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)
        }

        // g gives the value of the function presented by `data` and `generator`.
        let g = move |x : &Loc<F, N>| {
            let f = move |&d| generator.support_for(d).apply(x);
            data.iter().map(f).sum::<F>()
        };
        // Minimise it.
        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*/) {
            // The function isn't refined enough in `cube`, so return None
            // to indicate that further subdivision is required.
            RefinerResult::NeedRefinement
        } else {
            // The data is refined enough, so return new hopefully better bounds
            // and the minimiser.
            let res = (x, v);
            let l = aggregator.lower();
            let bounds = if l > v {
                eprintln!("imprecision!");
                Bounds(l, l)
            } else {
                Bounds(v, v)
            };
            RefinerResult::Uncertain(bounds, Some(res))
        }
    }

    fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) {
        match (*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> {
    /// The upper/lower bound to check for
    bound : F,
    /// Tolerance for function value estimation.
    tolerance : F,
    /// Maximum number of steps to execute the refiner for
    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,
}

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
    ) -> RefinerResult<Bounds<F>, Self::Result> {
        if aggregator.upper() <= self.bound + self.tolerance {
            // Below upper bound within tolerances. Indicate uncertain success.
            RefinerResult::Uncertain(*aggregator, true)
        } else if aggregator.lower() >= self.bound - self.tolerance {
            // Above upper bound within tolerances. Indicate certain failure.
            RefinerResult::Certain(false)
        } else if step < self.max_steps {
            // No decision possible, but within step bounds - further subdivision is required.
            RefinerResult::NeedRefinement
        } else {
            // No decision possible, but past step bounds
            RefinerResult::Uncertain(*aggregator, false)
        }
    }

    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> {
    type Result = bool;
    type Sorting = UpperBoundSorting<F>;

    fn refine(
        &self,
        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.
            RefinerResult::Uncertain(*aggregator, true)
        } else if aggregator.upper() <= self.bound + self.tolerance {
            // Below lower bound within tolerances. Indicate certain failure.
            RefinerResult::Certain(false)
        } else if step < self.max_steps {
            // No decision possible, but within step bounds - further subdivision is required.
            RefinerResult::NeedRefinement
        } else {
            // No decision possible, but past step bounds
            RefinerResult::Uncertain(*aggregator, false)
        }
    }

    fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result) {
        *r1 = *r1 && r2;
    }
}

// FIXME: The most likely reason for the “Refiner failure” expectation in the methods below
// is numerical inaccuracy: the `glb` maintained in `HeapContainer` (`refine.rs`) becomes bigger
// than the *upper bound* of nodes attempted to be inserted into the `heap` in the container.
// But the `glb` is there exactly to prevent that. Due to numerical inaccuracy, however, a
// newly subdivided node may have lower upper bound than the original lower bound that should
// have been above the `glb` since the node was picked from the queue. Due to the subdivision
// process, if a node whose lower bound is at the `glb` is picked, all of its refined subnodes
// should have lower bound at least the old `glb`, so in a single-threaded situation there should
// always be nodes above the `glb` in the queue. In a multi-threaded situation a node below the
// `glb` may be picked by some thread. When that happens, that thread inserts no new nodes into
// the queue. If the queue empties as a result of that, the thread goes to wait for other threads
// to produce results. Since some node had a node whose lower bound was above the `glb`, eventually
// 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()
    }

    /// 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.")
    }

    /// 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()
    }

    /// 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.")
    }

    /// 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.")
    }

    /// 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.")
    }
}