changeset
Improve documentation
Sat, 22 Oct 2022 14:40:28 +0300
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Sat, 22 Oct 2022 14:40:28 +0300
- changeset 5
- 59dc4c5883f4
- parent 4
- 61b068c50e25
- child 6
- d80b87b8acd0
Improve documentation
--- a/README.md Mon Oct 24 10:52:19 2022 +0300 +++ b/README.md Sat Oct 22 14:40:28 2022 +0300 @@ -3,22 +3,17 @@ Author: Tuomo Valkonen <tuomov@iki.fi> -**This crate is very much incomplete.** - This crate/repository/package contains some general utility routines and tools for implementing iterative algorithms and numerical computing in Rust. Former versions of the package were for Julia. They are no longer mintained. Included are: - * [Linear operator](crate::linops) abstraction, with support for matrices via `nalgebra` - or `ndarray` if this crate is compiled with the relevant features. - * [Iteration mechanism](crate::iterate) for generically implementing intermittent verbosity and stopping rules in algorithms. + * [Linear operator](crate::linops), [Euclidean space][crate::euclidean], and [norm][crate::norms] abstractions with support for matrices and vectors via `nalgebra`. + * [Iteration abstraction](crate::iterate) for generically implementing intermittent verbosity and stopping rules in iterative algorithms. * The iterators can employ a [logger](crate::logger) to store the intermittent results - (e.g. function values) in memory to later write into a `csv`-file. - * [Bisection trees](crate::bisection_tree) for branch-and-bound optimisation. + (e.g. function values) in memory to later [write][crate::tabledump] into a `csv`-file. + * Geometrical [bisection trees](crate::bisection_tree) for efficient representations of sums of functions and branch-and-bound optimisation. -<!--The code is used by [ImageTools][] and both, for example, by : ---> - -## Installation +<!--## Installation ??? +-->
--- a/src/bisection_tree.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,6 +1,52 @@ -/// -/// Geometrical bisection trees -/// +/*! +Geometrical bisection trees and efficient representations of sums of functions. + +[Bisection trees][BTImpl] subdivide at each branching node a [cubical][crate::sets::Cube] domain of +dimension $N$ into $P=2^N$ subcubes at a branching point. Objects implementing [`Support`] can be +[inserted][BTImpl::insert] into the leaves of the tree via a low storage requirement identifier. +The identifier is inserted into *every* leaf node of the tree intersected by the cube provided +by [`Support::support_hint`]. + +Bisection tree functions or [`BTFN`]s use bisection trees for efficient presentations of sums +$\sum_{k=1}^m f_k$ of mathematical functions $f_k$, where each $f_k$ typically has a small +support, i.e., set $\\{x \mid f(x) ≠ 0 \\}$, compared to the support of the entire sum. +To do so, the data type representing $f_k$ needs to implement [`Support`]. To evaluate the +value of the sum, each $f_k$ also needs to implement [`Mapping`][crate::mapping::Mapping]. +Moreover, the sum needs to be represented by a [`SupportGenerator`] that associates to a +low-storage-requirement identifier (typically `usize`) an object of the type that represents +$f_k$. [`BTFN`]s support basic vector space operations, and [minimisation][BTFN::minimise] and +[maximisation][BTFN::maximise] via a [branch-and-bound strategy][BTSearch::search_and_refine]. + +The nodes of a bisection tree also store aggregate information about the objects stored in the +tree via an [`Aggregator`]. This way, rough upper and lower [bound][Bounds] estimates on +$\sum_{k=1}^m f_k$ can be efficiently maintained in the tree, based on corresponding estimates +on each $f_k$. This is done [locally][LocalAnalysis] for every node and corresponding subcube +of the domain of the tree. + +## Type parameters + +The various types and traits herein depend on several type parameters: + +- `N` is the geometrical dimension and `F` the numerical type for coordinates. +- `P` is the branching factor and would generally equal $2^N$, but has to be explicitly specified + due to current restrictions in Rust's const generics. +- `D` is the data/identifier stored in the bisection tree, and `A` the [aggregate][Aggregator] + information of all the data in a branch. + +## Starting points + +[`BTImpl::new`] creates a new bisection tree. [`BTImpl::insert`] can be used to insert new data +into it. [`BTImpl::iter_at`] iterates over the data at a point `x` of the domain of the tree. + +A new [`BTFN`] can be constructed with [`BTFN::construct`] given an implementation of +a [`SupportGenerator`]. They can be summed and multipliced by a schalar using standard arithmetic +operations. The types of the objects in two summed `BTFN`s do not need to be the same. +To find an approximate minimum of a `BTFN` using a branch-and-bound strategy, +use [`BTFN::minimise`]. [`Bounded::bounds`] provides a shortcut to [`GlobalAnalysis`] with the +[`Bounds`] aggregator. If the rough bounds so obtained do not indicate that the `BTFN` is in some +given bounds, instead of doing a full minimisation and maximisation for higher quality bounds, +it is more efficient to use [`BTFN::has_upper_bound`] and [`BTFN::has_lower_bound`]. +*/ mod supportid; pub use supportid::*;
--- a/src/bisection_tree/aggregator.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree/aggregator.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,17 +1,33 @@ +/*! +Aggregation / summarisation of information in branches of bisection trees. +*/ + use crate::types::*; use crate::sets::Set; -/// Trait for aggregating information about a branch of a [`super::BT`] bisection tree. +/// Trait for aggregating information about a branch of a [bisection tree][super::BT]. +/// +/// Currently [`Bounds`] is the only provided aggregator. +/// It keeps track of upper and lower bounds of a function representeed by the `BT` by +/// summing [`Bounds`] produced by [`LocalAnalysis`][super::support::LocalAnalysis] of the +/// [`Support`][super::support::Support]s of the data stored in the tree. +/// For the `Bounds` aggregator: +/// * [`Self::aggregate`] sums input bounds to the current bound. This provides a conservative +/// estimate of the upper and lower bounds of a sum of functions. +/// * [`Self::summarise`] takes the maximum of the input bounds. This calculates the bounds +/// of a function on a greater domain from bounds on subdomains +/// (in practise [`Cube`][crate::sets::Cube]s). +/// pub trait Aggregator : Clone + std::fmt::Debug { - // Aggregate a new leaf data point to current state. + /// Aggregate a new data to current state. fn aggregate<I>(&mut self, aggregates : I) where I : Iterator<Item=Self>; - // Summarise several other aggregators, resetting current state. + /// Summarise several other aggregators, resetting current state. fn summarise<'a, I>(&'a mut self, aggregates : I) where I : Iterator<Item=&'a Self>; - /// Initialisation of aggregate data for an empty lower node. + /// Create a new “empty” aggregate data. fn new() -> Self; } @@ -32,11 +48,15 @@ } /// Upper and lower bounds on an `F`-valued function. -/// Returned by [`super::support::Bounded::bounds`]. #[derive(Copy,Clone,Debug)] -pub struct Bounds<F>(pub F, pub F); +pub struct Bounds<F>( + /// Lower bound + pub F, + /// Upper bound + pub F +); -impl<F : Num> Bounds<F> { +impl<F : Copy> Bounds<F> { /// Returns the lower bound #[inline] pub fn lower(&self) -> F { self.0 } @@ -47,8 +67,9 @@ } impl<F : Float> Bounds<F> { - /// Returns a uniform bound on the function (maximum over the absolute values of the - /// upper and lower bound). + /// 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; @@ -107,6 +128,7 @@ 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; @@ -114,7 +136,7 @@ l1 <= l2 && u2 <= u1 } - // Return the greatest bound contained by both argument bounds, if one exists. + /// 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;
--- a/src/bisection_tree/bt.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree/bt.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,9 +1,13 @@ + +/*! +Bisection tree basics, [`BT`] type and the [`BTImpl`] trait. +*/ use std::slice::{Iter,IterMut}; use std::iter::once; use std::rc::Rc; use serde::{Serialize, Deserialize}; -pub use nalgebra::Const; +pub(super) use nalgebra::Const; use itertools::izip; use crate::iter::{MapF,Mappable}; @@ -15,32 +19,43 @@ map2_indexed, collect_into_array_unchecked }; -pub use crate::sets::Cube; -pub use crate::loc::Loc; +use crate::sets::Cube; +use crate::loc::Loc; use super::support::*; use super::aggregator::*; +/// 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 enum NodeOption<F : Num, D, A : Aggregator, const N : usize, const P : usize> { +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. Uninitialised, - // TODO: replace with QuickVec fast and w/o allocs on single elements. + /// Indicates a leaf node containing a copy-on-write reference-counted vector + /// of data of type `D`. Leaf(Rc<Vec<D>>), + /// Indicates a branch node, cotaning a copy-on-write reference to the [`Branches`]. Branches(Rc<Branches<F, D, A, N, P>>), } /// Node of a [`BT`] bisection tree. +/// +/// 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> { + /// The data or branches under the node. pub(super) data : NodeOption<F, D, A, N, P>, /// Aggregator for `data`. pub(super) aggregator : A, } -/// Branch information of a [`Node`] of a [`BT`] bisection tree. +/// 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 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>, /// Subnodes @@ -79,14 +94,25 @@ } } +/// 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 + std::fmt::Debug { + /// Lower depth type. type Lower : Depth; + /// Returns a lower depth, if there still is one. fn lower(&self) -> Option<Self::Lower>; + /// Returns a lower depth or self if this is the lowest depth. fn lower_or(&self) -> Self::Lower; } +/// Dynamic (runtime) [`Depth`] for a [`BT`]. #[derive(Copy,Clone,Debug,Serialize,Deserialize)] -pub struct DynamicDepth(pub u8); +pub struct DynamicDepth( + /// The depth + pub u8 +); + impl Depth for DynamicDepth { type Lower = Self; #[inline] @@ -120,7 +146,11 @@ } impl_constdepth!(1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32); - +/// Trait for counting the branching factor of a [`BT`] of dimension `N`. +/// +/// 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> {} macro_rules! impl_branchcount { ($($n:literal)*) => { $( @@ -133,31 +163,53 @@ 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() } + /// Returns the node within `Self` containing the point `x`. + /// + /// 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] - pub fn get_node(&self, x : &Loc<F,N>) -> &Node<F,D,A,N,P> { + 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 data `D` in a [`BT`]. pub struct BTIter<'a, D> { iter : Iter<'a, D>, } -pub struct SubcubeIter<'b, F : Float, const N : usize, const P : usize> { +impl<'a, D> Iterator for BTIter<'a,D> { + type Item = &'a D; + #[inline] + fn next(&mut self) -> Option<&'a D> { + self.iter.next() + } +} + +/// 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, } +/// 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) -> Cube<F, N> { +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 { [branch, end] @@ -187,7 +239,9 @@ where Const<P> : BranchCount<N>, A : Aggregator { - pub fn new_with<S : LocalAnalysis <F, A, N>>( + /// 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 ) -> Self { @@ -201,14 +255,16 @@ } } - /// Get an iterator over the aggregators of the nodes of this branch head. + /// Returns an iterator over the aggregators of the nodes directly under this branch head. #[inline] - pub fn aggregators(&self) -> MapF<Iter<'_, Node<F,D,A,N,P>>, &'_ A> { + pub(super) fn aggregators(&self) -> MapF<Iter<'_, Node<F,D,A,N,P>>, &'_ A> { self.nodes.iter().mapF(Node::get_aggregator) } + /// Returns an iterator over the subcubes of `domain` subdivided at the branching point + /// of `self`. #[inline] - pub fn iter_subcubes<'b>(&self, domain : &'b Cube<F, N>) + pub(super) fn iter_subcubes<'b>(&self, domain : &'b Cube<F, N>) -> SubcubeIter<'b, F, N, P> { SubcubeIter { domain : domain, @@ -217,23 +273,33 @@ } } - /// Iterate over all nodes and corresponding subcubes of self. + /* + /// Returns an iterator over all nodes and corresponding subcubes of `self`. #[inline] - pub fn nodes_and_cubes<'a, 'b>(&'a self, domain : &'b Cube<F, N>) + pub(super) fn nodes_and_cubes<'a, 'b>(&'a self, domain : &'b Cube<F, N>) -> std::iter::Zip<Iter<'a, Node<F,D,A,N,P>>, SubcubeIter<'b, F, N, P>> { self.nodes.iter().zip(self.iter_subcubes(domain)) } + */ - /// Mutably iterate over all nodes and corresponding subcubes of self. + /// Mutably iterate over all nodes and corresponding subcubes of `self`. #[inline] - pub fn nodes_and_cubes_mut<'a, 'b>(&'a mut self, domain : &'b Cube<F, N>) + 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) } /// Insert data into the branch. - pub fn insert<M : Depth, S : LocalAnalysis<F, A, N>>( + /// + /// The parameters are as follows: + /// * `domain` is the cube corresponding to this branch. + /// * `d` is the data to be inserted + /// * `new_leaf_depth` is the depth relative to `self` at which the data is to be inserted. + /// * `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<M : Depth, S : LocalAnalysis<F, A, N>>( &mut self, domain : &Cube<F,N>, d : D, @@ -253,8 +319,13 @@ } } - /// Construct a new instance for a different aggregator - pub fn convert_aggregator<ANew, G>( + /// Construct a new instance of the branch for a different aggregator. + /// + /// The `generator` is used to convert the data of type `D` of the branch into corresponding + /// [`Support`]s. The `domain` is the cube corresponding to `self`. + /// The type parameter `ANew´ is the new aggregator, and needs to be implemented for the + /// generator's `SupportType`. + pub(super) fn convert_aggregator<ANew, G>( self, generator : &G, domain : &Cube<F, N> @@ -274,8 +345,11 @@ } } - /// Recalculate aggregator after changes to generator - pub fn refresh_aggregator<G>( + /// Recalculate aggregator after changes to generator. + /// + /// The `generator` is used to convert the data of type `D` of the branch into corresponding + /// [`Support`]s. The `domain` is the cube corresponding to `self`. + pub(super) fn refresh_aggregator<G>( &mut self, generator : &G, domain : &Cube<F, N> @@ -292,8 +366,9 @@ where Const<P> : BranchCount<N>, A : Aggregator { + /// Create a new node #[inline] - pub fn new() -> Self { + pub(super) fn new() -> Self { Node { data : NodeOption::Uninitialised, aggregator : A::new(), @@ -302,7 +377,7 @@ /// Get leaf data #[inline] - pub fn get_leaf_data(&self, x : &Loc<F, N>) -> Option<&Vec<D>> { + pub(super) fn get_leaf_data(&self, x : &Loc<F, N>) -> Option<&Vec<D>> { match self.data { NodeOption::Uninitialised => None, NodeOption::Leaf(ref data) => Some(data), @@ -312,13 +387,24 @@ /// Returns a reference to the aggregator of this node #[inline] - pub fn get_aggregator(&self) -> &A { + pub(super) fn get_aggregator(&self) -> &A { &self.aggregator } - /// Insert `d` into the tree. If an `Incomplete` node is encountered, a new - /// leaf is created at a minimum depth of `new_leaf_depth` - pub fn insert<M : Depth, S : LocalAnalysis <F, A, N>>( + /// Insert data under the node. + /// + /// The parameters are as follows: + /// * `domain` is the cube corresponding to this branch. + /// * `d` is the data to be inserted + /// * `new_leaf_depth` is the depth relative to `self` at which new leaves are created. + /// * `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. + /// + /// If `self` is already [`NodeOption::Leaf`], the data is inserted directly in this node. + /// 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<M : Depth, S : LocalAnalysis <F, A, N>>( &mut self, domain : &Cube<F,N>, d : D, @@ -361,8 +447,13 @@ } } - /// Construct a new instance for a different aggregator - pub fn convert_aggregator<ANew, G>( + /// Construct a new instance of the node for a different aggregator + /// + /// The `generator` is used to convert the data of type `D` of the node into corresponding + /// [`Support`]s. The `domain` is the cube corresponding to `self`. + /// The type parameter `ANew´ is the new aggregator, and needs to be implemented for the + /// generator's `SupportType`. + pub(super) fn convert_aggregator<ANew, G>( mut self, generator : &G, domain : &Cube<F, N> @@ -402,8 +493,11 @@ } } - /// Refresh aggregator after changes to generator - pub fn refresh_aggregator<G>( + /// Refresh aggregator after changes to generator. + /// + /// The `generator` is used to convert the data of type `D` of the node into corresponding + /// [`Support`]s. The `domain` is the cube corresponding to `self`. + pub(super) fn refresh_aggregator<G>( &mut self, generator : &G, domain : &Cube<F, N> @@ -427,21 +521,10 @@ } } -impl<'a, D> Iterator for BTIter<'a,D> { - type Item = &'a D; - #[inline] - fn next(&mut self) -> Option<&'a D> { - self.iter.next() - } -} - - -// -// BT -// - -/// Internal structure to hide the `const P : usize` parameter of [`Node`] until -/// const generics are flexible enough to fix `P=pow(2, N)`. +/// Helper trait for working with [`Node`]s without the knowledge of `P`. +/// +/// 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, @@ -449,24 +532,41 @@ type Node : Clone + std::fmt::Debug; } +/// Helper structure for looking up a [`Node`] without the knowledge of `P`. +/// +/// This can be removed once Rust's const generics are flexible enough to allow fixing +/// `P=pow(2, N)`. #[derive(Debug)] pub struct BTNodeLookup; -/// Interface to a [`BT`] bisection tree. +/// 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> { + /// The data type stored in the tree type Data : 'static + Copy; + /// The depth type of the tree type Depth : Depth; + /// The type for the [aggregate information][Aggregator] about the `Data` stored in each node + /// of the tree. 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; - /// Insert d into the `BisectionTree`. + /// 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>>( &mut self, d : Self::Data, support : &S ); - /// Construct a new instance for a different aggregator + /// 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, @@ -474,21 +574,26 @@ G::SupportType : LocalAnalysis<F, ANew, N>; - /// Refresh aggregator after changes to generator + /// 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>; - /// Iterarate items at x + /// Iterarate all [`Self::Data`] items at the point `x` of the domain. fn iter_at<'a>(&'a self, x : &'a Loc<F,N>) -> BTIter<'a, Self::Data>; - /// Create a new instance + /// Create a new tree on `domain` of indicated `depth`. fn new(domain : Cube<F, N>, depth : Self::Depth) -> Self; } -/// The main bisection tree structure. The interface operations are via [`BTImpl`] -/// to hide the `const P : usize` parameter until const generics are flexible enough -/// to fix `P=pow(2, N)`. +/// The main bisection tree structure. +/// +/// 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, @@ -497,8 +602,11 @@ A : Aggregator, const N : usize, > where BTNodeLookup : BTNode<F, D, A, N> { + /// The depth of the tree (initial, before refinement) pub(super) depth : M, + /// The domain of the toplevel node pub(super) domain : Cube<F, N>, + /// The toplevel node of the tree pub(super) topnode : <BTNodeLookup as BTNode<F, D, A, N>>::Node, } @@ -521,7 +629,6 @@ type Agg = A; type Converted<ANew> = BT<M,F,D,ANew,$n> where ANew : Aggregator; - /// Insert `d` into the tree. fn insert<S: LocalAnalysis<F, A, $n>>( &mut self, d : D, @@ -535,7 +642,6 @@ ); } - /// Construct a new instance for a different aggregator fn convert_aggregator<ANew, G>(self, generator : &G) -> Self::Converted<ANew> where ANew : Aggregator, G : SupportGenerator<F, $n, Id=D>, @@ -548,14 +654,12 @@ } } - /// Refresh aggregator after changes to generator fn refresh_aggregator<G>(&mut self, generator : &G) where G : SupportGenerator<F, $n, Id=Self::Data>, G::SupportType : LocalAnalysis<F, Self::Agg, $n> { self.topnode.refresh_aggregator(generator, &self.domain); } - /// Iterate elements at `x`. fn iter_at<'a>(&'a self, x : &'a Loc<F,$n>) -> BTIter<'a,D> { match self.topnode.get_leaf_data(x) { Some(data) => BTIter { iter : data.iter() },
--- a/src/bisection_tree/btfn.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree/btfn.rs Sat Oct 22 14:40:28 2022 +0300 @@ -2,13 +2,12 @@ use numeric_literals::replace_float_literals; use std::iter::Sum; use std::marker::PhantomData; -pub use nalgebra::Const; use crate::types::Float; use crate::mapping::Mapping; use crate::linops::Linear; use crate::sets::Set; -pub use crate::sets::Cube; -pub use crate::loc::Loc; +use crate::sets::Cube; +use crate::loc::Loc; use super::support::*; use super::bt::*; use super::refine::*; @@ -17,13 +16,16 @@ use crate::fe_model::base::RealLocalModel; use crate::fe_model::p2_local_model::*; -/// Presentation for functions constructed as a sum of components with limited support. -/// Lookup of relevant components at a specific point of a [`Loc<F,N>`] domain is done -/// with a bisection tree (BT). We do notn impose that the `BT` parameter would be an in -/// instance [`BTImpl`] to allow performance improvements via “pre-BTFNs” that do not -/// construct a BT until being added to another function. Addition and subtraction always -/// use a copy-on-write clone of the BT of the left-hand-side operand, discarding -/// the BT of the right-hand-side operand. +/// Presentation for (mathematical) functions constructed as a sum of components functions with +/// typically small support. +/// +/// The domain of the function is [`Loc`]`<F, N>`, where `F` is the type of floating point numbers, +/// and `N` the dimension. +/// +/// The `generator` lists the component functions that have to implement [`Support`]. +/// Identifiers of the components ([`SupportGenerator::Id`], usually `usize`) are stored stored +/// in a [bisection tree][BTImpl], when one is provided as `bt`. However `bt` may also be `()` +/// for a [`PreBTFN`] that is only useful for vector space operations with a full [`BTFN`]. #[derive(Clone,Debug)] pub struct BTFN< F : Float, @@ -42,16 +44,13 @@ G::SupportType : LocalAnalysis<F, BT::Agg, N>, BT : BTImpl<F, N> { - /// Returns the [`Support`] corresponding to the identifier `d`. - /// Simply wraps [`SupportGenerator::support_for`]. - #[inline] - fn support_for(&self, d : G::Id) -> G::SupportType { - self.generator.support_for(d) - } - - /// Create a new BTFN. The bisection tree `bt` should be pre-initialised to - /// correspond to the `generator`. Use [`Self::construct`] if no preinitialised tree - /// is available. Use [`Self::new_refresh`] when the aggregator may need updates. + /// Create a new BTFN from a support generator and a pre-initialised bisection tree. + /// + /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`. + /// Use [`Self::construct`] if no preinitialised tree is available. Use [`Self::new_refresh`] + /// when the aggregators of the tree may need updates. + /// + /// See the documentation for [`BTFN`] on the role of the `generator`. pub fn new(bt : BT, generator : G) -> Self { BTFN { bt : bt, @@ -60,8 +59,13 @@ } } - /// Create a new BTFN. The bisection tree `bt` should be pre-initialised to - /// correspond to the `generator`, but the aggregator may be out of date. + /// Create a new BTFN support generator and a pre-initialised bisection tree, + /// cloning the tree and refreshing aggregators. + /// + /// The bisection tree `bt` should be pre-initialised to correspond to the `generator`, but + /// the aggregator may be out of date. + /// + /// See the documentation for [`BTFN`] on the role of the `generator`. pub fn new_refresh(bt : &BT, generator : G) -> Self { // clone().refresh_aggregator(…) as opposed to convert_aggregator // ensures that type is maintained. Due to Rc-pointer copy-on-write, @@ -71,8 +75,11 @@ BTFN::new(btnew, generator) } - /// Create a new BTFN. Unlike [`Self::new`], this constructs the bisection tree based on - /// the generator. + /// Create a new BTFN from a support generator, domain, and depth for a new [`BT`]. + /// + /// The top node of the created [`BT`] will have the given `domain`. + /// + /// See the documentation for [`BTFN`] on the role of the `generator`. pub fn construct(domain : Cube<F, N>, depth : BT::Depth, generator : G) -> Self { let mut bt = BT::new(domain, depth); for (d, support) in generator.all_data() { @@ -81,7 +88,10 @@ Self::new(bt, generator) } - /// Construct a new instance for a different aggregator + /// Convert the aggregator of the [`BTFN`] to a different one. + /// + /// This will construct a [`BTFN`] with the same components and generator as the (consumed) + /// `self`, but a new `BT` with [`Aggregator`]s of type `ANew`. pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N> where ANew : Aggregator, G : SupportGenerator<F, N, Id=BT::Data>, @@ -89,16 +99,8 @@ BTFN::new(self.bt.convert_aggregator(&self.generator), self.generator) } - /// Change the [`BTImpl`]. - pub fn instantiate< - BTNew : BTImpl<F, N, Data=BT::Data>, - > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<F, G, BTNew, N> - where G::SupportType : LocalAnalysis<F, BTNew::Agg, N> { - BTFN::construct(domain, depth, self.generator) - } - /// Change the generator (after, e.g., a scaling of the latter). - pub fn new_generator(&self, generator : G) -> Self { + fn new_generator(&self, generator : G) -> Self { BTFN::new_refresh(&self.bt, generator) } @@ -109,13 +111,31 @@ } -/// A “pre-BTFN” with no bisection tree. +impl<F : Float, G, BT, const N : usize> +BTFN<F, G, BT, N> +where G : SupportGenerator<F, N> { + /// Change the [bisection tree][BTImpl] of the [`BTFN`] to a different one. + /// + /// This can be used to convert a [`PreBTFN`] to a full [`BTFN`], or the change + /// the aggreagator; see also [`self.convert_aggregator`]. + pub fn instantiate< + BTNew : BTImpl<F, N, Data=G::Id>, + > (self, domain : Cube<F, N>, depth : BTNew::Depth) -> BTFN<F, G, BTNew, N> + where G::SupportType : LocalAnalysis<F, BTNew::Agg, N> { + BTFN::construct(domain, depth, self.generator) + } +} + +/// A BTFN with no bisection tree. +/// +/// Most BTFN methods are not available, but if a BTFN is going to be summed with another +/// before other use, it will be more efficient to not construct an unnecessary bisection tree +/// that would be shortly dropped. pub type PreBTFN<F, G, const N : usize> = BTFN<F, G, (), N>; impl<F : Float, G, const N : usize> PreBTFN<F, G, N> where G : SupportGenerator<F, N> { - /// Create a new “pre-BTFN” with no bisection tree. This can feasibly only be - /// multiplied by a scalar or added to or subtracted from a proper [`BTFN`]. + /// Create a new [`PreBTFN`] with no bisection tree. pub fn new_pre(generator : G) -> Self { BTFN { bt : (), @@ -214,14 +234,12 @@ G2 : SupportGenerator<F, N, Id=usize> + Clone, G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, G2::SupportType : LocalAnalysis<F, BT1::Agg, N>, - /*&'b G2 : std::ops::Neg<Output=G2>*/ { + &'b G2 : std::ops::Neg<Output=G2> { type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; #[inline] fn sub(self, other : &'b BTFN<F, G2, BT2, N>) -> Self::Output { - // FIXME: the compiler crashes in an overflow on this. - //$preprocess(self).add_another(std::ops::Neg::neg(&other.generator)) - $preprocess(self).add_another(other.generator.clone().neg()) + $preprocess(self).add_another(-&other.generator) } } } @@ -371,7 +389,7 @@ type Codomain = V; fn value(&self, x : &'a Loc<F,N>) -> Self::Codomain { - self.bt.iter_at(x).map(|&d| self.support_for(d).value(x)).sum() + self.bt.iter_at(x).map(|&d| self.generator.support_for(d).value(x)).sum() } } @@ -385,7 +403,7 @@ type Codomain = V; fn value(&self, x : Loc<F,N>) -> Self::Codomain { - self.bt.iter_at(&x).map(|&d| self.support_for(d).value(x)).sum() + self.bt.iter_at(&x).map(|&d| self.generator.support_for(d).value(x)).sum() } } @@ -421,7 +439,13 @@ } /// Helper trait for performing approximate minimisation using P2 elements. +/// +/// `U` is the domain, generally [`Loc`]`<F, N>`, and `F` the type of floating point numbers. +/// `Self` is generally a set of `U`, for example, [`Cube`]`<F, N>`. pub trait P2Minimise<U, F : Float> : Set<U> { + /// Minimise `g` over the set presented by `Self`. + /// + /// The function returns `(x, v)` where `x` is the minimiser `v` an approximation of `g(x)`. fn p2_minimise<G : Fn(&U) -> F>(&self, g : G) -> (U, F); } @@ -450,7 +474,7 @@ 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, 3, 3>::new( + let m1 = P2LocalModel::<F, 2, 3>::new( &[a, b, c, ab, bc, ca], &[va, vb, vc, vab, vbc, vca] ); @@ -458,7 +482,7 @@ let r1@(_, v1) = m1.minimise(&s1); let s2 = Simplex([c, d, a]); - let m2 = P2LocalModel::<F, 2, 3, 3, 3>::new( + let m2 = P2LocalModel::<F, 2, 3>::new( &[c, d, a, cd, da, ca], &[vc, vd, va, vcd, vda, vca] ); @@ -489,7 +513,7 @@ }; let [ve, vf] = [g(&e), g(&f)]; - let m1 = P2LocalModel::<F, 2, 3, 3, 3>::new( + let m1 = P2LocalModel::<F, 2, 3>::new( &[a, b, c, d, e, f], &[va, vb, vc, vd, ve, vf], ); @@ -499,18 +523,24 @@ } } +/// Helper type to use [`P2Refiner`] for maximisation. struct RefineMax; + +/// Helper type to use [`P2Refiner`] for minimisation. struct RefineMin; -/// TODO: update doc. -/// Simple refiner that keeps the difference between upper and lower bounds -/// of [`Bounded`] functions within $ε$. The functions have to also be continuous -/// for this refiner to ever succeed. +/// A bisection tree [`Refiner`] for maximising or minimising a [`BTFN`]. +/// /// The type parameter `T` should be either [`RefineMax`] or [`RefineMin`]. struct P2Refiner<F : Float, T> { + /// The maximum / minimum should be above / below this threshold. + /// If the threshold cannot be satisfied, the refiner will return `None`. bound : Option<F>, + /// Tolerance for function value estimation. tolerance : F, + /// Maximum number of steps to execute the refiner for max_steps : usize, + /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes. #[allow(dead_code)] // `how` is just for type system purposes. how : T, } @@ -602,12 +632,21 @@ } +/// A bisection tree [`Refiner`] for checking that a [`BTFN`] is within a stated +//// upper or lower bound. +/// +/// The type parameter `T` should be either [`RefineMax`] for upper bound or [`RefineMin`] +/// for lower bound. struct BoundRefiner<F : Float, T> { + /// The upper/lower bound to check for bound : F, + /// Tolerance for function value estimation. tolerance : F, + /// Maximum number of steps to execute the refiner for max_steps : usize, #[allow(dead_code)] // `how` is just for type system purposes. + /// Either [`RefineMax`] or [`RefineMin`]. Used only for type system purposes. how : T, } @@ -641,6 +680,35 @@ } } +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 : &Vec<G::Id>, + _generator : &G, + step : usize + ) -> RefinerResult<Bounds<F>, Self::Result> { + if aggregator.lower() >= self.bound - self.tolerance { + // Above lower bound within tolerances. Indicate uncertain success. + RefinerResult::Uncertain(*aggregator, true) + } else if aggregator.upper() <= self.bound + self.tolerance { + // Below lower bound within tolerances. Indicate certain failure. + RefinerResult::Certain(false) + } else if step < self.max_steps { + // No decision possible, but within step bounds - further subdivision is required. + RefinerResult::NeedRefinement + } else { + // No decision possible, but past step bounds + RefinerResult::Uncertain(*aggregator, false) + } + } +} impl<F : Float, G, BT, const N : usize> BTFN<F, G, BT, N> where BT : BTSearch<F, N, Agg=Bounds<F>>, @@ -649,41 +717,59 @@ + LocalAnalysis<F, Bounds<F>, N>, Cube<F, N> : P2Minimise<Loc<F, N>, F> { - /// Try to maximise the BTFN within stated value `tolerance`, taking at most `max_steps` - /// refinement steps. Returns the maximiser and the maximum value. + /// 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() } - /// Try to maximise the BTFN within stated value `tolerance`, taking at most `max_steps` - /// refinement steps. Returns the maximiser and the maximum value when one is found above - /// the `bound` threshold, otherwise `None`. + /// 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.") } - /// Try to minimise the BTFN within stated value `tolerance`, taking at most `max_steps` - /// refinement steps. Returns the minimiser and the minimum value. + /// 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() } - /// Try to minimise the BTFN within stated value `tolerance`, taking at most `max_steps` - /// refinement steps. Returns the minimiser and the minimum value when one is found below - /// the `bound` threshold, otherwise `None`. + /// Minimise the `BTFN` within stated value `tolerance` subject to a lower bound. + /// + /// At most `max_steps` refinement steps are taken. + /// Returns the approximate minimiser and the corresponding function value when one is found + /// above the `bound` threshold, otherwise `None`. pub fn minimise_below(&mut self, bound : F, tolerance : F, max_steps : usize) -> Option<(Loc<F, N>, F)> { let refiner = P2Refiner{ tolerance, max_steps, how : RefineMin, bound : Some(bound) }; self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.") } - /// Verify that the BTFN is has stated upper `bound` within stated `tolerance, taking at most - /// `max_steps` refinment steps. + + /// Verify that the `BTFN` has a given upper `bound` within indicated `tolerance`. + /// + /// At most `max_steps` refinement steps are taken. pub fn has_upper_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool { let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMax }; self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.") } + + /// Verify that the `BTFN` has a given lower `bound` within indicated `tolerance`. + /// + /// At most `max_steps` refinement steps are taken. + pub fn has_lower_bound(&mut self, bound : F, tolerance : F, max_steps : usize) -> bool { + let refiner = BoundRefiner{ bound, tolerance, max_steps, how : RefineMin }; + self.bt.search_and_refine(&refiner, &self.generator).expect("Refiner failure.") + } }
--- a/src/bisection_tree/either.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree/either.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,18 +1,27 @@ + +use std::iter::Chain; use crate::types::*; +use crate::mapping::Mapping; +use crate::iter::{Mappable,MapF,MapZ}; +use crate::sets::Cube; +use crate::loc::Loc; + use super::support::*; -use crate::mapping::Mapping; -use std::iter::Chain; -use crate::iter::{Mappable,MapF,MapZ}; use super::aggregator::*; -// TODO: try to remove F and N +/// A structure for storing two [`SupportGenerator`]s summed/chain together. +/// +/// This is needed to work with sums of different types of [`Support`]s. #[derive(Debug,Clone)] pub struct BothGenerators<A, B>( pub(super) A, pub(super) B, ); +/// A structure for a [`Support`] that can be either `A` or `B`. +/// +/// This is needed to work with sums of different types of [`Support`]s. #[derive(Debug,Clone)] pub enum EitherSupport<A, B> { Left(A), @@ -32,7 +41,7 @@ >; impl<G1, G2> BothGenerators<G1, G2> { - /// Helper for [`all_data`]. + /// Helper for [`all_left_data`]. #[inline] fn map_left<F : Float, const N : usize>((d, support) : (G1::Id, G1::SupportType)) -> (usize, EitherSupport<G1::SupportType, G2::SupportType>) @@ -43,7 +52,7 @@ (id.into(), EitherSupport::Left(support)) } - /// Helper for [`all_data`]. + /// Helper for [`all_right_data`]. #[inline] fn map_right<F : Float, const N : usize>(n0 : &usize, (d, support) : (G2::Id, G2::SupportType)) -> (usize, EitherSupport<G1::SupportType, G2::SupportType>) @@ -54,6 +63,9 @@ ((n0+id).into(), EitherSupport::Right(support)) } + /// Calls [`SupportGenerator::all_data`] on the “left” support generator. + /// + /// Converts both the id and the [`Support`] into a form that corresponds to `BothGenerators`. #[inline] pub(super) fn all_left_data<F : Float, const N : usize>(&self) -> MapF<G1::AllDataIter<'_>, (usize, EitherSupport<G1::SupportType, G2::SupportType>)> @@ -62,6 +74,9 @@ self.0.all_data().mapF(Self::map_left) } + /// Calls [`SupportGenerator::all_data`] on the “right” support generator. + /// + /// Converts both the id and the [`Support`] into a form that corresponds to `BothGenerators`. #[inline] pub(super) fn all_right_data<F : Float, const N : usize>(&self) -> MapZ<G2::AllDataIter<'_>, usize, (usize, EitherSupport<G1::SupportType, G2::SupportType>)>
--- a/src/bisection_tree/refine.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree/refine.rs Sat Oct 22 14:40:28 2022 +0300 @@ -4,15 +4,17 @@ use std::rc::Rc; use std::marker::PhantomData; use crate::types::*; +use crate::nanleast::NaNLeast; +use crate::sets::Cube; use super::support::*; use super::bt::*; use super::aggregator::*; -use crate::nanleast::NaNLeast; /// Trait for sorting [`Aggregator`]s for [`BT`] refinement. /// -/// The sorting involves two sorting keys, the “upper” and the “lower” key. Any [`BT`] branches +/// The sorting involves two sorting keys, the “upper” and the “lower” key. Any [`BT`] nodes /// with upper key less the lower key of another are discarded from the refinement process. +/// Nodes with the highest upper sorting key are picked for refinement. pub trait AggregatorSorting { // Priority type Agg : Aggregator; @@ -24,7 +26,7 @@ /// Returns upper sorting key fn sort_upper(aggregator : &Self::Agg) -> Self::Sort; - /// Returns bottom sorting key. + /// Returns a sorting key that is less than any other sorting key. fn bottom() -> Self::Sort; } @@ -67,10 +69,11 @@ fn bottom() -> Self::Sort { NaNLeast(F::NEG_INFINITY) } } -/// Result type of [`Refiner::refine`] for a refiner producing a result of type `R` acting on -/// an [`Aggregator`] of type `A`. +/// Return type of [`Refiner::refine`]. +/// +/// The parameter `R` is the result type of the refiner acting on an [`Aggregator`] of type `A`. pub enum RefinerResult<A : Aggregator, R> { - /// Indicates in insufficiently refined state: the [`BT`] needs to be further refined. + /// Indicates an insufficiently refined state: the [`BT`] needs to be further refined. NeedRefinement, /// Indicates a certain result `R`, stop refinement immediately. Certain(R), @@ -81,20 +84,42 @@ use RefinerResult::*; -/// A `Refiner` is used to determine whether an [`Aggregator`] `A` is sufficiently refined within -/// a [`Cube`] of a [`BT`], and in such a case, produce a desired result (e.g. a maximum value of -/// a function). +/// A `Refiner` is used to search a [`BT`], refining the subdivision when necessary. +/// +/// The search is performed by [`BTSearch::search_and_refine`]. +/// The `Refiner` is used to determine whether an [`Aggregator`] `A` stored in the [`BT`] is +/// sufficiently refined within a [`Cube`], and in such a case, produce a desired result (e.g. +/// a maximum value of a function). pub trait Refiner<F : Float, A, G, const N : usize> where F : Num, A : Aggregator, G : SupportGenerator<F, N> { + /// The result type of the refiner type Result : std::fmt::Debug; + /// The sorting to be employed by [`BTSearch::search_and_refine`] on node aggregators + /// to detemrine node priority. type Sorting : AggregatorSorting<Agg = A>; - /// Determines whether `aggregator` is sufficiently refined within `cube`. - /// Should return a possibly refined version of the `aggregator` and an arbitrary value of - /// the result type of the refiner. + /// Determines whether `aggregator` is sufficiently refined within `domain`. + /// + /// If the aggregator is sufficiently refined that the desired `Self::Result` can be produced, + /// a [`RefinerResult`]`::Certain` or `Uncertain` should be returned, depending on + /// the confidence of the solution. In the uncertain case an improved aggregator should also + /// be included. If the result cannot be produced, `NeedRefinement` should be + /// returned. + /// + /// For example, if the refiner is used to minimise a function presented by the `BT`, + /// an `Uncertain` result can be used to return a local maximum of the function on `domain` + /// The result can be claimed `Certain` if it is a global maximum. In that case the + /// refinment will stop immediately. A `NeedRefinement` result indicates that the `aggregator` + /// and/or `domain` are not sufficiently refined to compute a lcoal maximum of sufficient + /// quality. + /// + /// The vector `data` stored all the data of the [`BT`] in the node corresponding to `domain`. + /// The `generator` can be used to convert `data` into [`Support`]s. The parameter `step` + /// counts the calls to `refine`, and can be used to stop the refinement when a maximum + /// number of steps is reached. fn refine( &self, aggregator : &A, @@ -191,7 +216,10 @@ } } -pub struct HeapContainer<'a, F, D, A, S, RResult, const N : usize, const P : usize> +/// This is a container for a [`BinaryHeap`] of [`RefinementInfo`]s together with tracking of +/// the greatest lower bound of the [`Aggregator`]s of the [`Node`]s therein accroding to +/// chosen [`AggregatorSorting`]. +struct HeapContainer<'a, F, D, A, S, RResult, const N : usize, const P : usize> where F : Float, D : 'static + Copy, Const<P> : BranchCount<N>, @@ -211,6 +239,7 @@ A : Aggregator, S : AggregatorSorting<Agg = A> { + /// Push `ri` into the [`BinaryHeap`]. Do greatest lower bound maintenance. fn push(&mut self, ri : RefinementInfo<'a, F, D, A, S, RResult, N, P>) { if ri.sort_upper() >= self.glb { let l = ri.sort_lower(); @@ -258,6 +287,8 @@ /// If [`RefinerResult::Uncertain`] is returned, the leaf is inserted back into the refinement /// queue `container`. If `self` is a branch, its subnodes are staged into `container` using /// [`Branches::stage_refine`]. + /// + /// `domain`, as usual, indicates the spatial area corresponding to `self`. fn search_and_refine<'a, 'b, R, G>( &'a mut self, domain : Cube<F,N>, @@ -329,16 +360,22 @@ } } -/// Helper trait for implementing a refining search on a [`BT`]. +/// Interface trait to a refining search on a [`BT`]. +/// +/// 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 BTSearch<F, const N : usize> : BTImpl<F, N> where F : Float { - /// Perform a refining search on [`Self`], as determined by `refiner`. Nodes are inserted - /// in a priority queue and processed in the order determined by the [`AggregatorSorting`] - /// [`Refiner::Sorting`]. Leaf nodes are subdivided until the refiner decides that a - /// sufficiently refined leaf node has been found, as determined by either the refiner - /// returning a [`RefinerResult::Certain`] result, or a previous [`RefinerResult::Uncertain`] - /// result is found again at the top of the priority queue. + /// Perform a search on [`Self`], as determined by `refiner`. + /// + /// Nodes are inserted in a priority queue and processed in the order determined by the + /// [`AggregatorSorting`] [`Refiner::Sorting`]. Leaf nodes are subdivided until the refiner + /// decides that a sufficiently refined leaf node has been found, as determined by either the + /// refiner returning a [`RefinerResult::Certain`] result, or a previous + /// [`RefinerResult::Uncertain`] result is found again at the top of the priority queue. + /// + /// The `generator` converts [`BTImpl::Data`] stored in the bisection tree into a [`Support`]. fn search_and_refine<'b, R, G>( &'b mut self, refiner : &R,
--- a/src/bisection_tree/support.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree/support.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,18 +1,22 @@ +/*! +Traits for representing the support of a [`Mapping`], and analysing the mapping on a [`Cube`]. +*/ use serde::Serialize; use std::ops::{MulAssign,DivAssign,Neg}; -pub use nalgebra::Const; -use crate::types::{Float,Num}; -use crate::maputil::{array_init,map2}; -use crate::mapping::{Mapping}; -pub use crate::sets::cube::Cube; -pub use crate::loc::Loc; +use crate::types::{Float, Num}; +use crate::maputil::map2; +use crate::mapping::Mapping; +use crate::sets::Cube; +use crate::loc::Loc; use super::aggregator::Bounds; use crate::norms::{Norm, L1, L2, Linfinity}; -/// A trait for encoding constant values +/// A trait for encoding constant [`Float`] values pub trait Constant : Copy + 'static + std::fmt::Debug + Into<Self::Type> { + /// The type of the value type Type : Float; + /// Returns the value of the constant fn value(&self) -> Self::Type; } @@ -24,41 +28,58 @@ /// 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 { - /// Return a cube containing the support of the function. - /// It may be larger than the actual support, but must contain it. + /// 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>; - /// Indicate whether `x` is in the support of the function. + /// Indicate whether `x` is in the support of the function represented by `self`. fn in_support(&self, x : &Loc<F,N>) -> bool; - // Indicate whether `cube` is fully in the support of the function. + // Indicate whether `cube` is fully in the support of the function represented by `self`. //fn fully_in_support(&self, cube : &Cube<F,N>) -> bool; - /// Return an optional hint for bisecting the support. Useful for nonsmooth - /// functions and finite element approximations to make the bisections compatible - /// with the grid or points of non-differentiability. + /// Return an optional hint for bisecting the support. + /// + /// The output along each axis a possible coordinate at which to bisect `cube`. + /// + /// This is useful for nonsmooth functions to make finite element models as used by + /// [`BTFN`][super::btfn::BTFN] minimisation/maximisation compatible with points of + /// non-differentiability. + /// + /// The default implementation returns `[None; N]`. #[inline] - fn bisection_hint(&self, _cube : &Cube<F, N>) -> [Option<F>; N] { - array_init(|| None) + #[allow(unused_variables)] + fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] { + [None; N] } - /// Shift the support by `x`. + /// Translate `self` by `x`. #[inline] fn shift(self, x : Loc<F, N>) -> Shift<Self, F, N> { Shift { shift : x, base_fn : self } } - /// Weight any corresponding [`Mapping`] + /// Multiply `self` by the scalar `a`. #[inline] fn weigh<C : Constant<Type=F>>(self, a : C) -> Weighted<Self, C> { Weighted { weight : a, base_fn : self } } } -/// Trait for globally analysing a certain property of a [`Mapping`]. +/// 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; } @@ -70,13 +91,22 @@ // } // } -/// Trait for locally analysing a certain property `A` of a [`Support`] or [`Mapping`] -/// within a [`Cube<F, N>`]. +/// 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 a [`Mapping`] on [`Loc<F, N>`]. +/// 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. @@ -90,7 +120,7 @@ impl<F : Float, T : GlobalAnalysis<F, Bounds<F>>> Bounded<F> for T { } -/// Shift of [`Support`] and [`Bounded`]. +/// 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>, @@ -170,10 +200,13 @@ impl_shift_norm!(L1 L2 Linfinity); -/// Weighting of [`Support`] and [`Bounded`]. +/// Weighting of a [`Support`] and [`Mapping`] by scalar multiplication; +/// output of [`Support::weigh`]. #[derive(Copy,Clone,Debug,Serialize)] pub struct Weighted<T, C : Constant> { + /// The weight pub weight : C, + /// The base [`Support`] or [`Mapping`] being weighted. pub base_fn : T, } @@ -296,10 +329,14 @@ impl_weighted_norm!(L1 L2 Linfinity); -/// Normalisation of [`Support`] and [`Bounded`] to L¹ norm 1. +/// 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>(pub T); +pub struct Normalised<T>( + /// The base [`Support`] or [`Mapping`]. + pub T +); impl<'a, T, F : Float, const N : usize> Mapping<&'a Loc<F, N>> for Normalised<T> where T : Norm<F, L1> + for<'b> Mapping<&'b Loc<F,N>, Codomain=F> { @@ -402,19 +439,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> { + /// The identification type type Id : 'static + Copy; + /// The type of the [`Support`] (often also a [`Mapping`]). 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; - /// Gives the [`Support`] for component function identified by `d`. - fn support_for(&self, d : Self::Id) -> Self::SupportType; + /// Returns the component identified by `id`. + /// + /// Panics if `id` is an invalid identifier. + fn support_for(&self, id : Self::Id) -> Self::SupportType; - /// Returns the number of different component functions in this generator. + /// Returns the number of different components in this generator. fn support_count(&self) -> usize; - /// Iterator over all pairs of (id, support). + /// Returns an iterator over all pairs of `(id, support)`. fn all_data(&self) -> Self::AllDataIter<'_>; }
--- a/src/bisection_tree/supportid.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/bisection_tree/supportid.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,5 +1,7 @@ use std::marker::PhantomData; +/// Standard identifier for a [`Support`][super::support::Support] stored in a +/// [`BTImpl`][super::bt::BTImpl] as [`BTImpl::Data`][super::bt::BTImpl::Data]. #[derive(Debug)] pub struct SupportId<S> { pub id : usize,
--- a/src/coefficients.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/coefficients.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,8 +1,10 @@ -///! This module implements some arithmetic functions that are valid in a ‘const context”, i.e., -///! can be computed at compile-time. +/*! +Arithmetic functions that are valid in a ‘const context”, i.e., can be computed at compile-time. +*/ -/// Calculate $n!$. -/// (The version in the crate `num_integer` deesn't currently work in a `const` context. +/// Calculate $n!$. Valid in a const-context. +/// +/// (The version in the crate `num_integer` deesn't currently work in a const-context.) pub const fn factorial(n : usize) -> usize { const fn f(i : usize, a : usize) -> usize { if i==0 { a } else { f(i-1, i*a) } @@ -10,14 +12,16 @@ f(n, 1) } -/// Calculate $\nchoosek{n}{k}$. -/// (The version in the crate `num_integer` deesn't currently work in a `const` context. +/// Calculate $n \choose k$. Valid in a const-context. +/// +/// (The version in the crate `num_integer` doesn't currently work in a const-context.) pub const fn binomial(n : usize, k : usize) -> usize { factorial(n)/(factorial(n-k)*factorial(k)) } -/// Calculate $n^k$. -/// (The version in the crate `num_integer` deesn't currently work in a `const` context. +/// Calculate $n^k$. Valid in a const-context. +/// +/// (The version in the crate `num_integer` doesn't currently work in a const-context.) pub const fn pow(n : usize, k : usize) -> usize { const fn pow_accum(n : usize, k : usize, accum : usize) -> usize { if k==0 { accum } else { pow_accum(n, k-1, accum*n) }
--- a/src/error.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/error.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,7 +1,13 @@ +/*! +Error passing helper types +*/ use std::error::Error; +/// A [`Result`] containing `T` or a dynamic error type pub type DynResult<T> = Result<T, Box<dyn Error>>; + +/// A [`Result`] containing `()` or a dynamic error type pub type DynError = DynResult<()>; #[derive(Clone, Debug)]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/euclidean.rs Sat Oct 22 14:40:28 2022 +0300 @@ -0,0 +1,85 @@ +/*! +Euclidean spaces. +*/ + +use crate::types::*; +use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg}; + +/// Space (type) with a defined dot product. +/// +/// `U` is the space of the multiplier, and `F` the space of scalars. +/// Since `U` ≠ `Self`, this trait can also implement dual products. +pub trait Dot<U, F> { + fn dot(&self, other : &U) -> F; +} + +/// Space (type) with Euclidean and vector space structure +/// +/// The type should implement vector space operations (addition, subtraction, scalar +/// multiplication and scalar division) along with their assignment versions, as well +/// as the [`Dot`] product with respect to `Self`. +pub trait Euclidean<F : Float> : Sized + Dot<Self,F> + + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F> + + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F> + + Add<Self, Output=<Self as Euclidean<F>>::Output> + + Sub<Self, Output=<Self as Euclidean<F>>::Output> + + for<'b> Add<&'b Self, Output=<Self as Euclidean<F>>::Output> + + for<'b> Sub<&'b Self, Output=<Self as Euclidean<F>>::Output> + + AddAssign<Self> + for<'b> AddAssign<&'b Self> + + SubAssign<Self> + for<'b> SubAssign<&'b Self> + + Neg<Output=<Self as Euclidean<F>>::Output> { + type Output : Euclidean<F>; + + /// Returns origin of same dimensions as `self`. + fn similar_origin(&self) -> <Self as Euclidean<F>>::Output; + + /// Calculate the square of the 2-norm, $\frac{1}{2}\\|x\\|_2^2$, where `self` is $x$. + #[inline] + fn norm2_squared(&self) -> F { + self.dot(self) + } + + /// Calculate the square of the 2-norm divided by 2, $\frac{1}{2}\\|x\\|_2^2$, + /// where `self` is $x$. + #[inline] + fn norm2_squared_div2(&self) -> F { + self.norm2_squared()/F::TWO + } + + /// Calculate the 2-norm $‖x‖_2$, where `self` is $x$. + #[inline] + fn norm2(&self) -> F { + self.norm2_squared().sqrt() + } + + /// Calculate the 2-distance squared $\\|x-y\\|_2^2$, where `self` is $x$. + fn dist2_squared(&self, y : &Self) -> F; + + /// Calculate the 2-distance $\\|x-y\\|_2$, where `self` is $x$. + #[inline] + fn dist2(&self, y : &Self) -> F { + self.dist2_squared(y).sqrt() + } + + /// Projection to the 2-ball. + #[inline] + fn proj_ball2(mut self, ρ : F) -> Self { + self.proj_ball2_mut(ρ); + self + } + + /// In-place projection to the 2-ball. + #[inline] + fn proj_ball2_mut(&mut self, ρ : F) { + let r = self.norm2(); + if r>ρ { + *self *= ρ/r + } + } +} + +/// Trait for [`Euclidean`] spaces with dimensions known at compile time. +pub trait StaticEuclidean<F : Float> : Euclidean<F> { + /// Returns the origin + fn origin() -> <Self as Euclidean<F>>::Output; +}
--- a/src/fe_model.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/fe_model.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,4 +1,6 @@ -/// Simple planar (1D and 2D) finite element discretisation on a box. +/*! +Simple planar (1D and 2D) finite element discretisation on a single simplex or box. +*/ pub mod base; pub mod p2_local_model;
--- a/src/fe_model/base.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/fe_model/base.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,6 +1,8 @@ -///! Base types for simple (finite) element models. +/*! +Base types for simple (finite) element models. +*/ -/// A local model can be evaluated for value and differential +/// A local function model that can be evaluated for value and differential. pub trait LocalModel<Domain, Codomain> { /// Get the value of the model at `x` fn value(&self, x : &Domain) -> Codomain;
--- a/src/fe_model/p2_local_model.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/fe_model/p2_local_model.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,16 +1,27 @@ +/*! +Second order polynomical (P2) models on real intervals and planar 2D simplices. +*/ + use crate::types::*; use crate::loc::Loc; use crate::sets::{Set,NPolygon,SpannedHalfspace}; use crate::linsolve::*; -use crate::norms::Dot; +use crate::euclidean::Dot; use super::base::{LocalModel,RealLocalModel}; use crate::sets::Cube; use numeric_literals::replace_float_literals; +/// Type for simplices of arbitrary dimension `N`. +/// +/// The type parameter `D` indicates the number of nodes. (Rust's const generics do not currently +/// allow its automatic calculation from `N`.) pub struct Simplex<F : Float, const N : usize, const D : usize>(pub [Loc<F, N>; D]); -pub type PlanarSimplex<F> = Simplex<F,2,3>; -pub type RealInterval<F> = Simplex<F,1,2>; +/// A two-dimensional planar simplex +pub type PlanarSimplex<F> = Simplex<F, 2, 3>; +/// A real interval +pub type RealInterval<F> = Simplex<F, 1, 2>; +/// Calculates (a+b)/2 #[inline] #[replace_float_literals(F::cast_from(literal))] pub(crate) fn midpoint<F : Float, const N : usize>(a : &Loc<F,N>, b : &Loc<F,N>) -> Loc<F, N> { @@ -94,12 +105,19 @@ } } +/// A trait for generating second order polynomial model of dimension `N` on `Self´. +/// +/// `Self` should present a subset aset of elements of the type [`Loc`]`<F, N>`. pub trait P2Model<F : Num, const N : usize> { + /// Implementation type of the second order polynomical model. + /// Typically a [`P2LocalModel`]. type Model : LocalModel<Loc<F,N>,F>; + /// Generates a second order polynomial model of the function `g` on `Self`. fn p2_model<G : Fn(&Loc<F, N>) -> F>(&self, g : G) -> Self::Model; } -pub struct P2LocalModel<F : Num, const N : usize, const E : usize, const V : usize, const Q : usize> { +/// A local second order polynomical model of dimension `N` with `E` edges +pub struct P2LocalModel<F : Num, const N : usize, const E : usize/*, const V : usize, const Q : usize*/> { a0 : F, a1 : Loc<F, N>, a2 : Loc<F, E>, @@ -120,7 +138,9 @@ } } -impl<F : Float> P2LocalModel<F, 1, 1, 2, 0> { +impl<F : Float> P2LocalModel<F, 1, 1/*, 2, 0*/> { + /// Creates a new 1D second order polynomical model based on three nodal coordinates and + /// corresponding function values. #[inline] pub fn new( &[n0, n1, n01] : &[Loc<F, 1>; 3], @@ -146,7 +166,7 @@ } impl<F : Float> P2Model<F,1> for RealInterval<F> { - type Model = P2LocalModel<F, 1, 1, 2, 0>; + type Model = P2LocalModel<F, 1, 1/*, 2, 0*/>; #[inline] fn p2_model<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> Self::Model { @@ -164,6 +184,7 @@ impl<F : Float> PlanarSimplex<F> { #[inline] + /// Returns the midpoints of all the edges of the simplex fn midpoints(&self) -> [Loc<F, 2>; 3] { let [ref n0, ref n1, ref n2] = &self.0; let n01 = midpoint(n0, n1); @@ -173,7 +194,9 @@ } } -impl<F : Float> P2LocalModel<F, 2, 3, 3, 3> { +impl<F : Float> P2LocalModel<F, 2, 3/*, 3, 3*/> { + /// Creates a new 2D second order polynomical model based on six nodal coordinates and + /// corresponding function values. #[inline] pub fn new( &[n0, n1, n2, n01, n12, n20] : &[Loc<F, 2>; 6], @@ -202,7 +225,7 @@ } impl<F : Float> P2Model<F,2> for PlanarSimplex<F> { - type Model = P2LocalModel<F, 2, 3, 3, 3>; + type Model = P2LocalModel<F, 2, 3/*, 3, 3*/>; #[inline] fn p2_model<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> Self::Model { @@ -217,7 +240,7 @@ macro_rules! impl_local_model { ($n:literal, $e:literal, $v:literal, $q:literal) => { - impl<F : Float> LocalModel<Loc<F, $n>, F> for P2LocalModel<F, $n, $e, $v, $q> { + impl<F : Float> LocalModel<Loc<F, $n>, F> for P2LocalModel<F, $n, $e/*, $v, $q*/> { #[inline] fn value(&self, x : &Loc<F,$n>) -> F { self.a0 + x.dot(&self.a1) + x.p2powers().dot(&self.a2) @@ -240,7 +263,8 @@ // #[replace_float_literals(F::cast_from(literal))] -impl<F : Float> P2LocalModel<F, 1, 1, 2, 0> { +impl<F : Float> P2LocalModel<F, 1, 1/*, 2, 0*/> { + /// Minimises the model along the edge `[x0, x1]`. #[inline] fn minimise_edge(&self, x0 : Loc<F, 1>, x1 : Loc<F,1>) -> (Loc<F,1>, F) { let &P2LocalModel{ @@ -271,7 +295,7 @@ } impl<'a, F : Float> RealLocalModel<RealInterval<F>,Loc<F,1>,F> -for P2LocalModel<F, 1, 1, 2, 0> { +for P2LocalModel<F, 1, 1/*, 2, 0*/> { #[inline] fn minimise(&self, &Simplex([x0, x1]) : &RealInterval<F>) -> (Loc<F,1>, F) { self.minimise_edge(x0, x1) @@ -279,8 +303,8 @@ } #[replace_float_literals(F::cast_from(literal))] -impl<F : Float> P2LocalModel<F, 2, 3, 3, 3> { - /// Minimise the 2D model on the edge {x0 + t(x1 - x0) | t ∈ [0, 1] } +impl<F : Float> P2LocalModel<F, 2, 3/*, 3, 3*/> { + /// Minimise the 2D model along the edge `[x0, x1] = {x0 + t(x1 - x0) | t ∈ [0, 1] }`. #[inline] fn minimise_edge(&self, x0 : &Loc<F,2>, x1 : &Loc<F,2>/*, v0 : F, v1 : F*/) -> (Loc<F,2>, F) { let &P2LocalModel { @@ -307,7 +331,7 @@ #[replace_float_literals(F::cast_from(literal))] impl<'a, F : Float> RealLocalModel<PlanarSimplex<F>,Loc<F,2>,F> -for P2LocalModel<F, 2, 3, 3, 3> { +for P2LocalModel<F, 2, 3/*, 3, 3*/> { #[inline] fn minimise(&self, el : &PlanarSimplex<F>) -> (Loc<F,2>, F) { let &P2LocalModel { @@ -347,7 +371,7 @@ #[replace_float_literals(F::cast_from(literal))] impl<'a, F : Float> RealLocalModel<Cube<F, 2>,Loc<F,2>,F> -for P2LocalModel<F, 2, 3, 3, 3> { +for P2LocalModel<F, 2, 3/*, 3, 3*/> { #[inline] fn minimise(&self, el : &Cube<F, 2>) -> (Loc<F,2>, F) { let &P2LocalModel { @@ -386,25 +410,6 @@ } } -// impl<F : Float> P2Model<F, 1> for Cube<F, 1> { -// type Model = CubeP2LocalModel<F,1,1,1>; - -// fn p2_model<G : FnMut(&Loc<F, 1>) -> F>(&self, fun : &G) -> Self::Model { -// CubeP2LocalModel([Simplex(self.corners()).p2_model(fun)]) -// } -// } - -// impl<F : Float> P2Model<F, 2> for Cube<F, 2> { -// type Model = CubeP2LocalModel<F,2,3,2>; - -// fn p2_model<G : FnMut(&Loc<F, 2>) -> F>(&self, fun : &G) -> Self::Model { -// let [a, b, c, d] = self.corners(); -// CubeP2LocalModel([Simplex([a, b, c]).p2_model(fun), -// Simplex([b, c, d]).p2_model(fun)]) -// } -// } - - #[cfg(test)] mod tests { use super::*;
--- a/src/iter.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/iter.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,5 +1,12 @@ -/// Iteration utilities; [stateful][StatefulIterator] and -/// [restartable][RestartableIterator] iterators. +/*! +Iteration utilities. + +This includes + * [Stateful][StatefulIterator] and [restartable][RestartableIterator] iterators. + * Variants of [`Iterator::map`] and [`Iterator::filter_map`], etc., that work on functions + instead of closures, forgoing exact closure type signature specification; see [`Mappable`]. + +*/ use crate::types::Integer; use std::marker::PhantomData; @@ -153,138 +160,6 @@ } } -/* -/// This is a helper struct returning iterators with a lifetime from trait functions. -/// `T` will typically be `&'a Self` or some other reference that can only be -/// specified in the trait function declaration without having to parametrise the entire -/// trait with `'a`. The parameter `ImplData` is supposed to be a unique type that idenfiers -/// the iterator for the implementation of the trait. Then `Iterator` is to be implemented -/// for `ProxyIterator<T, ImplState>`. For example: -/// TODO: not updated. -/// ``` -/// trait Test { -/// type ImplData; -/// fn returns_iterator<'a>(&'a self) -> ProxyIterator<&'a Self, Self::ImplData>; -/// } -/// struct MyIndex{ -/// index : usize, -/// } -/// impl<const N : usize> Test for [f64; N] { -/// type ImplData = MyIndex; -/// fn returns_iterator<'a>(&'a self) -> ProxyIterator<&'a Self, MyIndex> { -/// TraitIterator::new(self, MyIndex{ index : 0 }) -/// } -/// } -/// impl<'a, const N : usize> Iterator for `TraitIterator<&'a [f64; N], MyIndex> { -/// type Item = f64; -/// fn next(&mut self) -> Option<f64> { -/// if self.state.index < N { -/// let v = self.data[self.state.index]; -/// self.state.index += 1; -/// Some(v) -/// } else { -/// None -/// } -/// } -/// } -/// ``` -*/ - -macro_rules! impl_trait_iterator { - ($name:ident, $proxyname:ident, &'a $($mut:ident)?, $($itemref:tt)*) => { - pub trait $name { - type Item : 'static; - type Data : ?Sized; - fn next_with<'a>(&mut self, data : &'a $($mut)? Self::Data) -> Option<$($itemref)* Self::Item>; - } - - pub struct $proxyname<'a, Impl> - where Impl : $name { - pub data : &'a $($mut)? Impl::Data, - pub iter : Impl, - } - - impl<'a, Impl> $proxyname<'a, Impl> - where Impl : $name { - #[inline] - pub fn new(data : &'a $($mut)? Impl::Data, iter : Impl) -> $proxyname<'a, Impl> { - $proxyname{ data : data, iter : iter} - } - } - - impl<'a, Impl> Iterator for $proxyname<'a, Impl> - where Impl : $name { - type Item = $($itemref)* Impl::Item; - #[inline] - fn next(&mut self) -> Option<Self::Item> { - self.iter.next_with(self.data) - } - } - } -} - -impl_trait_iterator!(TraitIterator, ProxyIterator, &'a, ); -impl_trait_iterator!(TraitIteratorRef, ProxyIteratorRef, &'a, &'a); -// TODO: Lifetime errors. Probably due to potential overlapping mutable borrows. -// So should use some unsafe pointer hacks or so. -//impl_trait_iterator!(TraitIteratorMut, ProxyIteratorMut, &'a mut, &'a mut); - -/* -pub struct FlattenTI<SI> -where SI : TraitIterator, - SI::Item : Iterator { - iter_iter : SI, - element_iter : Option<SI::Item>, -} - -impl<SI> FlattenTI<SI> -where SI : TraitIterator, - SI::Item : Iterator { - fn new(iter : SI) -> Self { - FlattenTI{ iter_iter : iter, element_iter : None } - } -} - -impl<SI : TraitIterator> TraitIterator for FlattenTI<SI> -where SI::Item : Iterator { - type Item = <SI::Item as Iterator>::Item; - type Data = SI::Data; - - fn next_with<'a>(&mut self, data : &'a $($mut)? SI::Data) -> Option<Self::Item> { - loop { - if let Some(ei) = self.element_iter { - if let n @ Some(_) = ei.next() { - return n - } - } else { - self.element_iter = self.iter_iter.next_with(data); - if let None = self.element_iter { - return None - } - } - } - } -} -*/ - -/* -pub enum EitherIterator<T, L : Iterator<Item=T>, R : Iterator<Item=T>> { - Left(L), - Right(R) -} - -impl<T, L : Iterator<Item=T>, R : Iterator<Item=T>> Iterator for EitherIterator<T, L, R> { - type Item = T; - #[inline] - fn next(&mut self) -> Option<T> { - match self { - EitherIterator::Left(ref l) => { l.next() } - EitherIterator::Right(ref r) => { r.next() } - } - } -} -*/ - /// A [`RestartableIterator`] over the range ´[start, end)`. #[derive(Clone,Copy,Debug)] pub struct RangeIter<T> { @@ -369,8 +244,8 @@ } } - -pub struct RestartableIter<'a, T> { +/// A restartable slice iterator. +pub struct RestartableSliceIter<'a, T> { inner : RangeIter<*const T>, _phantom : PhantomData<&'a[T]> } @@ -385,7 +260,8 @@ } } -impl<'a, T : 'a> RestartableIter<'a, T> { +impl<'a, T : 'a> RestartableSliceIter<'a, T> { + /// Converts `Some` pointer to `Some` reference fn map_result(result : Option<*const T>) -> Option<&'a T> { match result { None => { None } @@ -393,14 +269,15 @@ } } + /// Creates a restartable iterator over `slice`. pub fn new(slice : &'a [T]) -> Self { let ptr_range = slice.as_ptr_range(); let inner = RangeIter::new(ptr_range.start, ptr_range.end); - RestartableIter{ inner : inner, _phantom : PhantomData } + RestartableSliceIter{ inner : inner, _phantom : PhantomData } } } -impl<'a, T : 'a> Iterator for RestartableIter<'a, T> { +impl<'a, T : 'a> Iterator for RestartableSliceIter<'a, T> { type Item = &'a T; #[inline] fn next(&mut self) -> Option<&'a T> { @@ -408,24 +285,25 @@ } } -impl<'a, T : 'a> ExactSizeIterator for RestartableIter<'a, T> { +impl<'a, T : 'a> ExactSizeIterator for RestartableSliceIter<'a, T> { #[inline] fn len(&self) -> usize { self.inner.len() } } -impl<'a, T : 'a> StatefulIterator for RestartableIter<'a, T> { +impl<'a, T : 'a> StatefulIterator for RestartableSliceIter<'a, T> { #[inline] fn current(&self) -> Option<Self::Item> { Self::map_result(self.inner.current()) } } -impl<'a, T : 'a> RestartableIterator for RestartableIter<'a, T> { +impl<'a, T : 'a> RestartableIterator for RestartableSliceIter<'a, T> { fn restart(&mut self) -> Option<Self::Item> { Self::map_result(self.inner.current()) } } +/// A restartable variant of [`std::iter::Cloned`]. pub struct RestartableCloned<I : Iterator> { inner : I, }
--- a/src/iterate.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/iterate.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,22 +1,23 @@ /*! +Tools for separating the computational steps of an iterative algorithm from stopping rules +and reporting. -This module provides tools for separating the computational steps of iterative -algorithms from visualisation and other reporting routines. The approach is heavily -indebted to functional programming. The computational step is to be implemented as a -closure. That closure gets passed a `state` parameter that implements an -[`if_verbose`][AlgIteratorState::if_verbose] method that is to be called to determine -whether function values or other potentially costly things need to be calculated on that -iteration. The parameter of [`if_verbose`][AlgIteratorState::if_verbose] is another -closure that does the necessary computation. +The computational step is to be implemented as a closure. That closure gets passed a `state` +parameter that implements an [`if_verbose`][AlgIteratorState::if_verbose] method that is to +be called to determine whether function values or other potentially costly things need to be +calculated on that iteration. The parameter of [`if_verbose`][AlgIteratorState::if_verbose] is +another closure that does the necessary computation. ## Simple example ```rust -use alg_tools::types::*; -use alg_tools::iterate::*; -let mut iter = AlgIteratorOptions{ max_iter : 100, - verbose_iter : Verbose::Every(10), - .. Default::default() }; +# use alg_tools::types::*; +# use alg_tools::iterate::*; +let mut iter = AlgIteratorOptions{ + max_iter : 100, + verbose_iter : Verbose::Every(10), + .. Default::default() +}; let mut x = 1 as float; iter.iterate(|state|{ // This is our computational step @@ -44,11 +45,8 @@ ``` */ -use colored::{Colorize,ColoredString}; -//use num_traits::Num; +use colored::{Colorize, ColoredString}; use core::fmt::Debug; -//use core::fmt::Display; -//use std::sync::mpsc::Receiver; use serde::{Serialize, Deserialize}; use cpu_time::ProcessTime; use std::marker::PhantomData; @@ -57,7 +55,7 @@ use crate::types::*; use crate::logger::*; -/// Create the displayed presentation for log items. Override for nice (condensed) presentation of rich log items, or define `show` if it is not defined for your type. +/// Create the displayed presentation for log items. pub trait LogRepr : Debug { fn logrepr(&self) -> ColoredString { format!("« {self:?} »").as_str().into() } } @@ -83,6 +81,9 @@ } } +/// Helper struct for returning results annotated with an additional string to +/// [`if_verbose`][AlgIteratorState::if_verbose]. The [`LogRepr`] implementation will +/// display that string when so decided by the specific [`AlgIterator`] in use. #[derive(Debug,Clone)] pub struct Annotated<F>(pub F, pub String); @@ -102,11 +103,17 @@ pub data : V } -fn make_logitem<V>(iter : usize, data : V) -> LogItem<V> { - LogItem{ iter, data } +impl<V> LogItem<V> { + /// Creates a new log item + fn new(iter : usize, data : V) -> Self { + LogItem{ iter, data } + } } -/// State of an [`AlgIterator`] +/// State of an [`AlgIterator`]. +/// +/// This is the parameter obtained by the closure passed to [`AlgIterator::iterate`] or +/// [`AlgIteratorFactory::iterate`]. pub trait AlgIteratorState { /// Call `call_objective` if this is a verbose iteration. /// @@ -117,6 +124,7 @@ /// [module documentation][self]. fn if_verbose<V, E : Error>(&self, calc_objective : impl FnMut() -> V) -> Step<V, E>; + /// Returns the current iteration count. fn iteration(&self) -> usize; } @@ -145,13 +153,16 @@ } } -/// An iterator for algorithms, produced by [`AlgIteratorFactory.prepare`]. +/// An iterator for algorithms, produced by [`AlgIteratorFactory::prepare`]. /// /// Typically not accessed directly, but transparently produced by an [`AlgIteratorFactory`]. /// Every [`AlgIteratorFactory`] has to implement a corresponding `AlgIterator`. pub trait AlgIterator : Sized { + /// The state type type State : AlgIteratorState; + /// The output type for [`Self::step`]. type Output; + /// The error type for [`Self::step`] and [`Self::iterate`]. type Err : Error; /// Advance the iterator. @@ -162,7 +173,7 @@ self.state().iteration() } - /// Return current iteration stats. + /// Return current state. fn state(&self) -> Self::State; /// Iterate the `AlgIterator` until termination. @@ -181,9 +192,7 @@ /// Converts the `AlgIterator` into a plain [`Iterator`]. /// - /// [`Step::Quiet`] results are discarded. - /// (This is to avoid having to manually implement [`IntoIterator`] or [`Iterator`] - /// for every single [`AlgIterator`].) + /// [`Step::Quiet`] results are discarded, and [`Step::Failure`] results **panic**. fn downcast(self) -> AlgIteratorI<Self> { AlgIteratorI(self) } @@ -191,7 +200,7 @@ /// Conversion of an `AlgIterator` into a plain [`Iterator`]. /// -/// The conversion discards [`Step::Quiet`] and panics on [`Step::Failure`]. +/// The conversion discards [`Step::Quiet`] and **panics** on [`Step::Failure`]. pub struct AlgIteratorI<A>(A); impl<A> Iterator for AlgIteratorI<A> @@ -210,24 +219,19 @@ } } -/// A factory for producing an [`AlgIterator`]. To use one: +/// A factory for producing an [`AlgIterator`]. /// -/// ```rust -/// use alg_tools::iterate::*; -/// let mut iter = AlgIteratorOptions::default(); -/// iter.iterate(|state|{ -/// // perform iterations -/// state.if_verbose(||{ -/// // calculate and return function value or other displayed data v -/// return 0 -/// }) -/// }) -/// ``` +/// For usage instructions see the [module documentation][self]. pub trait AlgIteratorFactory<V> : Sized { type Iter<F, E> : AlgIterator<State = Self::State, Output = Self::Output, Err = E> where F : FnMut(&Self::State) -> Step<V, E>, E : Error; + /// The state type of the corresponding [`AlgIterator`]. + /// A reference to this is passed to the closures passed to methods such as [`Self::iterate`]. type State : AlgIteratorState; + /// The output type of the corresponding [`AlgIterator`]. + /// This is the output of the closures passed to methods such as [`Self::iterate`] after + /// mappings performed by each [`AlgIterator`] implementation. type Output; /// Prepare an [`AlgIterator`], consuming the factory. @@ -240,10 +244,10 @@ /// Iterate the the closure `step`. /// - /// The closure should accept one `state` ([`AlgIteratorState`] parameter, - /// and return the output of `state.if_verbose`, [`Step::Terminated`] to indicate completion - /// for other reason, or [`Step::Failure`] for termination for failure. - /// For usage instructions see the [module documentation][self]. + /// The closure should accept a `state` parameter (satisfying the trait [`AlgIteratorState`]). + /// It should return the output of + /// `state.`[`if_verbose`][AlgIteratorState::if_verbose], [`Step::Terminated`] to indicate + /// completion for other reason, or [`Step::Failure`] for termination for failure. /// /// This method is equivalent to [`Self::prepare`] followed by [`AlgIterator::iterate`]. #[inline] @@ -255,9 +259,9 @@ /// Iterate the the closure `step`. /// - /// The closure should accept one `state` ([`AlgIteratorState`] parameter, - /// and return the output of `state.if_verbose`. For a fallible closure, - /// use [`Self::iterate_fallible`]. + /// The closure should accept a `state` parameter (satisfying the trait [`AlgIteratorState`]), + /// It should return the output of + /// `state.`[`if_verbose`][AlgIteratorState::if_verbose]. /// /// For usage instructions see the [module documentation][self]. /// @@ -271,9 +275,9 @@ /// Iterate the closure `step` with data produced by `datasource`. /// - /// The closure should accept a `state` ([`AlgIteratorState`] parameter, and a data parameter - /// taken from `datasource` It should return the output of - /// `state.[if_verbose][AlgIteratorState.if_verbose]`, [`Step::Terminated`] to indicate + /// The closure should accept a `state` parameter (satisfying the trait [`AlgIteratorState`]), + /// and a data parameter taken from `datasource`. It should return the output of + /// `state.`[`if_verbose`][AlgIteratorState::if_verbose], [`Step::Terminated`] to indicate /// completion for other reason, or [`Step::Failure`] for termination for failure. /// /// If the `datasource` runs out of data, the iterator is considered having terminated @@ -292,9 +296,9 @@ /// Iterate the closure `step` with data produced by `datasource`. /// - /// The closure should accept a `state` ([`AlgIteratorState`] parameter, and a data parameter - /// taken from `datasource` It should return the output of - /// `state.[if_verbose][AlgIteratorState.if_verbose]`. + /// The closure should accept a `state` parameter (satisfying the trait [`AlgIteratorState`]), + /// and a data parameter taken from `datasource`. It should return the output of + /// `state.`[`if_verbose`][AlgIteratorState::if_verbose]. /// /// If the `datasource` runs out of data, the iterator is considered having terminated /// successsfully. @@ -319,6 +323,24 @@ // } /// Add logging to the iterator produced by the factory. + /// + /// Returns a new factory whose corresponding [`AlgIterator`] only inserts that data into the + /// log without passing it onwards. + /// + /// Use as: + /// ```rust + /// # use alg_tools::iterate::*; + /// # use alg_tools::logger::*; + /// let iter = AlgIteratorOptions::default(); + /// let mut log = Logger::new(); + /// iter.into_log(&mut log).iterate(|state|{ + /// // perform iterations + /// state.if_verbose(||{ + /// // calculate and return function value or other displayed data v + /// return 0 + /// }) + /// }) + /// ``` fn into_log<'log>(self, logger : &'log mut Logger<Self::Output>) -> LoggingIteratorFactory<'log, Self::Output, Self> where Self : Sized { @@ -329,6 +351,8 @@ } /// Map the output of the iterator produced by the factory. + /// + /// Returns a new factory. fn mapped<U, G>(self, map : G) -> MappingIteratorFactory<G, Self> where Self : Sized, @@ -339,11 +363,14 @@ } } - /// Adds iteration number to the output. Typically followed by [`Self::into_log`]. + /// Adds iteration number to the output. + /// + /// Returns a new factory. + /// Typically followed by [`Self::into_log`]. fn with_iteration_number(self) -> MappingIteratorFactory<fn(usize, Self::Output) -> LogItem<Self::Output>, Self> where Self : Sized { - self.mapped(make_logitem) + self.mapped(LogItem::new) } /// Add timing to the iterator produced by the factory. @@ -370,7 +397,9 @@ fn is_quiet(&self) -> bool { false } } -/// Options for [`BasicAlgIteratorFactory`]. Use as: +/// Options for [`BasicAlgIteratorFactory`]. +/// +/// Use as: /// ``` /// # use alg_tools::iterate::*; /// let iter = AlgIteratorOptions{ max_iter : 10000, .. Default::default() }; @@ -431,8 +460,12 @@ /// State of a `BasicAlgIterator` #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] pub struct BasicState { + /// Current iteration iter : usize, + /// Whether the iteration is verbose, i.e., results should be displayed. + /// Requires `calc` to be `true`. verbose : bool, + /// Whether results should be calculated. calc : bool, } @@ -443,7 +476,7 @@ _phantoms : PhantomData<V>, } -/// The simplest [`AlgIterator`], reated by [`BasicAlgIteratorFactory`] +/// The simplest [`AlgIterator`], created by [`BasicAlgIteratorFactory`] #[derive(Clone,Debug)] pub struct BasicAlgIterator<F, V, E : Error> { options : AlgIteratorOptions, @@ -580,11 +613,13 @@ // Stall detecting iteration function. // -/// An [`AlgIteratorFactory`] for an [`AlgIterator`] that detects “stall” -/// $(v_{k+n}-v_k)/v_k ≤ θ$, where $n$ the distance between [`Step::Result`] iterations. +/// An [`AlgIteratorFactory`] for an [`AlgIterator`] that detects “stall”. +/// +/// We define stall as $(v_{k+n}-v_k)/v_k ≤ θ$, where $n$ the distance between +/// [`Step::Result`] iterations, and $θ$ is the provided `stall` parameter. #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] pub struct StallIteratorFactory<U : Num, BaseFactory> { - /// Basic options + /// An [`AlgIteratorFactory`] on which to build on pub base_options : BaseFactory, /// Stalling threshold $θ$. pub stall : U, @@ -660,9 +695,9 @@ /// return value is less than `target`, and terminates if it is. #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] pub struct ValueIteratorFactory<U : Num, BaseFactory> { - /// Basic options + /// An [`AlgIteratorFactory`] on which to build on pub base_options : BaseFactory, - /// Stalling threshold $θ$. + /// Target value pub target : U, } @@ -734,26 +769,15 @@ // /// [`AlgIteratorFactory`] for a logging [`AlgIterator`]. +/// +/// Typically produced with [`AlgIteratorFactory::into_log`]. /// The `Output` of the corresponding [`LoggingIterator`] is `()`: -/// The data is only inserted into the log, and not passed onwards. -/// Use as: -/// ```rust -/// use alg_tools::iterate::*; -/// use alg_tools::logger::*; -/// let iter = AlgIteratorOptions::default(); -/// let mut log = Logger::new(); -/// iter.into_log(&mut log).iterate(|state|{ -/// // perform iterations -/// state.if_verbose(||{ -/// // calculate and return function value or other displayed data v -/// return 0 -/// }) -/// }) -/// ``` #[derive(Debug)] pub struct LoggingIteratorFactory<'log, U, BaseFactory> { - pub base_options : BaseFactory, - pub logger : &'log mut Logger<U>, + /// Base [`AlgIteratorFactory`] on which to build + base_options : BaseFactory, + /// The `Logger` to use. + logger : &'log mut Logger<U>, } /// Iterator produced by `LoggingIteratorFactory`. @@ -821,13 +845,18 @@ } /// This [`AlgIteratorFactory`] allows output mapping. +/// +/// Typically produced with [`AlgIteratorFactory::mapped`]. #[derive(Debug)] pub struct MappingIteratorFactory<G, BaseFactory> { - pub base_options : BaseFactory, - pub map : G, + /// Base [`AlgIteratorFactory`] on which to build + base_options : BaseFactory, + /// A closure `G : Fn(usize, BaseFactory::Output) -> U` that gets the current iteration + /// and the output of the base factory as input, and produces a new output. + map : G, } -/// Iterator produced by `MappingIteratorFactory`. +/// [`AlgIterator`] produced by [`MappingIteratorFactory`]. pub struct MappingIterator<G, BaseIterator> { base_iterator : BaseIterator, map : G, @@ -892,7 +921,7 @@ // Timing iterator // -/// An [`AlgIteratorFactory`] for an [`AlgIterator`] that adds CPU time spent to verbose events. +/// An [`AlgIteratorFactory`] for an [`AlgIterator`] that adds spent CPU time to verbose events. #[derive(Debug)] pub struct TimingIteratorFactory<BaseFactory>(pub BaseFactory);
--- a/src/lib.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/lib.rs Sat Oct 22 14:40:28 2022 +0300 @@ -28,6 +28,7 @@ pub mod error; pub mod maputil; pub mod tuple; +pub mod euclidean; pub mod norms; #[macro_use] pub mod loc;
--- a/src/lingrid.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/lingrid.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,27 +1,42 @@ -/// Production of equally spaced nodes within intervals. +/*! +Linear grids. + +These are multi-dimensional intervals $\prod_{i=1}^N [a_i, b_i]$ divided along each dimension +into n_i equally-spaced nodes, with $a_i$ the first node and $b_i$ the last node along each +dimension. + +The [`LinSpace`]s provided by this module are similar to [`num::range_step_inclusive`], but as an +iterator they are [restartable][RestartableIterator] and parametrised by the number of nodes +instead of a step. This way it can be ensured that $a_i$ and $b_i$ are the last and the first node. + +The starting points for the use of this module are the [`linspace`], [`lingrid`], and +[`lingrid_centered`] functions. They return a [`LinSpace`]s that implements [`IntoIterator`] for +iteration over the grid. Additional utility functions are in the [`Grid`] trait. +*/ use crate::types::*; use crate::loc::Loc; use crate::sets::Cube; use crate::iter::{RestartableIterator, StatefulIterator}; use crate::maputil::{map2, map4}; -use serde::Serialize; +use serde::{Serialize, Deserialize}; // TODO: rewrite this using crate::sets::Cube. -/// The interval `[start, end]` divided into `count` nodes. -/// Implementes `IntoIterator` to iterate over the nodes of the `LinSpace`. -/// It is similar to [`num::range_step_inclusive`], but [restartable][RestartableIterator] and -/// parametrised by the number of nodes instead of a step. This way it can be ensured that the -/// last item produced is equal to `end`. To use other types than `float` and `f32`, create -/// the `LinSpace` struct directly. -#[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] -pub struct LinSpace<F, I> { - pub start : F, - pub end : F, +/// An abstraction of possibly multi-dimensional linear grids. +/// +/// `U` is typically a `F` for a `Float` `F` for one-dimensional grids created by `linspace`, +/// or [`Loc`]`<F, N>` for multi-dimensional grids created by `lingrid`. +/// In the first case `count` of nodes is `usize`, and in the second case `[usize; N]`. +#[derive(Clone, Copy, Debug, Serialize, Deserialize, Eq, PartialEq)] +pub struct LinSpace<U, I> { + pub start : U, + pub end : U, pub count : I, } +/// A `N`-dimensional interval divided into an indicated number of equally-spaced nodes along +/// each dimension. #[allow(type_alias_bounds)] // Need it to access F::CompatibleSize. pub type LinGrid<F : Float, const N : usize> = LinSpace<Loc<F, N>, [usize; N]>; @@ -30,7 +45,7 @@ LinSpace{ start : start, end : end, count : count } } -/// Create a multi-dimensional linear grid. +/// Creates a multi-dimensional linear grid. /// /// The first and last point in each dimension are the boundaries of the corresponding /// dimensions of `cube`, and there are `count` nodes along each dimension. @@ -65,13 +80,19 @@ current : Option<I>, } -pub trait Grid<F, I> { - fn entry_linear_unchecked(&self, i : usize) -> F; - //fn entry(&self, i : I) -> Option<F>; - fn entry_unchecked(&self, i : &I) -> F; +/// Abstraction of a linear grid over space `U` with multi-dimensional index set `I`. +pub trait Grid<U, I> { + /// Converts a linear index `i` into a grid point. + fn entry_linear_unchecked(&self, i : usize) -> U; + // Converts a multi-dimensional index `i` into a grid point. + fn entry_unchecked(&self, i : &I) -> U; + + // fn entry(&self, i : I) -> Option<F> } +/// Helper trait for iteration of [`Grid`]s. pub trait GridIteration<F, I> { + /// Returns the next multi-dimensional index (not yet converted into grid point). fn next_index(&mut self) -> Option<I>; } @@ -102,9 +123,12 @@ #[inline] fn next_index(&mut self) -> Option<I> { match self.current { - None if I::ZERO < self.lingrid.count => { self.current = Some(I::ZERO); self.current } - Some(v) if v+I::ONE < self.lingrid.count => { self.current = Some(v+I::ONE); self.current } - _ => { None } + None if I::ZERO < self.lingrid.count + => { self.current = Some(I::ZERO); self.current } + Some(v) if v+I::ONE < self.lingrid.count + => { self.current = Some(v+I::ONE); self.current } + _ + => { None } } } }
--- a/src/linops.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/linops.rs Sat Oct 22 14:40:28 2022 +0300 @@ -26,7 +26,7 @@ self.axpy(1.0, x, 0.0) } - /// Computes `y = αx`, where `y` is `Self`. + /// Computes `y = αx`, where `y` is `Self`. fn scale_from(&mut self, α : F, x : &X) { self.axpy(α, x, 0.0) } @@ -35,15 +35,15 @@ /// Efficient in-place application for [`Linear`] operators. #[replace_float_literals(F::cast_from(literal))] pub trait GEMV<F : Num, X, Y = <Self as Linear<X>>::Codomain> : Linear<X> { - // Computes `y = αAx + βy`, where `A` is `Self`. + /// Computes `y = αAx + βy`, where `A` is `Self`. fn gemv(&self, y : &mut Y, α : F, x : &X, β : F); - // Computes `y = Ax`, where `A` is `Self` + /// Computes `y = Ax`, where `A` is `Self` fn apply_mut(&self, y : &mut Y, x : &X){ self.gemv(y, 1.0, x, 0.0) } - // Computes `y += Ax`, where `A` is `Self` + /// Computes `y += Ax`, where `A` is `Self` fn apply_add(&self, y : &mut Y, x : &X){ self.gemv(y, 1.0, x, 1.0) } @@ -59,9 +59,9 @@ fn opnorm_bound(&self) -> Self::FloatType; } -/// Linear operator application into mutable target. The [`AsRef`] bound -/// is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; -/// the former is assumed to be e.g. a view into the latter. +// Linear operator application into mutable target. The [`AsRef`] bound +// is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; +// the former is assumed to be e.g. a view into the latter. /*impl<X,Y,T> Fn(&X) -> Y for T where T : Linear<X,Codomain=Y> { fn call(&self, x : &X) -> Y { @@ -69,7 +69,7 @@ } }*/ -/// Trait for forming the adjoint operator of an operator $A$=`Self`. +/// Trait for forming the adjoint operator of `Self`. pub trait Adjointable<X,Yʹ> : Linear<X> { type AdjointCodomain; type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a; @@ -82,7 +82,9 @@ }*/ } -/// Trait for forming a preadjoint of an operator $A$, i.e., an operator $A_*$ +/// Trait for forming a preadjoint of an operator. +/// +/// For an operator $A$ this is an operator $A_*$ /// such that its adjoint $(A_*)^*=A$. The space `X` is the domain of the `Self` /// operator. The space `Ypre` is the predual of its codomain, and should be the /// domain of the adjointed operator. `Self::Preadjoint` should be @@ -95,7 +97,7 @@ fn preadjoint(&self) -> Self::Preadjoint<'_>; } -/// Adjointable operators $A: X → Y$ on between reflexibe spaces $X$ and $Y$. +/// Adjointable operators $A: X → Y$ on between reflexive spaces $X$ and $Y$. pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {} impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
--- a/src/linsolve.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/linsolve.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,4 +1,6 @@ -/// Linear solvers for small problems. +/*! +Linear equation solvers for small problems stored in Rust arrays. +*/ use crate::types::Float; use std::mem::MaybeUninit;
--- a/src/loc.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/loc.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,18 +1,26 @@ /*! -Array containers with vectorspace operations on floats. -For working with vectors in $ℝ^2$ or $ℝ^3$. +Array containers that support vector space operations on floats. +For working with small vectors in $ℝ^2$ or $ℝ^3$. */ use std::ops::{Add,Sub,AddAssign,SubAssign,Mul,Div,MulAssign,DivAssign,Neg,Index,IndexMut}; use std::slice::{Iter,IterMut}; use crate::types::{Float,Num,SignedNum}; use crate::maputil::{FixedLength,FixedLengthMut,map1,map2,map1_mut,map2_mut}; +use crate::euclidean::*; use crate::norms::*; use crate::linops::AXPY; use serde::ser::{Serialize, Serializer, SerializeSeq}; +/// A container type for (short) `N`-dimensional vectors of element type `F`. +/// +/// Supports basic operations of an [`Euclidean`] space, several [`Norm`]s, and +/// fused [`AXPY`] operations, among others. #[derive(Copy,Clone,Debug,PartialEq,Eq)] -pub struct Loc<F, const N : usize>(pub [F; N]); +pub struct Loc<F, const N : usize>( + /// An array of the elements of the vector + pub [F; N] +); // Need to manually implement as [F; N] serialisation is provided only for some N. impl<F, const N : usize> Serialize for Loc<F, N> @@ -32,16 +40,19 @@ } impl<F, const N : usize> Loc<F, N> { + /// Creates a new `Loc` vector from an array. #[inline] pub fn new(arr : [F; N]) -> Self { Loc(arr) } - + + /// Returns an iterator over the elements of the vector #[inline] pub fn iter(&self) -> Iter<'_, F> { self.0.iter() } + /// Returns an iterator over mutable references to the elements of the vector #[inline] pub fn iter_mut(&mut self) -> IterMut<'_, F> { self.0.iter_mut() @@ -49,46 +60,59 @@ } impl<F : Copy, const N : usize> Loc<F, N> { + /// Maps `g` over the elements of the vector, returning a new [`Loc`] vector #[inline] - pub fn map<H, G : Fn(F) -> H>(&self, g : G) -> Loc<H, N> { + pub fn map<H>(&self, g : impl Fn(F) -> H) -> Loc<H, N> { Loc::new(map1(self, |u| g(*u))) } + /// Maps `g` over pairs of elements of two vectors, retuning a new one. #[inline] - pub fn map2<H, G : Fn(F, F) -> H>(&self, other : &Self, g : G) -> Loc<H, N> { + pub fn map2<H>(&self, other : &Self, g : impl Fn(F, F) -> H) -> Loc<H, N> { Loc::new(map2(self, other, |u, v| g(*u, *v))) } + /// Maps `g` over mutable references to elements of the vector. #[inline] - pub fn map_mut<G : Fn(&mut F)>(&mut self, g : G) { + pub fn map_mut(&mut self, g : impl Fn(&mut F)) { map1_mut(self, g) } + /// Maps `g` over pairs of mutable references to elements of `self, and elements + /// of `other` vector. #[inline] - pub fn map2_mut<G : Fn(&mut F, F)>(&mut self, other : &Self, g : G) { + pub fn map2_mut(&mut self, other : &Self, g : impl Fn(&mut F, F)) { map2_mut(self, other, |u, v| g(u, *v)) } + /// Maps `g` over the elements of `self` and returns the product of the results. #[inline] - pub fn product_map<G : Fn(F) -> A, A : Num>(&self, f : G) -> A { + pub fn product_map<A : Num>(&self, g : impl Fn(F) -> A) -> A { match N { - 1 => f(unsafe { *self.0.get_unchecked(0) }), - 2 => f(unsafe { *self.0.get_unchecked(0) }) * - f(unsafe { *self.0.get_unchecked(1) }), - 3 => f(unsafe { *self.0.get_unchecked(0) }) * - f(unsafe { *self.0.get_unchecked(1) }) * - f(unsafe { *self.0.get_unchecked(2) }), - _ => self.iter().fold(A::ONE, |m, &x| m*f(x)) + 1 => g(unsafe { *self.0.get_unchecked(0) }), + 2 => g(unsafe { *self.0.get_unchecked(0) }) * + g(unsafe { *self.0.get_unchecked(1) }), + 3 => g(unsafe { *self.0.get_unchecked(0) }) * + g(unsafe { *self.0.get_unchecked(1) }) * + g(unsafe { *self.0.get_unchecked(2) }), + _ => self.iter().fold(A::ONE, |m, &x| m * g(x)) } } } +/// Construct a [`Loc`]. +/// +/// Use as +/// ``` +/// # use alg_tools::loc::Loc; +/// # use alg_tools::loc; +/// let x = loc![1.0, 2.0]; +/// ``` #[macro_export] macro_rules! loc { ($($x:expr),+ $(,)?) => { Loc::new([$($x),+]) } } -// Conversions impl<F, const N : usize> From<[F; N]> for Loc<F, N> { #[inline]
--- a/src/logger.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/logger.rs Sat Oct 22 14:40:28 2022 +0300 @@ -19,12 +19,12 @@ Logger{ data : Vec::new() } } - /// Store the value `v` in the log at index `i`. + /// Store the value `v` in the log. pub fn log(&mut self, v : V) -> () { self.data.push(v); } - /// Get logged data. + /// Get the logged data as a [`Vec`] array. pub fn data(&self) -> &Vec<V> { &self.data }
--- a/src/mapping.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/mapping.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,43 +1,60 @@ -/// Traits for (mathematical) functions. +/*! +Traits for mathematical functions. +*/ use std::marker::PhantomData; use crate::types::{Float}; use serde::Serialize; use crate::loc::Loc; +/// A mapping from `Domain` to `Codomain`. pub trait Mapping<Domain> { type Codomain; + /// Calculate the value of the mapping at `x`. fn value(&self, x : Domain) -> Self::Codomain; } +/// A helper trait alias for referring to `Mapping`s from references to floats. pub trait RealRefMapping<F : Float, const N : usize> : for<'a> Mapping<&'a Loc<F, N>, Codomain=F> {} impl<F : Float, T, const N : usize> RealRefMapping<F, N> for T where T : for<'a> Mapping<&'a Loc<F, N>, Codomain=F> {} + +/// A differentiable mapping from `Domain` to [`Mapping::Codomain`], with differentials +/// `Differential`. pub trait DifferentiableMapping<Domain> : Mapping<Domain> { type Differential; + /// Calculate the differentialeof the mapping at `x`. fn differential(&self, x : Domain) -> Self::Differential; } +/// A `Mapping` whose minimum and maximum can be computed. pub trait RealMapping<Domain> : Mapping<Domain> where Self::Codomain : Float { + /// Calculate a minimum and a minimiser of the mapping. fn minimise(&self, tolerance : Self::Codomain) -> (Domain, Self::Codomain); + /// Calculate a maximum and a maximiser of the mapping. fn maximise(&self, tolerance : Self::Codomain) -> (Domain, Self::Codomain); } -// -// Sum -// - +/// A sum of [`Mapping`]s. #[derive(Serialize, Debug, Clone)] pub struct Sum<Domain, M : Mapping<Domain>> { components : Vec<M>, _domain : PhantomData<Domain>, } +impl<Domain, M : Mapping<Domain>> Sum<Domain, M> { + /// Construct from an iterator. + pub fn new<I : Iterator<Item = M>>(iter : I) -> Self { + Sum { components : iter.collect(), _domain : PhantomData } + } +} + + impl<Domain, M> Mapping<Domain> for Sum<Domain,M> where M : Mapping<Domain>, M :: Codomain : std::iter::Sum,
--- a/src/maputil.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/maputil.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,18 +1,30 @@ -/// -/// Mapping utilities -/// +/*! +Utilities for mapping over various container types. +*/ use std::mem::MaybeUninit; use itertools::izip; +/// Trait for a fixed-length container type. +/// +/// Implemented by [`Loc`][crate::loc::Loc] vectors, [`Cube`][crate::sets::Cube]s, +/// and basic arrays. pub trait FixedLength<const N : usize> { + /// Type of elements of the container. type Elem; - type Iter : Iterator<Item=Self::Elem>; + /// Type of iterators over the elements of the container. + type Iter : Iterator<Item = Self::Elem>; + + /// Returns an iteartor over the elements of the container. fn fl_iter(self) -> Self::Iter; } +/// Trait for a mutable fixed-length container type. pub trait FixedLengthMut<const N : usize> : FixedLength<N> { + /// Type of iterators over references to mutable elements of the container. type IterMut<'a> : Iterator<Item=&'a mut Self::Elem> where Self : 'a; + + /// Returns an iterator over mutable references to elements of the container. fn fl_iter_mut(&mut self) -> Self::IterMut<'_>; } @@ -50,7 +62,7 @@ macro_rules! make_mapmany { ($name:ident, $name_indexed:ident, $var0:ident $($var:ident)* ; $etype0:ident $($etype:ident)*, $ctype0:ident $($ctype:ident)*) => { - /// Map over multiple [`FixedLength`] containers, returning an array. + /// Map over [`FixedLength`] container(s), returning an array. #[inline] pub fn $name< $etype0, @@ -69,8 +81,7 @@ collect_into_array_unchecked(map) } - /// Map over multiple [`FixedLength`] containers, returning an array. - /// This version also passes the index to the mapping function. + /// Map over [`FixedLength`] containers(s) and element indices, returning an array. #[inline] pub fn $name_indexed< $etype0, @@ -161,7 +172,9 @@ // } + /// Iterator returned by [`foldmap`][FoldMappable::foldmap] applied to an iterator. + pub struct FoldMap<I : Iterator<Item=A>, A, B, J : Copy, F : Fn(J, A) -> (J, B)> { iter : I, f : F, @@ -260,9 +273,6 @@ } } - - - impl<A, J : Copy, const N : usize> FoldMappable<A, J> for [A; N] { type Output<B, F> = [B; N] where F : Fn(J, A) -> (J, B); #[inline] @@ -298,7 +308,6 @@ } } -/// /// This is taken and simplified from core::array to not involve `ControlFlow`, /// `Try` etc. (Pulling everything including `NeverShortCircuit` turned out /// too much to maintain here.)
--- a/src/nalgebra_support.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/nalgebra_support.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,4 +1,12 @@ -///! Integration with nalgebra +/*! +Integration with nalgebra. + +This module mainly implements [`Euclidean`], [`Norm`], [`Dot`], [`Linear`], etc. for [`nalgebra`] +matrices and vectors. +It also provides [`ToNalgebraRealField`] as a vomit-inducingly ugly workaround to nalgebra +force-feeding its own versions of the same basic mathematical methods on `f32` and `f64` as +[`num_traits`] does. +*/ use nalgebra::{ Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar, @@ -14,7 +22,7 @@ use std::ops::Mul; use num_traits::identities::{Zero, One}; use crate::linops::*; -use crate::norms::Dot; +use crate::euclidean::*; use crate::types::Float; use crate::norms::*; @@ -247,16 +255,34 @@ } } -/// Helper trait to hide the symbols of `nalgebra::RealField` -/// while allowing nalgebra to be used in subroutines. +/// Helper trait to hide the symbols of [`nalgebra::RealField`]. +/// +/// By assuming `ToNalgebraRealField` intead of `nalgebra::RealField` as a trait bound, +/// functions can piggyback `nalgebra::RealField` without exponsing themselves to it. +/// Thus methods from [`num_traits`] can be used directly without similarly named methods +/// from [`nalgebra`] conflicting with them. Only when absolutely necessary to work with +/// nalgebra, one can convert to the nalgebra view of the same type using the methods of +/// this trait. pub trait ToNalgebraRealField : Float { + /// The nalgebra type corresponding to this type. Usually same as `Self`. + /// + /// This type only carries `nalgebra` traits. type NalgebraType : RealField; + /// The “mixed” type corresponding to this type. Usually same as `Self`. + /// + /// This type carries both `num_traits` and `nalgebra` traits. type MixedType : RealField + Float; + /// Convert to the nalgebra view of `self`. fn to_nalgebra(self) -> Self::NalgebraType; + + /// Convert to the mixed (nalgebra and num_traits) view of `self`. fn to_nalgebra_mixed(self) -> Self::MixedType; + /// Convert from the nalgebra view of `self`. fn from_nalgebra(t : Self::NalgebraType) -> Self; + + /// Convert from the mixed (nalgebra and num_traits) view to `self`. fn from_nalgebra_mixed(t : Self::MixedType) -> Self; }
--- a/src/nanleast.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/nanleast.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,13 +1,21 @@ +/*! +This module provides an artificial total order of floating point numbers. + +The [`NaNLeast`]`<F>` container for `F` a [`Float`] puts `F` in an [`Ord`] total order with +`NaN` the least element. This allows the numbers to be sorted with `NaN` the “least important” +element when looking for the maximum. Thus erroneous computations producing `NaN` can be ignored +when there are good results. +*/ + use crate::types::Float; use std::cmp::{PartialOrd,Ord,Ordering,Ordering::*}; -/// A container for floating point numbers for sorting NaN as a least element -/// (to essentially ignore errorneous computations producing NaNs). +/// A container for floating point numbers. +/// +/// The implementation of [`Ord`] for this type type sorts `NaN` as the least element. #[derive(Debug, Clone, Copy)] pub struct NaNLeast<F : Float>(pub F); -/// Compare floating point numbers ordering nan as the least element. - impl<F : Float> Ord for NaNLeast<F> { #[inline] fn cmp(&self, NaNLeast(b) : &Self) -> Ordering {
--- a/src/norms.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/norms.rs Sat Oct 22 14:40:28 2022 +0300 @@ -2,120 +2,62 @@ Norms, projections, etc. */ +use serde::Serialize; use crate::types::*; -use std::ops::{Mul,MulAssign,Div,DivAssign,Add,Sub,AddAssign,SubAssign,Neg}; -use serde::Serialize; -//use std::iter::Sum; - -// -// Dot products -// - -/// Space with a defined dot product. -pub trait Dot<U,F> { - fn dot(&self, other : &U) -> F; -} - -//self.iter().zip(other.iter()).map(|(&x,&y)| x*y).sum() - -// -// Euclidean spaces -// - -pub trait Euclidean<F : Float> : Sized + Dot<Self,F> - + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F> - + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F> - + Add<Self, Output=<Self as Euclidean<F>>::Output> - + Sub<Self, Output=<Self as Euclidean<F>>::Output> - + for<'b> Add<&'b Self, Output=<Self as Euclidean<F>>::Output> - + for<'b> Sub<&'b Self, Output=<Self as Euclidean<F>>::Output> - + AddAssign<Self> + for<'b> AddAssign<&'b Self> - + SubAssign<Self> + for<'b> SubAssign<&'b Self> - + Neg<Output=<Self as Euclidean<F>>::Output> { - type Output : Euclidean<F>; - - /// Returns the origin of same dimensions as `self`. - fn similar_origin(&self) -> <Self as Euclidean<F>>::Output; - - /// Calculate the square of the 2-norm. - #[inline] - fn norm2_squared(&self) -> F { - self.dot(self) - } - - /// Calculate the square of the 2-norm divided by 2. - #[inline] - fn norm2_squared_div2(&self) -> F { - self.norm2_squared()/F::TWO - } - - /// Calculate the 2-norm. - #[inline] - fn norm2(&self) -> F { - self.norm2_squared().sqrt() - } - - /// Calculate the 2-distance squared. - fn dist2_squared(&self, other : &Self) -> F; - - /// Calculate the 2-distance. - #[inline] - fn dist2(&self, other : &Self) -> F { - self.dist2_squared(other).sqrt() - } - - /// Project to the 2-ball. - #[inline] - fn proj_ball2(mut self, ρ : F) -> Self { - self.proj_ball2_mut(ρ); - self - } - - /// Project to the 2-ball in-place. - #[inline] - fn proj_ball2_mut(&mut self, ρ : F) { - let r = self.norm2(); - if r>ρ { - *self *= ρ/r - } - } -} - -/// Trait for [`Euclidean`] spaces with dimensions known at compile time. -pub trait StaticEuclidean<F : Float> : Euclidean<F> { - /// Returns the origin - fn origin() -> <Self as Euclidean<F>>::Output; -} +use crate::euclidean::*; // // Abstract norms // +/// Exponent type for the 1-[`Norm`]. #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] pub struct L1; +/// Exponent type for the 2-[`Norm`]. #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] pub struct L2; +/// Exponent type for the ∞-[`Norm`]. #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] pub struct Linfinity; +/// Exponent type for 2,1-[`Norm`]. +/// (1-norm over a domain Ω, 2-norm of a vector at each point of the domain.) #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] pub struct L21; +/// A Huber/Moreau–Yosida smoothed [`L1`] norm. (Not a norm itself.) +/// +/// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher +/// values more smoothing. Behaviour with γ < 0 is undefined. #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] pub struct HuberL1<F : Float>(pub F); +/// A Huber/Moreau–Yosida smoothed [`L21`] norm. (Not a norm itself.) +/// +/// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher +/// values more smoothing. Behaviour with γ < 0 is undefined. #[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)] pub struct HuberL21<F : Float>(pub F); +/// A normed space (type) with exponent or other type `Exponent` for the norm. +/// +/// Use as +/// ``` +/// # use alg_tools::norms::{Norm, L1, L2, Linfinity}; +/// # use alg_tools::loc::Loc; +/// let x = Loc([1.0, 2.0, 3.0]); +/// +/// println!("{}, {} {}", x.norm(L1), x.norm(L2), x.norm(Linfinity)) +/// ``` pub trait Norm<F, Exponent> { - /// Calculate the norm + /// Calculate the norm. fn norm(&self, _p : Exponent) -> F; } -/// Indicates that a [`Norm`] is dominated by another norm (`Exponent`) on `Elem` with the -/// corresponding field `F`. +/// Indicates that the `Self`-[`Norm`] is dominated by the `Exponent`-`Norm` on the space +/// `Elem` with the corresponding field `F`. pub trait Dominated<F : Num, Exponent, Elem> { /// Indicates the factor $c$ for the inequality $‖x‖ ≤ C ‖x‖_p$. fn norm_factor(&self, p : Exponent) -> F; @@ -126,19 +68,30 @@ } } +/// Trait for distances with respect to a norm. pub trait Dist<F,Exponent> : Norm<F, Exponent> { /// Calculate the distance fn dist(&self, other : &Self, _p : Exponent) -> F; } +/// Trait for Euclidean projections to the `Exponent`-[`Norm`]-ball. +/// +/// Use as +/// ``` +/// # use alg_tools::norms::{Projection, L2, Linfinity}; +/// # use alg_tools::loc::Loc; +/// let x = Loc([1.0, 2.0, 3.0]); +/// +/// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity)); +/// ``` pub trait Projection<F, Exponent> : Norm<F, Exponent> + Euclidean<F> where F : Float { - /// Project to the norm-ball. + /// Projection of `self` to the `q`-norm-ball of radius ρ. fn proj_ball(mut self, ρ : F, q : Exponent) -> Self { self.proj_ball_mut(ρ, q); self } - /// Project to the norm-ball in-place. + /// In-place projection of `self` to the `q`-norm-ball of radius ρ. fn proj_ball_mut(&mut self, ρ : F, _q : Exponent); } @@ -165,11 +118,11 @@ xn } else { if xn > γ { - xn-γ/F::TWO + xn-γ / F::TWO } else if xn<(-γ) { - -xn-γ/F::TWO + -xn-γ / F::TWO } else { - xnsq/(F::TWO*γ) + xnsq / (F::TWO * γ) } } } @@ -187,47 +140,3 @@ } } -/* -#[inline] -pub fn mean<V>(x : V) -> V::Field where V : ValidArray { - x.iter().sum()/x.len() -} - -#[inline] -pub fn proj_nonneg_mut<V>(x : &mut V) -> &mut V where V : ValidArray { - x.iter_mut().for_each(|&mut p| if p < 0 { *p = 0 } ); - x -} -*/ - - -// -// 2,1-norm generic implementation -// - -/* -pub trait InnerVectors { - type Item; - type Iter : Iterator<Item=Self::Item>; - fn inner_vectors(&self) -> &Self::Iter; -} - -pub trait InnerVectorsMut : InnerVectors { - type IterMut : Iterator<Item=Self::Item>; - fn inner_vectors_mut(&self) -> &mut Self::Item; -} - -impl<F : Float + Sum, T : InnerVectors> Norm<F, L21> for T where T::Item : Norm<F, L2> { - fn norm(&self, _ : L21) -> F { - self.inner_vectors().map(|t| t.norm(L2)).sum() - } -} - -impl<F : Float + Sum, T : InnerVectorsMut + Euclidean<F>> Projection<F, L21> -for T where T::ItemMut : Projection<F, L2> { - fn proj_ball_mut(&mut self, _ : L21) { - self.inner_vectors_mut().for_each(|t| t.proj_ball_mut(L2)); - } -} -*/ -
--- a/src/sets.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/sets.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,20 +1,25 @@ -/// Various types of sets. +/*! +This module provides various sets and traits for them. +*/ + use std::ops::{RangeFull,RangeFrom,Range,RangeInclusive,RangeTo,RangeToInclusive}; use std::marker::PhantomData; use crate::types::*; -use crate::loc::{Loc}; -use crate::norms::Dot; +use crate::loc::Loc; +use crate::euclidean::Dot; use serde::Serialize; -pub mod cube; -pub use cube::*; +mod cube; +pub use cube::Cube; /// Trait for arbitrary sets. The parameter `U` is the element type. pub trait Set<U> where U : ?Sized { - /// Check for containment. + /// Check for element containment fn contains(&self, item : &U) -> bool; } +/// Additional ordering (besides [`PartialOrd`]) of a subfamily of sets: +/// greatest lower bound and least upper bound. pub trait SetOrd : Sized { /// Returns the smallest set of same class contain both parameters. fn common(&self, other : &Self) -> Self; @@ -48,10 +53,15 @@ impl_ranges!(RangeFrom,Range,RangeInclusive,RangeTo,RangeToInclusive); -/// Halfspace. The orthogonal (dual) vectors `A` are supposed to implement [`Dot`]`<U,F>` -/// for `U` the element type. +/// Halfspaces described by an orthogonal vector and an offset. +/// +/// The halfspace is $H = \\{ t v + a \mid a^⊤ v = 0 \\}$, where $v$ is the orthogonal +/// vector and $t$ the offset. +/// +/// `U` is the element type, `F` the floating point number type, and `A` the type of the +/// orthogonal (dual) vectors. They need implement [`Dot`]`<U, F>`. #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] -pub struct Halfspace<A, F, U> where A : Dot<U,F>, F : Float { +pub struct Halfspace<A, F, U> where A : Dot<U, F>, F : Float { pub orthogonal : A, pub offset : F, _phantom : PhantomData<U>, @@ -64,8 +74,11 @@ } } -pub trait SpannedHalfspace<F,U> where F : Float { - type A : Dot<U,F>; +/// Trait for generating a halfspace spanned by another set `Self` of elements of type `U`. +pub trait SpannedHalfspace<F, U> where F : Float { + /// Type of the orthogonal vector describing the halfspace. + type A : Dot<U, F>; + /// Returns the halfspace spanned by this set. fn spanned_halfspace(&self) -> Halfspace<Self::A, F, U>; }
--- a/src/sets/cube.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/sets/cube.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,4 +1,20 @@ -///! Presetnation of cubes $[a_1, b_1] × ⋯ × [a_n, b_n]$ +/*! +Multi-dimensional cubes. + +This module provides the [`Cube`] type for multi-dimensional cubes $∏_{i=1}^N [a_i, b_i)$. + +As an example, to create a the two-dimensional cube $[0, 1] × [-1, 1]$, you can +``` +# use alg_tools::sets::cube::Cube; +let cube = Cube::new([[0.0, 1.0], [-1.0, 1.0]]); +``` +or +``` +# use alg_tools::sets::cube::Cube; +# use alg_tools::types::float; +let cube : Cube<float, 2> = [[0.0, 1.0], [-1.0, 1.0]].into(); +``` +*/ use serde::ser::{Serialize, Serializer, SerializeTupleStruct}; use crate::types::*; @@ -12,7 +28,8 @@ map2, }; -/// A half-open `N`-cube of elements of type `U`. +/// A multi-dimensional cube $∏_{i=1}^N [a_i, b_i)$ with the starting and ending points +/// along $a_i$ and $b_i$ along each dimension of type `U`. #[derive(Copy, Clone, Debug, Eq, PartialEq)] pub struct Cube<U : Num, const N : usize>(pub(super) [[U; 2]; N]); @@ -82,51 +99,68 @@ } impl<U : Num, const N : usize> Cube<U, N> { + /// Maps `f` over the triples $\\{(i, a\_i, b\_i)\\}\_{i=1}^N$ + /// of the cube $∏_{i=1}^N [a_i, b_i)$. #[inline] pub fn map_indexed<T>(&self, f : impl Fn(usize, U, U) -> T) -> [T; N] { map1_indexed(self, |i, &[a, b]| f(i, a, b)) } + /// Maps `f` over the tuples $\\{(a\_i, b\_i)\\}\_{i=1}^N$ + /// of the cube $∏_{i=1}^N [a_i, b_i)$. #[inline] pub fn map<T>(&self, f : impl Fn(U, U) -> T) -> [T; N] { map1(self, |&[a, b]| f(a, b)) } + /// Iterates over the start and end coordinates $\{(a_i, b_i)\}_{i=1}^N$ of the cube along + /// each dimension. #[inline] pub fn iter_coords(&self) -> std::slice::Iter<'_, [U; 2]> { self.0.iter() } + /// Returns the “start” coordinate $a_i$ of the cube $∏_{i=1}^N [a_i, b_i)$. #[inline] pub fn start(&self, i : usize) -> U { self.0[i][0] } + /// Returns the end coordinate $a_i$ of the cube $∏_{i=1}^N [a_i, b_i)$. #[inline] pub fn end(&self, i : usize) -> U { self.0[i][1] } - + + /// Returns the “start” $(a_1, … ,a_N)$ of the cube $∏_{i=1}^N [a_i, b_i)$ + /// spanned between $(a_1, … ,a_N)$ and $(b_1, … ,b_N)$. #[inline] pub fn span_start(&self) -> Loc<U, N> { Loc::new(self.map(|a, _b| a)) } + /// Returns the end $(b_1, … ,b_N)$ of the cube $∏_{i=1}^N [a_i, b_i)$ + /// spanned between $(a_1, … ,a_N)$ and $(b_1, … ,b_N)$. #[inline] pub fn span_end(&self) -> Loc<U, N> { Loc::new(self.map(|_a, b| b)) } + /// Iterates over the corners $\{(c_1, … ,c_N) | c_i ∈ \{a_i, b_i\}\}$ of the cube + /// $∏_{i=1}^N [a_i, b_i)$. #[inline] pub fn iter_corners(&self) -> CubeCornersIter<'_, U, N> { CubeCornersIter{ index : 0, cube : self } } + /// Returns the width-`N`-tuple $(b_1-a_1, … ,b_N-a_N)$ of the cube $∏_{i=1}^N [a_i, b_i)$. #[inline] pub fn width(&self) -> Loc<U, N> { Loc::new(self.map(|a, b| b-a)) } + /// Translates the cube $∏_{i=1}^N [a_i, b_i)$ by the `shift` $(s_1, … , s_N)$ to + /// $∏_{i=1}^N [a_i+s_i, b_i+s_i)$. #[inline] pub fn shift(&self, shift : &Loc<U, N>) -> Self { let mut cube = self.clone(); @@ -137,6 +171,7 @@ cube } + /// Creates a new cube from an array. #[inline] pub fn new(data : [[U; 2]; N]) -> Self { Cube(data) @@ -152,6 +187,7 @@ impl<U : Num> Cube<U, 1> { /// Get the corners of the cube. + /// /// TODO: generic implementation once const-generics can be involved in /// calculations. #[inline] @@ -163,6 +199,7 @@ impl<U : Num> Cube<U, 2> { /// Get the corners of the cube in counter-clockwise order. + /// /// TODO: generic implementation once const-generics can be involved in /// calculations. #[inline] @@ -177,6 +214,7 @@ impl<U : Num> Cube<U, 3> { /// Get the corners of the cube. + /// /// TODO: generic implementation once const-generics can be involved in /// calculations. #[inline]
--- a/src/tabledump.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/tabledump.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,3 +1,7 @@ + +/*! +Helper traits and functions for dumping data into CSV or JSON files. +*/ use std::slice::Iter; use csv; @@ -5,14 +9,14 @@ use serde::Serialize; use crate::error::DynError; -/// Write a CSV from an iterator -pub fn write_csv<I, J>(mut iter : I, f : String) -> DynError +/// Write a CSV from an iterator of [serde-serializable][Serialize] items. +pub fn write_csv<I, J>(mut iter : I, filename : String) -> DynError where I : Iterator<Item=J>, J : Serialize { let wtr = csv::WriterBuilder::new() .has_headers(true) .delimiter(b'\t') - .from_path(f); + .from_path(filename); wtr.and_then(|mut w|{ //w.write_record(self.tabledump_headers())?; iter.try_for_each(|item| w.serialize(item)) @@ -20,25 +24,30 @@ Ok(()) } -/// Write a JSON from an iterator -pub fn write_json<I, J>(iter : I, f : String) -> DynError +/// Write a JSON from an iterator of [serde-serializable][Serialize] items. +pub fn write_json<I, J>(iter : I, filename : String) -> DynError where I : Iterator<Item=J>, J : Serialize { let v : Vec<J> = iter.collect(); - serde_json::to_writer_pretty(std::fs::File::create(f)?, &v)?; + serde_json::to_writer_pretty(std::fs::File::create(filename)?, &v)?; Ok(()) } -/// Helper trait for dumping data in a CSV or JSON file. +/// Trait that should be implemented by types that can be dumped into CSV (or JSON) files. +/// +/// The method that needs to be implemented is [`Self::tabledump_entries`] that prodvides an +/// iterator over the items to be dumped. Moreover, the items themselves need to be +/// [serde-serializable][`Serialize`]. Caveats about [`csv`] crate restrictions on writable +/// data apply for CSV writing. pub trait TableDump<'a> where <Self::Iter as Iterator>::Item : Serialize { - /// Iterator over the rows + /// Iterator over the items to be dumped. type Iter : Iterator; // Return the headers of the CSV file. //fn tabledump_headers(&'a self) -> Vec<String>; - /// Return an iterator over the rows of the CSV file. + /// Return an iterator over the rows that would be dumped into a CSV file. fn tabledump_entries(&'a self) -> Self::Iter; /// Write a CSV file.
--- a/src/tuple.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/tuple.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,19 +1,31 @@ -// -// Tuple tools -// +/*! +Helper traits for working with tuples. +*/ /// A trait for removing the first item from a tuple and returning the rest. /// Only implemented for a few low dimensions. pub trait TupleOps : Sized { + /// First element of the tuple. type Head; + /// Rest of the tuple. type Tail; + + /// Nested 2-tuple conversioon of an n-tuple. type Nested; - /// Remove the first item from a tuple and returning the rest. + /// Remove the first item from a tuple and return the rest. fn tail(self) -> Self::Tail { self.split_tail().1 } + + /// Return the first element and the rest of the tuple fn split_tail(self) -> (Self::Head, Self::Tail); + + /// Convert `(a,b,c,…)` into `(a, (b, (c, …)))`. fn nest(self) -> Self::Nested; + + /// Convert `(a, (b, (c, …)))` into `(a,b,c,…)`. fn from_nested(nested : Self::Nested) -> Self; + + // Convert `head` and `tail` into `(head, tail…)`. fn from_head_tail(head : Self::Head, tail : Self::Tail) -> Self; }
--- a/src/types.rs Mon Oct 24 10:52:19 2022 +0300 +++ b/src/types.rs Sat Oct 22 14:40:28 2022 +0300 @@ -1,6 +1,15 @@ -///! Some useful (numerical) types and type traits +/*! +Some useful (numerical) types and traits. + +The traits are based on corresponding ones in [`num_traits`], but try to fill some gaps in the +super-traits and available constants. -use trait_set::trait_set; +As [`nalgebra`] unnecessarily provides many of the same methods as [`num_traits`], to avoid having +to refer to the methods with the full path, it is often necesary to use [`ToNalgebraRealField`][crate::nalgebra_support::ToNalgebraRealField] to hide the nalgebra implementations until +absolutely necessary to use nalgebra. +*/ + +//use trait_set::trait_set; pub use num_traits::Float as NumTraitsFloat; // needed to re-export functions. pub use num_traits::cast::AsPrimitive; @@ -40,7 +49,7 @@ i8 i16 i32 i64 i128 isize f32 f64); -/// Trait for numeric types +/// Trait for general numeric types pub trait Num : 'static + Copy + num::Num + num_traits::NumAssign + std::iter::Sum + std::iter::Product + std::fmt::Debug + std::fmt::Display + serde::Serialize @@ -64,10 +73,10 @@ impl<U : Num + num::Signed + std::ops::Neg<Output=Self>> SignedNum for U { } /// Trait for floating point numbers -pub trait Float : SignedNum + num::Float + From<Self::CompatibleSize> { - /// An unsigned integer that can be used for indexing operations and - /// converted to F without loss. - type CompatibleSize : CompatibleUnsigned<Self>; +pub trait Float : SignedNum + num::Float /*+ From<Self::CompatibleSize>*/ { + // An unsigned integer that can be used for indexing operations and + // converted to F without loss. + //type CompatibleSize : CompatibleUnsigned<Self>; const PI : Self; const E : Self; @@ -113,10 +122,10 @@ impl_num_consts!(f32 f64); impl Float for f64 { - #[cfg(any(target_pointer_width = "128", target_pointer_width = "64"))] + /*#[cfg(any(target_pointer_width = "128", target_pointer_width = "64"))] type CompatibleSize = u32; #[cfg(any(target_pointer_width = "32", target_pointer_width = "16"))] - type CompatibleSize = usize; + type CompatibleSize = usize;*/ const PI : Self = std::f64::consts::PI; const E : Self = std::f64::consts::E; @@ -128,10 +137,12 @@ } impl Float for f32 { + /* #[cfg(any(target_pointer_width = "128", target_pointer_width = "64", target_pointer_width = "32"))] type CompatibleSize = u16; #[cfg(any(target_pointer_width = "16"))] type CompatibleSize = usize; + */ const PI : Self = std::f32::consts::PI; const E : Self = std::f32::consts::E; @@ -142,7 +153,9 @@ const FRAC_2_SQRT_PI : Self = std::f32::consts::FRAC_2_SQRT_PI; } +/* trait_set! { pub trait CompatibleUnsigned<F : Float> = Unsigned + Into<F>; pub trait CompatibleSigned<F : Float> = Signed + Into<F>; } +*/
