# HG changeset patch # User Tuomo Valkonen # Date 1746122818 18000 # Node ID 6aa955ad8122f2a1b9343b37e3adcb20891f3979 # Parent 495448cca603576ff48a3fc827e515bc0313589a Transpose loc parameters to allow f64 defaults diff -r 495448cca603 -r 6aa955ad8122 src/bisection_tree/bt.rs --- 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 { /// Point for subdivision of the (unstored) [`Cube`] corresponding to the node. - pub(super) branch_at: Loc, + pub(super) branch_at: Loc, /// Subnodes pub(super) nodes: [Node; 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) -> usize { + fn get_node_index(&self, x: &Loc) -> 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) -> &Node { + fn get_node(&self, x: &Loc) -> &Node { &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, - branch_at: Loc, + domain: &'b Cube, + branch_at: Loc, index: usize, } /// Returns the `i`:th subcube of `domain` subdivided at `branch_at`. #[inline] fn get_subcube( - branch_at: &Loc, - domain: &Cube, + branch_at: &Loc, + domain: &Cube, i: usize, -) -> Cube { +) -> Cube { 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; + type Item = Cube; #[inline] fn next(&mut self) -> Option { 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 + Support>( - domain: &Cube, + pub(super) fn new_with + Support>( + domain: &Cube, 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) -> SubcubeIter<'b, F, N, P> { + pub(super) fn iter_subcubes<'b>(&self, domain: &'b Cube) -> 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) + pub(super) fn nodes_and_cubes<'a, 'b>(&'a self, domain : &'b Cube) -> std::iter::Zip>, 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, + domain: &'b Cube, ) -> std::iter::Zip>, 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, + domain: &'smaller Cube, task_budget: TaskBudget<'scope, 'refs>, - guard: impl Fn(&Node, &Cube) -> bool + Send + 'smaller, - mut f: impl for<'a> FnMut(&mut Node, &Cube, TaskBudget<'smaller, 'a>) + guard: impl Fn(&Node, &Cube) -> bool + Send + 'smaller, + mut f: impl for<'a> FnMut(&mut Node, &Cube, 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 + Support>( + pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis + Support>( &mut self, - domain: &Cube, + domain: &Cube, d: D, new_leaf_depth: M, support: &S, @@ -360,11 +360,11 @@ pub(super) fn convert_aggregator( self, generator: &G, - domain: &Cube, + domain: &Cube, ) -> Branches where ANew: Aggregator, - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, { let branch_at = self.branch_at; @@ -387,10 +387,10 @@ pub(super) fn refresh_aggregator<'refs, 'scope, G>( &mut self, generator: &G, - domain: &Cube, + domain: &Cube, task_budget: TaskBudget<'scope, 'refs>, ) where - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, { self.recurse( @@ -422,7 +422,7 @@ /* /// Get leaf data #[inline] - pub(super) fn get_leaf_data(&self, x : &Loc) -> Option<&Vec> { + pub(super) fn get_leaf_data(&self, x : &Loc) -> Option<&Vec> { 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) -> Option> { + pub(super) fn get_leaf_data_iter(&self, x: &Loc) -> Option> { 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 + Support>( + pub(super) fn insert<'refs, 'scope, M: Depth, S: LocalAnalysis + Support>( &mut self, - domain: &Cube, + domain: &Cube, d: D, new_leaf_depth: M, support: &S, @@ -515,11 +515,11 @@ pub(super) fn convert_aggregator( mut self, generator: &G, - domain: &Cube, + domain: &Cube, ) -> Node where ANew: Aggregator, - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, { // 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, + domain: &Cube, task_budget: TaskBudget<'scope, 'refs>, ) where - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, { 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: +pub trait BTImpl: std::fmt::Debug + Clone + GlobalAnalysis { /// 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: BTImpl + type Converted: BTImpl 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 + Support>( + fn insert + Support>( &mut self, d: Self::Data, support: &S, @@ -642,7 +642,7 @@ fn convert_aggregator(self, generator: &G) -> Self::Converted where ANew: Aggregator, - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis; /// Refreshes the aggregator of the three after possible changes to the support generator. @@ -651,19 +651,19 @@ /// into corresponding [`Support`]s. fn refresh_aggregator(&mut self, generator: &G) where - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis; /// Returns an iterator over all [`Self::Data`] items at the point `x` of the domain. - fn iter_at(&self, x: &Loc) -> std::slice::Iter<'_, Self::Data>; + fn iter_at(&self, x: &Loc) -> std::slice::Iter<'_, Self::Data>; /* /// Returns all [`Self::Data`] items at the point `x` of the domain. - fn data_at(&self, x : &Loc) -> Arc>; + fn data_at(&self, x : &Loc) -> Arc>; */ /// Create a new tree on `domain` of indicated `depth`. - fn new(domain: Cube, depth: Self::Depth) -> Self; + fn new(domain: Cube, 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, + pub(super) domain: Cube, /// The toplevel node of the tree pub(super) topnode: >::Node, } @@ -693,7 +693,7 @@ type Node = Node; } - impl BTImpl for BT + impl BTImpl<$n, F> for BT where M : Depth, F : Float, D : 'static + Copy + Send + Sync + std::fmt::Debug, @@ -703,7 +703,7 @@ type Agg = A; type Converted = BT where ANew : Aggregator; - fn insert + Support>( + fn insert + Support< $n, F>>( &mut self, d : D, support : &S @@ -721,7 +721,7 @@ fn convert_aggregator(self, generator : &G) -> Self::Converted where ANew : Aggregator, - G : SupportGenerator, + G : SupportGenerator< $n, F, Id=D>, G::SupportType : LocalAnalysis { let topnode = self.topnode.convert_aggregator(generator, &self.domain); @@ -733,22 +733,22 @@ } fn refresh_aggregator(&mut self, generator : &G) - where G : SupportGenerator, + where G : SupportGenerator< $n, F, Id=Self::Data>, G::SupportType : LocalAnalysis { with_task_budget(|task_budget| self.topnode.refresh_aggregator(generator, &self.domain, task_budget) ) } - /*fn data_at(&self, x : &Loc) -> Arc> { + /*fn data_at(&self, x : &Loc<$n, F>) -> Arc> { self.topnode.get_leaf_data(x).unwrap_or_else(|| Arc::new(Vec::new())) }*/ - fn iter_at(&self, x : &Loc) -> 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, depth : M) -> Self { + fn new(domain : Cube<$n, F>, depth : M) -> Self { BT { depth : depth, domain : domain, diff -r 495448cca603 -r 6aa955ad8122 src/bisection_tree/btfn.rs --- 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`]``, where `F` is the type of floating point numbers, +/// The domain of the function is [`Loc`]``, 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, BT /*: BTImpl*/, const N: usize> /*where G::SupportType : LocalAnalysis*/ +pub struct BTFN, BT /*: BTImpl< N, F>*/, const N: usize> /*where G::SupportType : LocalAnalysis*/ { bt: BT, generator: Arc, @@ -39,18 +39,18 @@ impl Space for BTFN where - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, - BT: BTImpl, + BT: BTImpl, { type Decomp = BasicDecomposition; } impl BTFN where - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, - BT: BTImpl, + BT: BTImpl, { /// 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, depth: BT::Depth, generator: G) -> Self { + pub fn construct(domain: Cube, depth: BT::Depth, generator: G) -> Self { Self::construct_arc(domain, depth, Arc::new(generator)) } - fn construct_arc(domain: Cube, depth: BT::Depth, generator: Arc) -> Self { + fn construct_arc(domain: Cube, depth: BT::Depth, generator: Arc) -> 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(self) -> BTFN, N> where ANew: Aggregator, - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, { BTFN::new_arc(self.bt.convert_aggregator(&*self.generator), self.generator) @@ -130,15 +130,15 @@ impl BTFN where - G: SupportGenerator, + G: SupportGenerator, { /// 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>( + pub fn instantiate>( self, - domain: Cube, + domain: Cube, depth: BTNew::Depth, ) -> BTFN where @@ -157,7 +157,7 @@ impl PreBTFN where - G: SupportGenerator, + G: SupportGenerator, { /// Create a new [`PreBTFN`] with no bisection tree. pub fn new_pre(generator: G) -> Self { @@ -171,14 +171,14 @@ impl BTFN where - G: SupportGenerator, + G: SupportGenerator, G::SupportType: LocalAnalysis, - BT: BTImpl, + BT: BTImpl, { /// Helper function for implementing [`std::ops::Add`]. fn add_another(&self, g2: Arc) -> BTFN, BT, N> where - G2: SupportGenerator, + G2: SupportGenerator, G2::SupportType: LocalAnalysis, { let mut bt = self.bt.clone(); @@ -201,9 +201,9 @@ impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Add> for $lhs - where BT1 : BTImpl, - G1 : SupportGenerator + $($extra_trait)?, - G2 : SupportGenerator, + where BT1 : BTImpl< N, F, Data=usize>, + G1 : SupportGenerator< N, F, Id=usize> + $($extra_trait)?, + G2 : SupportGenerator< N, F, Id=usize>, G1::SupportType : LocalAnalysis, G2::SupportType : LocalAnalysis { type Output = BTFN, BT1, N>; @@ -216,9 +216,9 @@ impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Add<&'b BTFN> for $lhs - where BT1 : BTImpl, - G1 : SupportGenerator + $($extra_trait)?, - G2 : SupportGenerator, + where BT1 : BTImpl< N, F, Data=usize>, + G1 : SupportGenerator< N, F, Id=usize> + $($extra_trait)?, + G2 : SupportGenerator< N, F, Id=usize>, G1::SupportType : LocalAnalysis, G2::SupportType : LocalAnalysis { @@ -239,9 +239,9 @@ impl<'a, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Sub> for $lhs - where BT1 : BTImpl, - G1 : SupportGenerator + $($extra_trait)?, - G2 : SupportGenerator, + where BT1 : BTImpl< N, F, Data=usize>, + G1 : SupportGenerator< N, F, Id=usize> + $($extra_trait)?, + G2 : SupportGenerator< N, F, Id=usize>, G1::SupportType : LocalAnalysis, G2::SupportType : LocalAnalysis { type Output = BTFN, BT1, N>; @@ -258,9 +258,9 @@ impl<'a, 'b, F : Float, G1, G2, BT1, BT2, const N : usize> std::ops::Sub<&'b BTFN> for $lhs - where BT1 : BTImpl, - G1 : SupportGenerator + $($extra_trait)?, - G2 : SupportGenerator + 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, G2::SupportType : LocalAnalysis, &'b G2 : std::ops::Neg { @@ -281,8 +281,8 @@ ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl std::ops::$trait_assign for BTFN where - BT: BTImpl, - G: SupportGenerator, + BT: BTImpl, + G: SupportGenerator, G::SupportType: LocalAnalysis, { #[inline] @@ -294,8 +294,8 @@ impl std::ops::$trait for BTFN where - BT: BTImpl, - G: SupportGenerator, + BT: BTImpl, + G: SupportGenerator, G::SupportType: LocalAnalysis, { type Output = Self; @@ -309,8 +309,8 @@ impl<'a, F: Float, G, BT, const N: usize> std::ops::$trait for &'a BTFN where - BT: BTImpl, - G: SupportGenerator, + BT: BTImpl, + G: SupportGenerator, G::SupportType: LocalAnalysis, &'a G: std::ops::$trait, { @@ -331,8 +331,8 @@ impl std::ops::$trait> 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 std::ops::$trait for BTFN where - BT: BTImpl, - G: SupportGenerator, + BT: BTImpl, + G: SupportGenerator, G::SupportType: LocalAnalysis, { type Output = Self; @@ -387,8 +387,8 @@ /*impl<'a, F : Float, G, BT, const N : usize> std::ops::$trait for &'a BTFN - where BT : BTImpl, - G : SupportGenerator, + where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis, &'a G : std::ops::$trait { type Output = BTFN; @@ -406,16 +406,16 @@ // Apply, Mapping, Differentiate // -impl Mapping> for BTFN +impl Mapping> for BTFN where - BT: BTImpl, - G: SupportGenerator, - G::SupportType: LocalAnalysis + Mapping, Codomain = V>, + BT: BTImpl, + G: SupportGenerator, + G::SupportType: LocalAnalysis + Mapping, Codomain = V>, V: Sum + Space, { type Codomain = V; - fn apply>>(&self, x: I) -> Self::Codomain { + fn apply>>(&self, x: I) -> Self::Codomain { let xc = x.cow(); self.bt .iter_at(&*xc) @@ -424,17 +424,17 @@ } } -impl DifferentiableImpl> for BTFN +impl DifferentiableImpl> for BTFN where - BT: BTImpl, - G: SupportGenerator, + BT: BTImpl, + G: SupportGenerator, G::SupportType: - LocalAnalysis + DifferentiableMapping, DerivativeDomain = V>, + LocalAnalysis + DifferentiableMapping, DerivativeDomain = V>, V: Sum + Space, { type Derivative = V; - fn differential_impl>>(&self, x: I) -> Self::Derivative { + fn differential_impl>>(&self, x: I) -> Self::Derivative { let xc = x.cow(); self.bt .iter_at(&*xc) @@ -449,8 +449,8 @@ impl GlobalAnalysis for BTFN where - BT: BTImpl, - G: SupportGenerator, + BT: BTImpl, + G: SupportGenerator, G::SupportType: LocalAnalysis, { #[inline] @@ -467,8 +467,8 @@ /* impl<'b, X, F : Float, G, BT, const N : usize> Apply<&'b X, F> for BTFN -where BT : BTImpl, - G : SupportGenerator, +where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis, X : for<'a> Apply<&'a BTFN, F> { @@ -480,8 +480,8 @@ impl Apply for BTFN -where BT : BTImpl, - G : SupportGenerator, +where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis, X : for<'a> Apply<&'a BTFN, F> { @@ -493,8 +493,8 @@ impl Linear for BTFN -where BT : BTImpl, - G : SupportGenerator, +where BT : BTImpl< N, F>, + G : SupportGenerator< N, F, Id=BT::Data>, G::SupportType : LocalAnalysis, X : for<'a> Apply<&'a BTFN, F> { type Codomain = F; @@ -512,16 +512,16 @@ fn p2_minimise F>(&self, g: G) -> (U, F); } -impl P2Minimise, F> for Cube { - fn p2_minimise) -> F>(&self, g: G) -> (Loc, F) { +impl P2Minimise, F> for Cube<1, F> { + fn p2_minimise) -> 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 P2Minimise, F> for Cube { - fn p2_minimise) -> F>(&self, g: G) -> (Loc, F) { +impl P2Minimise, F> for Cube<2, F> { + fn p2_minimise) -> 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 Refiner, G, N> for P2Refiner where - Cube: P2Minimise, F>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + Cube: P2Minimise, F>, + G: SupportGenerator, + G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, { - type Result = Option<(Loc, F)>; + type Result = Option<(Loc, F)>; type Sorting = UpperBoundSorting; fn refine( &self, aggregator: &Bounds, - cube: &Cube, + cube: &Cube, 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| { + let g = move |x: &Loc| { let f = move |&d| generator.support_for(d).apply(x); -data.iter().map(f).sum::() }; @@ -669,17 +669,17 @@ impl Refiner, G, N> for P2Refiner where - Cube: P2Minimise, F>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + Cube: P2Minimise, F>, + G: SupportGenerator, + G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, { - type Result = Option<(Loc, F)>; + type Result = Option<(Loc, F)>; type Sorting = LowerBoundSorting; fn refine( &self, aggregator: &Bounds, - cube: &Cube, + cube: &Cube, 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| { + let g = move |x: &Loc| { let f = move |&d| generator.support_for(d).apply(x); data.iter().map(f).sum::() }; @@ -755,7 +755,7 @@ impl Refiner, G, N> for BoundRefiner where - G: SupportGenerator, + G: SupportGenerator, { type Result = bool; type Sorting = UpperBoundSorting; @@ -763,7 +763,7 @@ fn refine( &self, aggregator: &Bounds, - _cube: &Cube, + _cube: &Cube, _data: &[G::Id], _generator: &G, step: usize, @@ -790,7 +790,7 @@ impl Refiner, G, N> for BoundRefiner where - G: SupportGenerator, + G: SupportGenerator, { type Result = bool; type Sorting = UpperBoundSorting; @@ -798,7 +798,7 @@ fn refine( &self, aggregator: &Bounds, - _cube: &Cube, + _cube: &Cube, _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 MinMaxMapping for BTFN +impl MinMaxMapping for BTFN where - BT: BTSearch>, - G: SupportGenerator, - G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, - Cube: P2Minimise, F>, + BT: BTSearch>, + G: SupportGenerator, + G::SupportType: Mapping, Codomain = F> + LocalAnalysis, N>, + Cube: P2Minimise, F>, { - fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc, F) { + fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc, F) { let refiner = P2Refiner { tolerance, max_steps, @@ -863,7 +863,7 @@ bound: F, tolerance: F, max_steps: usize, - ) -> Option<(Loc, F)> { + ) -> Option<(Loc, 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) { + fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc, F) { let refiner = P2Refiner { tolerance, max_steps, @@ -893,7 +893,7 @@ bound: F, tolerance: F, max_steps: usize, - ) -> Option<(Loc, F)> { + ) -> Option<(Loc, F)> { let refiner = P2Refiner { tolerance, max_steps, diff -r 495448cca603 -r 6aa955ad8122 src/bisection_tree/either.rs --- 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( - pub(super) Arc, - pub(super) Arc, -); +#[derive(Debug, Clone)] +pub struct BothGenerators(pub(super) Arc, pub(super) Arc); /// 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 { +#[derive(Debug, Clone)] +pub enum EitherSupport { 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, - G2 : SupportGenerator, - const N : usize + 'a, + F, + G1: SupportGenerator, + G2: SupportGenerator, + const N: usize, > = Chain< - MapF, (usize, EitherSupport)>, - MapZ, usize, (usize, EitherSupport)>, + MapF, (usize, EitherSupport)>, + MapZ, usize, (usize, EitherSupport)>, >; impl BothGenerators { /// Helper for [`all_left_data`]. #[inline] - fn map_left((d, support) : (G1::Id, G1::SupportType)) - -> (usize, EitherSupport) - where G1 : SupportGenerator, - G2 : SupportGenerator { - - let id : usize = d.into(); + fn map_left( + (d, support): (G1::Id, G1::SupportType), + ) -> (usize, EitherSupport) + where + G1: SupportGenerator, + G2: SupportGenerator, + { + let id: usize = d.into(); (id.into(), EitherSupport::Left(support)) } /// Helper for [`all_right_data`]. #[inline] - fn map_right(n0 : &usize, (d, support) : (G2::Id, G2::SupportType)) - -> (usize, EitherSupport) - where G1 : SupportGenerator, - G2 : SupportGenerator { - - let id : usize = d.into(); - ((n0+id).into(), EitherSupport::Right(support)) + fn map_right( + n0: &usize, + (d, support): (G2::Id, G2::SupportType), + ) -> (usize, EitherSupport) + where + G1: SupportGenerator, + G2: SupportGenerator, + { + 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(&self) - -> MapF, (usize, EitherSupport)> - where G1 : SupportGenerator, - G2 : SupportGenerator { + pub(super) fn all_left_data( + &self, + ) -> MapF, (usize, EitherSupport)> + where + G1: SupportGenerator, + G2: SupportGenerator, + { 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(&self) - -> MapZ, usize, (usize, EitherSupport)> - where G1 : SupportGenerator, - G2 : SupportGenerator { + pub(super) fn all_right_data( + &self, + ) -> MapZ, usize, (usize, EitherSupport)> + where + G1: SupportGenerator, + G2: SupportGenerator, + { let n0 = self.0.support_count(); self.1.all_data().mapZ(n0, Self::map_right) } } -impl -SupportGenerator -for BothGenerators -where G1 : SupportGenerator, - G2 : SupportGenerator { - +impl SupportGenerator for BothGenerators +where + G1: SupportGenerator, + G2: SupportGenerator, +{ type Id = usize; - type SupportType = EitherSupport; - type AllDataIter<'a> = BothAllDataIter<'a, F, G1, G2, N> where G1 : 'a, G2 : 'a; + type SupportType = EitherSupport; + 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 Support for EitherSupport -where S1 : Support, - S2 : Support { - +impl Support for EitherSupport +where + S1: Support, + S2: Support, +{ #[inline] - fn support_hint(&self) -> Cube { + fn support_hint(&self) -> Cube { 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) -> bool { + fn in_support(&self, x: &Loc) -> 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) -> [Option; N] { + fn bisection_hint(&self, cube: &Cube) -> [Option; N] { match self { EitherSupport::Left(ref a) => a.bisection_hint(cube), EitherSupport::Right(ref b) => b.bisection_hint(cube), @@ -155,13 +160,14 @@ } } -impl LocalAnalysis for EitherSupport -where A : Aggregator, - S1 : LocalAnalysis, - S2 : LocalAnalysis, { - +impl LocalAnalysis for EitherSupport +where + A: Aggregator, + S1: LocalAnalysis, + S2: LocalAnalysis, +{ #[inline] - fn local_analysis(&self, cube : &Cube) -> A { + fn local_analysis(&self, cube: &Cube) -> A { match self { EitherSupport::Left(ref a) => a.local_analysis(cube), EitherSupport::Right(ref b) => b.local_analysis(cube), @@ -169,11 +175,12 @@ } } -impl GlobalAnalysis for EitherSupport -where A : Aggregator, - S1 : GlobalAnalysis, - S2 : GlobalAnalysis, { - +impl GlobalAnalysis for EitherSupport +where + A: Aggregator, + S1: GlobalAnalysis, + S2: GlobalAnalysis, +{ #[inline] fn global_analysis(&self) -> A { match self { @@ -183,17 +190,17 @@ } } -impl Mapping for EitherSupport +impl Mapping for EitherSupport where - F : Space, - X : Space, - S1 : Mapping, - S2 : Mapping, + F: Space, + X: Space, + S1: Mapping, + S2: Mapping, { type Codomain = F; #[inline] - fn apply>(&self, x : I) -> F { + fn apply>(&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 DifferentiableImpl for EitherSupport +impl DifferentiableImpl for EitherSupport where - O : Space, - X : Space, - S1 : DifferentiableMapping, - S2 : DifferentiableMapping, + O: Space, + X: Space, + S1: DifferentiableMapping, + S2: DifferentiableMapping, { type Derivative = O; #[inline] - fn differential_impl>(&self, x : I) -> O { + fn differential_impl>(&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 - std::ops::$trait_assign - for BothGenerators - where G1 : std::ops::$trait_assign + Clone, - G2 : std::ops::$trait_assign + Clone, { + impl std::ops::$trait_assign for BothGenerators + where + G1: std::ops::$trait_assign + Clone, + G2: std::ops::$trait_assign + 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 - for &'a BothGenerators - where &'a G1 : std::ops::$trait, - &'a G2 : std::ops::$trait { + impl<'a, F: Float, G1, G2> std::ops::$trait for &'a BothGenerators + where + &'a G1: std::ops::$trait, + &'a G2: std::ops::$trait, + { type Output = BothGenerators; #[inline] - fn $fn(self, t : F) -> BothGenerators { - BothGenerators(Arc::new(self.0.$fn(t)), - Arc::new(self.1.$fn(t))) + fn $fn(self, t: F) -> BothGenerators { + 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 std::ops::Neg for BothGenerators -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; #[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 +*/ diff -r 495448cca603 -r 6aa955ad8122 src/bisection_tree/refine.rs --- 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(PhantomData); +pub struct UpperBoundSorting(PhantomData); /// An [`AggregatorSorting`] for [`Bounds`], using the upper/lower bound as the lower/upper key. /// /// See [`UpperBoundSorting`] for the opposite ordering. -pub struct LowerBoundSorting(PhantomData); +pub struct LowerBoundSorting(PhantomData); -impl AggregatorSorting for UpperBoundSorting { +impl AggregatorSorting for UpperBoundSorting { type Agg = Bounds; type Sort = NaNLeast; #[inline] - fn sort_lower(aggregator : &Bounds) -> Self::Sort { NaNLeast(aggregator.lower()) } - - #[inline] - fn sort_upper(aggregator : &Bounds) -> Self::Sort { NaNLeast(aggregator.upper()) } + fn sort_lower(aggregator: &Bounds) -> Self::Sort { + NaNLeast(aggregator.lower()) + } #[inline] - fn bottom() -> Self::Sort { NaNLeast(F::NEG_INFINITY) } + fn sort_upper(aggregator: &Bounds) -> Self::Sort { + NaNLeast(aggregator.upper()) + } + + #[inline] + fn bottom() -> Self::Sort { + NaNLeast(F::NEG_INFINITY) + } } - -impl AggregatorSorting for LowerBoundSorting { +impl AggregatorSorting for LowerBoundSorting { type Agg = Bounds; type Sort = NaNLeast; #[inline] - fn sort_upper(aggregator : &Bounds) -> Self::Sort { NaNLeast(-aggregator.lower()) } + fn sort_upper(aggregator: &Bounds) -> Self::Sort { + NaNLeast(-aggregator.lower()) + } #[inline] - fn sort_lower(aggregator : &Bounds) -> Self::Sort { NaNLeast(-aggregator.upper()) } + fn sort_lower(aggregator: &Bounds) -> 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 { +pub enum RefinerResult { /// 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 : Sync + Send + 'static -where F : Num, - A : Aggregator, - G : SupportGenerator { - +pub trait Refiner: Sync + Send + 'static +where + F: Num, + A: Aggregator, + G: SupportGenerator, +{ /// 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; + type Sorting: AggregatorSorting; /// Determines whether `aggregator` is sufficiently refined within `domain`. /// @@ -124,42 +135,45 @@ /// number of steps is reached. fn refine( &self, - aggregator : &A, - domain : &Cube, - data : &[G::Id], - generator : &G, - step : usize, + aggregator: &A, + domain: &Cube, + data: &[G::Id], + generator: &G, + step: usize, ) -> RefinerResult; /// 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 { +struct RefinementInfo<'a, F, D, A, S, RResult, const N: usize, const P: usize> +where + F: Float, + D: 'static, + A: Aggregator, + S: AggregatorSorting, +{ /// Domain of `node` - cube : Cube, + cube: Cube, /// Node to be refined - node : &'a mut Node, + node: &'a mut Node, /// 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, + sorting: PhantomData, } -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 { - +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, +{ #[inline] - fn with_aggregator(&self, f : impl FnOnce(&A) -> U) -> U { + fn with_aggregator(&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 { - +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, +{ #[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 { - - #[inline] - fn partial_cmp(&self, other : &Self) -> Option { 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 { +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, +{ + #[inline] + fn partial_cmp(&self, other: &Self) -> Option { + 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 { +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, +{ +} +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, +{ #[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

