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