changeset
Transpose loc parameters to allow f64 defaults
Thu, 01 May 2025 13:06:58 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Thu, 01 May 2025 13:06:58 -0500
- branch
- dev
- changeset 124
- 6aa955ad8122
- parent 122
- 495448cca603
- child 125
- 25b6c8b20d1b
Transpose loc parameters to allow f64 defaults
--- a/src/bisection_tree/bt.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/bisection_tree/bt.rs Thu May 01 13:06:58 2025 -0500 @@ -51,7 +51,7 @@ #[derive(Clone, Debug)] pub(super) struct Branches<F: Num, D, A: Aggregator, const N: usize, const P: usize> { /// Point for subdivision of the (unstored) [`Cube`] corresponding to the node. - pub(super) branch_at: Loc<F, N>, + pub(super) branch_at: Loc<N, F>, /// Subnodes pub(super) nodes: [Node<F, D, A, N, P>; P], } @@ -184,7 +184,7 @@ /// /// 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 { + fn get_node_index(&self, x: &Loc<N, F>) -> 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() @@ -195,25 +195,25 @@ /// This only takes the branch subdivision point $d$ into account, so is always succesfull. /// Thus, for this point, each branch corresponds to a quadrant of $ℝ^N$ relative to $d$. #[inline] - fn get_node(&self, x: &Loc<F, N>) -> &Node<F, D, A, N, P> { + fn get_node(&self, x: &Loc<N, F>) -> &Node<F, D, A, N, P> { &self.nodes[self.get_node_index(x)] } } /// An iterator over the $P=2^N$ subcubes of a [`Cube`] subdivided at a point `d`. pub(super) struct SubcubeIter<'b, F: Float, const N: usize, const P: usize> { - domain: &'b Cube<F, N>, - branch_at: Loc<F, N>, + domain: &'b Cube<N, F>, + branch_at: Loc<N, F>, 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>, + branch_at: &Loc<N, F>, + domain: &Cube<N, F>, i: usize, -) -> Cube<F, N> { +) -> Cube<N, F> { map2_indexed(branch_at, domain, move |j, &branch, &[start, end]| { if i & (1 << j) != 0 { [branch, end] @@ -225,7 +225,7 @@ } impl<'a, 'b, F: Float, const N: usize, const P: usize> Iterator for SubcubeIter<'b, F, N, P> { - type Item = Cube<F, N>; + type Item = Cube<N, F>; #[inline] fn next(&mut self) -> Option<Self::Item> { if self.index < P { @@ -246,8 +246,8 @@ { /// 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> + Support<F, N>>( - domain: &Cube<F, N>, + pub(super) fn new_with<S: LocalAnalysis<F, A, N> + Support<N, F>>( + domain: &Cube<N, F>, support: &S, ) -> Self { let hint = support.bisection_hint(domain); @@ -272,7 +272,7 @@ /// Returns an iterator over the subcubes of `domain` subdivided at the branching point /// of `self`. #[inline] - pub(super) fn iter_subcubes<'b>(&self, domain: &'b Cube<F, N>) -> SubcubeIter<'b, F, N, P> { + pub(super) fn iter_subcubes<'b>(&self, domain: &'b Cube<N, F>) -> SubcubeIter<'b, F, N, P> { SubcubeIter { domain: domain, branch_at: self.branch_at, @@ -283,7 +283,7 @@ /* /// Returns an iterator over all nodes and corresponding subcubes of `self`. #[inline] - pub(super) 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<N, F>) -> std::iter::Zip<Iter<'a, Node<F,D,A,N,P>>, SubcubeIter<'b, F, N, P>> { self.nodes.iter().zip(self.iter_subcubes(domain)) } @@ -293,7 +293,7 @@ #[inline] pub(super) fn nodes_and_cubes_mut<'a, 'b>( &'a mut self, - domain: &'b Cube<F, N>, + domain: &'b Cube<N, F>, ) -> 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) @@ -303,10 +303,10 @@ #[inline] fn recurse<'scope, 'smaller, 'refs>( &'smaller mut self, - domain: &'smaller Cube<F, N>, + domain: &'smaller Cube<N, F>, task_budget: TaskBudget<'scope, 'refs>, - guard: impl Fn(&Node<F, D, A, N, P>, &Cube<F, N>) -> bool + Send + 'smaller, - mut f: impl for<'a> FnMut(&mut Node<F, D, A, N, P>, &Cube<F, N>, TaskBudget<'smaller, 'a>) + guard: impl Fn(&Node<F, D, A, N, P>, &Cube<N, F>) -> bool + Send + 'smaller, + mut f: impl for<'a> FnMut(&mut Node<F, D, A, N, P>, &Cube<N, F>, TaskBudget<'smaller, 'a>) + Send + Copy + 'smaller, @@ -332,9 +332,9 @@ /// * `support` is the [`Support`] that is used determine with which subcubes of `domain` /// (at subdivision depth `new_leaf_depth`) the data `d` is to be associated with. /// - pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis<F, A, N> + Support<F, N>>( + pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis<F, A, N> + Support<N, F>>( &mut self, - domain: &Cube<F, N>, + domain: &Cube<N, F>, d: D, new_leaf_depth: M, support: &S, @@ -360,11 +360,11 @@ pub(super) fn convert_aggregator<ANew, G>( self, generator: &G, - domain: &Cube<F, N>, + domain: &Cube<N, F>, ) -> Branches<F, D, ANew, N, P> where ANew: Aggregator, - G: SupportGenerator<F, N, Id = D>, + G: SupportGenerator<N, F, Id = D>, G::SupportType: LocalAnalysis<F, ANew, N>, { let branch_at = self.branch_at; @@ -387,10 +387,10 @@ pub(super) fn refresh_aggregator<'refs, 'scope, G>( &mut self, generator: &G, - domain: &Cube<F, N>, + domain: &Cube<N, F>, task_budget: TaskBudget<'scope, 'refs>, ) where - G: SupportGenerator<F, N, Id = D>, + G: SupportGenerator<N, F, Id = D>, G::SupportType: LocalAnalysis<F, A, N>, { self.recurse( @@ -422,7 +422,7 @@ /* /// Get leaf data #[inline] - pub(super) fn get_leaf_data(&self, x : &Loc<F, N>) -> Option<&Vec<D>> { + pub(super) fn get_leaf_data(&self, x : &Loc<N, F>) -> Option<&Vec<D>> { match self.data { NodeOption::Uninitialised => None, NodeOption::Leaf(ref data) => Some(data), @@ -432,7 +432,7 @@ /// Get leaf data iterator #[inline] - pub(super) fn get_leaf_data_iter(&self, x: &Loc<F, N>) -> Option<std::slice::Iter<'_, D>> { + pub(super) fn get_leaf_data_iter(&self, x: &Loc<N, F>) -> Option<std::slice::Iter<'_, D>> { match self.data { NodeOption::Uninitialised => None, NodeOption::Leaf(ref data) => Some(data.iter()), @@ -459,9 +459,9 @@ /// If `self` is a [`NodeOption::Branches`], the data is passed to branches whose subcubes /// `support` intersects. If an [`NodeOption::Uninitialised`] node is encountered, a new leaf is /// created at a minimum depth of `new_leaf_depth`. - pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis<F, A, N> + Support<F, N>>( + pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis<F, A, N> + Support<N, F>>( &mut self, - domain: &Cube<F, N>, + domain: &Cube<N, F>, d: D, new_leaf_depth: M, support: &S, @@ -515,11 +515,11 @@ pub(super) fn convert_aggregator<ANew, G>( mut self, generator: &G, - domain: &Cube<F, N>, + domain: &Cube<N, F>, ) -> Node<F, D, ANew, N, P> where ANew: Aggregator, - G: SupportGenerator<F, N, Id = D>, + G: SupportGenerator<N, F, Id = D>, G::SupportType: LocalAnalysis<F, ANew, N>, { // The mem::replace is needed due to the [`Drop`] implementation to extract self.data. @@ -561,10 +561,10 @@ pub(super) fn refresh_aggregator<'refs, 'scope, G>( &mut self, generator: &G, - domain: &Cube<F, N>, + domain: &Cube<N, F>, task_budget: TaskBudget<'scope, 'refs>, ) where - G: SupportGenerator<F, N, Id = D>, + G: SupportGenerator<N, F, Id = D>, G::SupportType: LocalAnalysis<F, A, N>, { match &mut self.data { @@ -610,7 +610,7 @@ /// 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>: +pub trait BTImpl<const N: usize, F: Float = f64>: std::fmt::Debug + Clone + GlobalAnalysis<F, Self::Agg> { /// The data type stored in the tree @@ -621,7 +621,7 @@ /// 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> + type Converted<ANew>: BTImpl<N, F, Data = Self::Data, Agg = ANew> where ANew: Aggregator; @@ -629,7 +629,7 @@ /// /// Every leaf node of the tree that intersects the `support` will contain a copy of /// `d`. - fn insert<S: LocalAnalysis<F, Self::Agg, N> + Support<F, N>>( + fn insert<S: LocalAnalysis<F, Self::Agg, N> + Support<N, F>>( &mut self, d: Self::Data, support: &S, @@ -642,7 +642,7 @@ fn convert_aggregator<ANew, G>(self, generator: &G) -> Self::Converted<ANew> where ANew: Aggregator, - G: SupportGenerator<F, N, Id = Self::Data>, + G: SupportGenerator<N, F, Id = Self::Data>, G::SupportType: LocalAnalysis<F, ANew, N>; /// Refreshes the aggregator of the three after possible changes to the support generator. @@ -651,19 +651,19 @@ /// into corresponding [`Support`]s. fn refresh_aggregator<G>(&mut self, generator: &G) where - G: SupportGenerator<F, N, Id = Self::Data>, + G: SupportGenerator<N, F, Id = Self::Data>, G::SupportType: LocalAnalysis<F, Self::Agg, N>; /// Returns an iterator over all [`Self::Data`] items at the point `x` of the domain. - fn iter_at(&self, x: &Loc<F, N>) -> std::slice::Iter<'_, Self::Data>; + fn iter_at(&self, x: &Loc<N, F>) -> std::slice::Iter<'_, Self::Data>; /* /// Returns all [`Self::Data`] items at the point `x` of the domain. - fn data_at(&self, x : &Loc<F,N>) -> Arc<Vec<Self::Data>>; + fn data_at(&self, x : &Loc<N, F>) -> Arc<Vec<Self::Data>>; */ /// Create a new tree on `domain` of indicated `depth`. - fn new(domain: Cube<F, N>, depth: Self::Depth) -> Self; + fn new(domain: Cube<N, F>, depth: Self::Depth) -> Self; } /// The main bisection tree structure. @@ -679,7 +679,7 @@ /// The depth of the tree (initial, before refinement) pub(super) depth: M, /// The domain of the toplevel node - pub(super) domain: Cube<F, N>, + pub(super) domain: Cube<N, F>, /// The toplevel node of the tree pub(super) topnode: <BTNodeLookup as BTNode<F, D, A, N>>::Node, } @@ -693,7 +693,7 @@ type Node = Node<F,D,A,$n,{pow(2, $n)}>; } - impl<M,F,D,A> BTImpl<F,$n> for BT<M,F,D,A,$n> + impl<M,F,D,A> BTImpl<$n, F> for BT<M,F,D,A,$n> where M : Depth, F : Float, D : 'static + Copy + Send + Sync + std::fmt::Debug, @@ -703,7 +703,7 @@ type Agg = A; type Converted<ANew> = BT<M,F,D,ANew,$n> where ANew : Aggregator; - fn insert<S: LocalAnalysis<F, A, $n> + Support<F, $n>>( + fn insert<S: LocalAnalysis<F, A, $n> + Support< $n, F>>( &mut self, d : D, support : &S @@ -721,7 +721,7 @@ fn convert_aggregator<ANew, G>(self, generator : &G) -> Self::Converted<ANew> where ANew : Aggregator, - G : SupportGenerator<F, $n, Id=D>, + G : SupportGenerator< $n, F, Id=D>, G::SupportType : LocalAnalysis<F, ANew, $n> { let topnode = self.topnode.convert_aggregator(generator, &self.domain); @@ -733,22 +733,22 @@ } fn refresh_aggregator<G>(&mut self, generator : &G) - where G : SupportGenerator<F, $n, Id=Self::Data>, + where G : SupportGenerator< $n, F, Id=Self::Data>, G::SupportType : LocalAnalysis<F, Self::Agg, $n> { with_task_budget(|task_budget| self.topnode.refresh_aggregator(generator, &self.domain, task_budget) ) } - /*fn data_at(&self, x : &Loc<F,$n>) -> Arc<Vec<D>> { + /*fn data_at(&self, x : &Loc<$n, F>) -> Arc<Vec<D>> { self.topnode.get_leaf_data(x).unwrap_or_else(|| Arc::new(Vec::new())) }*/ - fn iter_at(&self, x : &Loc<F,$n>) -> std::slice::Iter<'_, D> { + fn iter_at(&self, x : &Loc<$n, F>) -> std::slice::Iter<'_, D> { self.topnode.get_leaf_data_iter(x).unwrap_or_else(|| [].iter()) } - fn new(domain : Cube<F, $n>, depth : M) -> Self { + fn new(domain : Cube<$n, F>, depth : M) -> Self { BT { depth : depth, domain : domain,
--- a/src/bisection_tree/btfn.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/bisection_tree/btfn.rs Thu May 01 13:06:58 2025 -0500 @@ -22,7 +22,7 @@ /// 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, +/// The domain of the function is [`Loc`]`<N, F>`, where `F` is the type of floating point numbers, /// and `N` the dimension. /// /// The `generator` lists the component functions that have to implement [`Support`]. @@ -30,7 +30,7 @@ /// in a [bisection tree][BTImpl], when one is provided as `bt`. However `bt` may also be `()` /// for a [`PreBTFN`] that is only useful for vector space operations with a full [`BTFN`]. #[derive(Clone, Debug)] -pub struct BTFN<F: Float, G: SupportGenerator<F, N>, BT /*: BTImpl<F, N>*/, const N: usize> /*where G::SupportType : LocalAnalysis<F, A, N>*/ +pub struct BTFN<F: Float, G: SupportGenerator<N, F>, BT /*: BTImpl< N, F>*/, const N: usize> /*where G::SupportType : LocalAnalysis<F, A, N>*/ { bt: BT, generator: Arc<G>, @@ -39,18 +39,18 @@ impl<F: Float, G, BT, const N: usize> Space for BTFN<F, G, BT, N> where - G: SupportGenerator<F, N, Id = BT::Data>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, - BT: BTImpl<F, N>, + BT: BTImpl<N, F>, { type Decomp = BasicDecomposition; } impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N> where - G: SupportGenerator<F, N, Id = BT::Data>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, - BT: BTImpl<F, N>, + BT: BTImpl<N, F>, { /// Create a new BTFN from a support generator and a pre-initialised bisection tree. /// @@ -92,11 +92,11 @@ /// The top node of the created [`BT`] will have the given `domain`. /// /// See the documentation for [`BTFN`] on the role of the `generator`. - pub fn construct(domain: Cube<F, N>, depth: BT::Depth, generator: G) -> Self { + pub fn construct(domain: Cube<N, F>, depth: BT::Depth, generator: G) -> Self { Self::construct_arc(domain, depth, Arc::new(generator)) } - fn construct_arc(domain: Cube<F, N>, depth: BT::Depth, generator: Arc<G>) -> Self { + fn construct_arc(domain: Cube<N, F>, depth: BT::Depth, generator: Arc<G>) -> Self { let mut bt = BT::new(domain, depth); for (d, support) in generator.all_data() { bt.insert(d, &support); @@ -111,7 +111,7 @@ pub fn convert_aggregator<ANew>(self) -> BTFN<F, G, BT::Converted<ANew>, N> where ANew: Aggregator, - G: SupportGenerator<F, N, Id = BT::Data>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, ANew, N>, { BTFN::new_arc(self.bt.convert_aggregator(&*self.generator), self.generator) @@ -130,15 +130,15 @@ impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N> where - G: SupportGenerator<F, N>, + G: SupportGenerator<N, F>, { /// 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>>( + pub fn instantiate<BTNew: BTImpl<N, F, Data = G::Id>>( self, - domain: Cube<F, N>, + domain: Cube<N, F>, depth: BTNew::Depth, ) -> BTFN<F, G, BTNew, N> where @@ -157,7 +157,7 @@ impl<F: Float, G, const N: usize> PreBTFN<F, G, N> where - G: SupportGenerator<F, N>, + G: SupportGenerator<N, F>, { /// Create a new [`PreBTFN`] with no bisection tree. pub fn new_pre(generator: G) -> Self { @@ -171,14 +171,14 @@ impl<F: Float, G, BT, const N: usize> BTFN<F, G, BT, N> where - G: SupportGenerator<F, N, Id = usize>, + G: SupportGenerator<N, F, Id = usize>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, - BT: BTImpl<F, N, Data = usize>, + BT: BTImpl<N, F, Data = usize>, { /// Helper function for implementing [`std::ops::Add`]. fn add_another<G2>(&self, g2: Arc<G2>) -> BTFN<F, BothGenerators<G, G2>, BT, N> where - G2: SupportGenerator<F, N, Id = usize>, + G2: SupportGenerator<N, F, Id = usize>, G2::SupportType: LocalAnalysis<F, BT::Agg, N>, { let mut bt = self.bt.clone(); @@ -201,9 +201,9 @@ impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Add<BTFN<F, G2, BT2, N>> for $lhs - where BT1 : BTImpl<F, N, Data=usize>, - G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, - G2 : SupportGenerator<F, N, Id=usize>, + where BT1 : BTImpl< N, F, Data=usize>, + G1 : SupportGenerator< N, F, Id=usize> + $($extra_trait)?, + G2 : SupportGenerator< N, F, Id=usize>, G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; @@ -216,9 +216,9 @@ impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Add<&'b BTFN<F, G2, BT2, N>> for $lhs - where BT1 : BTImpl<F, N, Data=usize>, - G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, - G2 : SupportGenerator<F, N, Id=usize>, + where BT1 : BTImpl< N, F, Data=usize>, + G1 : SupportGenerator< N, F, Id=usize> + $($extra_trait)?, + G2 : SupportGenerator< N, F, Id=usize>, G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { @@ -239,9 +239,9 @@ impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Sub<BTFN<F, G2, BT2, N>> for $lhs - where BT1 : BTImpl<F, N, Data=usize>, - G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, - G2 : SupportGenerator<F, N, Id=usize>, + where BT1 : BTImpl< N, F, Data=usize>, + G1 : SupportGenerator< N, F, Id=usize> + $($extra_trait)?, + G2 : SupportGenerator< N, F, Id=usize>, G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, G2::SupportType : LocalAnalysis<F, BT1::Agg, N> { type Output = BTFN<F, BothGenerators<G1, G2>, BT1, N>; @@ -258,9 +258,9 @@ impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Sub<&'b BTFN<F, G2, BT2, N>> for $lhs - where BT1 : BTImpl<F, N, Data=usize>, - G1 : SupportGenerator<F, N, Id=usize> + $($extra_trait)?, - G2 : SupportGenerator<F, N, Id=usize> + Clone, + where BT1 : BTImpl< N, F, Data=usize>, + G1 : SupportGenerator< N, F, Id=usize> + $($extra_trait)?, + G2 : SupportGenerator< N, F, Id=usize> + Clone, G1::SupportType : LocalAnalysis<F, BT1::Agg, N>, G2::SupportType : LocalAnalysis<F, BT1::Agg, N>, &'b G2 : std::ops::Neg<Output=G2> { @@ -281,8 +281,8 @@ ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl<F: Float, G, BT, const N: usize> std::ops::$trait_assign<F> for BTFN<F, G, BT, N> where - BT: BTImpl<F, N>, - G: SupportGenerator<F, N, Id = BT::Data>, + BT: BTImpl<N, F>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, { #[inline] @@ -294,8 +294,8 @@ impl<F: Float, G, BT, const N: usize> std::ops::$trait<F> for BTFN<F, G, BT, N> where - BT: BTImpl<F, N>, - G: SupportGenerator<F, N, Id = BT::Data>, + BT: BTImpl<N, F>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, { type Output = Self; @@ -309,8 +309,8 @@ impl<'a, F: Float, G, BT, const N: usize> std::ops::$trait<F> for &'a BTFN<F, G, BT, N> where - BT: BTImpl<F, N>, - G: SupportGenerator<F, N, Id = BT::Data>, + BT: BTImpl<N, F>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, &'a G: std::ops::$trait<F, Output = G>, { @@ -331,8 +331,8 @@ impl<G, BT, const N : usize> std::ops::$trait<BTFN<$f, G, BT, N>> for $f - where BT : BTImpl<$f, N>, - G : SupportGenerator<$f, N, Id=BT::Data>, + where BT : BTImpl< N, $f>, + G : SupportGenerator< N, $f, Id=BT::Data>, G::SupportType : LocalAnalysis<$f, BT::Agg, N> { type Output = BTFN<$f, G, BT, N>; #[inline] @@ -346,8 +346,8 @@ impl<'a, G, BT, const N : usize> std::ops::$trait<&'a BTFN<$f, G, BT, N>> for $f - where BT : BTImpl<$f, N>, - G : SupportGenerator<$f, N, Id=BT::Data> + Clone, + where BT : BTImpl< N, $f>, + G : SupportGenerator< N, $f, Id=BT::Data> + Clone, G::SupportType : LocalAnalysis<$f, BT::Agg, N>, // FIXME: This causes compiler overflow /*&'a G : std::ops::$trait<$f,Output=G>*/ { @@ -371,8 +371,8 @@ ($trait:ident, $fn:ident) => { impl<F: Float, G, BT, const N: usize> std::ops::$trait for BTFN<F, G, BT, N> where - BT: BTImpl<F, N>, - G: SupportGenerator<F, N, Id = BT::Data>, + BT: BTImpl<N, F>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, { type Output = Self; @@ -387,8 +387,8 @@ /*impl<'a, F : Float, G, BT, const N : usize> std::ops::$trait for &'a BTFN<F, G, BT, N> - where BT : BTImpl<F, N>, - G : SupportGenerator<F, N, Id=BT::Data>, + where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis<F, BT::Agg, N>, &'a G : std::ops::$trait<Output=G> { type Output = BTFN<F, G, BT, N>; @@ -406,16 +406,16 @@ // Apply, Mapping, Differentiate // -impl<F: Float, G, BT, V, const N: usize> Mapping<Loc<F, N>> for BTFN<F, G, BT, N> +impl<F: Float, G, BT, V, const N: usize> Mapping<Loc<N, F>> for BTFN<F, G, BT, N> where - BT: BTImpl<F, N>, - G: SupportGenerator<F, N, Id = BT::Data>, - G::SupportType: LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<F, N>, Codomain = V>, + BT: BTImpl<N, F>, + G: SupportGenerator<N, F, Id = BT::Data>, + G::SupportType: LocalAnalysis<F, BT::Agg, N> + Mapping<Loc<N, F>, Codomain = V>, V: Sum + Space, { type Codomain = V; - fn apply<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Codomain { + fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain { let xc = x.cow(); self.bt .iter_at(&*xc) @@ -424,17 +424,17 @@ } } -impl<F: Float, G, BT, V, const N: usize> DifferentiableImpl<Loc<F, N>> for BTFN<F, G, BT, N> +impl<F: Float, G, BT, V, const N: usize> DifferentiableImpl<Loc<N, F>> for BTFN<F, G, BT, N> where - BT: BTImpl<F, N>, - G: SupportGenerator<F, N, Id = BT::Data>, + BT: BTImpl<N, F>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: - LocalAnalysis<F, BT::Agg, N> + DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>, + LocalAnalysis<F, BT::Agg, N> + DifferentiableMapping<Loc<N, F>, DerivativeDomain = V>, V: Sum + Space, { type Derivative = V; - fn differential_impl<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Derivative { + fn differential_impl<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Derivative { let xc = x.cow(); self.bt .iter_at(&*xc) @@ -449,8 +449,8 @@ impl<F: Float, G, BT, const N: usize> GlobalAnalysis<F, BT::Agg> for BTFN<F, G, BT, N> where - BT: BTImpl<F, N>, - G: SupportGenerator<F, N, Id = BT::Data>, + BT: BTImpl<N, F>, + G: SupportGenerator<N, F, Id = BT::Data>, G::SupportType: LocalAnalysis<F, BT::Agg, N>, { #[inline] @@ -467,8 +467,8 @@ /* impl<'b, X, F : Float, G, BT, const N : usize> Apply<&'b X, F> for BTFN<F, G, BT, N> -where BT : BTImpl<F, N>, - G : SupportGenerator<F, N, Id=BT::Data>, +where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis<F, BT::Agg, N>, X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> { @@ -480,8 +480,8 @@ impl<X, F : Float, G, BT, const N : usize> Apply<X, F> for BTFN<F, G, BT, N> -where BT : BTImpl<F, N>, - G : SupportGenerator<F, N, Id=BT::Data>, +where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis<F, BT::Agg, N>, X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> { @@ -493,8 +493,8 @@ impl<X, F : Float, G, BT, const N : usize> Linear<X> for BTFN<F, G, BT, N> -where BT : BTImpl<F, N>, - G : SupportGenerator<F, N, Id=BT::Data>, +where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis<F, BT::Agg, N>, X : for<'a> Apply<&'a BTFN<F, G, BT, N>, F> { type Codomain = F; @@ -512,16 +512,16 @@ fn p2_minimise<G: Fn(&U) -> F>(&self, g: G) -> (U, F); } -impl<F: Float> P2Minimise<Loc<F, 1>, F> for Cube<F, 1> { - fn p2_minimise<G: Fn(&Loc<F, 1>) -> F>(&self, g: G) -> (Loc<F, 1>, F) { +impl<F: Float> P2Minimise<Loc<1, F>, F> for Cube<1, F> { + fn p2_minimise<G: Fn(&Loc<1, F>) -> F>(&self, g: G) -> (Loc<1, F>, F) { let interval = Simplex(self.corners()); interval.p2_model(&g).minimise(&interval) } } #[replace_float_literals(F::cast_from(literal))] -impl<F: Float> P2Minimise<Loc<F, 2>, F> for Cube<F, 2> { - fn p2_minimise<G: Fn(&Loc<F, 2>) -> F>(&self, g: G) -> (Loc<F, 2>, F) { +impl<F: Float> P2Minimise<Loc<2, F>, F> for Cube<2, F> { + fn p2_minimise<G: Fn(&Loc<2, F>) -> F>(&self, g: G) -> (Loc<2, F>, F) { if false { // Split into two triangle (simplex) with separate P2 model in each. // The six nodes of each triangle are the corners and the edges. @@ -606,17 +606,17 @@ impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for P2Refiner<F, RefineMax> where - Cube<F, N>: P2Minimise<Loc<F, N>, F>, - G: SupportGenerator<F, N>, - G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, + Cube<N, F>: P2Minimise<Loc<N, F>, F>, + G: SupportGenerator<N, F>, + G::SupportType: Mapping<Loc<N, F>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, { - type Result = Option<(Loc<F, N>, F)>; + type Result = Option<(Loc<N, F>, F)>; type Sorting = UpperBoundSorting<F>; fn refine( &self, aggregator: &Bounds<F>, - cube: &Cube<F, N>, + cube: &Cube<N, F>, data: &[G::Id], generator: &G, step: usize, @@ -630,7 +630,7 @@ } // g gives the negative of the value of the function presented by `data` and `generator`. - let g = move |x: &Loc<F, N>| { + let g = move |x: &Loc<N, F>| { let f = move |&d| generator.support_for(d).apply(x); -data.iter().map(f).sum::<F>() }; @@ -669,17 +669,17 @@ impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for P2Refiner<F, RefineMin> where - Cube<F, N>: P2Minimise<Loc<F, N>, F>, - G: SupportGenerator<F, N>, - G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, + Cube<N, F>: P2Minimise<Loc<N, F>, F>, + G: SupportGenerator<N, F>, + G::SupportType: Mapping<Loc<N, F>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, { - type Result = Option<(Loc<F, N>, F)>; + type Result = Option<(Loc<N, F>, F)>; type Sorting = LowerBoundSorting<F>; fn refine( &self, aggregator: &Bounds<F>, - cube: &Cube<F, N>, + cube: &Cube<N, F>, data: &[G::Id], generator: &G, step: usize, @@ -693,7 +693,7 @@ } // g gives the value of the function presented by `data` and `generator`. - let g = move |x: &Loc<F, N>| { + let g = move |x: &Loc<N, F>| { let f = move |&d| generator.support_for(d).apply(x); data.iter().map(f).sum::<F>() }; @@ -755,7 +755,7 @@ impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for BoundRefiner<F, RefineMax> where - G: SupportGenerator<F, N>, + G: SupportGenerator<N, F>, { type Result = bool; type Sorting = UpperBoundSorting<F>; @@ -763,7 +763,7 @@ fn refine( &self, aggregator: &Bounds<F>, - _cube: &Cube<F, N>, + _cube: &Cube<N, F>, _data: &[G::Id], _generator: &G, step: usize, @@ -790,7 +790,7 @@ impl<F: Float, G, const N: usize> Refiner<F, Bounds<F>, G, N> for BoundRefiner<F, RefineMin> where - G: SupportGenerator<F, N>, + G: SupportGenerator<N, F>, { type Result = bool; type Sorting = UpperBoundSorting<F>; @@ -798,7 +798,7 @@ fn refine( &self, aggregator: &Bounds<F>, - _cube: &Cube<F, N>, + _cube: &Cube<N, F>, _data: &[G::Id], _generator: &G, step: usize, @@ -838,14 +838,14 @@ // there should be a result, or new nodes above the `glb` inserted into the queue. Then the waiting // threads can also continue processing. If, however, numerical inaccuracy destroyes the `glb`, // the queue may run out, and we get “Refiner failure”. -impl<F: Float, G, BT, const N: usize> MinMaxMapping<F, N> for BTFN<F, G, BT, N> +impl<F: Float, G, BT, const N: usize> MinMaxMapping<N, F> for BTFN<F, G, BT, N> where - BT: BTSearch<F, N, Agg = Bounds<F>>, - G: SupportGenerator<F, N, Id = BT::Data>, - G::SupportType: Mapping<Loc<F, N>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, - Cube<F, N>: P2Minimise<Loc<F, N>, F>, + BT: BTSearch<N, F, Agg = Bounds<F>>, + G: SupportGenerator<N, F, Id = BT::Data>, + G::SupportType: Mapping<Loc<N, F>, Codomain = F> + LocalAnalysis<F, Bounds<F>, N>, + Cube<N, F>: P2Minimise<Loc<N, F>, F>, { - fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F) { + fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc<N, F>, F) { let refiner = P2Refiner { tolerance, max_steps, @@ -863,7 +863,7 @@ bound: F, tolerance: F, max_steps: usize, - ) -> Option<(Loc<F, N>, F)> { + ) -> Option<(Loc<N, F>, F)> { let refiner = P2Refiner { tolerance, max_steps, @@ -875,7 +875,7 @@ .expect("Refiner failure.") } - fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F) { + fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc<N, F>, F) { let refiner = P2Refiner { tolerance, max_steps, @@ -893,7 +893,7 @@ bound: F, tolerance: F, max_steps: usize, - ) -> Option<(Loc<F, N>, F)> { + ) -> Option<(Loc<N, F>, F)> { let refiner = P2Refiner { tolerance, max_steps,
--- a/src/bisection_tree/either.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/bisection_tree/either.rs Thu May 01 13:06:58 2025 -0500 @@ -1,36 +1,26 @@ - use std::iter::Chain; use std::sync::Arc; +use crate::iter::{MapF, MapZ, Mappable}; +use crate::loc::Loc; +use crate::mapping::{DifferentiableImpl, DifferentiableMapping, Instance, Mapping, Space}; +use crate::sets::Cube; use crate::types::*; -use crate::mapping::{ - Instance, - Mapping, - DifferentiableImpl, - DifferentiableMapping, - Space, -}; -use crate::iter::{Mappable, MapF, MapZ}; -use crate::sets::Cube; -use crate::loc::Loc; +use super::aggregator::*; use super::support::*; -use super::aggregator::*; /// 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) Arc<A>, - pub(super) Arc<B>, -); +#[derive(Debug, Clone)] +pub struct BothGenerators<A, B>(pub(super) Arc<A>, pub(super) Arc<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> { +#[derive(Debug, Clone)] +pub enum EitherSupport<B, A> { Left(A), Right(B), } @@ -38,46 +28,55 @@ // We need type alias bounds to access associate types. #[allow(type_alias_bounds)] type BothAllDataIter< - 'a, F, - G1 : SupportGenerator<F, N>, - G2 : SupportGenerator<F, N>, - const N : usize + 'a, + F, + G1: SupportGenerator<N, F>, + G2: SupportGenerator<N, F>, + const N: usize, > = Chain< - MapF<G1::AllDataIter<'a>, (usize, EitherSupport<G1::SupportType, G2::SupportType>)>, - MapZ<G2::AllDataIter<'a>, usize, (usize, EitherSupport<G1::SupportType, G2::SupportType>)>, + MapF<G1::AllDataIter<'a>, (usize, EitherSupport<G2::SupportType, G1::SupportType>)>, + MapZ<G2::AllDataIter<'a>, usize, (usize, EitherSupport<G2::SupportType, G1::SupportType>)>, >; impl<G1, G2> BothGenerators<G1, G2> { /// 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>) - where G1 : SupportGenerator<F, N, Id=usize>, - G2 : SupportGenerator<F, N, Id=usize> { - - let id : usize = d.into(); + fn map_left<F: Float, const N: usize>( + (d, support): (G1::Id, G1::SupportType), + ) -> (usize, EitherSupport<G2::SupportType, G1::SupportType>) + where + G1: SupportGenerator<N, F, Id = usize>, + G2: SupportGenerator<N, F, Id = usize>, + { + let id: usize = d.into(); (id.into(), EitherSupport::Left(support)) } /// 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>) - where G1 : SupportGenerator<F, N, Id=usize>, - G2 : SupportGenerator<F, N, Id=usize> { - - let id : usize = d.into(); - ((n0+id).into(), EitherSupport::Right(support)) + fn map_right<F: Float, const N: usize>( + n0: &usize, + (d, support): (G2::Id, G2::SupportType), + ) -> (usize, EitherSupport<G2::SupportType, G1::SupportType>) + where + G1: SupportGenerator<N, F, Id = usize>, + G2: SupportGenerator<N, F, Id = usize>, + { + let id: usize = d.into(); + ((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>)> - where G1 : SupportGenerator<F, N, Id=usize>, - G2 : SupportGenerator<F, N, Id=usize> { + pub(super) fn all_left_data<F: Float, const N: usize>( + &self, + ) -> MapF<G1::AllDataIter<'_>, (usize, EitherSupport<G2::SupportType, G1::SupportType>)> + where + G1: SupportGenerator<N, F, Id = usize>, + G2: SupportGenerator<N, F, Id = usize>, + { self.0.all_data().mapF(Self::map_left) } @@ -85,33 +84,38 @@ /// /// 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>)> - where G1 : SupportGenerator<F, N, Id=usize>, - G2 : SupportGenerator<F, N, Id=usize> { + pub(super) fn all_right_data<F: Float, const N: usize>( + &self, + ) -> MapZ<G2::AllDataIter<'_>, usize, (usize, EitherSupport<G2::SupportType, G1::SupportType>)> + where + G1: SupportGenerator<N, F, Id = usize>, + G2: SupportGenerator<N, F, Id = usize>, + { let n0 = self.0.support_count(); self.1.all_data().mapZ(n0, Self::map_right) } } -impl<F : Float, G1, G2, const N : usize> -SupportGenerator<F, N> -for BothGenerators<G1, G2> -where G1 : SupportGenerator<F, N, Id=usize>, - G2 : SupportGenerator<F, N, Id=usize> { - +impl<F: Float, G1, G2, const N: usize> SupportGenerator<N, F> for BothGenerators<G1, G2> +where + G1: SupportGenerator<N, F, Id = usize>, + G2: SupportGenerator<N, F, Id = usize>, +{ type Id = usize; - type SupportType = EitherSupport<G1::SupportType, G2::SupportType>; - type AllDataIter<'a> = BothAllDataIter<'a, F, G1, G2, N> where G1 : 'a, G2 : 'a; + type SupportType = EitherSupport<G2::SupportType, G1::SupportType>; + type AllDataIter<'a> + = BothAllDataIter<'a, F, G1, G2, N> + where + G1: 'a, + G2: 'a; #[inline] - fn support_for(&self, id : Self::Id) - -> Self::SupportType { + fn support_for(&self, id: Self::Id) -> Self::SupportType { let n0 = self.0.support_count(); if id < n0 { EitherSupport::Left(self.0.support_for(id.into())) } else { - EitherSupport::Right(self.1.support_for((id-n0).into())) + EitherSupport::Right(self.1.support_for((id - n0).into())) } } @@ -126,12 +130,13 @@ } } -impl<F: Float, S1, S2, const N : usize> Support<F, N> for EitherSupport<S1, S2> -where S1 : Support<F, N>, - S2 : Support<F, N> { - +impl<F: Float, S1, S2, const N: usize> Support<N, F> for EitherSupport<S2, S1> +where + S1: Support<N, F>, + S2: Support<N, F>, +{ #[inline] - fn support_hint(&self) -> Cube<F,N> { + fn support_hint(&self) -> Cube<N, F> { match self { EitherSupport::Left(ref a) => a.support_hint(), EitherSupport::Right(ref b) => b.support_hint(), @@ -139,7 +144,7 @@ } #[inline] - fn in_support(&self, x : &Loc<F,N>) -> bool { + fn in_support(&self, x: &Loc<N, F>) -> bool { match self { EitherSupport::Left(ref a) => a.in_support(x), EitherSupport::Right(ref b) => b.in_support(x), @@ -147,7 +152,7 @@ } #[inline] - fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] { + fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] { match self { EitherSupport::Left(ref a) => a.bisection_hint(cube), EitherSupport::Right(ref b) => b.bisection_hint(cube), @@ -155,13 +160,14 @@ } } -impl<F : Float, A, S1, S2, const N : usize> LocalAnalysis<F, A, N> for EitherSupport<S1, S2> -where A : Aggregator, - S1 : LocalAnalysis<F, A, N>, - S2 : LocalAnalysis<F, A, N>, { - +impl<F: Float, A, S1, S2, const N: usize> LocalAnalysis<F, A, N> for EitherSupport<S2, S1> +where + A: Aggregator, + S1: LocalAnalysis<F, A, N>, + S2: LocalAnalysis<F, A, N>, +{ #[inline] - fn local_analysis(&self, cube : &Cube<F, N>) -> A { + fn local_analysis(&self, cube: &Cube<N, F>) -> A { match self { EitherSupport::Left(ref a) => a.local_analysis(cube), EitherSupport::Right(ref b) => b.local_analysis(cube), @@ -169,11 +175,12 @@ } } -impl<F : Float, A, S1, S2> GlobalAnalysis<F, A> for EitherSupport<S1, S2> -where A : Aggregator, - S1 : GlobalAnalysis<F, A>, - S2 : GlobalAnalysis<F, A>, { - +impl<F: Float, A, S1, S2> GlobalAnalysis<F, A> for EitherSupport<S2, S1> +where + A: Aggregator, + S1: GlobalAnalysis<F, A>, + S2: GlobalAnalysis<F, A>, +{ #[inline] fn global_analysis(&self) -> A { match self { @@ -183,17 +190,17 @@ } } -impl<F, S1, S2, X> Mapping<X> for EitherSupport<S1, S2> +impl<F, S1, S2, X> Mapping<X> for EitherSupport<S2, S1> where - F : Space, - X : Space, - S1 : Mapping<X, Codomain=F>, - S2 : Mapping<X, Codomain=F>, + F: Space, + X: Space, + S1: Mapping<X, Codomain = F>, + S2: Mapping<X, Codomain = F>, { type Codomain = F; #[inline] - fn apply<I : Instance<X>>(&self, x : I) -> F { + fn apply<I: Instance<X>>(&self, x: I) -> F { match self { EitherSupport::Left(ref a) => a.apply(x), EitherSupport::Right(ref b) => b.apply(x), @@ -201,17 +208,17 @@ } } -impl<X, S1, S2, O> DifferentiableImpl<X> for EitherSupport<S1, S2> +impl<X, S1, S2, O> DifferentiableImpl<X> for EitherSupport<S2, S1> where - O : Space, - X : Space, - S1 : DifferentiableMapping<X, DerivativeDomain=O>, - S2 : DifferentiableMapping<X, DerivativeDomain=O>, + O: Space, + X: Space, + S1: DifferentiableMapping<X, DerivativeDomain = O>, + S2: DifferentiableMapping<X, DerivativeDomain = O>, { type Derivative = O; #[inline] - fn differential_impl<I : Instance<X>>(&self, x : I) -> O { + fn differential_impl<I: Instance<X>>(&self, x: I) -> O { match self { EitherSupport::Left(ref a) => a.differential(x), EitherSupport::Right(ref b) => b.differential(x), @@ -221,44 +228,47 @@ macro_rules! make_either_scalarop_rhs { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { - impl<F : Float, G1, G2> - std::ops::$trait_assign<F> - for BothGenerators<G1, G2> - where G1 : std::ops::$trait_assign<F> + Clone, - G2 : std::ops::$trait_assign<F> + Clone, { + impl<F: Float, G1, G2> std::ops::$trait_assign<F> for BothGenerators<G1, G2> + where + G1: std::ops::$trait_assign<F> + Clone, + G2: std::ops::$trait_assign<F> + Clone, + { #[inline] - fn $fn_assign(&mut self, t : F) { + fn $fn_assign(&mut self, t: F) { Arc::make_mut(&mut self.0).$fn_assign(t); Arc::make_mut(&mut self.1).$fn_assign(t); } } - impl<'a, F : Float, G1, G2> - std::ops::$trait<F> - for &'a BothGenerators<G1, G2> - where &'a G1 : std::ops::$trait<F,Output=G1>, - &'a G2 : std::ops::$trait<F,Output=G2> { + impl<'a, F: Float, G1, G2> std::ops::$trait<F> for &'a BothGenerators<G1, G2> + where + &'a G1: std::ops::$trait<F, Output = G1>, + &'a G2: std::ops::$trait<F, Output = G2>, + { type Output = BothGenerators<G1, G2>; #[inline] - fn $fn(self, t : F) -> BothGenerators<G1, G2> { - BothGenerators(Arc::new(self.0.$fn(t)), - Arc::new(self.1.$fn(t))) + fn $fn(self, t: F) -> BothGenerators<G1, G2> { + BothGenerators(Arc::new(self.0.$fn(t)), Arc::new(self.1.$fn(t))) } } - } + }; } make_either_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); make_either_scalarop_rhs!(Div, div, DivAssign, div_assign); impl<G1, G2> std::ops::Neg for BothGenerators<G1, G2> -where G1 : std::ops::Neg + Clone, - G2 : std::ops::Neg + Clone, { +where + G1: std::ops::Neg + Clone, + G2: std::ops::Neg + Clone, +{ type Output = BothGenerators<G1::Output, G2::Output>; #[inline] fn neg(self) -> Self::Output { - BothGenerators(Arc::new(Arc::unwrap_or_clone(self.0).neg()), - Arc::new(Arc::unwrap_or_clone(self.1).neg())) + BothGenerators( + Arc::new(Arc::unwrap_or_clone(self.0).neg()), + Arc::new(Arc::unwrap_or_clone(self.1).neg()), + ) } } /* @@ -270,4 +280,4 @@ BothGenerators(self.0.neg(), self.1.neg()) } } -*/ \ No newline at end of file +*/
--- a/src/bisection_tree/refine.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/bisection_tree/refine.rs Thu May 01 13:06:58 2025 -0500 @@ -1,32 +1,31 @@ - -use std::collections::BinaryHeap; -use std::cmp::{PartialOrd, Ord, Ordering, Ordering::*, max}; -use std::marker::PhantomData; -use std::sync::{Arc, Mutex, MutexGuard, Condvar}; -use crate::types::*; +use super::aggregator::*; +use super::bt::*; +use super::support::*; use crate::nanleast::NaNLeast; +use crate::parallelism::TaskBudget; +use crate::parallelism::{thread_pool, thread_pool_size}; use crate::sets::Cube; -use crate::parallelism::{thread_pool_size, thread_pool}; -use super::support::*; -use super::bt::*; -use super::aggregator::*; -use crate::parallelism::TaskBudget; +use crate::types::*; +use std::cmp::{max, Ord, Ordering, Ordering::*, PartialOrd}; +use std::collections::BinaryHeap; +use std::marker::PhantomData; +use std::sync::{Arc, Condvar, Mutex, MutexGuard}; /// Trait for sorting [`Aggregator`]s for [`BT`] refinement. /// /// 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 : Sync + Send + 'static { +pub trait AggregatorSorting: Sync + Send + 'static { // Priority - type Agg : Aggregator; - type Sort : Ord + Copy + std::fmt::Debug + Sync + Send; + type Agg: Aggregator; + type Sort: Ord + Copy + std::fmt::Debug + Sync + Send; /// Returns lower sorting key - fn sort_lower(aggregator : &Self::Agg) -> Self::Sort; + fn sort_lower(aggregator: &Self::Agg) -> Self::Sort; /// Returns upper sorting key - fn sort_upper(aggregator : &Self::Agg) -> Self::Sort; + fn sort_upper(aggregator: &Self::Agg) -> Self::Sort; /// Returns a sorting key that is less than any other sorting key. fn bottom() -> Self::Sort; @@ -35,53 +34,64 @@ /// An [`AggregatorSorting`] for [`Bounds`], using the upper/lower bound as the upper/lower key. /// /// See [`LowerBoundSorting`] for the opposite ordering. -pub struct UpperBoundSorting<F : Float>(PhantomData<F>); +pub struct UpperBoundSorting<F: Float>(PhantomData<F>); /// An [`AggregatorSorting`] for [`Bounds`], using the upper/lower bound as the lower/upper key. /// /// See [`UpperBoundSorting`] for the opposite ordering. -pub struct LowerBoundSorting<F : Float>(PhantomData<F>); +pub struct LowerBoundSorting<F: Float>(PhantomData<F>); -impl<F : Float> AggregatorSorting for UpperBoundSorting<F> { +impl<F: Float> AggregatorSorting for UpperBoundSorting<F> { type Agg = Bounds<F>; type Sort = NaNLeast<F>; #[inline] - fn sort_lower(aggregator : &Bounds<F>) -> Self::Sort { NaNLeast(aggregator.lower()) } - - #[inline] - fn sort_upper(aggregator : &Bounds<F>) -> Self::Sort { NaNLeast(aggregator.upper()) } + fn sort_lower(aggregator: &Bounds<F>) -> Self::Sort { + NaNLeast(aggregator.lower()) + } #[inline] - fn bottom() -> Self::Sort { NaNLeast(F::NEG_INFINITY) } + fn sort_upper(aggregator: &Bounds<F>) -> Self::Sort { + NaNLeast(aggregator.upper()) + } + + #[inline] + fn bottom() -> Self::Sort { + NaNLeast(F::NEG_INFINITY) + } } - -impl<F : Float> AggregatorSorting for LowerBoundSorting<F> { +impl<F: Float> AggregatorSorting for LowerBoundSorting<F> { type Agg = Bounds<F>; type Sort = NaNLeast<F>; #[inline] - fn sort_upper(aggregator : &Bounds<F>) -> Self::Sort { NaNLeast(-aggregator.lower()) } + fn sort_upper(aggregator: &Bounds<F>) -> Self::Sort { + NaNLeast(-aggregator.lower()) + } #[inline] - fn sort_lower(aggregator : &Bounds<F>) -> Self::Sort { NaNLeast(-aggregator.upper()) } + fn sort_lower(aggregator: &Bounds<F>) -> Self::Sort { + NaNLeast(-aggregator.upper()) + } #[inline] - fn bottom() -> Self::Sort { NaNLeast(F::NEG_INFINITY) } + fn bottom() -> Self::Sort { + NaNLeast(F::NEG_INFINITY) + } } /// 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> { +pub enum RefinerResult<A: Aggregator, R> { /// Indicates an insufficiently refined state: the [`BT`] needs to be further refined. NeedRefinement, /// Indicates a certain result `R`, stop refinement immediately. Certain(R), /// Indicates an uncertain result: continue refinement until candidates have been exhausted /// or a certain result found. - Uncertain(A, R) + Uncertain(A, R), } use RefinerResult::*; @@ -92,16 +102,17 @@ /// 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> : Sync + Send + 'static -where F : Num, - A : Aggregator, - G : SupportGenerator<F, N> { - +pub trait Refiner<F: Float, A, G, const N: usize>: Sync + Send + 'static +where + F: Num, + A: Aggregator, + G: SupportGenerator<N, F>, +{ /// The result type of the refiner - type Result : std::fmt::Debug + Sync + Send + 'static; + type Result: std::fmt::Debug + Sync + Send + 'static; /// The sorting to be employed by [`BTSearch::search_and_refine`] on node aggregators /// to detemrine node priority. - type Sorting : AggregatorSorting<Agg = A>; + type Sorting: AggregatorSorting<Agg = A>; /// Determines whether `aggregator` is sufficiently refined within `domain`. /// @@ -124,42 +135,45 @@ /// number of steps is reached. fn refine( &self, - aggregator : &A, - domain : &Cube<F, N>, - data : &[G::Id], - generator : &G, - step : usize, + aggregator: &A, + domain: &Cube<N, F>, + data: &[G::Id], + generator: &G, + step: usize, ) -> RefinerResult<A, Self::Result>; /// Fuse two [`Self::Result`]s (needed in threaded refinement). - fn fuse_results(r1 : &mut Self::Result, r2 : Self::Result); + fn fuse_results(r1: &mut Self::Result, r2: Self::Result); } /// Structure for tracking the refinement process in a [`BinaryHeap`]. -struct RefinementInfo<'a, F, D, A, S, RResult, const N : usize, const P : usize> -where F : Float, - D : 'static, - A : Aggregator, - S : AggregatorSorting<Agg = A> { +struct RefinementInfo<'a, F, D, A, S, RResult, const N: usize, const P: usize> +where + F: Float, + D: 'static, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ /// Domain of `node` - cube : Cube<F, N>, + cube: Cube<N, F>, /// Node to be refined - node : &'a mut Node<F, D, A, N, P>, + node: &'a mut Node<F, D, A, N, P>, /// Result and improve aggregator for the [`Refiner`] - refiner_info : Option<(A, RResult)>, + refiner_info: Option<(A, RResult)>, /// For [`AggregatorSorting`] being used for the type system - sorting : PhantomData<S>, + sorting: PhantomData<S>, } -impl<'a, F, D, A, S, RResult, const N : usize, const P : usize> -RefinementInfo<'a, F, D, A, S, RResult, N, P> -where F : Float, - D : 'static, - A : Aggregator, - S : AggregatorSorting<Agg = A> { - +impl<'a, F, D, A, S, RResult, const N: usize, const P: usize> + RefinementInfo<'a, F, D, A, S, RResult, N, P> +where + F: Float, + D: 'static, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ #[inline] - fn with_aggregator<U>(&self, f : impl FnOnce(&A) -> U) -> U { + fn with_aggregator<U>(&self, f: impl FnOnce(&A) -> U) -> U { match self.refiner_info { Some((ref agg, _)) => f(agg), None => f(&self.node.aggregator), @@ -177,93 +191,105 @@ } } -impl<'a, F, D, A, S, RResult, const N : usize, const P : usize> PartialEq -for RefinementInfo<'a, F, D, A, S, RResult, N, P> -where F : Float, - D : 'static, - A : Aggregator, - S : AggregatorSorting<Agg = A> { - +impl<'a, F, D, A, S, RResult, const N: usize, const P: usize> PartialEq + for RefinementInfo<'a, F, D, A, S, RResult, N, P> +where + F: Float, + D: 'static, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ #[inline] - fn eq(&self, other : &Self) -> bool { self.cmp(other) == Equal } -} - -impl<'a, F, D, A, S, RResult, const N : usize, const P : usize> PartialOrd -for RefinementInfo<'a, F, D, A, S, RResult, N, P> -where F : Float, - D : 'static, - A : Aggregator, - S : AggregatorSorting<Agg = A> { - - #[inline] - fn partial_cmp(&self, other : &Self) -> Option<Ordering> { Some(self.cmp(other)) } + fn eq(&self, other: &Self) -> bool { + self.cmp(other) == Equal + } } -impl<'a, F, D, A, S, RResult, const N : usize, const P : usize> Eq -for RefinementInfo<'a, F, D, A, S, RResult, N, P> -where F : Float, - D : 'static, - A : Aggregator, - S : AggregatorSorting<Agg = A> { +impl<'a, F, D, A, S, RResult, const N: usize, const P: usize> PartialOrd + for RefinementInfo<'a, F, D, A, S, RResult, N, P> +where + F: Float, + D: 'static, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ + #[inline] + fn partial_cmp(&self, other: &Self) -> Option<Ordering> { + Some(self.cmp(other)) + } } -impl<'a, F, D, A, S, RResult, const N : usize, const P : usize> Ord -for RefinementInfo<'a, F, D, A, S, RResult, N, P> -where F : Float, - D : 'static, - A : Aggregator, - S : AggregatorSorting<Agg = A> { +impl<'a, F, D, A, S, RResult, const N: usize, const P: usize> Eq + for RefinementInfo<'a, F, D, A, S, RResult, N, P> +where + F: Float, + D: 'static, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ +} +impl<'a, F, D, A, S, RResult, const N: usize, const P: usize> Ord + for RefinementInfo<'a, F, D, A, S, RResult, N, P> +where + F: Float, + D: 'static, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ #[inline] - fn cmp(&self, other : &Self) -> Ordering { - self.with_aggregator(|agg1| other.with_aggregator(|agg2| { - match S::sort_upper(agg1).cmp(&S::sort_upper(agg2)) { + fn cmp(&self, other: &Self) -> Ordering { + self.with_aggregator(|agg1| { + other.with_aggregator(|agg2| match S::sort_upper(agg1).cmp(&S::sort_upper(agg2)) { Equal => S::sort_lower(agg1).cmp(&S::sort_lower(agg2)), order => order, - } - })) + }) + }) } } /// 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>, - A : Aggregator, - S : AggregatorSorting<Agg = A> { +struct HeapContainer<'a, F, D, A, S, RResult, const N: usize, const P: usize> +where + F: Float, + D: 'static + Copy, + Const<P>: BranchCount<N>, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ /// Priority queue of nodes to be refined - heap : BinaryHeap<RefinementInfo<'a, F, D, A, S, RResult, N, P>>, + heap: BinaryHeap<RefinementInfo<'a, F, D, A, S, RResult, N, P>>, /// Maximum of node sorting lower bounds seen in the heap - glb : S::Sort, + glb: S::Sort, /// Number of insertions in the heap since previous prune - insert_counter : usize, + insert_counter: usize, /// If a result has been found by some refinment threat, it is stored here - result : Option<RResult>, + result: Option<RResult>, /// Refinement step counter - step : usize, + step: usize, /// Number of threads currently processing (not sleeping) - n_processing : usize, + n_processing: usize, /// Threshold for heap pruning - heap_prune_threshold : usize, + heap_prune_threshold: usize, } -impl<'a, F, D, A, S, RResult, const N : usize, const P : usize> -HeapContainer<'a, F, D, A, S, RResult, N, P> -where F : Float, - D : 'static + Copy, - Const<P> : BranchCount<N>, - A : Aggregator, - S : AggregatorSorting<Agg = A> { - +impl<'a, F, D, A, S, RResult, const N: usize, const P: usize> + HeapContainer<'a, F, D, A, S, RResult, N, P> +where + F: Float, + D: 'static + Copy, + Const<P>: BranchCount<N>, + A: Aggregator, + S: AggregatorSorting<Agg = A>, +{ /// Push `ri` into the [`BinaryHeap`]. Do greatest lower bound maintenance. /// /// Returns a boolean indicating whether the push was actually performed due to glb /// filtering or not. #[inline] - fn push(&mut self, ri : RefinementInfo<'a, F, D, A, S, RResult, N, P>) -> bool { + fn push(&mut self, ri: RefinementInfo<'a, F, D, A, S, RResult, N, P>) -> bool { if ri.sort_upper() >= self.glb { let l = ri.sort_lower(); self.heap.push(ri); @@ -276,37 +302,38 @@ } } -impl<F : Float, D, A, const N : usize, const P : usize> -Branches<F,D,A,N,P> -where Const<P> : BranchCount<N>, - A : Aggregator, - D : 'static + Copy + Send + Sync { - +impl<F: Float, D, A, const N: usize, const P: usize> Branches<F, D, A, N, P> +where + Const<P>: BranchCount<N>, + A: Aggregator, + D: 'static + Copy + Send + Sync, +{ /// Stage all subnodes of `self` into the refinement queue `container`. fn stage_refine<'a, S, RResult>( &'a mut self, - domain : Cube<F,N>, - container : &mut HeapContainer<'a, F, D, A, S, RResult, N, P>, - ) where S : AggregatorSorting<Agg = A> { + domain: Cube<N, F>, + container: &mut HeapContainer<'a, F, D, A, S, RResult, N, P>, + ) where + S: AggregatorSorting<Agg = A>, + { // Insert all subnodes into the refinement heap. for (node, cube) in self.nodes_and_cubes_mut(&domain) { container.push(RefinementInfo { cube, node, - refiner_info : None, - sorting : PhantomData, + refiner_info: None, + sorting: PhantomData, }); } } } - -impl<F : Float, D, A, const N : usize, const P : usize> -Node<F,D,A,N,P> -where Const<P> : BranchCount<N>, - A : Aggregator, - D : 'static + Copy + Send + Sync { - +impl<F: Float, D, A, const N: usize, const P: usize> Node<F, D, A, N, P> +where + Const<P>: BranchCount<N>, + A: Aggregator, + D: 'static + Copy + Send + Sync, +{ /// If `self` is a leaf node, uses the `refiner` to determine whether further subdivision /// is required to get a sufficiently refined solution for the problem the refiner is used /// to solve. If the refiner returns [`RefinerResult::Certain`] result, it is returned. @@ -316,17 +343,18 @@ /// /// `domain`, as usual, indicates the spatial area corresponding to `self`. fn search_and_refine<'a, 'b, 'c, R, G>( - self : &'a mut Self, - domain : Cube<F,N>, - refiner : &R, - generator : &G, - container_arc : &'c Arc<Mutex<HeapContainer<'a, F, D, A, R::Sorting, R::Result, N, P>>>, - step : usize + self: &'a mut Self, + domain: Cube<N, F>, + refiner: &R, + generator: &G, + container_arc: &'c Arc<Mutex<HeapContainer<'a, F, D, A, R::Sorting, R::Result, N, P>>>, + step: usize, ) -> Result<R::Result, MutexGuard<'c, HeapContainer<'a, F, D, A, R::Sorting, R::Result, N, P>>> - where R : Refiner<F, A, G, N>, - G : SupportGenerator<F, N, Id=D>, - G::SupportType : LocalAnalysis<F, A, N> { - + where + R: Refiner<F, A, G, N>, + G: SupportGenerator<N, F, Id = D>, + G::SupportType: LocalAnalysis<F, A, N>, + { //drop(container); // Refine a leaf. @@ -368,26 +396,27 @@ unsafe { Arc::get_mut_unchecked(arc_b) } .stage_refine(domain, &mut *container); #[cfg(not(nightly))] - Arc::get_mut(arc_b).unwrap() + Arc::get_mut(arc_b) + .unwrap() .stage_refine(domain, &mut *container); - - return Err(container) - }, + + return Err(container); + } _ => unreachable!("This cannot happen"), } } } res - }, + } NodeOption::Branches(ref mut b) => { // Insert branches into refinement priority queue. let mut container = container_arc.lock().unwrap(); Arc::make_mut(b).stage_refine(domain, &mut *container); - return Err(container) - }, + return Err(container); + } NodeOption::Uninitialised => { refiner.refine(&self.aggregator, &domain, &[], generator, step) - }, + } }; match res { @@ -399,17 +428,17 @@ // aggregator. let mut container = container_arc.lock().unwrap(); container.push(RefinementInfo { - cube : domain, - node : self, - refiner_info : Some((agg, val)), - sorting : PhantomData, + cube: domain, + node: self, + refiner_info: Some((agg, val)), + sorting: PhantomData, }); Err(container) - }, + } Certain(val) => { // The refiner gave a certain result so return it to allow early termination Ok(val) - }, + } NeedRefinement => { // This should only happen when we run into NodeOption::Uninitialised above. // There's really nothing to do. @@ -423,9 +452,10 @@ /// /// 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 { - +pub trait BTSearch<const N: usize, F = f64>: BTImpl<N, F> +where + F: Float, +{ /// Perform a search on on `Self`, as determined by `refiner`. /// /// Nodes are inserted in a priority queue and processed in the order determined by the @@ -437,26 +467,28 @@ /// 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, - generator : &Arc<G>, + refiner: R, + generator: &Arc<G>, ) -> Option<R::Result> - where R : Refiner<F, Self::Agg, G, N> + Sync + Send + 'static, - G : SupportGenerator<F, N, Id=Self::Data> + Sync + Send + 'static, - G::SupportType : LocalAnalysis<F, Self::Agg, N>; + where + R: Refiner<F, Self::Agg, G, N> + Sync + Send + 'static, + G: SupportGenerator<N, F, Id = Self::Data> + Sync + Send + 'static, + G::SupportType: LocalAnalysis<F, Self::Agg, N>; } -fn refinement_loop<F : Float, D, A, R, G, const N : usize, const P : usize> ( - wakeup : Option<Arc<Condvar>>, - refiner : &R, - generator_arc : &Arc<G>, - container_arc : &Arc<Mutex<HeapContainer<F, D, A, R::Sorting, R::Result, N, P>>>, -) where A : Aggregator, - R : Refiner<F, A, G, N>, - G : SupportGenerator<F, N, Id=D>, - G::SupportType : LocalAnalysis<F, A, N>, - Const<P> : BranchCount<N>, - D : 'static + Copy + Sync + Send + std::fmt::Debug { - +fn refinement_loop<F: Float, D, A, R, G, const N: usize, const P: usize>( + wakeup: Option<Arc<Condvar>>, + refiner: &R, + generator_arc: &Arc<G>, + container_arc: &Arc<Mutex<HeapContainer<F, D, A, R::Sorting, R::Result, N, P>>>, +) where + A: Aggregator, + R: Refiner<F, A, G, N>, + G: SupportGenerator<N, F, Id = D>, + G::SupportType: LocalAnalysis<F, A, N>, + Const<P>: BranchCount<N>, + D: 'static + Copy + Sync + Send + std::fmt::Debug, +{ let mut did_park = true; let mut container = container_arc.lock().unwrap(); @@ -471,7 +503,7 @@ // Some refinement task/thread has found a result, return if container.result.is_some() { container.n_processing -= 1; - break 'main + break 'main; } match container.heap.pop() { @@ -489,7 +521,7 @@ container = c.wait(container).unwrap(); continue 'get_next; } else { - break 'main + break 'main; } } }; @@ -502,7 +534,7 @@ // Terminate based on a “best possible” result. container.result = Some(result); container.n_processing -= 1; - break 'main + break 'main; } // Do priority queue maintenance @@ -513,10 +545,8 @@ container.glb = glb; // Prune container.heap.retain(|ri| ri.sort_upper() >= glb); - }, - None => { - container.glb = R::Sorting::bottom() } + None => container.glb = R::Sorting::bottom(), } container.insert_counter = 0; } @@ -525,8 +555,14 @@ drop(container); // … and process the node. We may get returned an already unlocked mutex. - match Node::search_and_refine(ri.node, ri.cube, refiner, &**generator_arc, - &container_arc, step) { + match Node::search_and_refine( + ri.node, + ri.cube, + refiner, + &**generator_arc, + &container_arc, + step, + ) { Ok(r) => { let mut container = container_arc.lock().unwrap(); // Terminate based on a certain result from the refiner @@ -534,8 +570,8 @@ Some(ref mut r_prev) => R::fuse_results(r_prev, r), None => container.result = Some(r), } - break 'main - }, + break 'main; + } Err(cnt) => { container = cnt; // Wake up another thread if one is sleeping; there should be now work in the @@ -545,7 +581,6 @@ } } } - } // Make sure no task is sleeping @@ -558,9 +593,9 @@ macro_rules! impl_btsearch { ($($n:literal)*) => { $( impl<'a, M, F, D, A> - BTSearch<F, $n> + BTSearch<$n, F> for BT<M,F,D,A,$n> - where //Self : BTImpl<F,$n,Data=D,Agg=A, Depth=M>, // <== automatically deduced + where //Self : BTImpl<$n, F, Data=D,Agg=A, Depth=M>, // <== automatically deduced M : Depth, F : Float + Send, A : Aggregator, @@ -571,7 +606,7 @@ generator : &Arc<G>, ) -> Option<R::Result> where R : Refiner<F, A, G, $n>, - G : SupportGenerator<F, $n, Id=D>, + G : SupportGenerator< $n, F, Id=D>, G::SupportType : LocalAnalysis<F, A, $n> { let mut init_container = HeapContainer { heap : BinaryHeap::new(), @@ -620,4 +655,3 @@ } impl_btsearch!(1 2 3 4); -
--- a/src/bisection_tree/support.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/bisection_tree/support.rs Thu May 01 13:06:58 2025 -0500 @@ -16,14 +16,14 @@ /// A trait for working with the supports of [`Mapping`]s. /// /// `Mapping` is not a super-trait to allow more general use. -pub trait Support<F: Num, const N: usize>: Sized + Sync + Send + 'static { +pub trait Support<const N: usize, F: Num>: Sized + Sync + Send + 'static { /// Return a cube containing the support of the function represented by `self`. /// /// The hint may be larger than the actual support, but must contain it. - fn support_hint(&self) -> Cube<F, N>; + fn support_hint(&self) -> Cube<N, F>; /// Indicate whether `x` is in the support of the function represented by `self`. - fn in_support(&self, x: &Loc<F, N>) -> bool; + fn in_support(&self, x: &Loc<N, F>) -> bool; // Indicate whether `cube` is fully in the support of the function represented by `self`. //fn fully_in_support(&self, cube : &Cube<F,N>) -> bool; @@ -39,13 +39,13 @@ /// The default implementation returns `[None; N]`. #[inline] #[allow(unused_variables)] - fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] { + fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] { [None; N] } /// Translate `self` by `x`. #[inline] - fn shift(self, x: Loc<F, N>) -> Shift<Self, F, N> { + fn shift(self, x: Loc<N, F>) -> Shift<Self, N, F> { Shift { shift: x, base_fn: self, @@ -55,46 +55,46 @@ /// 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>, +pub struct Shift<T, const N: usize, F = f64> { + shift: Loc<N, F>, base_fn: T, } -impl<'a, T, V: Space, F: Float, const N: usize> Mapping<Loc<F, N>> for Shift<T, F, N> +impl<'a, T, V: Space, F: Float, const N: usize> Mapping<Loc<N, F>> for Shift<T, N, F> where - T: Mapping<Loc<F, N>, Codomain = V>, + T: Mapping<Loc<N, F>, Codomain = V>, { type Codomain = V; #[inline] - fn apply<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Codomain { + fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain { self.base_fn.apply(x.own() - &self.shift) } } -impl<'a, T, V: Space, F: Float, const N: usize> DifferentiableImpl<Loc<F, N>> for Shift<T, F, N> +impl<'a, T, V: Space, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for Shift<T, N, F> where - T: DifferentiableMapping<Loc<F, N>, DerivativeDomain = V>, + T: DifferentiableMapping<Loc<N, F>, DerivativeDomain = V>, { type Derivative = V; #[inline] - fn differential_impl<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Derivative { + fn differential_impl<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Derivative { self.base_fn.differential(x.own() - &self.shift) } } -impl<'a, T, F: Float, const N: usize> Support<F, N> for Shift<T, F, N> +impl<'a, T, F: Float, const N: usize> Support<N, F> for Shift<T, N, F> where - T: Support<F, N>, + T: Support<N, F>, { #[inline] - fn support_hint(&self) -> Cube<F, N> { + fn support_hint(&self) -> Cube<N, F> { self.base_fn.support_hint().shift(&self.shift) } #[inline] - fn in_support(&self, x: &Loc<F, N>) -> bool { + fn in_support(&self, x: &Loc<N, F>) -> bool { self.base_fn.in_support(&(x - &self.shift)) } @@ -104,13 +104,13 @@ // } #[inline] - fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] { + fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] { let base_hint = self.base_fn.bisection_hint(cube); map2(base_hint, &self.shift, |h, s| h.map(|z| z + *s)) } } -impl<'a, T, F: Float, const N: usize> GlobalAnalysis<F, Bounds<F>> for Shift<T, F, N> +impl<'a, T, F: Float, const N: usize> GlobalAnalysis<F, Bounds<F>> for Shift<T, N, F> where T: LocalAnalysis<F, Bounds<F>, N>, { @@ -120,20 +120,20 @@ } } -impl<'a, T, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Shift<T, F, N> +impl<'a, T, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Shift<T, N, F> where T: LocalAnalysis<F, Bounds<F>, N>, { #[inline] - fn local_analysis(&self, cube: &Cube<F, N>) -> Bounds<F> { + fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> { self.base_fn.local_analysis(&cube.shift(&(-self.shift))) } } macro_rules! impl_shift_norm { ($($norm:ident)*) => { $( - impl<'a, T, F : Float, const N : usize> Norm<F, $norm> for Shift<T,F,N> - where T : Norm<F, $norm> { + impl<'a, T, F : Float, const N : usize> Norm<$norm, F> for Shift<T, N, F> + where T : Norm<$norm, F> { #[inline] fn norm(&self, n : $norm) -> F { self.base_fn.norm(n) @@ -144,18 +144,18 @@ impl_shift_norm!(L1 L2 Linfinity); -impl<'a, T, F: Float, C, const N: usize> Support<F, N> for Weighted<T, C> +impl<'a, T, F: Float, C, const N: usize> Support<N, F> for Weighted<T, C> where - T: Support<F, N>, + T: Support<N, F>, C: Constant<Type = F>, { #[inline] - fn support_hint(&self) -> Cube<F, N> { + fn support_hint(&self) -> Cube<N, F> { self.base_fn.support_hint() } #[inline] - fn in_support(&self, x: &Loc<F, N>) -> bool { + fn in_support(&self, x: &Loc<N, F>) -> bool { self.base_fn.in_support(x) } @@ -164,7 +164,7 @@ // } #[inline] - fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] { + fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] { self.base_fn.bisection_hint(cube) } } @@ -191,7 +191,7 @@ C: Constant<Type = F>, { #[inline] - fn local_analysis(&self, cube: &Cube<F, N>) -> Bounds<F> { + fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> { let Bounds(lower, upper) = self.base_fn.local_analysis(cube); debug_assert!(lower <= upper); match self.weight.value() { @@ -240,8 +240,8 @@ macro_rules! impl_weighted_norm { ($($norm:ident)*) => { $( - impl<'a, T, F : Float> Norm<F, $norm> for Weighted<T,F> - where T : Norm<F, $norm> { + impl<'a, T, F : Float> Norm<$norm, F> for Weighted<T,F> + where T : Norm<$norm, F> { #[inline] fn norm(&self, n : $norm) -> F { self.base_fn.norm(n) * self.weight.abs() @@ -261,14 +261,14 @@ pub T, ); -impl<'a, T, F: Float, const N: usize> Mapping<Loc<F, N>> for Normalised<T> +impl<'a, T, F: Float, const N: usize> Mapping<Loc<N, F>> for Normalised<T> where - T: Norm<F, L1> + Mapping<Loc<F, N>, Codomain = F>, + T: Norm<L1, F> + Mapping<Loc<N, F>, Codomain = F>, { type Codomain = F; #[inline] - fn apply<I: Instance<Loc<F, N>>>(&self, x: I) -> Self::Codomain { + fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain { let w = self.0.norm(L1); if w == F::ZERO { F::ZERO @@ -278,17 +278,17 @@ } } -impl<'a, T, F: Float, const N: usize> Support<F, N> for Normalised<T> +impl<'a, T, F: Float, const N: usize> Support<N, F> for Normalised<T> where - T: Norm<F, L1> + Support<F, N>, + T: Norm<L1, F> + Support<N, F>, { #[inline] - fn support_hint(&self) -> Cube<F, N> { + fn support_hint(&self) -> Cube<N, F> { self.0.support_hint() } #[inline] - fn in_support(&self, x: &Loc<F, N>) -> bool { + fn in_support(&self, x: &Loc<N, F>) -> bool { self.0.in_support(x) } @@ -297,14 +297,14 @@ // } #[inline] - fn bisection_hint(&self, cube: &Cube<F, N>) -> [Option<F>; N] { + fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] { self.0.bisection_hint(cube) } } impl<'a, T, F: Float> GlobalAnalysis<F, Bounds<F>> for Normalised<T> where - T: Norm<F, L1> + GlobalAnalysis<F, Bounds<F>>, + T: Norm<L1, F> + GlobalAnalysis<F, Bounds<F>>, { #[inline] fn global_analysis(&self) -> Bounds<F> { @@ -318,10 +318,10 @@ impl<'a, T, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Normalised<T> where - T: Norm<F, L1> + LocalAnalysis<F, Bounds<F>, N>, + T: Norm<L1, F> + LocalAnalysis<F, Bounds<F>, N>, { #[inline] - fn local_analysis(&self, cube: &Cube<F, N>) -> Bounds<F> { + fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> { let Bounds(lower, upper) = self.0.local_analysis(cube); debug_assert!(lower <= upper); let w = self.0.norm(L1); @@ -330,9 +330,9 @@ } } -impl<'a, T, F: Float> Norm<F, L1> for Normalised<T> +impl<'a, T, F: Float> Norm<L1, F> for Normalised<T> where - T: Norm<F, L1>, + T: Norm<L1, F>, { #[inline] fn norm(&self, _: L1) -> F { @@ -347,8 +347,8 @@ macro_rules! impl_normalised_norm { ($($norm:ident)*) => { $( - impl<'a, T, F : Float> Norm<F, $norm> for Normalised<T> - where T : Norm<F, $norm> + Norm<F, L1> { + impl<'a, T, F : Float> Norm<$norm, F> for Normalised<T> + where T : Norm<$norm, F> + Norm<L1, F> { #[inline] fn norm(&self, n : $norm) -> F { let w = self.0.norm(L1); @@ -361,26 +361,26 @@ impl_normalised_norm!(L2 Linfinity); /* -impl<F : Num, S : Support<F, N>, const N : usize> LocalAnalysis<F, NullAggregator, N> for S { - fn local_analysis(&self, _cube : &Cube<F, N>) -> NullAggregator { NullAggregator } +impl<F : Num, S : Support< N, F>, const N : usize> LocalAnalysis<F, NullAggregator, N> for S { + fn local_analysis(&self, _cube : &Cube<N, F>) -> NullAggregator { NullAggregator } } impl<F : Float, S : Bounded<F>, const N : usize> LocalAnalysis<F, Bounds<F>, N> for S { #[inline] - fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> { + fn local_analysis(&self, cube : &Cube<N, F>) -> Bounds<F> { self.bounds(cube) } }*/ /// Generator of [`Support`]-implementing component functions based on low storage requirement /// [ids][`Self::Id`]. -pub trait SupportGenerator<F: Float, const N: usize>: +pub trait SupportGenerator<const N: usize, F: Float = f64>: MulAssign<F> + DivAssign<F> + Neg<Output = Self> + Clone + Sync + Send + 'static { /// The identification type type Id: 'static + Copy; /// The type of the [`Support`] (often also a [`Mapping`]). - type SupportType: 'static + Support<F, N>; + type SupportType: 'static + Support<N, F>; /// An iterator over all the [`Support`]s of the generator. type AllDataIter<'a>: Iterator<Item = (Self::Id, Self::SupportType)> where
--- a/src/bounds.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/bounds.rs Thu May 01 13:06:58 2025 -0500 @@ -40,7 +40,7 @@ /// 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; + fn local_analysis(&self, cube: &Cube<N, F>) -> A; } /// Trait for determining the upper and lower bounds of an float-valued [`Mapping`]. @@ -59,12 +59,12 @@ impl<F: Float, T: GlobalAnalysis<F, Bounds<F>>> Bounded<F> for T {} /// A [`RealMapping`] that provides rough bounds as well as minimisation and maximisation. -pub trait MinMaxMapping<F: Float, const N: usize>: RealMapping<F, N> + Bounded<F> { +pub trait MinMaxMapping<const N: usize, F: Float = f64>: RealMapping<N, F> + Bounded<F> { /// Maximise the mapping within stated value `tolerance`. /// /// At most `max_steps` refinement steps are taken. /// Returns the approximate maximiser and the corresponding function value. - fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F); + fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc<N, F>, F); /// Maximise the mapping within stated value `tolerance` subject to a lower bound. /// @@ -76,13 +76,13 @@ bound: F, tolerance: F, max_steps: usize, - ) -> Option<(Loc<F, N>, F)>; + ) -> Option<(Loc<N, F>, F)>; /// Minimise the mapping within stated value `tolerance`. /// /// At most `max_steps` refinement steps are taken. /// Returns the approximate minimiser and the corresponding function value. - fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc<F, N>, F); + fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc<N, F>, F); /// Minimise the mapping within stated value `tolerance` subject to a lower bound. /// @@ -94,7 +94,7 @@ bound: F, tolerance: F, max_steps: usize, - ) -> Option<(Loc<F, N>, F)>; + ) -> Option<(Loc<N, F>, F)>; /// Verify that the mapping has a given upper `bound` within indicated `tolerance`. ///
--- a/src/collection.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/collection.rs Thu May 01 13:06:58 2025 -0500 @@ -5,20 +5,24 @@ use crate::loc::Loc; /// An abstract collection of elements. -pub trait Collection : IntoIterator<Item = Self::Element> { +pub trait Collection: IntoIterator<Item = Self::Element> { /// Type of elements of the collection type Element; /// Iterator over references to elements of the collection - type RefsIter<'a> : Iterator<Item=&'a Self::Element> where Self : 'a; + type RefsIter<'a>: Iterator<Item = &'a Self::Element> + where + Self: 'a; /// Returns an iterator over references to elements of the collection. fn iter_refs(&self) -> Self::RefsIter<'_>; } /// An abstract collection of mutable elements. -pub trait CollectionMut : Collection { +pub trait CollectionMut: Collection { /// Iterator over references to elements of the collection - type RefsIterMut<'a> : Iterator<Item=&'a mut Self::Element> where Self : 'a; + type RefsIterMut<'a>: Iterator<Item = &'a mut Self::Element> + where + Self: 'a; /// Returns an iterator over references to elements of the collection. fn iter_refs_mut(&mut self) -> Self::RefsIterMut<'_>; @@ -51,4 +55,4 @@ slice_like_collection!(Vec<E> where E); slice_like_collection!([E; N] where E, const N : usize); -slice_like_collection!(Loc<E, N> where E, const N : usize); +slice_like_collection!(Loc<N, E> where E, const N : usize);
--- a/src/convex.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/convex.rs Thu May 01 13:06:58 2025 -0500 @@ -98,7 +98,7 @@ impl<F, E, Domain> Mapping<Domain> for NormConstraint<F, E> where - Domain: Space + Norm<F, E>, + Domain: Space + Norm<E, F>, F: Float, E: NormExponent, { @@ -126,8 +126,8 @@ where E: HasDualExponent, F: Float, - Domain: HasDual<F> + Norm<F, E> + Normed<F>, - <Domain as HasDual<F>>::DualSpace: Norm<F, E::DualExp>, + Domain: HasDual<F> + Norm<E, F> + Normed<F>, + <Domain as HasDual<F>>::DualSpace: Norm<E::DualExp, F>, { type Conjugate<'a> = NormConstraint<F, E::DualExp> @@ -147,8 +147,8 @@ C: Constant<Type = F>, E: HasDualExponent, F: Float, - Domain: HasDual<F> + Norm<F, E> + Space, - <Domain as HasDual<F>>::DualSpace: Norm<F, E::DualExp>, + Domain: HasDual<F> + Norm<E, F> + Space, + <Domain as HasDual<F>>::DualSpace: Norm<E::DualExp, F>, { type Conjugate<'a> = NormConstraint<F, E::DualExp> @@ -165,7 +165,7 @@ impl<Domain, E, F> Prox<Domain> for NormConstraint<F, E> where - Domain: Space + Norm<F, E>, + Domain: Space + Norm<E, F>, E: NormExponent, F: Float, NormProjection<F, E>: Mapping<Domain, Codomain = Domain>,
--- a/src/direct_product.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/direct_product.rs Thu May 01 13:06:58 2025 -0500 @@ -511,15 +511,15 @@ } } -impl<F, A, B, ExpA, ExpB, ExpJ> Norm<F, PairNorm<ExpA, ExpB, ExpJ>> for Pair<A, B> +impl<F, A, B, ExpA, ExpB, ExpJ> Norm<PairNorm<ExpA, ExpB, ExpJ>, F> for Pair<A, B> where F: Num, ExpA: NormExponent, ExpB: NormExponent, ExpJ: NormExponent, - A: Norm<F, ExpA>, - B: Norm<F, ExpB>, - Loc<F, 2>: Norm<F, ExpJ>, + A: Norm<ExpA, F>, + B: Norm<ExpB, F>, + Loc<2, F>: Norm<ExpJ, F>, { fn norm(&self, PairNorm(expa, expb, expj): PairNorm<ExpA, ExpB, ExpJ>) -> F { Loc([self.0.norm(expa), self.1.norm(expb)]).norm(expj)
--- a/src/discrete_gradient.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/discrete_gradient.rs Thu May 01 13:06:58 2025 -0500 @@ -392,7 +392,7 @@ where B: Discretisation<F>, F: Float + nalgebra::RealField, - DVector<F>: Norm<F, L2>, + DVector<F>: Norm<L2, F>, { fn opnorm_bound(&self, _: L2, _: L2) -> DynResult<F> { // Fuck nalgebra. @@ -405,7 +405,7 @@ where B: Discretisation<F>, F: Float + nalgebra::RealField, - DVector<F>: Norm<F, L2>, + DVector<F>: Norm<L2, F>, { fn opnorm_bound(&self, _: L2, _: L2) -> DynResult<F> { // Fuck nalgebra.
--- a/src/euclidean.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/euclidean.rs Thu May 01 13:06:58 2025 -0500 @@ -2,31 +2,37 @@ Euclidean spaces. */ -use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg}; -use crate::types::*; use crate::instance::Instance; use crate::norms::{HasDual, Reflexive}; +use crate::types::*; +use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign}; /// 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 an inner product. -pub trait Euclidean<F : Float> : HasDual<F, DualSpace=Self> + Reflexive<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> +pub trait Euclidean<F: Float = f64>: + HasDual<F, DualSpace = Self> + + Reflexive<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>; + type Output: Euclidean<F>; // Inner product - fn dot<I : Instance<Self>>(&self, other : I) -> F; + fn dot<I: Instance<Self>>(&self, other: I) -> F; /// Calculate the square of the 2-norm, $\frac{1}{2}\\|x\\|_2^2$, where `self` is $x$. /// @@ -38,7 +44,7 @@ /// where `self` is $x$. #[inline] fn norm2_squared_div2(&self) -> F { - self.norm2_squared()/F::TWO + self.norm2_squared() / F::TWO } /// Calculate the 2-norm $‖x‖_2$, where `self` is $x$. @@ -48,33 +54,33 @@ } /// Calculate the 2-distance squared $\\|x-y\\|_2^2$, where `self` is $x$. - fn dist2_squared<I : Instance<Self>>(&self, y : I) -> F; + fn dist2_squared<I: Instance<Self>>(&self, y: I) -> F; /// Calculate the 2-distance $\\|x-y\\|_2$, where `self` is $x$. #[inline] - fn dist2<I : Instance<Self>>(&self, y : I) -> F { + fn dist2<I: Instance<Self>>(&self, y: I) -> F { self.dist2_squared(y).sqrt() } /// Projection to the 2-ball. #[inline] - fn proj_ball2(mut self, ρ : F) -> Self { + 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) { + fn proj_ball2_mut(&mut self, ρ: F) { let r = self.norm2(); - if r>ρ { - *self *= ρ/r + if r > ρ { + *self *= ρ / r } } } /// Trait for [`Euclidean`] spaces with dimensions known at compile time. -pub trait StaticEuclidean<F : Float> : Euclidean<F> { +pub trait StaticEuclidean<F: Float = f64>: Euclidean<F> { /// Returns the origin fn origin() -> <Self as Euclidean<F>>::Output; }
--- a/src/fe_model/p2_local_model.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/fe_model/p2_local_model.rs Thu May 01 13:06:58 2025 -0500 @@ -2,21 +2,21 @@ 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 super::base::{LocalModel, RealLocalModel}; use crate::euclidean::Euclidean; use crate::instance::Instance; -use super::base::{LocalModel,RealLocalModel}; +use crate::linsolve::*; +use crate::loc::Loc; use crate::sets::Cube; +use crate::sets::{NPolygon, Set, SpannedHalfspace}; +use crate::types::*; 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 struct Simplex<F: Float, const N: usize, const D: usize>(pub [Loc<N, F>; D]); /// A two-dimensional planar simplex pub type PlanarSimplex<F> = Simplex<F, 2, 3>; /// A real interval @@ -25,26 +25,29 @@ /// 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> { - (a+b)/2.0 +pub(crate) fn midpoint<F: Float, const N: usize>(a: &Loc<N, F>, b: &Loc<N, F>) -> Loc<N, F> { + (a + b) / 2.0 } -impl<'a, F : Float> Set<Loc<F,1>> for RealInterval<F> { +impl<'a, F: Float> Set<Loc<1, F>> for RealInterval<F> { #[inline] - fn contains<I : Instance<Loc<F, 1>>>(&self, z : I) -> bool { + fn contains<I: Instance<Loc<1, F>>>(&self, z: I) -> bool { let &Loc([x]) = z.ref_instance(); let &[Loc([x0]), Loc([x1])] = &self.0; (x0 < x && x < x1) || (x1 < x && x < x0) } } -impl<'a, F : Float> Set<Loc<F,2>> for PlanarSimplex<F> { +impl<'a, F: Float> Set<Loc<2, F>> for PlanarSimplex<F> { #[inline] - fn contains<I : Instance<Loc<F, 2>>>(&self, z : I) -> bool { + fn contains<I: Instance<Loc<2, F>>>(&self, z: I) -> bool { let &[x0, x1, x2] = &self.0; - NPolygon([[x0, x1].spanned_halfspace(), - [x1, x2].spanned_halfspace(), - [x2, x0].spanned_halfspace()]).contains(z) + NPolygon([ + [x0, x1].spanned_halfspace(), + [x1, x2].spanned_halfspace(), + [x2, x0].spanned_halfspace(), + ]) + .contains(z) } } @@ -58,121 +61,118 @@ } #[replace_float_literals(F::cast_from(literal))] -impl<F : Float> P2Powers for Loc<F, 1> { - type Output = Loc<F, 1>; - type Full = Loc<F, 3>; - type Diff = Loc<Loc<F, 1>, 1>; +impl<F: Float> P2Powers for Loc<1, F> { + type Output = Loc<1, F>; + type Full = Loc<3, F>; + type Diff = Loc<1, Loc<1, F>>; #[inline] fn p2powers(&self) -> Self::Output { let &Loc([x0]) = self; - [x0*x0].into() + [x0 * x0].into() } #[inline] fn p2powers_full(&self) -> Self::Full { let &Loc([x0]) = self; - [1.0, x0, x0*x0].into() + [1.0, x0, x0 * x0].into() } #[inline] fn p2powers_diff(&self) -> Self::Diff { let &Loc([x0]) = self; - [[x0+x0].into()].into() + [[x0 + x0].into()].into() } } #[replace_float_literals(F::cast_from(literal))] -impl<F : Float> P2Powers for Loc<F, 2> { - type Output = Loc<F, 3>; - type Full = Loc<F, 6>; - type Diff = Loc<Loc<F, 3>, 2>; +impl<F: Float> P2Powers for Loc<2, F> { + type Output = Loc<3, F>; + type Full = Loc<6, F>; + type Diff = Loc<2, Loc<3, F>>; #[inline] fn p2powers(&self) -> Self::Output { let &Loc([x0, x1]) = self; - [x0*x0, x0*x1, x1*x1].into() + [x0 * x0, x0 * x1, x1 * x1].into() } #[inline] fn p2powers_full(&self) -> Self::Full { let &Loc([x0, x1]) = self; - [1.0, x0, x1, x0*x0, x0*x1, x1*x1].into() + [1.0, x0, x1, x0 * x0, x0 * x1, x1 * x1].into() } #[inline] fn p2powers_diff(&self) -> Self::Diff { let &Loc([x0, x1]) = self; - [[x0+x0, x1, 0.0].into(), [0.0, x0, x1+x1].into()].into() + [[x0 + x0, x1, 0.0].into(), [0.0, x0, x1 + x1].into()].into() } } /// 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> { +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>; + type Model: LocalModel<Loc<N, F>, 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; + fn p2_model<G: Fn(&Loc<N, F>) -> F>(&self, g: G) -> Self::Model; } /// 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>, - //node_values : Loc<F, V>, - //edge_values : Loc<F, Q>, +pub struct P2LocalModel< + F: Num, + const N: usize, + const E: usize, /*, const V : usize, const Q : usize*/ +> { + a0: F, + a1: Loc<N, F>, + a2: Loc<E, F>, + //node_values : Loc<V, F>, + //edge_values : Loc<Q, F>, } // // 1D planar model construction // -impl<F : Float> RealInterval<F> { +impl<F: Float> RealInterval<F> { #[inline] - fn midpoints(&self) -> [Loc<F, 1>; 1] { + fn midpoints(&self) -> [Loc<1, F>; 1] { let [ref n0, ref n1] = &self.0; let n01 = midpoint(n0, n1); [n01] } } -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], - &[v0, v1, v01] : &[F; 3], - ) -> Self { - let p = move |x : &Loc<F, 1>, v : F| { + pub fn new(&[n0, n1, n01]: &[Loc<1, F>; 3], &[v0, v1, v01]: &[F; 3]) -> Self { + let p = move |x: &Loc<1, F>, v: F| { let Loc([c, d, e]) = x.p2powers_full(); [c, d, e, v] }; - let [a0, a1, a11] = linsolve([ - p(&n0, v0), - p(&n1, v1), - p(&n01, v01) - ]); + let [a0, a1, a11] = linsolve([p(&n0, v0), p(&n1, v1), p(&n01, v01)]); P2LocalModel { - a0 : a0, - a1 : [a1].into(), - a2 : [a11].into(), + a0: a0, + a1: [a1].into(), + a2: [a11].into(), //node_values : [v0, v1].into(), //edge_values: [].into(), } } } -impl<F : Float> P2Model<F,1> for RealInterval<F> { - type Model = P2LocalModel<F, 1, 1/*, 2, 0*/>; +impl<F: Float> P2Model<F, 1> for RealInterval<F> { + type Model = P2LocalModel<F, 1, 1 /*, 2, 0*/>; #[inline] - fn p2_model<G : Fn(&Loc<F, 1>) -> F>(&self, g : G) -> Self::Model { - let [n01] = self.midpoints(); + fn p2_model<G: Fn(&Loc<1, F>) -> F>(&self, g: G) -> Self::Model { + let [n01] = self.midpoints(); let [n0, n1] = self.0; let vals = [g(&n0), g(&n1), g(&n01)]; let nodes = [n0, n1, n01]; @@ -184,10 +184,10 @@ // 2D planar model construction // -impl<F : Float> PlanarSimplex<F> { +impl<F: Float> PlanarSimplex<F> { #[inline] /// Returns the midpoints of all the edges of the simplex - fn midpoints(&self) -> [Loc<F, 2>; 3] { + fn midpoints(&self) -> [Loc<2, F>; 3] { let [ref n0, ref n1, ref n2] = &self.0; let n01 = midpoint(n0, n1); let n12 = midpoint(n1, n2); @@ -196,15 +196,15 @@ } } -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], - &[v0, v1, v2, v01, v12, v20] : &[F; 6], + &[n0, n1, n2, n01, n12, n20]: &[Loc<2, F>; 6], + &[v0, v1, v2, v01, v12, v20]: &[F; 6], ) -> Self { - let p = move |x : &Loc<F,2>, v :F| { + let p = move |x: &Loc<2, F>, v: F| { let Loc([c, d, e, f, g, h]) = x.p2powers_full(); [c, d, e, f, g, h, v] }; @@ -217,20 +217,20 @@ p(&n20, v20), ]); P2LocalModel { - a0 : a0, - a1 : [a1, a2].into(), - a2 : [a11, a12, a22].into(), + a0: a0, + a1: [a1, a2].into(), + a2: [a11, a12, a22].into(), //node_values : [v0, v1, v2].into(), //edge_values: [v01, v12, v20].into(), } } } -impl<F : Float> P2Model<F,2> for PlanarSimplex<F> { - type Model = P2LocalModel<F, 2, 3/*, 3, 3*/>; +impl<F: Float> P2Model<F, 2> for PlanarSimplex<F> { + type Model = P2LocalModel<F, 2, 3 /*, 3, 3*/>; #[inline] - fn p2_model<G : Fn(&Loc<F, 2>) -> F>(&self, g : G) -> Self::Model { + fn p2_model<G: Fn(&Loc<2, F>) -> F>(&self, g: G) -> Self::Model { let midpoints = self.midpoints(); let [ref n0, ref n1, ref n2] = self.0; let [ref n01, ref n12, ref n20] = midpoints; @@ -242,125 +242,137 @@ 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<$n, F>, F> for P2LocalModel<F, $n, /*, $v, $q*/ $e> { #[inline] - fn value(&self, x : &Loc<F,$n>) -> F { + fn value(&self, x: &Loc<$n, F>) -> F { self.a0 + x.dot(&self.a1) + x.p2powers().dot(&self.a2) } #[inline] - fn differential(&self, x : &Loc<F,$n>) -> Loc<F,$n> { + fn differential(&self, x: &Loc<$n, F>) -> Loc<$n, F> { self.a1 + x.p2powers_diff().map(|di| di.dot(&self.a2)) } } - } + }; } impl_local_model!(1, 1, 2, 0); impl_local_model!(2, 3, 3, 3); - // // Minimisation // #[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{ - a1 : Loc([a1]), - a2 : Loc([a11]), + fn minimise_edge(&self, x0: Loc<1, F>, x1: Loc<1, F>) -> (Loc<1, F>, F) { + let &P2LocalModel { + a1: Loc([a1]), + a2: Loc([a11]), //node_values : Loc([v0, v1]), .. - } = self; + } = self; // We do this in cases, first trying for an interior solution, then edges. // For interior solution, first check determinant; no point trying if non-positive if a11 > 0.0 { // An interior solution x[1] has to satisfy // 2a₁₁*x[1] + a₁ =0 // This gives - let t = -a1/(2.0*a11); + let t = -a1 / (2.0 * a11); let (Loc([t0]), Loc([t1])) = (x0, x1); if (t0 <= t && t <= t1) || (t1 <= t && t <= t0) { let x = [t].into(); let v = self.value(&x); - return (x, v) + return (x, v); } } let v0 = self.value(&x0); let v1 = self.value(&x1); - if v0 < v1 { (x0, v0) } else { (x1, v1) } + if v0 < v1 { + (x0, v0) + } else { + (x1, v1) + } } } -impl<'a, F : Float> RealLocalModel<RealInterval<F>,Loc<F,1>,F> -for P2LocalModel<F, 1, 1/*, 2, 0*/> { +impl<'a, F: Float> RealLocalModel<RealInterval<F>, Loc<1, F>, F> + for P2LocalModel<F, 1, 1 /*, 2, 0*/> +{ #[inline] - fn minimise(&self, &Simplex([x0, x1]) : &RealInterval<F>) -> (Loc<F,1>, F) { + fn minimise(&self, &Simplex([x0, x1]): &RealInterval<F>) -> (Loc<1, F>, F) { self.minimise_edge(x0, x1) } } #[replace_float_literals(F::cast_from(literal))] -impl<F : Float> P2LocalModel<F, 2, 3/*, 3, 3*/> { +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) { + fn minimise_edge( + &self, + x0: &Loc<2, F>, + x1: &Loc<2, F>, /*, v0 : F, v1 : F*/ + ) -> (Loc<2, F>, F) { let &P2LocalModel { a0, - a1 : Loc([a1, a2]), - a2 : Loc([a11, a12, a22]), + a1: Loc([a1, a2]), + a2: Loc([a11, a12, a22]), .. - } = self; + } = self; let &Loc([x00, x01]) = x0; - let d@Loc([d0, d1]) = x1 - x0; - let b0 = a0 + a1*x00 + a2*x01 + a11*x00*x00 + a12*x00*x01 + a22*x01*x01; - let b1 = a1*d0 + a2*d1 + 2.0*a11*d0*x00 + a12*(d0*x01 + d1*x00) + 2.0*a22*d1*x01; - let b11 = a11*d0*d0 + a12*d0*d1 + a22*d1*d1; + let d @ Loc([d0, d1]) = x1 - x0; + let b0 = a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + a12 * x00 * x01 + a22 * x01 * x01; + let b1 = a1 * d0 + + a2 * d1 + + 2.0 * a11 * d0 * x00 + + a12 * (d0 * x01 + d1 * x00) + + 2.0 * a22 * d1 * x01; + let b11 = a11 * d0 * d0 + a12 * d0 * d1 + a22 * d1 * d1; let edge_1d_model = P2LocalModel { - a0 : b0, - a1 : Loc([b1]), - a2 : Loc([b11]), + a0: b0, + a1: Loc([b1]), + a2: Loc([b11]), //node_values : Loc([v0, v1]), }; let (Loc([t]), v) = edge_1d_model.minimise_edge(0.0.into(), 1.0.into()); - (x0 + d*t, v) + (x0 + d * t, v) } } #[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*/> { +impl<'a, F: Float> RealLocalModel<PlanarSimplex<F>, Loc<2, F>, F> + for P2LocalModel<F, 2, 3 /*, 3, 3*/> +{ #[inline] - fn minimise(&self, el : &PlanarSimplex<F>) -> (Loc<F,2>, F) { + fn minimise(&self, el: &PlanarSimplex<F>) -> (Loc<2, F>, F) { let &P2LocalModel { - a1 : Loc([a1, a2]), - a2 : Loc([a11, a12, a22]), + a1: Loc([a1, a2]), + a2: Loc([a11, a12, a22]), //node_values : Loc([v0, v1, v2]), .. - } = self; + } = self; // We do this in cases, first trying for an interior solution, then edges. // For interior solution, first check determinant; no point trying if non-positive - let r = 2.0*(a11*a22-a12*a12); + let r = 2.0 * (a11 * a22 - a12 * a12); if r > 0.0 { // An interior solution (x[1], x[2]) has to satisfy // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0 // This gives - let x = [(a22*a1-a12*a2)/r, (a12*a1-a11*a2)/r].into(); + let x = [(a22 * a1 - a12 * a2) / r, (a12 * a1 - a11 * a2) / r].into(); if el.contains(&x) { - return (x, self.value(&x)) + return (x, self.value(&x)); } } let &[ref x0, ref x1, ref x2] = &el.0; let mut min_edge = self.minimise_edge(x0, x1); - let more_edge = [self.minimise_edge(x1, x2), - self.minimise_edge(x2, x0)]; - + let more_edge = [self.minimise_edge(x1, x2), self.minimise_edge(x2, x0)]; + for edge in more_edge { if edge.1 < min_edge.1 { min_edge = edge; @@ -372,35 +384,38 @@ } #[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*/> { +impl<'a, F: Float> RealLocalModel<Cube<2, F>, Loc<2, F>, F> + for P2LocalModel<F, 2, 3 /*, 3, 3*/> +{ #[inline] - fn minimise(&self, el : &Cube<F, 2>) -> (Loc<F,2>, F) { + fn minimise(&self, el: &Cube<2, F>) -> (Loc<2, F>, F) { let &P2LocalModel { - a1 : Loc([a1, a2]), - a2 : Loc([a11, a12, a22]), + a1: Loc([a1, a2]), + a2: Loc([a11, a12, a22]), //node_values : Loc([v0, v1, v2]), .. - } = self; + } = self; // We do this in cases, first trying for an interior solution, then edges. // For interior solution, first check determinant; no point trying if non-positive - let r = 2.0*(a11*a22-a12*a12); + let r = 2.0 * (a11 * a22 - a12 * a12); if r > 0.0 { // An interior solution (x[1], x[2]) has to satisfy // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0 // This gives - let x = [(a22*a1-a12*a2)/r, (a12*a1-a11*a2)/r].into(); + let x = [(a22 * a1 - a12 * a2) / r, (a12 * a1 - a11 * a2) / r].into(); if el.contains(&x) { - return (x, self.value(&x)) + return (x, self.value(&x)); } } let [x0, x1, x2, x3] = el.corners(); let mut min_edge = self.minimise_edge(&x0, &x1); - let more_edge = [self.minimise_edge(&x1, &x2), - self.minimise_edge(&x2, &x3), - self.minimise_edge(&x3, &x0)]; + let more_edge = [ + self.minimise_edge(&x1, &x2), + self.minimise_edge(&x2, &x3), + self.minimise_edge(&x3, &x0), + ]; for edge in more_edge { if edge.1 < min_edge.1 { @@ -422,7 +437,7 @@ let domain = Simplex(vertices); // A simple quadratic function for which the approximation is exact on reals, // and appears exact on f64 as well. - let f = |&Loc([x]) : &Loc<f64, 1>| x*x + x + 1.0; + let f = |&Loc([x]): &Loc<1, f64>| x * x + x + 1.0; let model = domain.p2_model(f); let xs = [Loc([0.5]), Loc([0.25])]; @@ -439,7 +454,7 @@ let domain = Simplex(vertices); // A simple quadratic function for which the approximation is exact on reals, // and appears exact on f64 as well. - let f = |&Loc([x, y]) : &Loc<f64, 2>| - (x*x + x*y + x - 2.0 * y + 1.0); + let f = |&Loc([x, y]): &Loc<2, f64>| -(x * x + x * y + x - 2.0 * y + 1.0); let model = domain.p2_model(f); let xs = [Loc([0.5, 0.5]), Loc([0.25, 0.25])];
--- a/src/lingrid.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/lingrid.rs Thu May 01 13:06:58 2025 -0500 @@ -15,46 +15,54 @@ iteration over the grid. Additional utility functions are in the [`Grid`] trait. */ -use crate::types::*; +use crate::iter::{RestartableIterator, StatefulIterator}; use crate::loc::Loc; +use crate::maputil::{map2, map4}; use crate::sets::Cube; -use crate::iter::{RestartableIterator, StatefulIterator}; -use crate::maputil::{map2, map4}; -use serde::{Serialize, Deserialize}; +use crate::types::*; +use serde::{Deserialize, Serialize}; // TODO: rewrite this using crate::sets::Cube. /// 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`. +/// or [`Loc`]`<N, F>` 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, + 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]>; +pub type LinGrid<const N: usize, F: Float = f64> = LinSpace<Loc<N, F>, [usize; N]>; /// Creates a [`LinSpace`] on the real line. -pub fn linspace<F : Float>(start : F, end : F, count : usize) -> LinSpace<F, usize> { - LinSpace{ start : start, end : end, count : count } +pub fn linspace<F: Float>(start: F, end: F, count: usize) -> LinSpace<F, usize> { + LinSpace { + start: start, + end: end, + count: count, + } } /// 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. -pub fn lingrid<F : Float, const N : usize>( - cube : &Cube<F, N>, - count : &[usize; N] -) -> LinSpace<Loc<F, N>, [usize; N]> { - LinSpace{ start : cube.span_start(), end : cube.span_end(), count : *count } +pub fn lingrid<F: Float, const N: usize>( + cube: &Cube<N, F>, + count: &[usize; N], +) -> LinSpace<Loc<N, F>, [usize; N]> { + LinSpace { + start: cube.span_start(), + end: cube.span_end(), + count: *count, + } } /// Create a multi-dimensional linear grid with centered nodes. @@ -63,30 +71,33 @@ /// inside `cube`. Thus, if $w\_i$ is the width of the cube along dimension $i$, and $n_i$ the number /// of nodes, the width of the subcube along this dimension is $h_i = w\_i/(n\_i+1)$, and the first /// and last nodes are at a distance $h\_i/2$ from the closest boundary. -pub fn lingrid_centered<F : Float, const N : usize>( - cube : &Cube<F, N>, - count : &[usize; N] -) -> LinSpace<Loc<F, N>, [usize; N]> { +pub fn lingrid_centered<F: Float, const N: usize>( + cube: &Cube<N, F>, + count: &[usize; N], +) -> LinSpace<Loc<N, F>, [usize; N]> { let h_div_2 = map2(cube.width(), count, |w, &n| w / F::cast_from(2 * (n + 1))); let span_start = map2(cube.span_start(), &h_div_2, |a, &t| a + t).into(); - let span_end = map2(cube.span_end(), &h_div_2, |b, &t| b - t).into(); - LinSpace{ start : span_start, end : span_end, count : *count } + let span_end = map2(cube.span_end(), &h_div_2, |b, &t| b - t).into(); + LinSpace { + start: span_start, + end: span_end, + count: *count, + } } - /// Iterator over a `LinSpace`. #[derive(Clone, Debug)] pub struct LinSpaceIterator<F, I> { - lingrid : LinSpace<F,I>, - current : Option<I>, + lingrid: LinSpace<F, I>, + current: Option<I>, } /// 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; + 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_unchecked(&self, i: &I) -> U; // fn entry(&self, i : I) -> Option<F> } @@ -97,7 +108,7 @@ fn next_index(&mut self) -> Option<I>; } -impl<F : Float + CastFrom<I>, I : Unsigned> Grid<F, I> for LinSpace<F, I> { +impl<F: Float + CastFrom<I>, I: Unsigned> Grid<F, I> for LinSpace<F, I> { /*fn entry(&self, i : I) -> Option<F> { if i < self.count { Some(self.entry_unchecked(i)) @@ -107,37 +118,40 @@ }*/ #[inline] - fn entry_linear_unchecked(&self, i : usize) -> F { + fn entry_linear_unchecked(&self, i: usize) -> F { self.entry_unchecked(&I::cast_from(i)) } #[inline] - fn entry_unchecked(&self, i : &I) -> F { + fn entry_unchecked(&self, i: &I) -> F { let idx = F::cast_from(*i); - let scale = F::cast_from(self.count-I::ONE); - self.start + (self.end-self.start)*idx/scale + let scale = F::cast_from(self.count - I::ONE); + self.start + (self.end - self.start) * idx / scale } } -impl<F : Float + CastFrom<I>, I : Unsigned> GridIteration<F, I> -for LinSpaceIterator<F, I> { +impl<F: Float + CastFrom<I>, I: Unsigned> GridIteration<F, I> for LinSpaceIterator<F, I> { #[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, } } } -impl<F : Float + CastFrom<I>, I : Unsigned, const N : usize> Grid<Loc<F,N>, [I; N]> -for LinSpace<Loc<F,N>, [I; N]> { +impl<F: Float + CastFrom<I>, I: Unsigned, const N: usize> Grid<Loc<N, F>, [I; N]> + for LinSpace<Loc<N, F>, [I; N]> +{ #[inline] - fn entry_linear_unchecked(&self, i_ : usize) -> Loc<F, N> { + fn entry_linear_unchecked(&self, i_: usize) -> Loc<N, F> { let mut i = I::cast_from(i_); let mut tmp = [I::ZERO; N]; for k in 0..N { @@ -148,58 +162,66 @@ } #[inline] - fn entry_unchecked(&self, i : &[I; N]) -> Loc<F, N> { - let LinSpace{ start, end, count } = self; + fn entry_unchecked(&self, i: &[I; N]) -> Loc<N, F> { + let LinSpace { start, end, count } = self; map4(i, start, end, count, |&ik, &sk, &ek, &ck| { let idx = F::cast_from(ik); - let scale = F::cast_from(ck-I::ONE); + let scale = F::cast_from(ck - I::ONE); sk + (ek - sk) * idx / scale - }).into() + }) + .into() } } -impl<F : Float + CastFrom<I>, I : Unsigned, const N : usize> GridIteration<Loc<F,N>, [I; N]> -for LinSpaceIterator<Loc<F,N>, [I; N]> { - +impl<F: Float + CastFrom<I>, I: Unsigned, const N: usize> GridIteration<Loc<N, F>, [I; N]> + for LinSpaceIterator<Loc<N, F>, [I; N]> +{ #[inline] fn next_index(&mut self) -> Option<[I; N]> { match self.current { - None if self.lingrid.count.iter().all(|v| I::ZERO < *v) => { + None if self.lingrid.count.iter().all(|v| I::ZERO < *v) => { self.current = Some([I::ZERO; N]); self.current - }, + } Some(ref mut v) => { for k in 0..N { let a = v[k] + I::ONE; if a < self.lingrid.count[k] { v[k] = a; - return self.current + return self.current; } else { v[k] = I::ZERO; } } None - }, - _ => None + } + _ => None, } } } -impl<F, I> IntoIterator for LinSpace<F,I> -where LinSpace<F, I> : Grid<F, I>, - LinSpaceIterator<F, I> : GridIteration<F, I> { +impl<F, I> IntoIterator for LinSpace<F, I> +where + LinSpace<F, I>: Grid<F, I>, + LinSpaceIterator<F, I>: GridIteration<F, I>, +{ type Item = F; - type IntoIter = LinSpaceIterator<F,I>; + type IntoIter = LinSpaceIterator<F, I>; #[inline] fn into_iter(self) -> Self::IntoIter { - LinSpaceIterator { lingrid : self, current : None } + LinSpaceIterator { + lingrid: self, + current: None, + } } } -impl<F, I> Iterator for LinSpaceIterator<F,I> -where LinSpace<F, I> : Grid<F, I>, - LinSpaceIterator<F, I> : GridIteration<F, I> { +impl<F, I> Iterator for LinSpaceIterator<F, I> +where + LinSpace<F, I>: Grid<F, I>, + LinSpaceIterator<F, I>: GridIteration<F, I>, +{ type Item = F; #[inline] fn next(&mut self) -> Option<F> { @@ -207,19 +229,24 @@ } } -impl<F, I> StatefulIterator for LinSpaceIterator<F,I> -where LinSpace<F, I> : Grid<F, I>, - LinSpaceIterator<F, I> : GridIteration<F, I> { +impl<F, I> StatefulIterator for LinSpaceIterator<F, I> +where + LinSpace<F, I>: Grid<F, I>, + LinSpaceIterator<F, I>: GridIteration<F, I>, +{ #[inline] fn current(&self) -> Option<F> { - self.current.as_ref().map(|c| self.lingrid.entry_unchecked(c)) + self.current + .as_ref() + .map(|c| self.lingrid.entry_unchecked(c)) } } - -impl<F, I> RestartableIterator for LinSpaceIterator<F,I> -where LinSpace<F, I> : Grid<F, I>, - LinSpaceIterator<F, I> : GridIteration<F, I> { +impl<F, I> RestartableIterator for LinSpaceIterator<F, I> +where + LinSpace<F, I>: Grid<F, I>, + LinSpaceIterator<F, I>: GridIteration<F, I>, +{ #[inline] fn restart(&mut self) -> Option<F> { self.current = None;
--- a/src/linops.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/linops.rs Thu May 01 13:06:58 2025 -0500 @@ -77,7 +77,7 @@ pub trait BoundedLinear<X, XExp, CodExp, F = f64>: Linear<X> where F: Num, - X: Space + Norm<F, XExp>, + X: Space + Norm<XExp, F>, XExp: NormExponent, CodExp: NormExponent, { @@ -182,7 +182,7 @@ impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X> where - X: Space + Clone + Norm<F, E>, + X: Space + Clone + Norm<E, F>, F: Num, E: NormExponent, { @@ -264,8 +264,8 @@ impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F> where - X: Space + Norm<F, E1>, - Y: AXPY<F> + Clone + Norm<F, E2>, + X: Space + Norm<E1, F>, + Y: AXPY<F> + Clone + Norm<E2, F>, F: Num, E1: NormExponent, E2: NormExponent, @@ -346,8 +346,8 @@ impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp> where F: Num, - X: Space + Norm<F, Xexp>, - Z: Space + Norm<F, Zexp>, + X: Space + Norm<Xexp, F>, + Z: Space + Norm<Zexp, F>, Xexp: NormExponent, Yexp: NormExponent, Zexp: NormExponent, @@ -680,8 +680,8 @@ BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> for RowOp<S, T> where F: Float, - A: Space + Norm<F, ExpA>, - B: Space + Norm<F, ExpB>, + A: Space + Norm<ExpA, F>, + B: Space + Norm<ExpB, F>, S: BoundedLinear<A, ExpA, ExpR, F>, T: BoundedLinear<B, ExpB, ExpR, F>, S::Codomain: Add<T::Codomain>, @@ -707,7 +707,7 @@ for ColOp<S, T> where F: Float, - A: Space + Norm<F, ExpA>, + A: Space + Norm<ExpA, F>, S: BoundedLinear<A, ExpA, ExpS, F>, T: BoundedLinear<A, ExpA, ExpT, F>, ExpA: NormExponent,
--- a/src/loc.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/loc.rs Thu May 01 13:06:58 2025 -0500 @@ -3,44 +3,45 @@ 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::euclidean::*; +use crate::instance::{BasicDecomposition, Instance}; +use crate::linops::{Linear, Mapping, AXPY}; +use crate::mapping::Space; +use crate::maputil::{map1, map1_mut, map2, map2_mut, FixedLength, FixedLengthMut}; +use crate::norms::*; +use crate::types::{Float, Num, SignedNum}; +use serde::ser::{Serialize, SerializeSeq, Serializer}; use std::fmt::{Display, Formatter}; -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, Mapping, Linear}; -use crate::instance::{Instance, BasicDecomposition}; -use crate::mapping::Space; -use serde::ser::{Serialize, Serializer, SerializeSeq}; - +use std::ops::{ + Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign, +}; +use std::slice::{Iter, IterMut}; /// 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>( +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +pub struct Loc<const N: usize, F = f64>( /// An array of the elements of the vector - pub [F; N] + pub [F; N], ); -impl<F : Display, const N : usize> Display for Loc<F, N>{ +impl<F: Display, const N: usize> Display for Loc<N, F> { // Required method fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { - write!(f, "[")?; - let mut comma = ""; - for e in self.iter() { + write!(f, "[")?; + let mut comma = ""; + for e in self.iter() { write!(f, "{comma}{e}")?; comma = ", "; - } - write!(f, "]") + } + write!(f, "]") } } // Need to manually implement as [F; N] serialisation is provided only for some N. -impl<F, const N : usize> Serialize for Loc<F, N> +impl<F, const N: usize> Serialize for Loc<N, F> where F: Serialize, { @@ -56,10 +57,10 @@ } } -impl<F, const N : usize> Loc<F, N> { +impl<F, const N: usize> Loc<N, F> { /// Creates a new `Loc` vector from an array. #[inline] - pub fn new(arr : [F; N]) -> Self { + pub fn new(arr: [F; N]) -> Self { Loc(arr) } @@ -76,43 +77,44 @@ } } -impl<F : Copy, const N : usize> Loc<F, N> { +impl<F: Copy, const N: usize> Loc<N, F> { /// Maps `g` over the elements of the vector, returning a new [`Loc`] vector #[inline] - pub fn map<H>(&self, g : impl Fn(F) -> H) -> Loc<H, N> { + pub fn map<H>(&self, g: impl Fn(F) -> H) -> Loc<N, H> { 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>(&self, other : &Self, g : impl Fn(F, F) -> H) -> Loc<H, N> { + pub fn map2<H>(&self, other: &Self, g: impl Fn(F, F) -> H) -> Loc<N, H> { 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(&mut self, g : impl Fn(&mut F)) { + 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(&mut self, other : &Self, g : impl Fn(&mut F, F)) { + 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<A : Num>(&self, g : impl Fn(F) -> A) -> A { + pub fn product_map<A: Num>(&self, g: impl Fn(F) -> A) -> A { match N { 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)) + 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)), } } } @@ -130,29 +132,28 @@ ($($x:expr),+ $(,)?) => { Loc::new([$($x),+]) } } - -impl<F, const N : usize> From<[F; N]> for Loc<F, N> { +impl<F, const N: usize> From<[F; N]> for Loc<N, F> { #[inline] - fn from(other: [F; N]) -> Loc<F, N> { + fn from(other: [F; N]) -> Loc<N, F> { Loc(other) } } -/*impl<F : Copy, const N : usize> From<&[F; N]> for Loc<F, N> { +/*impl<F : Copy, const N : usize> From<&[F; N]> for Loc<N, F> { #[inline] - fn from(other: &[F; N]) -> Loc<F, N> { + fn from(other: &[F; N]) -> Loc<N, F> { Loc(*other) } }*/ -impl<F> From<F> for Loc<F, 1> { +impl<F> From<F> for Loc<1, F> { #[inline] - fn from(other: F) -> Loc<F, 1> { + fn from(other: F) -> Loc<1, F> { Loc([other]) } } -impl<F> Loc<F, 1> { +impl<F> Loc<1, F> { #[inline] pub fn flatten1d(self) -> F { let Loc([v]) = self; @@ -160,22 +161,21 @@ } } -impl<F, const N : usize> From<Loc<F, N>> for [F; N] { +impl<F, const N: usize> From<Loc<N, F>> for [F; N] { #[inline] - fn from(other : Loc<F, N>) -> [F; N] { + fn from(other: Loc<N, F>) -> [F; N] { other.0 } } -/*impl<F : Copy, const N : usize> From<&Loc<F, N>> for [F; N] { +/*impl<F : Copy, const N : usize> From<&Loc<N, F>> for [F; N] { #[inline] - fn from(other : &Loc<F, N>) -> [F; N] { + fn from(other : &Loc<N, F>) -> [F; N] { other.0 } }*/ - -impl<F, const N : usize> IntoIterator for Loc<F, N> { +impl<F, const N: usize> IntoIterator for Loc<N, F> { type Item = <[F; N] as IntoIterator>::Item; type IntoIter = <[F; N] as IntoIterator>::IntoIter; @@ -187,20 +187,24 @@ // Indexing -impl<F, Ix, const N : usize> Index<Ix> for Loc<F,N> -where [F; N] : Index<Ix> { +impl<F, Ix, const N: usize> Index<Ix> for Loc<N, F> +where + [F; N]: Index<Ix>, +{ type Output = <[F; N] as Index<Ix>>::Output; #[inline] - fn index(&self, ix : Ix) -> &Self::Output { + fn index(&self, ix: Ix) -> &Self::Output { self.0.index(ix) } } -impl<F, Ix, const N : usize> IndexMut<Ix> for Loc<F,N> -where [F; N] : IndexMut<Ix> { +impl<F, Ix, const N: usize> IndexMut<Ix> for Loc<N, F> +where + [F; N]: IndexMut<Ix>, +{ #[inline] - fn index_mut(&mut self, ix : Ix) -> &mut Self::Output { + fn index_mut(&mut self, ix: Ix) -> &mut Self::Output { self.0.index_mut(ix) } } @@ -209,61 +213,61 @@ macro_rules! make_binop { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { - impl<F : Num, const N : usize> $trait<Loc<F,N>> for Loc<F, N> { - type Output = Loc<F, N>; + impl<F: Num, const N: usize> $trait<Loc<N, F>> for Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(mut self, other : Loc<F, N>) -> Self::Output { + fn $fn(mut self, other: Loc<N, F>) -> Self::Output { self.$fn_assign(other); self } } - impl<'a, F : Num, const N : usize> $trait<&'a Loc<F,N>> for Loc<F, N> { - type Output = Loc<F, N>; + impl<'a, F: Num, const N: usize> $trait<&'a Loc<N, F>> for Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(mut self, other : &'a Loc<F, N>) -> Self::Output { + fn $fn(mut self, other: &'a Loc<N, F>) -> Self::Output { self.$fn_assign(other); self } } - impl<'b, F : Num, const N : usize> $trait<Loc<F,N>> for &'b Loc<F, N> { - type Output = Loc<F, N>; + impl<'b, F: Num, const N: usize> $trait<Loc<N, F>> for &'b Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(self, other : Loc<F, N>) -> Self::Output { + fn $fn(self, other: Loc<N, F>) -> Self::Output { self.map2(&other, |a, b| a.$fn(b)) } } - impl<'a, 'b, F : Num, const N : usize> $trait<&'a Loc<F,N>> for &'b Loc<F, N> { - type Output = Loc<F, N>; + impl<'a, 'b, F: Num, const N: usize> $trait<&'a Loc<N, F>> for &'b Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(self, other : &'a Loc<F, N>) -> Self::Output { + fn $fn(self, other: &'a Loc<N, F>) -> Self::Output { self.map2(other, |a, b| a.$fn(b)) } } - impl<F : Num, const N : usize> $trait_assign<Loc<F,N>> for Loc<F, N> { + impl<F: Num, const N: usize> $trait_assign<Loc<N, F>> for Loc<N, F> { #[inline] - fn $fn_assign(&mut self, other : Loc<F, N>) { + fn $fn_assign(&mut self, other: Loc<N, F>) { self.map2_mut(&other, |a, b| a.$fn_assign(b)) } } - impl<'a, F : Num, const N : usize> $trait_assign<&'a Loc<F,N>> for Loc<F, N> { + impl<'a, F: Num, const N: usize> $trait_assign<&'a Loc<N, F>> for Loc<N, F> { #[inline] - fn $fn_assign(&mut self, other : &'a Loc<F, N>) { + fn $fn_assign(&mut self, other: &'a Loc<N, F>) { self.map2_mut(other, |a, b| a.$fn_assign(b)) } } - } + }; } make_binop!(Add, add, AddAssign, add_assign); make_binop!(Sub, sub, SubAssign, sub_assign); -impl<F : Float, const N : usize> std::iter::Sum for Loc<F, N> { - fn sum<I: Iterator<Item = Loc<F, N>>>(mut iter: I) -> Self { +impl<F: Float, const N: usize> std::iter::Sum for Loc<N, F> { + fn sum<I: Iterator<Item = Loc<N, F>>>(mut iter: I) -> Self { match iter.next() { None => Self::ORIGIN, Some(mut v) => { @@ -278,62 +282,61 @@ macro_rules! make_scalarop_rhs { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { - impl<F : Num, const N : usize> $trait<F> for Loc<F, N> { - type Output = Loc<F, N>; + impl<F: Num, const N: usize> $trait<F> for Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(self, b : F) -> Self::Output { + fn $fn(self, b: F) -> Self::Output { self.map(|a| a.$fn(b)) } } - impl<'a, F : Num, const N : usize> $trait<&'a F> for Loc<F, N> { - type Output = Loc<F, N>; + impl<'a, F: Num, const N: usize> $trait<&'a F> for Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(self, b : &'a F) -> Self::Output { + fn $fn(self, b: &'a F) -> Self::Output { self.map(|a| a.$fn(*b)) } } - impl<'b, F : Num, const N : usize> $trait<F> for &'b Loc<F, N> { - type Output = Loc<F, N>; + impl<'b, F: Num, const N: usize> $trait<F> for &'b Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(self, b : F) -> Self::Output { + fn $fn(self, b: F) -> Self::Output { self.map(|a| a.$fn(b)) } } - impl<'a, 'b, F : Float, const N : usize> $trait<&'a F> for &'b Loc<F, N> { - type Output = Loc<F, N>; + impl<'a, 'b, F: Float, const N: usize> $trait<&'a F> for &'b Loc<N, F> { + type Output = Loc<N, F>; #[inline] - fn $fn(self, b : &'a F) -> Self::Output { + fn $fn(self, b: &'a F) -> Self::Output { self.map(|a| a.$fn(*b)) } } - impl<F : Num, const N : usize> $trait_assign<F> for Loc<F, N> { + impl<F: Num, const N: usize> $trait_assign<F> for Loc<N, F> { #[inline] - fn $fn_assign(&mut self, b : F) { + fn $fn_assign(&mut self, b: F) { self.map_mut(|a| a.$fn_assign(b)); } } - impl<'a, F : Num, const N : usize> $trait_assign<&'a F> for Loc<F, N> { + impl<'a, F: Num, const N: usize> $trait_assign<&'a F> for Loc<N, F> { #[inline] - fn $fn_assign(&mut self, b : &'a F) { + fn $fn_assign(&mut self, b: &'a F) { self.map_mut(|a| a.$fn_assign(*b)); } } - } + }; } - make_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); make_scalarop_rhs!(Div, div, DivAssign, div_assign); macro_rules! make_unaryop { ($trait:ident, $fn:ident) => { - impl<F : SignedNum, const N : usize> $trait for Loc<F, N> { - type Output = Loc<F, N>; + impl<F: SignedNum, const N: usize> $trait for Loc<N, F> { + type Output = Loc<N, F>; #[inline] fn $fn(mut self) -> Self::Output { self.map_mut(|a| *a = (*a).$fn()); @@ -341,48 +344,48 @@ } } - impl<'a, F : SignedNum, const N : usize> $trait for &'a Loc<F, N> { - type Output = Loc<F, N>; + impl<'a, F: SignedNum, const N: usize> $trait for &'a Loc<N, F> { + type Output = Loc<N, F>; #[inline] fn $fn(self) -> Self::Output { self.map(|a| a.$fn()) } } - } + }; } make_unaryop!(Neg, neg); macro_rules! make_scalarop_lhs { ($trait:ident, $fn:ident; $($f:ident)+) => { $( - impl<const N : usize> $trait<Loc<$f,N>> for $f { - type Output = Loc<$f, N>; + impl<const N : usize> $trait<Loc<N, $f>> for $f { + type Output = Loc<N, $f>; #[inline] - fn $fn(self, v : Loc<$f,N>) -> Self::Output { + fn $fn(self, v : Loc<N, $f>) -> Self::Output { v.map(|b| self.$fn(b)) } } - impl<'a, const N : usize> $trait<&'a Loc<$f,N>> for $f { - type Output = Loc<$f, N>; + impl<'a, const N : usize> $trait<&'a Loc<N, $f>> for $f { + type Output = Loc<N, $f>; #[inline] - fn $fn(self, v : &'a Loc<$f,N>) -> Self::Output { + fn $fn(self, v : &'a Loc<N, $f>) -> Self::Output { v.map(|b| self.$fn(b)) } } - impl<'b, const N : usize> $trait<Loc<$f,N>> for &'b $f { - type Output = Loc<$f, N>; + impl<'b, const N : usize> $trait<Loc<N, $f>> for &'b $f { + type Output = Loc<N, $f>; #[inline] - fn $fn(self, v : Loc<$f,N>) -> Self::Output { + fn $fn(self, v : Loc<N, $f>) -> Self::Output { v.map(|b| self.$fn(b)) } } - impl<'a, 'b, const N : usize> $trait<&'a Loc<$f,N>> for &'b $f { - type Output = Loc<$f, N>; + impl<'a, 'b, const N : usize> $trait<&'a Loc<N, $f>> for &'b $f { + type Output = Loc<N, $f>; #[inline] - fn $fn(self, v : &'a Loc<$f, N>) -> Self::Output { + fn $fn(self, v : &'a Loc<N, $f>) -> Self::Output { v.map(|b| self.$fn(b)) } } @@ -396,21 +399,21 @@ macro_rules! domination { ($norm:ident, $dominates:ident) => { - impl<F : Float, const N : usize> Dominated<F, $dominates, Loc<F, N>> for $norm { + impl<F: Float, const N: usize> Dominated<F, $dominates, Loc<N, F>> for $norm { #[inline] - fn norm_factor(&self, _p : $dominates) -> F { + fn norm_factor(&self, _p: $dominates) -> F { F::ONE } #[inline] - fn from_norm(&self, p_norm : F, _p : $dominates) -> F { + fn from_norm(&self, p_norm: F, _p: $dominates) -> F { p_norm } } }; ($norm:ident, $dominates:ident, $fn:path) => { - impl<F : Float, const N : usize> Dominated<F, $dominates, Loc<F, N>> for $norm { + impl<F: Float, const N: usize> Dominated<F, $dominates, Loc<N, F>> for $norm { #[inline] - fn norm_factor(&self, _p : $dominates) -> F { + fn norm_factor(&self, _p: $dominates) -> F { $fn(F::cast_from(N)) } } @@ -429,16 +432,17 @@ domination!(Linfinity, L2); domination!(L2, L1); -impl<F : Float,const N : usize> Euclidean<F> for Loc<F, N> { +impl<F: Float, const N: usize> Euclidean<F> for Loc<N, F> { type Output = Self; /// This implementation is not stabilised as it's meant to be used for very small vectors. /// Use [`nalgebra`] for larger vectors. #[inline] - fn dot<I : Instance<Self>>(&self, other : I) -> F { - self.0.iter() - .zip(other.ref_instance().0.iter()) - .fold(F::ZERO, |m, (&v, &w)| m + v * w) + fn dot<I: Instance<Self>>(&self, other: I) -> F { + self.0 + .iter() + .zip(other.ref_instance().0.iter()) + .fold(F::ZERO, |m, (&v, &w)| m + v * w) } /// This implementation is not stabilised as it's meant to be used for very small vectors. @@ -448,16 +452,19 @@ self.iter().fold(F::ZERO, |m, &v| m + v * v) } - fn dist2_squared<I : Instance<Self>>(&self, other : I) -> F { + fn dist2_squared<I: Instance<Self>>(&self, other: I) -> F { self.iter() .zip(other.ref_instance().iter()) - .fold(F::ZERO, |m, (&v, &w)| { let d = v - w; m + d * d }) + .fold(F::ZERO, |m, (&v, &w)| { + let d = v - w; + m + d * d + }) } #[inline] fn norm2(&self) -> F { // Optimisation for N==1 that avoids squaring and square rooting. - if N==1 { + if N == 1 { unsafe { self.0.get_unchecked(0) }.abs() } else { self.norm2_squared().sqrt() @@ -465,10 +472,10 @@ } #[inline] - fn dist2<I : Instance<Self>>(&self, other : I) -> F { + fn dist2<I: Instance<Self>>(&self, other: I) -> F { // Optimisation for N==1 that avoids squaring and square rooting. let otherr = other.ref_instance(); - if N==1 { + if N == 1 { unsafe { *self.0.get_unchecked(0) - *otherr.0.get_unchecked(0) }.abs() } else { self.dist2_squared(other).sqrt() @@ -476,28 +483,28 @@ } } -impl<F : Num, const N : usize> Loc<F, N> { +impl<F: Num, const N: usize> Loc<N, F> { /// Origin point - pub const ORIGIN : Self = Loc([F::ZERO; N]); + pub const ORIGIN: Self = Loc([F::ZERO; N]); } -impl<F : Num + std::ops::Neg<Output=F>, const N : usize> Loc<F, N> { +impl<F: Num + std::ops::Neg<Output = F>, const N: usize> Loc<N, F> { /// Reflects along the given coordinate - pub fn reflect(mut self, i : usize) -> Self { + pub fn reflect(mut self, i: usize) -> Self { self[i] = -self[i]; self } /// Reflects along the given coordinate sequences - pub fn reflect_many<I : IntoIterator<Item=usize>>(mut self, idxs : I) -> Self { + pub fn reflect_many<I: IntoIterator<Item = usize>>(mut self, idxs: I) -> Self { for i in idxs { - self[i]=-self[i] + self[i] = -self[i] } self } } -impl<F : std::ops::Neg<Output=F>> Loc<F, 2> { +impl<F: std::ops::Neg<Output = F>> Loc<2, F> { /// Reflect x coordinate pub fn reflect_x(self) -> Self { let Loc([x, y]) = self; @@ -511,18 +518,17 @@ } } -impl<F : Float> Loc<F, 2> { +impl<F: Float> Loc<2, F> { /// Rotate by angle φ - pub fn rotate(self, φ : F) -> Self { + pub fn rotate(self, φ: F) -> Self { let Loc([x, y]) = self; let sin_φ = φ.sin(); let cos_φ = φ.cos(); - [cos_φ * x - sin_φ * y, - sin_φ * x + cos_φ * y].into() + [cos_φ * x - sin_φ * y, sin_φ * x + cos_φ * y].into() } } -impl<F : std::ops::Neg<Output=F>> Loc<F, 3> { +impl<F: std::ops::Neg<Output = F>> Loc<3, F> { /// Reflect x coordinate pub fn reflect_x(self) -> Self { let Loc([x, y, z]) = self; @@ -542,39 +548,32 @@ } } -impl<F : Float> Loc<F, 3> { +impl<F: Float> Loc<3, F> { /// Rotate by angles (yaw, pitch, roll) - pub fn rotate(self, Loc([φ, ψ, θ]) : Self) -> Self { + pub fn rotate(self, Loc([φ, ψ, θ]): Self) -> Self { let Loc([mut x, mut y, mut z]) = self; let sin_φ = φ.sin(); let cos_φ = φ.cos(); - [x, y, z] = [cos_φ * x - sin_φ *y, - sin_φ * x + cos_φ * y, - z]; + [x, y, z] = [cos_φ * x - sin_φ * y, sin_φ * x + cos_φ * y, z]; let sin_ψ = ψ.sin(); let cos_ψ = ψ.cos(); - [x, y, z] = [cos_ψ * x + sin_ψ * z, - y, - -sin_ψ * x + cos_ψ * z]; + [x, y, z] = [cos_ψ * x + sin_ψ * z, y, -sin_ψ * x + cos_ψ * z]; let sin_θ = θ.sin(); let cos_θ = θ.cos(); - [x, y, z] = [x, - cos_θ * y - sin_θ * z, - sin_θ * y + cos_θ * z]; + [x, y, z] = [x, cos_θ * y - sin_θ * z, sin_θ * y + cos_θ * z]; [x, y, z].into() } } -impl<F : Float,const N : usize> StaticEuclidean<F> for Loc<F, N> { +impl<F: Float, const N: usize> StaticEuclidean<F> for Loc<N, F> { #[inline] fn origin() -> Self { Self::ORIGIN } } - /// The default norm for `Loc` is [`L2`]. -impl<F : Float, const N : usize> Normed<F> for Loc<F, N> { +impl<F: Float, const N: usize> Normed<F> for Loc<N, F> { type NormExp = L2; #[inline] @@ -593,22 +592,26 @@ } } -impl<F : Float, const N : usize> HasDual<F> for Loc<F, N> { +impl<F: Float, const N: usize> HasDual<F> for Loc<N, F> { type DualSpace = Self; } -impl<F : Float, const N : usize> Norm<F, L2> for Loc<F, N> { +impl<F: Float, const N: usize> Norm<L2, F> for Loc<N, F> { #[inline] - fn norm(&self, _ : L2) -> F { self.norm2() } + fn norm(&self, _: L2) -> F { + self.norm2() + } } -impl<F : Float, const N : usize> Dist<F, L2> for Loc<F, N> { +impl<F: Float, const N: usize> Dist<F, L2> for Loc<N, F> { #[inline] - fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> F { self.dist2(other) } + fn dist<I: Instance<Self>>(&self, other: I, _: L2) -> F { + self.dist2(other) + } } /* Implemented automatically as Euclidean. -impl<F : Float, const N : usize> Projection<F, L2> for Loc<F, N> { +impl<F : Float, const N : usize> Projection<F, L2> for Loc<N, F> { #[inline] fn proj_ball_mut(&mut self, ρ : F, nrm : L2) { let n = self.norm(nrm); @@ -618,53 +621,53 @@ } }*/ -impl<F : Float, const N : usize> Norm<F, L1> for Loc<F, N> { +impl<F: Float, const N: usize> Norm<L1, F> for Loc<N, F> { /// This implementation is not stabilised as it's meant to be used for very small vectors. /// Use [`nalgebra`] for larger vectors. #[inline] - fn norm(&self, _ : L1) -> F { + fn norm(&self, _: L1) -> F { self.iter().fold(F::ZERO, |m, v| m + v.abs()) } } -impl<F : Float, const N : usize> Dist<F, L1> for Loc<F, N> { +impl<F: Float, const N: usize> Dist<F, L1> for Loc<N, F> { #[inline] - fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> F { + fn dist<I: Instance<Self>>(&self, other: I, _: L1) -> F { self.iter() .zip(other.ref_instance().iter()) - .fold(F::ZERO, |m, (&v, &w)| m + (v-w).abs() ) + .fold(F::ZERO, |m, (&v, &w)| m + (v - w).abs()) } } -impl<F : Float, const N : usize> Projection<F, Linfinity> for Loc<F, N> { +impl<F: Float, const N: usize> Projection<F, Linfinity> for Loc<N, F> { #[inline] - fn proj_ball_mut(&mut self, ρ : F, _ : Linfinity) { - self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) + fn proj_ball_mut(&mut self, ρ: F, _: Linfinity) { + self.iter_mut() + .for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) } } -impl<F : Float, const N : usize> Norm<F, Linfinity> for Loc<F, N> { +impl<F: Float, const N: usize> Norm<Linfinity, F> for Loc<N, F> { /// This implementation is not stabilised as it's meant to be used for very small vectors. /// Use [`nalgebra`] for larger vectors. #[inline] - fn norm(&self, _ : Linfinity) -> F { + fn norm(&self, _: Linfinity) -> F { self.iter().fold(F::ZERO, |m, v| m.max(v.abs())) } } -impl<F : Float, const N : usize> Dist<F, Linfinity> for Loc<F, N> { +impl<F: Float, const N: usize> Dist<F, Linfinity> for Loc<N, F> { #[inline] - fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> F { + fn dist<I: Instance<Self>>(&self, other: I, _: Linfinity) -> F { self.iter() .zip(other.ref_instance().iter()) - .fold(F::ZERO, |m, (&v, &w)| m.max((v-w).abs())) + .fold(F::ZERO, |m, (&v, &w)| m.max((v - w).abs())) } } - // Misc. -impl<A, const N : usize> FixedLength<N> for Loc<A,N> { +impl<A, const N: usize> FixedLength<N> for Loc<N, A> { type Iter = std::array::IntoIter<A, N>; type Elem = A; #[inline] @@ -673,15 +676,18 @@ } } -impl<A, const N : usize> FixedLengthMut<N> for Loc<A,N> { - type IterMut<'a> = std::slice::IterMut<'a, A> where A : 'a; +impl<A, const N: usize> FixedLengthMut<N> for Loc<N, A> { + type IterMut<'a> + = std::slice::IterMut<'a, A> + where + A: 'a; #[inline] fn fl_iter_mut(&mut self) -> Self::IterMut<'_> { self.iter_mut() } } -impl<'a, A, const N : usize> FixedLength<N> for &'a Loc<A,N> { +impl<'a, A, const N: usize> FixedLength<N> for &'a Loc<N, A> { type Iter = std::slice::Iter<'a, A>; type Elem = &'a A; #[inline] @@ -690,38 +696,37 @@ } } - -impl<F : Num, const N : usize> Space for Loc<F, N> { +impl<F: Num, const N: usize> Space for Loc<N, F> { type Decomp = BasicDecomposition; } -impl<F : Float, const N : usize> Mapping<Loc<F, N>> for Loc<F, N> { +impl<F: Float, const N: usize> Mapping<Loc<N, F>> for Loc<N, F> { type Codomain = F; - fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain { + fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain { x.eval(|x̃| self.dot(x̃)) } } -impl<F : Float, const N : usize> Linear<Loc<F, N>> for Loc<F, N> { } +impl<F: Float, const N: usize> Linear<Loc<N, F>> for Loc<N, F> {} -impl<F : Float, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> { +impl<F: Float, const N: usize> AXPY<F, Loc<N, F>> for Loc<N, F> { type Owned = Self; #[inline] - fn axpy<I : Instance<Loc<F, N>>>(&mut self, α : F, x : I, β : F) { + fn axpy<I: Instance<Loc<N, F>>>(&mut self, α: F, x: I, β: F) { x.eval(|x̃| { if β == F::ZERO { - map2_mut(self, x̃, |yi, xi| { *yi = α * (*xi) }) + map2_mut(self, x̃, |yi, xi| *yi = α * (*xi)) } else { - map2_mut(self, x̃, |yi, xi| { *yi = β * (*yi) + α * (*xi) }) + map2_mut(self, x̃, |yi, xi| *yi = β * (*yi) + α * (*xi)) } }) } #[inline] - fn copy_from<I : Instance<Loc<F, N>>>(&mut self, x : I) { - x.eval(|x̃| map2_mut(self, x̃, |yi, xi| *yi = *xi )) + fn copy_from<I: Instance<Loc<N, F>>>(&mut self, x: I) { + x.eval(|x̃| map2_mut(self, x̃, |yi, xi| *yi = *xi)) } #[inline] @@ -734,4 +739,3 @@ *self = Self::ORIGIN; } } -
--- a/src/mapping.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/mapping.rs Thu May 01 13:06:58 2025 -0500 @@ -40,7 +40,7 @@ X: Space, T: Mapping<X, Codomain = Domain>, E: NormExponent, - Domain: Norm<F, E>, + Domain: Norm<E, F>, F: Num, { Composition { @@ -65,19 +65,19 @@ } } -/// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc<F, N>`] to `F`. -pub trait RealMapping<F: Float, const N: usize>: Mapping<Loc<F, N>, Codomain = F> {} +/// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc<N, F>`] to `F`. +pub trait RealMapping<const N: usize, F: Float = f64>: Mapping<Loc<N, F>, Codomain = F> {} -impl<F: Float, T, const N: usize> RealMapping<F, N> for T where T: Mapping<Loc<F, N>, Codomain = F> {} +impl<F: Float, T, const N: usize> RealMapping<N, F> for T where T: Mapping<Loc<N, F>, Codomain = F> {} -/// A helper trait alias for referring to [`Mapping`]s from [`Loc<F, N>`] to [`Loc<F, M>`]. -pub trait RealVectorField<F: Float, const N: usize, const M: usize>: - Mapping<Loc<F, N>, Codomain = Loc<F, M>> +/// A helper trait alias for referring to [`Mapping`]s from [`Loc<N, F>`] to [`Loc<M, F>`]. +pub trait RealVectorField<const N: usize, const M: usize, F: Float = f64>: + Mapping<Loc<N, F>, Codomain = Loc<M, F>> { } -impl<F: Float, T, const N: usize, const M: usize> RealVectorField<F, N, M> for T where - T: Mapping<Loc<F, N>, Codomain = Loc<F, M>> +impl<F: Float, T, const N: usize, const M: usize> RealVectorField<N, M, F> for T where + T: Mapping<Loc<N, F>, Codomain = Loc<M, F>> { } @@ -102,14 +102,14 @@ } /// Automatically implemented shorthand for referring to differentiable [`Mapping`]s from -/// [`Loc<F, N>`] to `F`. -pub trait DifferentiableRealMapping<F: Float, const N: usize>: - DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain = Loc<F, N>> +/// [`Loc<N, F>`] to `F`. +pub trait DifferentiableRealMapping<const N: usize, F: Float>: + DifferentiableMapping<Loc<N, F>, Codomain = F, DerivativeDomain = Loc<N, F>> { } -impl<F: Float, T, const N: usize> DifferentiableRealMapping<F, N> for T where - T: DifferentiableMapping<Loc<F, N>, Codomain = F, DerivativeDomain = Loc<F, N>> +impl<F: Float, T, const N: usize> DifferentiableRealMapping<N, F> for T where + T: DifferentiableMapping<Loc<N, F>, Codomain = F, DerivativeDomain = Loc<N, F>> { } @@ -186,7 +186,7 @@ impl<F: Space, X, G> Mapping<X> for FlattenedCodomain<X, F, G> where X: Space, - G: Mapping<X, Codomain = Loc<F, 1>>, + G: Mapping<X, Codomain = Loc<1, F>>, { type Codomain = F; @@ -198,7 +198,7 @@ /// An auto-trait for constructing a [`FlattenCodomain`] structure for /// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`. -pub trait FlattenCodomain<X: Space, F>: Mapping<X, Codomain = Loc<F, 1>> + Sized { +pub trait FlattenCodomain<X: Space, F>: Mapping<X, Codomain = Loc<1, F>> + Sized { /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`. fn flatten_codomain(self) -> FlattenedCodomain<X, F, Self> { FlattenedCodomain { @@ -208,9 +208,9 @@ } } -impl<X: Space, F, G: Sized + Mapping<X, Codomain = Loc<F, 1>>> FlattenCodomain<X, F> for G {} +impl<X: Space, F, G: Sized + Mapping<X, Codomain = Loc<1, F>>> FlattenCodomain<X, F> for G {} -/// Container for dimensional slicing [`Loc`]`<F, N>` codomain of a [`Mapping`] to `F`. +/// Container for dimensional slicing [`Loc`]`<N, F>` codomain of a [`Mapping`] to `F`. pub struct SlicedCodomain<'a, X, F, G: Clone, const N: usize> { g: Cow<'a, G>, slice: usize, @@ -221,7 +221,7 @@ where X: Space, F: Copy + Space, - G: Mapping<X, Codomain = Loc<F, N>> + Clone, + G: Mapping<X, Codomain = Loc<N, F>> + Clone, { type Codomain = F; @@ -235,8 +235,8 @@ /// An auto-trait for constructing a [`FlattenCodomain`] structure for /// flattening the codomain of a [`Mapping`] from [`Loc`]`<F, 1>` to `F`. -pub trait SliceCodomain<X: Space, F: Copy, const N: usize>: - Mapping<X, Codomain = Loc<F, N>> + Clone + Sized +pub trait SliceCodomain<X: Space, const N: usize, F: Copy = f64>: + Mapping<X, Codomain = Loc<N, F>> + Clone + Sized { /// Flatten the codomain from [`Loc`]`<F, 1>` to `F`. fn slice_codomain(self, slice: usize) -> SlicedCodomain<'static, X, F, Self, N> { @@ -259,8 +259,8 @@ } } -impl<X: Space, F: Copy, G: Sized + Mapping<X, Codomain = Loc<F, N>> + Clone, const N: usize> - SliceCodomain<X, F, N> for G +impl<X: Space, F: Copy, G: Sized + Mapping<X, Codomain = Loc<N, F>> + Clone, const N: usize> + SliceCodomain<X, N, F> for G { }
--- a/src/mapping/dataterm.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/mapping/dataterm.rs Thu May 01 13:06:58 2025 -0500 @@ -70,7 +70,7 @@ } } -//+ AdjointProductBoundedBy<RNDM<F, N>, P, FloatType = F>, +//+ AdjointProductBoundedBy<RNDM<N, F>, P, FloatType = F>, impl<F, X, A, G> Mapping<X> for DataTerm<F, X, A, G> where
--- a/src/nalgebra_support.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/nalgebra_support.rs Thu May 01 13:06:58 2025 -0500 @@ -8,95 +8,118 @@ [`num_traits`] does. */ +use crate::euclidean::*; +use crate::instance::Instance; +use crate::linops::*; +use crate::mapping::{BasicDecomposition, Space}; +use crate::norms::*; +use crate::types::Float; +use nalgebra::base::allocator::Allocator; +use nalgebra::base::constraint::{SameNumberOfColumns, SameNumberOfRows, ShapeConstraint}; +use nalgebra::base::dimension::*; use nalgebra::{ - Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar, - ClosedAddAssign, ClosedMulAssign, SimdComplexField, Vector, OVector, RealField, - LpNorm, UniformNorm + ClosedAddAssign, ClosedMulAssign, DefaultAllocator, Dim, LpNorm, Matrix, OMatrix, OVector, + RealField, Scalar, SimdComplexField, Storage, StorageMut, UniformNorm, Vector, }; -use nalgebra::base::constraint::{ - ShapeConstraint, SameNumberOfRows, SameNumberOfColumns -}; -use nalgebra::base::dimension::*; -use nalgebra::base::allocator::Allocator; +use num_traits::identities::{One, Zero}; use std::ops::Mul; -use num_traits::identities::{Zero, One}; -use crate::linops::*; -use crate::euclidean::*; -use crate::mapping::{Space, BasicDecomposition}; -use crate::types::Float; -use crate::norms::*; -use crate::instance::Instance; -impl<SM,N,M,E> Space for Matrix<E,N,M,SM> +impl<SM, N, M, E> Space for Matrix<E, N, M, SM> where - SM: Storage<E,N,M> + Clone, - N : Dim, M : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, - DefaultAllocator : Allocator<N,M>, + SM: Storage<E, N, M> + Clone, + N: Dim, + M: Dim, + E: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator: Allocator<N, M>, { type Decomp = BasicDecomposition; } -impl<SM,SV,N,M,K,E> Mapping<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, - DefaultAllocator : Allocator<N,K>, - DefaultAllocator : Allocator<M,K>, - DefaultAllocator : Allocator<N,M>, - DefaultAllocator : Allocator<M,N> { - type Codomain = OMatrix<E,N,K>; +impl<SM, SV, N, M, K, E> Mapping<Matrix<E, M, K, SV>> for Matrix<E, N, M, SM> +where + SM: Storage<E, N, M>, + SV: Storage<E, M, K> + Clone, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator: Allocator<N, K>, + DefaultAllocator: Allocator<M, K>, + DefaultAllocator: Allocator<N, M>, + DefaultAllocator: Allocator<M, N>, +{ + type Codomain = OMatrix<E, N, K>; #[inline] - fn apply<I : Instance<Matrix<E,M,K,SV>>>( - &self, x : I - ) -> Self::Codomain { + fn apply<I: Instance<Matrix<E, M, K, SV>>>(&self, x: I) -> Self::Codomain { x.either(|owned| self.mul(owned), |refr| self.mul(refr)) } } - -impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, - DefaultAllocator : Allocator<N,K>, - DefaultAllocator : Allocator<M,K>, - DefaultAllocator : Allocator<N,M>, - DefaultAllocator : Allocator<M,N> { +impl<'a, SM, SV, N, M, K, E> Linear<Matrix<E, M, K, SV>> for Matrix<E, N, M, SM> +where + SM: Storage<E, N, M>, + SV: Storage<E, M, K> + Clone, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator: Allocator<N, K>, + DefaultAllocator: Allocator<M, K>, + DefaultAllocator: Allocator<N, M>, + DefaultAllocator: Allocator<M, N>, +{ } -impl<SM,SV1,SV2,N,M,K,E> GEMV<E, Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: StorageMut<E,N,K>, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float, - DefaultAllocator : Allocator<N,K>, - DefaultAllocator : Allocator<M,K>, - DefaultAllocator : Allocator<N,M>, - DefaultAllocator : Allocator<M,N> { - +impl<SM, SV1, SV2, N, M, K, E> GEMV<E, Matrix<E, M, K, SV1>, Matrix<E, N, K, SV2>> + for Matrix<E, N, M, SM> +where + SM: Storage<E, N, M>, + SV1: Storage<E, M, K> + Clone, + SV2: StorageMut<E, N, K>, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + Float, + DefaultAllocator: Allocator<N, K>, + DefaultAllocator: Allocator<M, K>, + DefaultAllocator: Allocator<N, M>, + DefaultAllocator: Allocator<M, N>, +{ #[inline] - fn gemv<I : Instance<Matrix<E,M,K,SV1>>>( - &self, y : &mut Matrix<E,N,K,SV2>, α : E, x : I, β : E + fn gemv<I: Instance<Matrix<E, M, K, SV1>>>( + &self, + y: &mut Matrix<E, N, K, SV2>, + α: E, + x: I, + β: E, ) { x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β)) } #[inline] - fn apply_mut<'a, I : Instance<Matrix<E,M,K,SV1>>>(&self, y : &mut Matrix<E,N,K,SV2>, x : I) { + fn apply_mut<'a, I: Instance<Matrix<E, M, K, SV1>>>(&self, y: &mut Matrix<E, N, K, SV2>, x: I) { x.eval(|x̃| self.mul_to(x̃, y)) } } -impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM> -where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone, - M : Dim, E : Scalar + Zero + One + Float, - DefaultAllocator : Allocator<M> { +impl<SM, SV1, M, E> AXPY<E, Vector<E, M, SV1>> for Vector<E, M, SM> +where + SM: StorageMut<E, M> + Clone, + SV1: Storage<E, M> + Clone, + M: Dim, + E: Scalar + Zero + One + Float, + DefaultAllocator: Allocator<M>, +{ type Owned = OVector<E, M>; #[inline] - fn axpy<I : Instance<Vector<E,M,SV1>>>(&mut self, α : E, x : I, β : E) { + fn axpy<I: Instance<Vector<E, M, SV1>>>(&mut self, α: E, x: I, β: E) { x.eval(|x̃| Matrix::axpy(self, α, x̃, β)) } #[inline] - fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) { + fn copy_from<I: Instance<Vector<E, M, SV1>>>(&mut self, y: I) { y.eval(|ỹ| Matrix::copy_from(self, ỹ)) } @@ -125,26 +148,40 @@ } }*/ -impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM> -where SM: StorageMut<E,M> + Clone, - M : Dim, E : Scalar + Zero + One + Float + RealField, - DefaultAllocator : Allocator<M> { +impl<SM, M, E> Projection<E, Linfinity> for Vector<E, M, SM> +where + SM: StorageMut<E, M> + Clone, + M: Dim, + E: Scalar + Zero + One + Float + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] - fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) { - self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) + fn proj_ball_mut(&mut self, ρ: E, _: Linfinity) { + self.iter_mut() + .for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ)) } } -impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> -for Matrix<E,N,M,SM> -where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: Storage<E,N,K> + Clone, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + SimdComplexField, - DefaultAllocator : Allocator<N,K>, - DefaultAllocator : Allocator<M,K>, - DefaultAllocator : Allocator<N,M>, - DefaultAllocator : Allocator<M,N> { - type AdjointCodomain = OMatrix<E,M,K>; - type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a; +impl<'own, SV1, SV2, SM, N, M, K, E> Adjointable<Matrix<E, M, K, SV1>, Matrix<E, N, K, SV2>> + for Matrix<E, N, M, SM> +where + SM: Storage<E, N, M>, + SV1: Storage<E, M, K> + Clone, + SV2: Storage<E, N, K> + Clone, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + SimdComplexField, + DefaultAllocator: Allocator<N, K>, + DefaultAllocator: Allocator<M, K>, + DefaultAllocator: Allocator<N, M>, + DefaultAllocator: Allocator<M, N>, +{ + type AdjointCodomain = OMatrix<E, M, K>; + type Adjoint<'a> + = OMatrix<E, M, N> + where + SM: 'a; #[inline] fn adjoint(&self) -> Self::Adjoint<'_> { @@ -160,7 +197,7 @@ m2: &Matrix<T, R2, C2, S2>, ) -> T::SimdRealField where - T: SimdComplexField, + T: SimdComplexField, R1: Dim, C1: Dim, S1: Storage<T, R1, C1>, @@ -177,38 +214,38 @@ // TODO: should allow different input storages in `Euclidean`. -impl<E,M,S> Euclidean<E> -for Vector<E,M,S> -where M : Dim, - S : StorageMut<E,M> + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { - +impl<E, M, S> Euclidean<E> for Vector<E, M, S> +where + M: Dim, + S: StorageMut<E, M> + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ type Output = OVector<E, M>; #[inline] - fn dot<I : Instance<Self>>(&self, other : I) -> E { - Vector::<E,M,S>::dot(self, other.ref_instance()) + fn dot<I: Instance<Self>>(&self, other: I) -> E { + Vector::<E, M, S>::dot(self, other.ref_instance()) } #[inline] fn norm2_squared(&self) -> E { - Vector::<E,M,S>::norm_squared(self) + Vector::<E, M, S>::norm_squared(self) } #[inline] - fn dist2_squared<I : Instance<Self>>(&self, other : I) -> E { + fn dist2_squared<I: Instance<Self>>(&self, other: I) -> E { metric_distance_squared(self, other.ref_instance()) } } -impl<E,M,S> StaticEuclidean<E> -for Vector<E,M,S> -where M : DimName, - S : StorageMut<E,M> + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { - +impl<E, M, S> StaticEuclidean<E> for Vector<E, M, S> +where + M: DimName, + S: StorageMut<E, M> + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] fn origin() -> OVector<E, M> { OVector::zeros() @@ -216,13 +253,13 @@ } /// The default norm for `Vector` is [`L2`]. -impl<E,M,S> Normed<E> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M> + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { - +impl<E, M, S> Normed<E> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M> + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ type NormExp = L2; #[inline] @@ -232,91 +269,95 @@ #[inline] fn is_zero(&self) -> bool { - Vector::<E,M,S>::norm_squared(self) == E::ZERO + Vector::<E, M, S>::norm_squared(self) == E::ZERO } } -impl<E,M,S> HasDual<E> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M> + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { +impl<E, M, S> HasDual<E> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M> + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ // TODO: Doesn't work with different storage formats. type DualSpace = Self; } -impl<E,M,S> Norm<E, L1> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M>, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { - +impl<E, M, S> Norm<L1, E> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M>, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] - fn norm(&self, _ : L1) -> E { + fn norm(&self, _: L1) -> E { nalgebra::Norm::norm(&LpNorm(1), self) } } -impl<E,M,S> Dist<E, L1> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M> + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { +impl<E, M, S> Dist<E, L1> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M> + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] - fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> E { + fn dist<I: Instance<Self>>(&self, other: I, _: L1) -> E { nalgebra::Norm::metric_distance(&LpNorm(1), self, other.ref_instance()) } } -impl<E,M,S> Norm<E, L2> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M>, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { - +impl<E, M, S> Norm<L2, E> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M>, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] - fn norm(&self, _ : L2) -> E { + fn norm(&self, _: L2) -> E { nalgebra::Norm::norm(&LpNorm(2), self) } } -impl<E,M,S> Dist<E, L2> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M> + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { +impl<E, M, S> Dist<E, L2> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M> + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] - fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> E { + fn dist<I: Instance<Self>>(&self, other: I, _: L2) -> E { nalgebra::Norm::metric_distance(&LpNorm(2), self, other.ref_instance()) } } -impl<E,M,S> Norm<E, Linfinity> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M>, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { - +impl<E, M, S> Norm<Linfinity, E> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M>, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] - fn norm(&self, _ : Linfinity) -> E { + fn norm(&self, _: Linfinity) -> E { nalgebra::Norm::norm(&UniformNorm, self) } } -impl<E,M,S> Dist<E, Linfinity> -for Vector<E,M,S> -where M : Dim, - S : Storage<E,M> + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator<M> { +impl<E, M, S> Dist<E, Linfinity> for Vector<E, M, S> +where + M: Dim, + S: Storage<E, M> + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator<M>, +{ #[inline] - fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> E { + fn dist<I: Instance<Self>>(&self, other: I, _: Linfinity) -> E { nalgebra::Norm::metric_distance(&UniformNorm, self, other.ref_instance()) } } @@ -329,15 +370,15 @@ /// 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 { +pub trait ToNalgebraRealField: Float { /// The nalgebra type corresponding to this type. Usually same as `Self`. /// /// This type only carries `nalgebra` traits. - type NalgebraType : RealField; + 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; + type MixedType: RealField + Float; /// Convert to the nalgebra view of `self`. fn to_nalgebra(self) -> Self::NalgebraType; @@ -346,10 +387,10 @@ fn to_nalgebra_mixed(self) -> Self::MixedType; /// Convert from the nalgebra view of `self`. - fn from_nalgebra(t : Self::NalgebraType) -> 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; + fn from_nalgebra_mixed(t: Self::MixedType) -> Self; } impl ToNalgebraRealField for f32 { @@ -357,17 +398,24 @@ type MixedType = f32; #[inline] - fn to_nalgebra(self) -> Self::NalgebraType { self } - - #[inline] - fn to_nalgebra_mixed(self) -> Self::MixedType { self } + fn to_nalgebra(self) -> Self::NalgebraType { + self + } #[inline] - fn from_nalgebra(t : Self::NalgebraType) -> Self { t } + fn to_nalgebra_mixed(self) -> Self::MixedType { + self + } #[inline] - fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } + fn from_nalgebra(t: Self::NalgebraType) -> Self { + t + } + #[inline] + fn from_nalgebra_mixed(t: Self::MixedType) -> Self { + t + } } impl ToNalgebraRealField for f64 { @@ -375,15 +423,22 @@ type MixedType = f64; #[inline] - fn to_nalgebra(self) -> Self::NalgebraType { self } + fn to_nalgebra(self) -> Self::NalgebraType { + self + } #[inline] - fn to_nalgebra_mixed(self) -> Self::MixedType { self } + fn to_nalgebra_mixed(self) -> Self::MixedType { + self + } #[inline] - fn from_nalgebra(t : Self::NalgebraType) -> Self { t } + fn from_nalgebra(t: Self::NalgebraType) -> Self { + t + } #[inline] - fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t } + fn from_nalgebra_mixed(t: Self::MixedType) -> Self { + t + } } -
--- a/src/norms.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/norms.rs Thu May 01 13:06:58 2025 -0500 @@ -92,7 +92,7 @@ /// /// println!("{}, {} {}", x.norm(L1), x.norm(L2), x.norm(Linfinity)) /// ``` -pub trait Norm<F: Num, Exponent: NormExponent> { +pub trait Norm<Exponent: NormExponent, F: Num = f64> { /// Calculate the norm. fn norm(&self, _p: Exponent) -> F; } @@ -110,7 +110,7 @@ } /// Trait for distances with respect to a norm. -pub trait Dist<F: Num, Exponent: NormExponent>: Norm<F, Exponent> + Space { +pub trait Dist<F: Num, Exponent: NormExponent>: Norm<Exponent, F> + Space { /// Calculate the distance fn dist<I: Instance<Self>>(&self, other: I, _p: Exponent) -> F; } @@ -125,7 +125,7 @@ /// /// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity)); /// ``` -pub trait Projection<F: Num, Exponent: NormExponent>: Norm<F, Exponent> + Sized +pub trait Projection<F: Num, Exponent: NormExponent>: Norm<Exponent, F> + Sized where F: Float, { @@ -139,14 +139,14 @@ fn proj_ball_mut(&mut self, ρ: F, q: Exponent); } -/*impl<F : Float, E : Euclidean<F>> Norm<F, L2> for E { +/*impl<F : Float, E : Euclidean<F>> Norm<L2, F> for E { #[inline] fn norm(&self, _p : L2) -> F { self.norm2() } fn dist(&self, other : &Self, _p : L2) -> F { self.dist2(other) } }*/ -impl<F: Float, E: Euclidean<F> + Norm<F, L2>> Projection<F, L2> for E { +impl<F: Float, E: Euclidean<F> + Norm<L2, F>> Projection<F, L2> for E { #[inline] fn proj_ball(self, ρ: F, _p: L2) -> Self { self.proj_ball2(ρ) @@ -176,7 +176,7 @@ } } -impl<F: Float, E: Euclidean<F> + Normed<F, NormExp = L2>> Norm<F, HuberL1<F>> for E { +impl<F: Float, E: Euclidean<F> + Normed<F, NormExp = L2>> Norm<HuberL1<F>, F> for E { fn norm(&self, huber: HuberL1<F>) -> F { huber.apply(self.norm2_squared()) } @@ -188,7 +188,7 @@ } } -// impl<F : Float, E : Norm<F, L2>> Norm<F, L21> for Vec<E> { +// impl<F : Float, E : Norm<L2, F>> Norm<L21, F> for Vec<E> { // fn norm(&self, _l21 : L21) -> F { // self.iter().map(|e| e.norm(L2)).sum() // } @@ -204,7 +204,7 @@ where F: Float, E: NormExponent, - Domain: Space + Norm<F, E>, + Domain: Space + Norm<E, F>, { type Codomain = F; @@ -214,7 +214,7 @@ } } -pub trait Normed<F: Num = f64>: Space + Norm<F, Self::NormExp> { +pub trait Normed<F: Num = f64>: Space + Norm<Self::NormExp, F> { type NormExp: NormExponent; fn norm_exponent(&self) -> Self::NormExp; @@ -283,7 +283,7 @@ impl<C, F, D> Norm<F, Weighted<$exponent, C>> for D where F: Float, - D: Norm<F, $exponent>, + D: Norm<$exponent, F>, C: Constant<Type = F>, { fn norm(&self, e: Weighted<$exponent, C>) -> F {
--- a/src/sets.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/sets.rs Thu May 01 13:06:58 2025 -0500 @@ -2,42 +2,51 @@ This module provides various sets and traits for them. */ -use std::ops::{RangeFull,RangeFrom,Range,RangeInclusive,RangeTo,RangeToInclusive}; -use crate::types::*; +use crate::euclidean::Euclidean; +use crate::instance::{Instance, Space}; use crate::loc::Loc; -use crate::euclidean::Euclidean; -use crate::instance::{Space, Instance}; +use crate::types::*; use serde::Serialize; +use std::ops::{Range, RangeFrom, RangeFull, RangeInclusive, RangeTo, RangeToInclusive}; pub mod cube; pub use cube::Cube; /// Trait for arbitrary sets. The parameter `U` is the element type. -pub trait Set<U> where U : Space { +pub trait Set<U> +where + U: Space, +{ /// Check for element containment - fn contains<I : Instance<U>>(&self, item : I) -> bool; + fn contains<I: Instance<U>>(&self, item: I) -> 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; +pub trait SetOrd: Sized { + /// Returns the smallest set of same class contain both parameters. + fn common(&self, other: &Self) -> Self; - /// Returns the greatest set of same class contaied by n both parameter sets. - fn intersect(&self, other : &Self) -> Option<Self>; + /// Returns the greatest set of same class contaied by n both parameter sets. + fn intersect(&self, other: &Self) -> Option<Self>; } -impl<U, const N : usize> Set<Loc<U, N>> -for Cube<U,N> -where U : Num + PartialOrd + Sized { - fn contains<I : Instance<Loc<U, N>>>(&self, item : I) -> bool { - self.0.iter().zip(item.ref_instance().iter()).all(|(s, x)| s.contains(x)) +impl<U, const N: usize> Set<Loc<N, U>> for Cube<N, U> +where + U: Num + PartialOrd + Sized, +{ + fn contains<I: Instance<Loc<N, U>>>(&self, item: I) -> bool { + self.0 + .iter() + .zip(item.ref_instance().iter()) + .all(|(s, x)| s.contains(x)) } } -impl<U : Space> Set<U> for RangeFull { - fn contains<I : Instance<U>>(&self, _item : I) -> bool { true } +impl<U: Space> Set<U> for RangeFull { + fn contains<I: Instance<U>>(&self, _item: I) -> bool { + true + } } macro_rules! impl_ranges { @@ -56,45 +65,59 @@ )* } } -impl_ranges!(RangeFrom,Range,RangeInclusive,RangeTo,RangeToInclusive); +impl_ranges!(RangeFrom, Range, RangeInclusive, RangeTo, RangeToInclusive); /// 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. -#[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] -pub struct Halfspace<A, F> where A : Euclidean<F>, F : Float { - pub orthogonal : A, - pub offset : F, +#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)] +pub struct Halfspace<A, F> +where + A: Euclidean<F>, + F: Float, +{ + pub orthogonal: A, + pub offset: F, } -impl<A,F> Halfspace<A,F> where A : Euclidean<F>, F : Float { +impl<A, F> Halfspace<A, F> +where + A: Euclidean<F>, + F: Float, +{ #[inline] - pub fn new(orthogonal : A, offset : F) -> Self { - Halfspace{ orthogonal : orthogonal, offset : offset } + pub fn new(orthogonal: A, offset: F) -> Self { + Halfspace { + orthogonal: orthogonal, + offset: offset, + } } } /// Trait for generating a halfspace spanned by another set `Self` of elements of type `U`. -pub trait SpannedHalfspace<F> where F : Float { +pub trait SpannedHalfspace<F> +where + F: Float, +{ /// Type of the orthogonal vector describing the halfspace. - type A : Euclidean<F>; + type A: Euclidean<F>; /// Returns the halfspace spanned by this set. fn spanned_halfspace(&self) -> Halfspace<Self::A, F>; } // TODO: Gram-Schmidt for higher N. -impl<F : Float> SpannedHalfspace<F> for [Loc<F, 1>; 2] { - type A = Loc<F, 1>; +impl<F: Float> SpannedHalfspace<F> for [Loc<1, F>; 2] { + type A = Loc<1, F>; fn spanned_halfspace(&self) -> Halfspace<Self::A, F> { let (x0, x1) = (self[0], self[1]); - Halfspace::new(x1-x0, x0[0]) + Halfspace::new(x1 - x0, x0[0]) } } // TODO: Gram-Schmidt for higher N. -impl<F : Float> SpannedHalfspace<F> for [Loc<F, 2>; 2] { - type A = Loc<F, 2>; +impl<F: Float> SpannedHalfspace<F> for [Loc<2, F>; 2] { + type A = Loc<2, F>; fn spanned_halfspace(&self) -> Halfspace<Self::A, F> { let (x0, x1) = (&self[0], &self[1]); let d = x1 - x0; @@ -104,28 +127,30 @@ } } -impl<A,F> Set<A> for Halfspace<A,F> +impl<A, F> Set<A> for Halfspace<A, F> where - A : Euclidean<F>, - F : Float, + A: Euclidean<F>, + F: Float, { #[inline] - fn contains<I : Instance<A>>(&self, item : I) -> bool { + fn contains<I: Instance<A>>(&self, item: I) -> bool { self.orthogonal.dot(item) >= self.offset } } /// Polygons defined by `N` `Halfspace`s. -#[derive(Clone,Copy,Debug,Eq,PartialEq)] -pub struct NPolygon<A, F, const N : usize>(pub [Halfspace<A,F>; N]) -where A : Euclidean<F>, F : Float; +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub struct NPolygon<A, const N: usize, F = f64>(pub [Halfspace<A, F>; N]) +where + A: Euclidean<F>, + F: Float; -impl<A,F,const N : usize> Set<A> for NPolygon<A,F,N> +impl<A, F, const N: usize> Set<A> for NPolygon<A, N, F> where - A : Euclidean<F>, - F : Float, + A: Euclidean<F>, + F: Float, { - fn contains<I : Instance<A>>(&self, item : I) -> bool { + fn contains<I: Instance<A>>(&self, item: I) -> bool { let r = item.ref_instance(); self.0.iter().all(|halfspace| halfspace.contains(r)) }
--- a/src/sets/cube.rs Thu May 01 08:40:33 2025 -0500 +++ b/src/sets/cube.rs Thu May 01 13:06:58 2025 -0500 @@ -16,25 +16,19 @@ ``` */ -use serde::ser::{Serialize, Serializer, SerializeTupleStruct}; -use crate::types::*; use crate::loc::Loc; +use crate::maputil::{map1, map1_indexed, map2, FixedLength, FixedLengthMut}; use crate::sets::SetOrd; -use crate::maputil::{ - FixedLength, - FixedLengthMut, - map1, - map1_indexed, - map2, -}; +use crate::types::*; +use serde::ser::{Serialize, SerializeTupleStruct, Serializer}; /// 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]); +pub struct Cube<const N: usize, U: Num = f64>(pub(super) [[U; 2]; N]); // Need to manually implement as [F; N] serialisation is provided only for some N. -impl<F : Num + Serialize, const N : usize> Serialize for Cube<F, N> +impl<F: Num + Serialize, const N: usize> Serialize for Cube<N, F> where F: Serialize, { @@ -50,7 +44,7 @@ } } -impl<A : Num, const N : usize> FixedLength<N> for Cube<A,N> { +impl<A: Num, const N: usize> FixedLength<N> for Cube<N, A> { type Iter = std::array::IntoIter<[A; 2], N>; type Elem = [A; 2]; #[inline] @@ -59,7 +53,7 @@ } } -impl<A : Num, const N : usize> FixedLengthMut<N> for Cube<A,N> { +impl<A: Num, const N: usize> FixedLengthMut<N> for Cube<N, A> { type IterMut<'a> = std::slice::IterMut<'a, [A; 2]>; #[inline] fn fl_iter_mut(&mut self) -> Self::IterMut<'_> { @@ -67,7 +61,7 @@ } } -impl<'a, A : Num, const N : usize> FixedLength<N> for &'a Cube<A,N> { +impl<'a, A: Num, const N: usize> FixedLength<N> for &'a Cube<N, A> { type Iter = std::slice::Iter<'a, [A; 2]>; type Elem = &'a [A; 2]; #[inline] @@ -76,15 +70,14 @@ } } - /// Iterator for [`Cube`] corners. -pub struct CubeCornersIter<'a, U : Num, const N : usize> { - index : usize, - cube : &'a Cube<U, N>, +pub struct CubeCornersIter<'a, U: Num, const N: usize> { + index: usize, + cube: &'a Cube<N, U>, } -impl<'a, U : Num, const N : usize> Iterator for CubeCornersIter<'a, U, N> { - type Item = Loc<U, N>; +impl<'a, U: Num, const N: usize> Iterator for CubeCornersIter<'a, U, N> { + type Item = Loc<N, U>; #[inline] fn next(&mut self) -> Option<Self::Item> { if self.index >= N { @@ -92,28 +85,30 @@ } else { let i = self.index; self.index += 1; - let arr = self.cube.map_indexed(|k, a, b| if (i>>k)&1 == 0 { a } else { b }); + let arr = self + .cube + .map_indexed(|k, a, b| if (i >> k) & 1 == 0 { a } else { b }); Some(arr.into()) } } } -impl<U : Num, const N : usize> Cube<U, N> { +impl<U: Num, const N: usize> Cube<N, U> { /// 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] { + 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] { + 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 + /// 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]> { @@ -122,27 +117,27 @@ /// Returns the “start” coordinate $a_i$ of the cube $∏_{i=1}^N [a_i, b_i)$. #[inline] - pub fn start(&self, i : usize) -> U { + 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 { + 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> { + pub fn span_start(&self) -> Loc<N, U> { 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> { + pub fn span_end(&self) -> Loc<N, U> { Loc::new(self.map(|_a, b| b)) } @@ -150,19 +145,22 @@ /// $∏_{i=1}^N [a_i, b_i)$. #[inline] pub fn iter_corners(&self) -> CubeCornersIter<'_, U, N> { - CubeCornersIter{ index : 0, cube : self } + 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)) + pub fn width(&self) -> Loc<N, U> { + 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 { + pub fn shift(&self, shift: &Loc<N, U>) -> Self { let mut cube = self.clone(); for i in 0..N { cube.0[i][0] += shift[i]; @@ -173,144 +171,158 @@ /// Creates a new cube from an array. #[inline] - pub fn new(data : [[U; 2]; N]) -> Self { + pub fn new(data: [[U; 2]; N]) -> Self { Cube(data) } } -impl<F : Float, const N : usize> Cube<F, N> { +impl<F: Float, const N: usize> Cube<N, F> { /// Returns the centre of the cube - pub fn center(&self) -> Loc<F, N> { + pub fn center(&self) -> Loc<N, F> { map1(self, |&[a, b]| (a + b) / F::TWO).into() } } -impl<U : Num> Cube<U, 1> { +impl<U: Num> Cube<1, U> { /// Get the corners of the cube. /// /// TODO: generic implementation once const-generics can be involved in /// calculations. #[inline] - pub fn corners(&self) -> [Loc<U, 1>; 2] { + pub fn corners(&self) -> [Loc<1, U>; 2] { let [[a, b]] = self.0; [a.into(), b.into()] } } -impl<U : Num> Cube<U, 2> { +impl<U: Num> Cube<2, U> { /// Get the corners of the cube in counter-clockwise order. /// /// TODO: generic implementation once const-generics can be involved in /// calculations. #[inline] - pub fn corners(&self) -> [Loc<U, 2>; 4] { - let [[a1, b1], [a2, b2]]=self.0; - [[a1, a2].into(), - [b1, a2].into(), - [b1, b2].into(), - [a1, b2].into()] + pub fn corners(&self) -> [Loc<2, U>; 4] { + let [[a1, b1], [a2, b2]] = self.0; + [ + [a1, a2].into(), + [b1, a2].into(), + [b1, b2].into(), + [a1, b2].into(), + ] } } -impl<U : Num> Cube<U, 3> { +impl<U: Num> Cube<3, U> { /// Get the corners of the cube. /// /// TODO: generic implementation once const-generics can be involved in /// calculations. #[inline] - pub fn corners(&self) -> [Loc<U, 3>; 8] { - let [[a1, b1], [a2, b2], [a3, b3]]=self.0; - [[a1, a2, a3].into(), - [b1, a2, a3].into(), - [b1, b2, a3].into(), - [a1, b2, a3].into(), - [a1, b2, b3].into(), - [b1, b2, b3].into(), - [b1, a2, b3].into(), - [a1, a2, b3].into()] + pub fn corners(&self) -> [Loc<3, U>; 8] { + let [[a1, b1], [a2, b2], [a3, b3]] = self.0; + [ + [a1, a2, a3].into(), + [b1, a2, a3].into(), + [b1, b2, a3].into(), + [a1, b2, a3].into(), + [a1, b2, b3].into(), + [b1, b2, b3].into(), + [b1, a2, b3].into(), + [a1, a2, b3].into(), + ] } } // TODO: Implement Add and Sub of Loc to Cube, and Mul and Div by U : Num. -impl<U : Num, const N : usize> From<[[U; 2]; N]> for Cube<U, N> { +impl<U: Num, const N: usize> From<[[U; 2]; N]> for Cube<N, U> { #[inline] - fn from(data : [[U; 2]; N]) -> Self { + fn from(data: [[U; 2]; N]) -> Self { Cube(data) } } -impl<U : Num, const N : usize> From<Cube<U, N>> for [[U; 2]; N] { +impl<U: Num, const N: usize> From<Cube<N, U>> for [[U; 2]; N] { #[inline] - fn from(Cube(data) : Cube<U, N>) -> Self { + fn from(Cube(data): Cube<N, U>) -> Self { data } } - -impl<U, const N : usize> Cube<U, N> where U : Num + PartialOrd { +impl<U, const N: usize> Cube<N, U> +where + U: Num + PartialOrd, +{ /// Checks whether the cube is non-degenerate, i.e., the start coordinate /// of each axis is strictly less than the end coordinate. #[inline] pub fn nondegenerate(&self) -> bool { self.0.iter().all(|range| range[0] < range[1]) } - + /// Checks whether the cube intersects some `other` cube. /// Matching boundary points are not counted, so `U` is ideally a [`Float`]. #[inline] - pub fn intersects(&self, other : &Cube<U, N>) -> bool { - self.iter_coords().zip(other.iter_coords()).all(|([a1, b1], [a2, b2])| { - a1 < b2 && a2 < b1 - }) + pub fn intersects(&self, other: &Cube<N, U>) -> bool { + self.iter_coords() + .zip(other.iter_coords()) + .all(|([a1, b1], [a2, b2])| a1 < b2 && a2 < b1) } /// Checks whether the cube contains some `other` cube. - pub fn contains_set(&self, other : &Cube<U, N>) -> bool { - self.iter_coords().zip(other.iter_coords()).all(|([a1, b1], [a2, b2])| { - a1 <= a2 && b1 >= b2 - }) + pub fn contains_set(&self, other: &Cube<N, U>) -> bool { + self.iter_coords() + .zip(other.iter_coords()) + .all(|([a1, b1], [a2, b2])| a1 <= a2 && b1 >= b2) } /// Produces the point of minimum $ℓ^p$-norm within the cube `self` for any $p$-norm. /// This is the point where each coordinate is closest to zero. #[inline] - pub fn minnorm_point(&self) -> Loc<U, N> { + pub fn minnorm_point(&self) -> Loc<N, U> { let z = U::ZERO; // As always, we assume that a ≤ b. self.map(|a, b| { debug_assert!(a <= b); match (a < z, z < b) { - (false, _) => a, - (_, false) => b, - (true, true) => z + (false, _) => a, + (_, false) => b, + (true, true) => z, } - }).into() + }) + .into() } /// Produces the point of maximum $ℓ^p$-norm within the cube `self` for any $p$-norm. /// This is the point where each coordinate is furthest from zero. #[inline] - pub fn maxnorm_point(&self) -> Loc<U, N> { + pub fn maxnorm_point(&self) -> Loc<N, U> { let z = U::ZERO; // As always, we assume that a ≤ b. self.map(|a, b| { debug_assert!(a <= b); match (a < z, z < b) { - (false, _) => b, - (_, false) => a, + (false, _) => b, + (_, false) => a, // A this stage we must have a < 0 (so U must be signed), and want to check // whether |a| > |b|. We can do this without assuming U to actually implement // `Neg` by comparing whether 0 > a + b. - (true, true) => if z > a + b { a } else { b } + (true, true) => { + if z > a + b { + a + } else { + b + } + } } - }).into() + }) + .into() } } macro_rules! impl_common { ($($t:ty)*, $min:ident, $max:ident) => { $( - impl<const N : usize> SetOrd for Cube<$t, N> { + impl<const N : usize> SetOrd for Cube<N, $t> { #[inline] fn common(&self, other : &Self) -> Self { map2(self, other, |&[a1, b1], &[a2, b2]| { @@ -338,7 +350,7 @@ #[cfg(feature = "nightly")] impl_common!(f32 f64, minimum, maximum); -impl<U : Num, const N : usize> std::ops::Index<usize> for Cube<U, N> { +impl<U: Num, const N: usize> std::ops::Index<usize> for Cube<N, U> { type Output = [U; 2]; #[inline] fn index(&self, index: usize) -> &Self::Output { @@ -346,7 +358,7 @@ } } -impl<U : Num, const N : usize> std::ops::IndexMut<usize> for Cube<U, N> { +impl<U: Num, const N: usize> std::ops::IndexMut<usize> for Cube<N, U> { #[inline] fn index_mut(&mut self, index: usize) -> &mut Self::Output { &mut self.0[index]