: BranchCount, - A : Aggregator, - S : AggregatorSorting { +struct HeapContainer<'a, F, D, A, S, RResult, const N: usize, const P: usize> +where + F: Float, + D: 'static + Copy, + Const

: BranchCount, + A: Aggregator, + S: AggregatorSorting, +{ /// Priority queue of nodes to be refined - heap : BinaryHeap>, + heap: BinaryHeap>, /// 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, + result: Option, /// 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

: BranchCount, - A : Aggregator, - S : AggregatorSorting { - +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

: BranchCount, + A: Aggregator, + S: AggregatorSorting, +{ /// 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 -Branches -where Const

: BranchCount, - A : Aggregator, - D : 'static + Copy + Send + Sync { - +impl Branches +where + Const

: BranchCount, + 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, - container : &mut HeapContainer<'a, F, D, A, S, RResult, N, P>, - ) where S : AggregatorSorting { + domain: Cube, + container: &mut HeapContainer<'a, F, D, A, S, RResult, N, P>, + ) where + S: AggregatorSorting, + { // 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 -Node -where Const

: BranchCount, - A : Aggregator, - D : 'static + Copy + Send + Sync { - +impl Node +where + Const

: BranchCount, + 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, - refiner : &R, - generator : &G, - container_arc : &'c Arc>>, - step : usize + self: &'a mut Self, + domain: Cube, + refiner: &R, + generator: &G, + container_arc: &'c Arc>>, + step: usize, ) -> Result>> - where R : Refiner, - G : SupportGenerator, - G::SupportType : LocalAnalysis { - + where + R: Refiner, + G: SupportGenerator, + G::SupportType: LocalAnalysis, + { //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 : BTImpl -where F : Float { - +pub trait BTSearch: BTImpl +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, + refiner: R, + generator: &Arc, ) -> Option - where R : Refiner + Sync + Send + 'static, - G : SupportGenerator + Sync + Send + 'static, - G::SupportType : LocalAnalysis; + where + R: Refiner + Sync + Send + 'static, + G: SupportGenerator + Sync + Send + 'static, + G::SupportType: LocalAnalysis; } -fn refinement_loop ( - wakeup : Option>, - refiner : &R, - generator_arc : &Arc, - container_arc : &Arc>>, -) where A : Aggregator, - R : Refiner, - G : SupportGenerator, - G::SupportType : LocalAnalysis, - Const

: BranchCount, - D : 'static + Copy + Sync + Send + std::fmt::Debug { - +fn refinement_loop( + wakeup: Option>, + refiner: &R, + generator_arc: &Arc, + container_arc: &Arc>>, +) where + A: Aggregator, + R: Refiner, + G: SupportGenerator, + G::SupportType: LocalAnalysis, + Const

: BranchCount, + 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 + BTSearch<$n, F> for BT - where //Self : BTImpl, // <== 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, ) -> Option where R : Refiner, - G : SupportGenerator, + G : SupportGenerator< $n, F, Id=D>, G::SupportType : LocalAnalysis { let mut init_container = HeapContainer { heap : BinaryHeap::new(), @@ -620,4 +655,3 @@ } impl_btsearch!(1 2 3 4); - diff -r 495448cca603 -r 6aa955ad8122 src/bisection_tree/support.rs --- 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: Sized + Sync + Send + 'static { +pub trait Support: 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; + fn support_hint(&self) -> Cube; /// Indicate whether `x` is in the support of the function represented by `self`. - fn in_support(&self, x: &Loc) -> bool; + fn in_support(&self, x: &Loc) -> bool; // Indicate whether `cube` is fully in the support of the function represented by `self`. //fn fully_in_support(&self, cube : &Cube) -> bool; @@ -39,13 +39,13 @@ /// The default implementation returns `[None; N]`. #[inline] #[allow(unused_variables)] - fn bisection_hint(&self, cube: &Cube) -> [Option; N] { + fn bisection_hint(&self, cube: &Cube) -> [Option; N] { [None; N] } /// Translate `self` by `x`. #[inline] - fn shift(self, x: Loc) -> Shift { + fn shift(self, x: Loc) -> Shift { 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 { - shift: Loc, +pub struct Shift { + shift: Loc, base_fn: T, } -impl<'a, T, V: Space, F: Float, const N: usize> Mapping> for Shift +impl<'a, T, V: Space, F: Float, const N: usize> Mapping> for Shift where - T: Mapping, Codomain = V>, + T: Mapping, Codomain = V>, { type Codomain = V; #[inline] - fn apply>>(&self, x: I) -> Self::Codomain { + fn apply>>(&self, x: I) -> Self::Codomain { self.base_fn.apply(x.own() - &self.shift) } } -impl<'a, T, V: Space, F: Float, const N: usize> DifferentiableImpl> for Shift +impl<'a, T, V: Space, F: Float, const N: usize> DifferentiableImpl> for Shift where - T: DifferentiableMapping, DerivativeDomain = V>, + T: DifferentiableMapping, DerivativeDomain = V>, { type Derivative = V; #[inline] - fn differential_impl>>(&self, x: I) -> Self::Derivative { + fn differential_impl>>(&self, x: I) -> Self::Derivative { self.base_fn.differential(x.own() - &self.shift) } } -impl<'a, T, F: Float, const N: usize> Support for Shift +impl<'a, T, F: Float, const N: usize> Support for Shift where - T: Support, + T: Support, { #[inline] - fn support_hint(&self) -> Cube { + fn support_hint(&self) -> Cube { self.base_fn.support_hint().shift(&self.shift) } #[inline] - fn in_support(&self, x: &Loc) -> bool { + fn in_support(&self, x: &Loc) -> bool { self.base_fn.in_support(&(x - &self.shift)) } @@ -104,13 +104,13 @@ // } #[inline] - fn bisection_hint(&self, cube: &Cube) -> [Option; N] { + fn bisection_hint(&self, cube: &Cube) -> [Option; 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> for Shift +impl<'a, T, F: Float, const N: usize> GlobalAnalysis> for Shift where T: LocalAnalysis, N>, { @@ -120,20 +120,20 @@ } } -impl<'a, T, F: Float, const N: usize> LocalAnalysis, N> for Shift +impl<'a, T, F: Float, const N: usize> LocalAnalysis, N> for Shift where T: LocalAnalysis, N>, { #[inline] - fn local_analysis(&self, cube: &Cube) -> Bounds { + fn local_analysis(&self, cube: &Cube) -> Bounds { 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 for Shift - where T : Norm { + impl<'a, T, F : Float, const N : usize> Norm<$norm, F> for Shift + 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 for Weighted +impl<'a, T, F: Float, C, const N: usize> Support for Weighted where - T: Support, + T: Support, C: Constant, { #[inline] - fn support_hint(&self) -> Cube { + fn support_hint(&self) -> Cube { self.base_fn.support_hint() } #[inline] - fn in_support(&self, x: &Loc) -> bool { + fn in_support(&self, x: &Loc) -> bool { self.base_fn.in_support(x) } @@ -164,7 +164,7 @@ // } #[inline] - fn bisection_hint(&self, cube: &Cube) -> [Option; N] { + fn bisection_hint(&self, cube: &Cube) -> [Option; N] { self.base_fn.bisection_hint(cube) } } @@ -191,7 +191,7 @@ C: Constant, { #[inline] - fn local_analysis(&self, cube: &Cube) -> Bounds { + fn local_analysis(&self, cube: &Cube) -> Bounds { 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 for Weighted - where T : Norm { + impl<'a, T, F : Float> Norm<$norm, F> for Weighted + 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> for Normalised +impl<'a, T, F: Float, const N: usize> Mapping> for Normalised where - T: Norm + Mapping, Codomain = F>, + T: Norm + Mapping, Codomain = F>, { type Codomain = F; #[inline] - fn apply>>(&self, x: I) -> Self::Codomain { + fn apply>>(&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 for Normalised +impl<'a, T, F: Float, const N: usize> Support for Normalised where - T: Norm + Support, + T: Norm + Support, { #[inline] - fn support_hint(&self) -> Cube { + fn support_hint(&self) -> Cube { self.0.support_hint() } #[inline] - fn in_support(&self, x: &Loc) -> bool { + fn in_support(&self, x: &Loc) -> bool { self.0.in_support(x) } @@ -297,14 +297,14 @@ // } #[inline] - fn bisection_hint(&self, cube: &Cube) -> [Option; N] { + fn bisection_hint(&self, cube: &Cube) -> [Option; N] { self.0.bisection_hint(cube) } } impl<'a, T, F: Float> GlobalAnalysis> for Normalised where - T: Norm + GlobalAnalysis>, + T: Norm + GlobalAnalysis>, { #[inline] fn global_analysis(&self) -> Bounds { @@ -318,10 +318,10 @@ impl<'a, T, F: Float, const N: usize> LocalAnalysis, N> for Normalised where - T: Norm + LocalAnalysis, N>, + T: Norm + LocalAnalysis, N>, { #[inline] - fn local_analysis(&self, cube: &Cube) -> Bounds { + fn local_analysis(&self, cube: &Cube) -> Bounds { 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 for Normalised +impl<'a, T, F: Float> Norm for Normalised where - T: Norm, + T: Norm, { #[inline] fn norm(&self, _: L1) -> F { @@ -347,8 +347,8 @@ macro_rules! impl_normalised_norm { ($($norm:ident)*) => { $( - impl<'a, T, F : Float> Norm for Normalised - where T : Norm + Norm { + impl<'a, T, F : Float> Norm<$norm, F> for Normalised + where T : Norm<$norm, F> + Norm { #[inline] fn norm(&self, n : $norm) -> F { let w = self.0.norm(L1); @@ -361,26 +361,26 @@ impl_normalised_norm!(L2 Linfinity); /* -impl, const N : usize> LocalAnalysis for S { - fn local_analysis(&self, _cube : &Cube) -> NullAggregator { NullAggregator } +impl, const N : usize> LocalAnalysis for S { + fn local_analysis(&self, _cube : &Cube) -> NullAggregator { NullAggregator } } impl, const N : usize> LocalAnalysis, N> for S { #[inline] - fn local_analysis(&self, cube : &Cube) -> Bounds { + fn local_analysis(&self, cube : &Cube) -> Bounds { self.bounds(cube) } }*/ /// Generator of [`Support`]-implementing component functions based on low storage requirement /// [ids][`Self::Id`]. -pub trait SupportGenerator: +pub trait SupportGenerator: MulAssign + DivAssign + Neg + Clone + Sync + Send + 'static { /// The identification type type Id: 'static + Copy; /// The type of the [`Support`] (often also a [`Mapping`]). - type SupportType: 'static + Support; + type SupportType: 'static + Support; /// An iterator over all the [`Support`]s of the generator. type AllDataIter<'a>: Iterator where diff -r 495448cca603 -r 6aa955ad8122 src/bounds.rs --- 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) -> A; + fn local_analysis(&self, cube: &Cube) -> A; } /// Trait for determining the upper and lower bounds of an float-valued [`Mapping`]. @@ -59,12 +59,12 @@ impl>> Bounded for T {} /// A [`RealMapping`] that provides rough bounds as well as minimisation and maximisation. -pub trait MinMaxMapping: RealMapping + Bounded { +pub trait MinMaxMapping: RealMapping + Bounded { /// 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); + fn maximise(&mut self, tolerance: F, max_steps: usize) -> (Loc, 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)>; + ) -> Option<(Loc, 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); + fn minimise(&mut self, tolerance: F, max_steps: usize) -> (Loc, 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)>; + ) -> Option<(Loc, F)>; /// Verify that the mapping has a given upper `bound` within indicated `tolerance`. /// diff -r 495448cca603 -r 6aa955ad8122 src/collection.rs --- 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 { +pub trait Collection: IntoIterator { /// Type of elements of the collection type Element; /// Iterator over references to elements of the collection - type RefsIter<'a> : Iterator where Self : 'a; + type RefsIter<'a>: Iterator + 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 where Self : 'a; + type RefsIterMut<'a>: Iterator + 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 where E); slice_like_collection!([E; N] where E, const N : usize); -slice_like_collection!(Loc where E, const N : usize); +slice_like_collection!(Loc where E, const N : usize); diff -r 495448cca603 -r 6aa955ad8122 src/convex.rs --- 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 Mapping for NormConstraint where - Domain: Space + Norm, + Domain: Space + Norm, F: Float, E: NormExponent, { @@ -126,8 +126,8 @@ where E: HasDualExponent, F: Float, - Domain: HasDual + Norm + Normed, - >::DualSpace: Norm, + Domain: HasDual + Norm + Normed, + >::DualSpace: Norm, { type Conjugate<'a> = NormConstraint @@ -147,8 +147,8 @@ C: Constant, E: HasDualExponent, F: Float, - Domain: HasDual + Norm + Space, - >::DualSpace: Norm, + Domain: HasDual + Norm + Space, + >::DualSpace: Norm, { type Conjugate<'a> = NormConstraint @@ -165,7 +165,7 @@ impl Prox for NormConstraint where - Domain: Space + Norm, + Domain: Space + Norm, E: NormExponent, F: Float, NormProjection: Mapping, diff -r 495448cca603 -r 6aa955ad8122 src/direct_product.rs --- 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 Norm> for Pair +impl Norm, F> for Pair where F: Num, ExpA: NormExponent, ExpB: NormExponent, ExpJ: NormExponent, - A: Norm, - B: Norm, - Loc: Norm, + A: Norm, + B: Norm, + Loc<2, F>: Norm, { fn norm(&self, PairNorm(expa, expb, expj): PairNorm) -> F { Loc([self.0.norm(expa), self.1.norm(expb)]).norm(expj) diff -r 495448cca603 -r 6aa955ad8122 src/discrete_gradient.rs --- 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: Float + nalgebra::RealField, - DVector: Norm, + DVector: Norm, { fn opnorm_bound(&self, _: L2, _: L2) -> DynResult { // Fuck nalgebra. @@ -405,7 +405,7 @@ where B: Discretisation, F: Float + nalgebra::RealField, - DVector: Norm, + DVector: Norm, { fn opnorm_bound(&self, _: L2, _: L2) -> DynResult { // Fuck nalgebra. diff -r 495448cca603 -r 6aa955ad8122 src/euclidean.rs --- 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 : HasDual + Reflexive - + Mul>::Output> + MulAssign - + Div>::Output> + DivAssign - + Add>::Output> - + Sub>::Output> - + for<'b> Add<&'b Self, Output=>::Output> - + for<'b> Sub<&'b Self, Output=>::Output> - + AddAssign + for<'b> AddAssign<&'b Self> - + SubAssign + for<'b> SubAssign<&'b Self> - + Neg>::Output> +pub trait Euclidean: + HasDual + + Reflexive + + Mul>::Output> + + MulAssign + + Div>::Output> + + DivAssign + + Add>::Output> + + Sub>::Output> + + for<'b> Add<&'b Self, Output = >::Output> + + for<'b> Sub<&'b Self, Output = >::Output> + + AddAssign + + for<'b> AddAssign<&'b Self> + + SubAssign + + for<'b> SubAssign<&'b Self> + + Neg>::Output> { - type Output : Euclidean; + type Output: Euclidean; // Inner product - fn dot>(&self, other : I) -> F; + fn dot>(&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>(&self, y : I) -> F; + fn dist2_squared>(&self, y: I) -> F; /// Calculate the 2-distance $\\|x-y\\|_2$, where `self` is $x$. #[inline] - fn dist2>(&self, y : I) -> F { + fn dist2>(&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 : Euclidean { +pub trait StaticEuclidean: Euclidean { /// Returns the origin fn origin() -> >::Output; } diff -r 495448cca603 -r 6aa955ad8122 src/fe_model/p2_local_model.rs --- 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(pub [Loc; D]); +pub struct Simplex(pub [Loc; D]); /// A two-dimensional planar simplex pub type PlanarSimplex = Simplex; /// A real interval @@ -25,26 +25,29 @@ /// Calculates (a+b)/2 #[inline] #[replace_float_literals(F::cast_from(literal))] -pub(crate) fn midpoint(a : &Loc, b : &Loc) -> Loc { - (a+b)/2.0 +pub(crate) fn midpoint(a: &Loc, b: &Loc) -> Loc { + (a + b) / 2.0 } -impl<'a, F : Float> Set> for RealInterval { +impl<'a, F: Float> Set> for RealInterval { #[inline] - fn contains>>(&self, z : I) -> bool { + fn contains>>(&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> for PlanarSimplex { +impl<'a, F: Float> Set> for PlanarSimplex { #[inline] - fn contains>>(&self, z : I) -> bool { + fn contains>>(&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 P2Powers for Loc { - type Output = Loc; - type Full = Loc; - type Diff = Loc, 1>; +impl 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 P2Powers for Loc { - type Output = Loc; - type Full = Loc; - type Diff = Loc, 2>; +impl 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`]``. -pub trait P2Model { +pub trait P2Model { /// Implementation type of the second order polynomical model. /// Typically a [`P2LocalModel`]. - type Model : LocalModel,F>; + type Model: LocalModel, F>; /// Generates a second order polynomial model of the function `g` on `Self`. - fn p2_model) -> F>(&self, g : G) -> Self::Model; + fn p2_model) -> F>(&self, g: G) -> Self::Model; } /// A local second order polynomical model of dimension `N` with `E` edges -pub struct P2LocalModel { - a0 : F, - a1 : Loc, - a2 : Loc, - //node_values : Loc, - //edge_values : Loc, +pub struct P2LocalModel< + F: Num, + const N: usize, + const E: usize, /*, const V : usize, const Q : usize*/ +> { + a0: F, + a1: Loc, + a2: Loc, + //node_values : Loc, + //edge_values : Loc, } // // 1D planar model construction // -impl RealInterval { +impl RealInterval { #[inline] - fn midpoints(&self) -> [Loc; 1] { + fn midpoints(&self) -> [Loc<1, F>; 1] { let [ref n0, ref n1] = &self.0; let n01 = midpoint(n0, n1); [n01] } } -impl P2LocalModel { +impl P2LocalModel { /// 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; 3], - &[v0, v1, v01] : &[F; 3], - ) -> Self { - let p = move |x : &Loc, 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 P2Model for RealInterval { - type Model = P2LocalModel; +impl P2Model for RealInterval { + type Model = P2LocalModel; #[inline] - fn p2_model) -> F>(&self, g : G) -> Self::Model { - let [n01] = self.midpoints(); + fn p2_model) -> 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 PlanarSimplex { +impl PlanarSimplex { #[inline] /// Returns the midpoints of all the edges of the simplex - fn midpoints(&self) -> [Loc; 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 P2LocalModel { +impl P2LocalModel { /// 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; 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, 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 P2Model for PlanarSimplex { - type Model = P2LocalModel; +impl P2Model for PlanarSimplex { + type Model = P2LocalModel; #[inline] - fn p2_model) -> F>(&self, g : G) -> Self::Model { + fn p2_model) -> 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 LocalModel, F> for P2LocalModel { + impl LocalModel, F> for P2LocalModel { #[inline] - fn value(&self, x : &Loc) -> 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) -> Loc { + 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 P2LocalModel { +impl P2LocalModel { /// Minimises the model along the edge `[x0, x1]`. #[inline] - fn minimise_edge(&self, x0 : Loc, x1 : Loc) -> (Loc, 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,Loc,F> -for P2LocalModel { +impl<'a, F: Float> RealLocalModel, Loc<1, F>, F> + for P2LocalModel +{ #[inline] - fn minimise(&self, &Simplex([x0, x1]) : &RealInterval) -> (Loc, F) { + fn minimise(&self, &Simplex([x0, x1]): &RealInterval) -> (Loc<1, F>, F) { self.minimise_edge(x0, x1) } } #[replace_float_literals(F::cast_from(literal))] -impl P2LocalModel { +impl P2LocalModel { /// Minimise the 2D model along the edge `[x0, x1] = {x0 + t(x1 - x0) | t ∈ [0, 1] }`. #[inline] - fn minimise_edge(&self, x0 : &Loc, x1 : &Loc/*, v0 : F, v1 : F*/) -> (Loc, 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,Loc,F> -for P2LocalModel { +impl<'a, F: Float> RealLocalModel, Loc<2, F>, F> + for P2LocalModel +{ #[inline] - fn minimise(&self, el : &PlanarSimplex) -> (Loc, F) { + fn minimise(&self, el: &PlanarSimplex) -> (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,Loc,F> -for P2LocalModel { +impl<'a, F: Float> RealLocalModel, Loc<2, F>, F> + for P2LocalModel +{ #[inline] - fn minimise(&self, el : &Cube) -> (Loc, 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| 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| - (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])]; diff -r 495448cca603 -r 6aa955ad8122 src/lingrid.rs --- 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`]`` for multi-dimensional grids created by `lingrid`. +/// or [`Loc`]`` 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 { - 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 = LinSpace, [usize; N]>; +pub type LinGrid = LinSpace, [usize; N]>; /// Creates a [`LinSpace`] on the real line. -pub fn linspace(start : F, end : F, count : usize) -> LinSpace { - LinSpace{ start : start, end : end, count : count } +pub fn linspace(start: F, end: F, count: usize) -> LinSpace { + 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( - cube : &Cube, - count : &[usize; N] -) -> LinSpace, [usize; N]> { - LinSpace{ start : cube.span_start(), end : cube.span_end(), count : *count } +pub fn lingrid( + cube: &Cube, + count: &[usize; N], +) -> LinSpace, [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( - cube : &Cube, - count : &[usize; N] -) -> LinSpace, [usize; N]> { +pub fn lingrid_centered( + cube: &Cube, + count: &[usize; N], +) -> LinSpace, [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 { - lingrid : LinSpace, - current : Option, + lingrid: LinSpace, + current: Option, } /// Abstraction of a linear grid over space `U` with multi-dimensional index set `I`. pub trait Grid { /// 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 } @@ -97,7 +108,7 @@ fn next_index(&mut self) -> Option; } -impl, I : Unsigned> Grid for LinSpace { +impl, I: Unsigned> Grid for LinSpace { /*fn entry(&self, i : I) -> Option { 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, I : Unsigned> GridIteration -for LinSpaceIterator { +impl, I: Unsigned> GridIteration for LinSpaceIterator { #[inline] fn next_index(&mut self) -> Option { 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, I : Unsigned, const N : usize> Grid, [I; N]> -for LinSpace, [I; N]> { +impl, I: Unsigned, const N: usize> Grid, [I; N]> + for LinSpace, [I; N]> +{ #[inline] - fn entry_linear_unchecked(&self, i_ : usize) -> Loc { + fn entry_linear_unchecked(&self, i_: usize) -> Loc { 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 { - let LinSpace{ start, end, count } = self; + fn entry_unchecked(&self, i: &[I; N]) -> Loc { + 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, I : Unsigned, const N : usize> GridIteration, [I; N]> -for LinSpaceIterator, [I; N]> { - +impl, I: Unsigned, const N: usize> GridIteration, [I; N]> + for LinSpaceIterator, [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 IntoIterator for LinSpace -where LinSpace : Grid, - LinSpaceIterator : GridIteration { +impl IntoIterator for LinSpace +where + LinSpace: Grid, + LinSpaceIterator: GridIteration, +{ type Item = F; - type IntoIter = LinSpaceIterator; + type IntoIter = LinSpaceIterator; #[inline] fn into_iter(self) -> Self::IntoIter { - LinSpaceIterator { lingrid : self, current : None } + LinSpaceIterator { + lingrid: self, + current: None, + } } } -impl Iterator for LinSpaceIterator -where LinSpace : Grid, - LinSpaceIterator : GridIteration { +impl Iterator for LinSpaceIterator +where + LinSpace: Grid, + LinSpaceIterator: GridIteration, +{ type Item = F; #[inline] fn next(&mut self) -> Option { @@ -207,19 +229,24 @@ } } -impl StatefulIterator for LinSpaceIterator -where LinSpace : Grid, - LinSpaceIterator : GridIteration { +impl StatefulIterator for LinSpaceIterator +where + LinSpace: Grid, + LinSpaceIterator: GridIteration, +{ #[inline] fn current(&self) -> Option { - self.current.as_ref().map(|c| self.lingrid.entry_unchecked(c)) + self.current + .as_ref() + .map(|c| self.lingrid.entry_unchecked(c)) } } - -impl RestartableIterator for LinSpaceIterator -where LinSpace : Grid, - LinSpaceIterator : GridIteration { +impl RestartableIterator for LinSpaceIterator +where + LinSpace: Grid, + LinSpaceIterator: GridIteration, +{ #[inline] fn restart(&mut self) -> Option { self.current = None; diff -r 495448cca603 -r 6aa955ad8122 src/linops.rs --- 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: Linear where F: Num, - X: Space + Norm, + X: Space + Norm, XExp: NormExponent, CodExp: NormExponent, { @@ -182,7 +182,7 @@ impl BoundedLinear for IdOp where - X: Space + Clone + Norm, + X: Space + Clone + Norm, F: Num, E: NormExponent, { @@ -264,8 +264,8 @@ impl<'a, F, X, XD, Y, E1, E2> BoundedLinear for ZeroOp<'a, X, XD, Y, F> where - X: Space + Norm, - Y: AXPY + Clone + Norm, + X: Space + Norm, + Y: AXPY + Clone + Norm, F: Num, E1: NormExponent, E2: NormExponent, @@ -346,8 +346,8 @@ impl BoundedLinear for Composition where F: Num, - X: Space + Norm, - Z: Space + Norm, + X: Space + Norm, + Z: Space + Norm, Xexp: NormExponent, Yexp: NormExponent, Zexp: NormExponent, @@ -680,8 +680,8 @@ BoundedLinear, PairNorm, ExpR, F> for RowOp where F: Float, - A: Space + Norm, - B: Space + Norm, + A: Space + Norm, + B: Space + Norm, S: BoundedLinear, T: BoundedLinear, S::Codomain: Add, @@ -707,7 +707,7 @@ for ColOp where F: Float, - A: Space + Norm, + A: Space + Norm, S: BoundedLinear, T: BoundedLinear, ExpA: NormExponent, diff -r 495448cca603 -r 6aa955ad8122 src/loc.rs --- 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( +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +pub struct Loc( /// An array of the elements of the vector - pub [F; N] + pub [F; N], ); -impl Display for Loc{ +impl Display for Loc { // 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 Serialize for Loc +impl Serialize for Loc where F: Serialize, { @@ -56,10 +57,10 @@ } } -impl Loc { +impl Loc { /// 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 Loc { +impl Loc { /// Maps `g` over the elements of the vector, returning a new [`Loc`] vector #[inline] - pub fn map(&self, g : impl Fn(F) -> H) -> Loc { + pub fn map(&self, g: impl Fn(F) -> H) -> Loc { Loc::new(map1(self, |u| g(*u))) } /// Maps `g` over pairs of elements of two vectors, retuning a new one. #[inline] - pub fn map2(&self, other : &Self, g : impl Fn(F, F) -> H) -> Loc { + pub fn map2(&self, other: &Self, g: impl Fn(F, F) -> H) -> Loc { 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(&self, g : impl Fn(F) -> A) -> A { + pub fn product_map(&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 From<[F; N]> for Loc { +impl From<[F; N]> for Loc { #[inline] - fn from(other: [F; N]) -> Loc { + fn from(other: [F; N]) -> Loc { Loc(other) } } -/*impl From<&[F; N]> for Loc { +/*impl From<&[F; N]> for Loc { #[inline] - fn from(other: &[F; N]) -> Loc { + fn from(other: &[F; N]) -> Loc { Loc(*other) } }*/ -impl From for Loc { +impl From for Loc<1, F> { #[inline] - fn from(other: F) -> Loc { + fn from(other: F) -> Loc<1, F> { Loc([other]) } } -impl Loc { +impl Loc<1, F> { #[inline] pub fn flatten1d(self) -> F { let Loc([v]) = self; @@ -160,22 +161,21 @@ } } -impl From> for [F; N] { +impl From> for [F; N] { #[inline] - fn from(other : Loc) -> [F; N] { + fn from(other: Loc) -> [F; N] { other.0 } } -/*impl From<&Loc> for [F; N] { +/*impl From<&Loc> for [F; N] { #[inline] - fn from(other : &Loc) -> [F; N] { + fn from(other : &Loc) -> [F; N] { other.0 } }*/ - -impl IntoIterator for Loc { +impl IntoIterator for Loc { type Item = <[F; N] as IntoIterator>::Item; type IntoIter = <[F; N] as IntoIterator>::IntoIter; @@ -187,20 +187,24 @@ // Indexing -impl Index for Loc -where [F; N] : Index { +impl Index for Loc +where + [F; N]: Index, +{ type Output = <[F; N] as Index>::Output; #[inline] - fn index(&self, ix : Ix) -> &Self::Output { + fn index(&self, ix: Ix) -> &Self::Output { self.0.index(ix) } } -impl IndexMut for Loc -where [F; N] : IndexMut { +impl IndexMut for Loc +where + [F; N]: IndexMut, +{ #[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 $trait> for Loc { - type Output = Loc; + impl $trait> for Loc { + type Output = Loc; #[inline] - fn $fn(mut self, other : Loc) -> Self::Output { + fn $fn(mut self, other: Loc) -> Self::Output { self.$fn_assign(other); self } } - impl<'a, F : Num, const N : usize> $trait<&'a Loc> for Loc { - type Output = Loc; + impl<'a, F: Num, const N: usize> $trait<&'a Loc> for Loc { + type Output = Loc; #[inline] - fn $fn(mut self, other : &'a Loc) -> Self::Output { + fn $fn(mut self, other: &'a Loc) -> Self::Output { self.$fn_assign(other); self } } - impl<'b, F : Num, const N : usize> $trait> for &'b Loc { - type Output = Loc; + impl<'b, F: Num, const N: usize> $trait> for &'b Loc { + type Output = Loc; #[inline] - fn $fn(self, other : Loc) -> Self::Output { + fn $fn(self, other: Loc) -> Self::Output { self.map2(&other, |a, b| a.$fn(b)) } } - impl<'a, 'b, F : Num, const N : usize> $trait<&'a Loc> for &'b Loc { - type Output = Loc; + impl<'a, 'b, F: Num, const N: usize> $trait<&'a Loc> for &'b Loc { + type Output = Loc; #[inline] - fn $fn(self, other : &'a Loc) -> Self::Output { + fn $fn(self, other: &'a Loc) -> Self::Output { self.map2(other, |a, b| a.$fn(b)) } } - impl $trait_assign> for Loc { + impl $trait_assign> for Loc { #[inline] - fn $fn_assign(&mut self, other : Loc) { + fn $fn_assign(&mut self, other: Loc) { self.map2_mut(&other, |a, b| a.$fn_assign(b)) } } - impl<'a, F : Num, const N : usize> $trait_assign<&'a Loc> for Loc { + impl<'a, F: Num, const N: usize> $trait_assign<&'a Loc> for Loc { #[inline] - fn $fn_assign(&mut self, other : &'a Loc) { + fn $fn_assign(&mut self, other: &'a Loc) { 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 std::iter::Sum for Loc { - fn sum>>(mut iter: I) -> Self { +impl std::iter::Sum for Loc { + fn sum>>(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 $trait for Loc { - type Output = Loc; + impl $trait for Loc { + type Output = Loc; #[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 { - type Output = Loc; + impl<'a, F: Num, const N: usize> $trait<&'a F> for Loc { + type Output = Loc; #[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 for &'b Loc { - type Output = Loc; + impl<'b, F: Num, const N: usize> $trait for &'b Loc { + type Output = Loc; #[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 { - type Output = Loc; + impl<'a, 'b, F: Float, const N: usize> $trait<&'a F> for &'b Loc { + type Output = Loc; #[inline] - fn $fn(self, b : &'a F) -> Self::Output { + fn $fn(self, b: &'a F) -> Self::Output { self.map(|a| a.$fn(*b)) } } - impl $trait_assign for Loc { + impl $trait_assign for Loc { #[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 { + impl<'a, F: Num, const N: usize> $trait_assign<&'a F> for Loc { #[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 $trait for Loc { - type Output = Loc; + impl $trait for Loc { + type Output = Loc; #[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 { - type Output = Loc; + impl<'a, F: SignedNum, const N: usize> $trait for &'a Loc { + type Output = Loc; #[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 $trait> for $f { - type Output = Loc<$f, N>; + impl $trait> for $f { + type Output = Loc; #[inline] - fn $fn(self, v : Loc<$f,N>) -> Self::Output { + fn $fn(self, v : Loc) -> 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> for $f { + type Output = Loc; #[inline] - fn $fn(self, v : &'a Loc<$f,N>) -> Self::Output { + fn $fn(self, v : &'a Loc) -> Self::Output { v.map(|b| self.$fn(b)) } } - impl<'b, const N : usize> $trait> for &'b $f { - type Output = Loc<$f, N>; + impl<'b, const N : usize> $trait> for &'b $f { + type Output = Loc; #[inline] - fn $fn(self, v : Loc<$f,N>) -> Self::Output { + fn $fn(self, v : Loc) -> 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> for &'b $f { + type Output = Loc; #[inline] - fn $fn(self, v : &'a Loc<$f, N>) -> Self::Output { + fn $fn(self, v : &'a Loc) -> Self::Output { v.map(|b| self.$fn(b)) } } @@ -396,21 +399,21 @@ macro_rules! domination { ($norm:ident, $dominates:ident) => { - impl Dominated> for $norm { + impl Dominated> 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 Dominated> for $norm { + impl Dominated> 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 Euclidean for Loc { +impl Euclidean for Loc { 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>(&self, other : I) -> F { - self.0.iter() - .zip(other.ref_instance().0.iter()) - .fold(F::ZERO, |m, (&v, &w)| m + v * w) + fn dot>(&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>(&self, other : I) -> F { + fn dist2_squared>(&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>(&self, other : I) -> F { + fn dist2>(&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 Loc { +impl Loc { /// Origin point - pub const ORIGIN : Self = Loc([F::ZERO; N]); + pub const ORIGIN: Self = Loc([F::ZERO; N]); } -impl, const N : usize> Loc { +impl, const N: usize> Loc { /// 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>(mut self, idxs : I) -> Self { + pub fn reflect_many>(mut self, idxs: I) -> Self { for i in idxs { - self[i]=-self[i] + self[i] = -self[i] } self } } -impl> Loc { +impl> Loc<2, F> { /// Reflect x coordinate pub fn reflect_x(self) -> Self { let Loc([x, y]) = self; @@ -511,18 +518,17 @@ } } -impl Loc { +impl 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> Loc { +impl> Loc<3, F> { /// Reflect x coordinate pub fn reflect_x(self) -> Self { let Loc([x, y, z]) = self; @@ -542,39 +548,32 @@ } } -impl Loc { +impl 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 StaticEuclidean for Loc { +impl StaticEuclidean for Loc { #[inline] fn origin() -> Self { Self::ORIGIN } } - /// The default norm for `Loc` is [`L2`]. -impl Normed for Loc { +impl Normed for Loc { type NormExp = L2; #[inline] @@ -593,22 +592,26 @@ } } -impl HasDual for Loc { +impl HasDual for Loc { type DualSpace = Self; } -impl Norm for Loc { +impl Norm for Loc { #[inline] - fn norm(&self, _ : L2) -> F { self.norm2() } + fn norm(&self, _: L2) -> F { + self.norm2() + } } -impl Dist for Loc { +impl Dist for Loc { #[inline] - fn dist>(&self, other : I, _ : L2) -> F { self.dist2(other) } + fn dist>(&self, other: I, _: L2) -> F { + self.dist2(other) + } } /* Implemented automatically as Euclidean. -impl Projection for Loc { +impl Projection for Loc { #[inline] fn proj_ball_mut(&mut self, ρ : F, nrm : L2) { let n = self.norm(nrm); @@ -618,53 +621,53 @@ } }*/ -impl Norm for Loc { +impl Norm for Loc { /// 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 Dist for Loc { +impl Dist for Loc { #[inline] - fn dist>(&self, other : I, _ : L1) -> F { + fn dist>(&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 Projection for Loc { +impl Projection for Loc { #[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 Norm for Loc { +impl Norm for Loc { /// 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 Dist for Loc { +impl Dist for Loc { #[inline] - fn dist>(&self, other : I, _ : Linfinity) -> F { + fn dist>(&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 FixedLength for Loc { +impl FixedLength for Loc { type Iter = std::array::IntoIter; type Elem = A; #[inline] @@ -673,15 +676,18 @@ } } -impl FixedLengthMut for Loc { - type IterMut<'a> = std::slice::IterMut<'a, A> where A : 'a; +impl FixedLengthMut for Loc { + 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 for &'a Loc { +impl<'a, A, const N: usize> FixedLength for &'a Loc { type Iter = std::slice::Iter<'a, A>; type Elem = &'a A; #[inline] @@ -690,38 +696,37 @@ } } - -impl Space for Loc { +impl Space for Loc { type Decomp = BasicDecomposition; } -impl Mapping> for Loc { +impl Mapping> for Loc { type Codomain = F; - fn apply>>(&self, x : I) -> Self::Codomain { + fn apply>>(&self, x: I) -> Self::Codomain { x.eval(|x̃| self.dot(x̃)) } } -impl Linear> for Loc { } +impl Linear> for Loc {} -impl AXPY> for Loc { +impl AXPY> for Loc { type Owned = Self; #[inline] - fn axpy>>(&mut self, α : F, x : I, β : F) { + fn axpy>>(&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>>(&mut self, x : I) { - x.eval(|x̃| map2_mut(self, x̃, |yi, xi| *yi = *xi )) + fn copy_from>>(&mut self, x: I) { + x.eval(|x̃| map2_mut(self, x̃, |yi, xi| *yi = *xi)) } #[inline] @@ -734,4 +739,3 @@ *self = Self::ORIGIN; } } - diff -r 495448cca603 -r 6aa955ad8122 src/mapping.rs --- 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, E: NormExponent, - Domain: Norm, + Domain: Norm, F: Num, { Composition { @@ -65,19 +65,19 @@ } } -/// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc`] to `F`. -pub trait RealMapping: Mapping, Codomain = F> {} +/// Automatically implemented shorthand for referring to [`Mapping`]s from [`Loc`] to `F`. +pub trait RealMapping: Mapping, Codomain = F> {} -impl RealMapping for T where T: Mapping, Codomain = F> {} +impl RealMapping for T where T: Mapping, Codomain = F> {} -/// A helper trait alias for referring to [`Mapping`]s from [`Loc`] to [`Loc`]. -pub trait RealVectorField: - Mapping, Codomain = Loc> +/// A helper trait alias for referring to [`Mapping`]s from [`Loc`] to [`Loc`]. +pub trait RealVectorField: + Mapping, Codomain = Loc> { } -impl RealVectorField for T where - T: Mapping, Codomain = Loc> +impl RealVectorField for T where + T: Mapping, Codomain = Loc> { } @@ -102,14 +102,14 @@ } /// Automatically implemented shorthand for referring to differentiable [`Mapping`]s from -/// [`Loc`] to `F`. -pub trait DifferentiableRealMapping: - DifferentiableMapping, Codomain = F, DerivativeDomain = Loc> +/// [`Loc`] to `F`. +pub trait DifferentiableRealMapping: + DifferentiableMapping, Codomain = F, DerivativeDomain = Loc> { } -impl DifferentiableRealMapping for T where - T: DifferentiableMapping, Codomain = F, DerivativeDomain = Loc> +impl DifferentiableRealMapping for T where + T: DifferentiableMapping, Codomain = F, DerivativeDomain = Loc> { } @@ -186,7 +186,7 @@ impl Mapping for FlattenedCodomain where X: Space, - G: Mapping>, + G: Mapping>, { type Codomain = F; @@ -198,7 +198,7 @@ /// An auto-trait for constructing a [`FlattenCodomain`] structure for /// flattening the codomain of a [`Mapping`] from [`Loc`]`` to `F`. -pub trait FlattenCodomain: Mapping> + Sized { +pub trait FlattenCodomain: Mapping> + Sized { /// Flatten the codomain from [`Loc`]`` to `F`. fn flatten_codomain(self) -> FlattenedCodomain { FlattenedCodomain { @@ -208,9 +208,9 @@ } } -impl>> FlattenCodomain for G {} +impl>> FlattenCodomain for G {} -/// Container for dimensional slicing [`Loc`]`` codomain of a [`Mapping`] to `F`. +/// Container for dimensional slicing [`Loc`]`` 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> + Clone, + G: Mapping> + 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`]`` to `F`. -pub trait SliceCodomain: - Mapping> + Clone + Sized +pub trait SliceCodomain: + Mapping> + Clone + Sized { /// Flatten the codomain from [`Loc`]`` to `F`. fn slice_codomain(self, slice: usize) -> SlicedCodomain<'static, X, F, Self, N> { @@ -259,8 +259,8 @@ } } -impl> + Clone, const N: usize> - SliceCodomain for G +impl> + Clone, const N: usize> + SliceCodomain for G { } diff -r 495448cca603 -r 6aa955ad8122 src/mapping/dataterm.rs --- 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, P, FloatType = F>, +//+ AdjointProductBoundedBy, P, FloatType = F>, impl Mapping for DataTerm where diff -r 495448cca603 -r 6aa955ad8122 src/nalgebra_support.rs --- 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 Space for Matrix +impl Space for Matrix where - SM: Storage + Clone, - N : Dim, M : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, - DefaultAllocator : Allocator, + SM: Storage + Clone, + N: Dim, + M: Dim, + E: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator: Allocator, { type Decomp = BasicDecomposition; } -impl Mapping> for Matrix -where SM: Storage, SV: Storage + Clone, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator { - type Codomain = OMatrix; +impl Mapping> for Matrix +where + SM: Storage, + SV: Storage + Clone, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, +{ + type Codomain = OMatrix; #[inline] - fn apply>>( - &self, x : I - ) -> Self::Codomain { + fn apply>>(&self, x: I) -> Self::Codomain { x.either(|owned| self.mul(owned), |refr| self.mul(refr)) } } - -impl<'a, SM,SV,N,M,K,E> Linear> for Matrix -where SM: Storage, SV: Storage + Clone, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator { +impl<'a, SM, SV, N, M, K, E> Linear> for Matrix +where + SM: Storage, + SV: Storage + Clone, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, +{ } -impl GEMV, Matrix> for Matrix -where SM: Storage, SV1: Storage + Clone, SV2: StorageMut, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator { - +impl GEMV, Matrix> + for Matrix +where + SM: Storage, + SV1: Storage + Clone, + SV2: StorageMut, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + Float, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, +{ #[inline] - fn gemv>>( - &self, y : &mut Matrix, α : E, x : I, β : E + fn gemv>>( + &self, + y: &mut Matrix, + α: E, + x: I, + β: E, ) { x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β)) } #[inline] - fn apply_mut<'a, I : Instance>>(&self, y : &mut Matrix, x : I) { + fn apply_mut<'a, I: Instance>>(&self, y: &mut Matrix, x: I) { x.eval(|x̃| self.mul_to(x̃, y)) } } -impl AXPY> for Vector -where SM: StorageMut + Clone, SV1: Storage + Clone, - M : Dim, E : Scalar + Zero + One + Float, - DefaultAllocator : Allocator { +impl AXPY> for Vector +where + SM: StorageMut + Clone, + SV1: Storage + Clone, + M: Dim, + E: Scalar + Zero + One + Float, + DefaultAllocator: Allocator, +{ type Owned = OVector; #[inline] - fn axpy>>(&mut self, α : E, x : I, β : E) { + fn axpy>>(&mut self, α: E, x: I, β: E) { x.eval(|x̃| Matrix::axpy(self, α, x̃, β)) } #[inline] - fn copy_from>>(&mut self, y : I) { + fn copy_from>>(&mut self, y: I) { y.eval(|ỹ| Matrix::copy_from(self, ỹ)) } @@ -125,26 +148,40 @@ } }*/ -impl Projection for Vector -where SM: StorageMut + Clone, - M : Dim, E : Scalar + Zero + One + Float + RealField, - DefaultAllocator : Allocator { +impl Projection for Vector +where + SM: StorageMut + Clone, + M: Dim, + E: Scalar + Zero + One + Float + RealField, + DefaultAllocator: Allocator, +{ #[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> -for Matrix -where SM: Storage, SV1: Storage + Clone, SV2: Storage + Clone, - N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + SimdComplexField, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator, - DefaultAllocator : Allocator { - type AdjointCodomain = OMatrix; - type Adjoint<'a> = OMatrix where SM : 'a; +impl<'own, SV1, SV2, SM, N, M, K, E> Adjointable, Matrix> + for Matrix +where + SM: Storage, + SV1: Storage + Clone, + SV2: Storage + Clone, + N: Dim, + M: Dim, + K: Dim, + E: Scalar + Zero + One + SimdComplexField, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, + DefaultAllocator: Allocator, +{ + type AdjointCodomain = OMatrix; + type Adjoint<'a> + = OMatrix + where + SM: 'a; #[inline] fn adjoint(&self) -> Self::Adjoint<'_> { @@ -160,7 +197,7 @@ m2: &Matrix, ) -> T::SimdRealField where - T: SimdComplexField, + T: SimdComplexField, R1: Dim, C1: Dim, S1: Storage, @@ -177,38 +214,38 @@ // TODO: should allow different input storages in `Euclidean`. -impl Euclidean -for Vector -where M : Dim, - S : StorageMut + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { - +impl Euclidean for Vector +where + M: Dim, + S: StorageMut + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ type Output = OVector; #[inline] - fn dot>(&self, other : I) -> E { - Vector::::dot(self, other.ref_instance()) + fn dot>(&self, other: I) -> E { + Vector::::dot(self, other.ref_instance()) } #[inline] fn norm2_squared(&self) -> E { - Vector::::norm_squared(self) + Vector::::norm_squared(self) } #[inline] - fn dist2_squared>(&self, other : I) -> E { + fn dist2_squared>(&self, other: I) -> E { metric_distance_squared(self, other.ref_instance()) } } -impl StaticEuclidean -for Vector -where M : DimName, - S : StorageMut + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { - +impl StaticEuclidean for Vector +where + M: DimName, + S: StorageMut + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ #[inline] fn origin() -> OVector { OVector::zeros() @@ -216,13 +253,13 @@ } /// The default norm for `Vector` is [`L2`]. -impl Normed -for Vector -where M : Dim, - S : Storage + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { - +impl Normed for Vector +where + M: Dim, + S: Storage + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ type NormExp = L2; #[inline] @@ -232,91 +269,95 @@ #[inline] fn is_zero(&self) -> bool { - Vector::::norm_squared(self) == E::ZERO + Vector::::norm_squared(self) == E::ZERO } } -impl HasDual -for Vector -where M : Dim, - S : Storage + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { +impl HasDual for Vector +where + M: Dim, + S: Storage + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ // TODO: Doesn't work with different storage formats. type DualSpace = Self; } -impl Norm -for Vector -where M : Dim, - S : Storage, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { - +impl Norm for Vector +where + M: Dim, + S: Storage, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ #[inline] - fn norm(&self, _ : L1) -> E { + fn norm(&self, _: L1) -> E { nalgebra::Norm::norm(&LpNorm(1), self) } } -impl Dist -for Vector -where M : Dim, - S : Storage + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { +impl Dist for Vector +where + M: Dim, + S: Storage + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ #[inline] - fn dist>(&self, other : I, _ : L1) -> E { + fn dist>(&self, other: I, _: L1) -> E { nalgebra::Norm::metric_distance(&LpNorm(1), self, other.ref_instance()) } } -impl Norm -for Vector -where M : Dim, - S : Storage, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { - +impl Norm for Vector +where + M: Dim, + S: Storage, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ #[inline] - fn norm(&self, _ : L2) -> E { + fn norm(&self, _: L2) -> E { nalgebra::Norm::norm(&LpNorm(2), self) } } -impl Dist -for Vector -where M : Dim, - S : Storage + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { +impl Dist for Vector +where + M: Dim, + S: Storage + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ #[inline] - fn dist>(&self, other : I, _ : L2) -> E { + fn dist>(&self, other: I, _: L2) -> E { nalgebra::Norm::metric_distance(&LpNorm(2), self, other.ref_instance()) } } -impl Norm -for Vector -where M : Dim, - S : Storage, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { - +impl Norm for Vector +where + M: Dim, + S: Storage, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ #[inline] - fn norm(&self, _ : Linfinity) -> E { + fn norm(&self, _: Linfinity) -> E { nalgebra::Norm::norm(&UniformNorm, self) } } -impl Dist -for Vector -where M : Dim, - S : Storage + Clone, - E : Float + Scalar + Zero + One + RealField, - DefaultAllocator : Allocator { +impl Dist for Vector +where + M: Dim, + S: Storage + Clone, + E: Float + Scalar + Zero + One + RealField, + DefaultAllocator: Allocator, +{ #[inline] - fn dist>(&self, other : I, _ : Linfinity) -> E { + fn dist>(&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 + } } - diff -r 495448cca603 -r 6aa955ad8122 src/norms.rs --- 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 { +pub trait Norm { /// Calculate the norm. fn norm(&self, _p: Exponent) -> F; } @@ -110,7 +110,7 @@ } /// Trait for distances with respect to a norm. -pub trait Dist: Norm + Space { +pub trait Dist: Norm + Space { /// Calculate the distance fn dist>(&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: Norm + Sized +pub trait Projection: Norm + Sized where F: Float, { @@ -139,14 +139,14 @@ fn proj_ball_mut(&mut self, ρ: F, q: Exponent); } -/*impl> Norm for E { +/*impl> Norm for E { #[inline] fn norm(&self, _p : L2) -> F { self.norm2() } fn dist(&self, other : &Self, _p : L2) -> F { self.dist2(other) } }*/ -impl + Norm> Projection for E { +impl + Norm> Projection for E { #[inline] fn proj_ball(self, ρ: F, _p: L2) -> Self { self.proj_ball2(ρ) @@ -176,7 +176,7 @@ } } -impl + Normed> Norm> for E { +impl + Normed> Norm, F> for E { fn norm(&self, huber: HuberL1) -> F { huber.apply(self.norm2_squared()) } @@ -188,7 +188,7 @@ } } -// impl> Norm for Vec { +// impl> Norm for Vec { // 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, + Domain: Space + Norm, { type Codomain = F; @@ -214,7 +214,7 @@ } } -pub trait Normed: Space + Norm { +pub trait Normed: Space + Norm { type NormExp: NormExponent; fn norm_exponent(&self) -> Self::NormExp; @@ -283,7 +283,7 @@ impl Norm> for D where F: Float, - D: Norm, + D: Norm<$exponent, F>, C: Constant, { fn norm(&self, e: Weighted<$exponent, C>) -> F { diff -r 495448cca603 -r 6aa955ad8122 src/sets.rs --- 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 where U : Space { +pub trait Set +where + U: Space, +{ /// Check for element containment - fn contains>(&self, item : I) -> bool; + fn contains>(&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; + /// Returns the greatest set of same class contaied by n both parameter sets. + fn intersect(&self, other: &Self) -> Option; } -impl Set> -for Cube -where U : Num + PartialOrd + Sized { - fn contains>>(&self, item : I) -> bool { - self.0.iter().zip(item.ref_instance().iter()).all(|(s, x)| s.contains(x)) +impl Set> for Cube +where + U: Num + PartialOrd + Sized, +{ + fn contains>>(&self, item: I) -> bool { + self.0 + .iter() + .zip(item.ref_instance().iter()) + .all(|(s, x)| s.contains(x)) } } -impl Set for RangeFull { - fn contains>(&self, _item : I) -> bool { true } +impl Set for RangeFull { + fn contains>(&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 where A : Euclidean, F : Float { - pub orthogonal : A, - pub offset : F, +#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)] +pub struct Halfspace +where + A: Euclidean, + F: Float, +{ + pub orthogonal: A, + pub offset: F, } -impl Halfspace where A : Euclidean, F : Float { +impl Halfspace +where + A: Euclidean, + 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 where F : Float { +pub trait SpannedHalfspace +where + F: Float, +{ /// Type of the orthogonal vector describing the halfspace. - type A : Euclidean; + type A: Euclidean; /// Returns the halfspace spanned by this set. fn spanned_halfspace(&self) -> Halfspace; } // TODO: Gram-Schmidt for higher N. -impl SpannedHalfspace for [Loc; 2] { - type A = Loc; +impl SpannedHalfspace for [Loc<1, F>; 2] { + type A = Loc<1, F>; fn spanned_halfspace(&self) -> Halfspace { 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 SpannedHalfspace for [Loc; 2] { - type A = Loc; +impl SpannedHalfspace for [Loc<2, F>; 2] { + type A = Loc<2, F>; fn spanned_halfspace(&self) -> Halfspace { let (x0, x1) = (&self[0], &self[1]); let d = x1 - x0; @@ -104,28 +127,30 @@ } } -impl Set for Halfspace +impl Set for Halfspace where - A : Euclidean, - F : Float, + A: Euclidean, + F: Float, { #[inline] - fn contains>(&self, item : I) -> bool { + fn contains>(&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(pub [Halfspace; N]) -where A : Euclidean, F : Float; +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub struct NPolygon(pub [Halfspace; N]) +where + A: Euclidean, + F: Float; -impl Set for NPolygon +impl Set for NPolygon where - A : Euclidean, - F : Float, + A: Euclidean, + F: Float, { - fn contains>(&self, item : I) -> bool { + fn contains>(&self, item: I) -> bool { let r = item.ref_instance(); self.0.iter().all(|halfspace| halfspace.contains(r)) } diff -r 495448cca603 -r 6aa955ad8122 src/sets/cube.rs --- 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(pub(super) [[U; 2]; N]); +pub struct Cube(pub(super) [[U; 2]; N]); // Need to manually implement as [F; N] serialisation is provided only for some N. -impl Serialize for Cube +impl Serialize for Cube where F: Serialize, { @@ -50,7 +44,7 @@ } } -impl FixedLength for Cube { +impl FixedLength for Cube { type Iter = std::array::IntoIter<[A; 2], N>; type Elem = [A; 2]; #[inline] @@ -59,7 +53,7 @@ } } -impl FixedLengthMut for Cube { +impl FixedLengthMut for Cube { 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 for &'a Cube { +impl<'a, A: Num, const N: usize> FixedLength for &'a Cube { 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, +pub struct CubeCornersIter<'a, U: Num, const N: usize> { + index: usize, + cube: &'a Cube, } -impl<'a, U : Num, const N : usize> Iterator for CubeCornersIter<'a, U, N> { - type Item = Loc; +impl<'a, U: Num, const N: usize> Iterator for CubeCornersIter<'a, U, N> { + type Item = Loc; #[inline] fn next(&mut self) -> Option { 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 Cube { +impl Cube { /// 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(&self, f : impl Fn(usize, U, U) -> T) -> [T; N] { + pub fn map_indexed(&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(&self, f : impl Fn(U, U) -> T) -> [T; N] { + pub fn map(&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 { + pub fn span_start(&self) -> Loc { 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 { + pub fn span_end(&self) -> Loc { 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 { - Loc::new(self.map(|a, b| b-a)) + pub fn width(&self) -> Loc { + 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) -> Self { + pub fn shift(&self, shift: &Loc) -> 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 Cube { +impl Cube { /// Returns the centre of the cube - pub fn center(&self) -> Loc { + pub fn center(&self) -> Loc { map1(self, |&[a, b]| (a + b) / F::TWO).into() } } -impl Cube { +impl 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; 2] { + pub fn corners(&self) -> [Loc<1, U>; 2] { let [[a, b]] = self.0; [a.into(), b.into()] } } -impl Cube { +impl 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; 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 Cube { +impl 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; 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 From<[[U; 2]; N]> for Cube { +impl From<[[U; 2]; N]> for Cube { #[inline] - fn from(data : [[U; 2]; N]) -> Self { + fn from(data: [[U; 2]; N]) -> Self { Cube(data) } } -impl From> for [[U; 2]; N] { +impl From> for [[U; 2]; N] { #[inline] - fn from(Cube(data) : Cube) -> Self { + fn from(Cube(data): Cube) -> Self { data } } - -impl Cube where U : Num + PartialOrd { +impl Cube +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) -> bool { - self.iter_coords().zip(other.iter_coords()).all(|([a1, b1], [a2, b2])| { - a1 < b2 && a2 < b1 - }) + pub fn intersects(&self, other: &Cube) -> 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) -> bool { - self.iter_coords().zip(other.iter_coords()).all(|([a1, b1], [a2, b2])| { - a1 <= a2 && b1 >= b2 - }) + pub fn contains_set(&self, other: &Cube) -> 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 { + pub fn minnorm_point(&self) -> Loc { 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 { + pub fn maxnorm_point(&self) -> Loc { 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 SetOrd for Cube<$t, N> { + impl SetOrd for Cube { #[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 std::ops::Index for Cube { +impl std::ops::Index for Cube { type Output = [U; 2]; #[inline] fn index(&self, index: usize) -> &Self::Output { @@ -346,7 +358,7 @@ } } -impl std::ops::IndexMut for Cube { +impl std::ops::IndexMut for Cube { #[inline] fn index_mut(&mut self, index: usize) -> &mut Self::Output { &mut self.0[index]