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