Mercurial > repos > alg_tools / changeset

--- a/src/bisection_tree/aggregator.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/bisection_tree/aggregator.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -4,6 +4,7 @@
 
 use crate::types::*;
 use crate::sets::Set;
+use crate::instance::Instance;
 
 /// Trait for aggregating information about a branch of a [bisection tree][super::BT].
 ///
@@ -136,10 +137,11 @@
 }
 
 impl<F : Float> Set<F> for Bounds<F> {
-    fn contains(&self, item : &F) -> bool {
+    fn contains<I : Instance<F>>(&self, item : I) -> bool {
+        let v = item.own();
         let &Bounds(l, u) = self;
         debug_assert!(l <= u);
-        l <= *item && *item <= u
+        l <= v && v <= u
     }
 }
 

--- a/src/bisection_tree/btfn.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/bisection_tree/btfn.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -505,7 +505,7 @@
 ///
 /// `U` is the domain, generally [`Loc`]`<F, N>`, and `F` the type of floating point numbers.
 /// `Self` is generally a set of `U`, for example, [`Cube`]`<F, N>`.
-pub trait P2Minimise<U, F : Float> : Set<U> {
+pub trait P2Minimise<U : Space, F : Float> : Set<U> {
     /// Minimise `g` over the set presented by `Self`.
     ///
     /// The function returns `(x, v)` where `x` is the minimiser `v` an approximation of `g(x)`.

--- a/src/direct_product.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/direct_product.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -10,7 +10,7 @@
 use serde::{Serialize, Deserialize};
 use crate::types::{Num, Float};
 use crate::{maybe_lifetime, maybe_ref};
-use crate::euclidean::{Dot, Euclidean};
+use crate::euclidean::Euclidean;
 use crate::instance::{Instance, InstanceMut, Decomposition, DecompositionMut, MyCow};
 use crate::mapping::Space;
 use crate::linops::AXPY;
@@ -237,18 +237,6 @@
 impl_pair_vectorspace_ops!((f32, f32), f32);
 impl_pair_vectorspace_ops!((f64, f64), f64);
 
-impl<A, B, U, V, F> Dot<Pair<U, V>, F> for Pair<A, B>
-where
-    A : Dot<U, F>,
-    B : Dot<V, F>,
-    F : Num
-{
-
-    fn dot(&self, Pair(ref u, ref v) : &Pair<U, V>) -> F {
-        self.0.dot(u) + self.1.dot(v)
-    }
-}
-
 type PairOutput<F, A, B> = Pair<<A as Euclidean<F>>::Output, <B as Euclidean<F>>::Output>;
 
 impl<A, B, F> Euclidean<F> for Pair<A, B>
@@ -257,7 +245,7 @@
     B : Euclidean<F>,
     F : Float,
     PairOutput<F, A, B> : Euclidean<F>,
-    Self : Sized + Dot<Self,F>
+    Self : Sized
           + Mul<F, Output=PairOutput<F, A, B>> + MulAssign<F>
           + Div<F, Output=PairOutput<F, A, B>> + DivAssign<F>
           + Add<Self, Output=PairOutput<F, A, B>>
@@ -270,6 +258,15 @@
 {
     type Output = PairOutput<F, A, B>;
 
+    fn dot<I : Instance<Self>>(&self, other : I) -> F {
+        let Pair(u, v) = other.decompose();
+        self.0.dot(u) + self.1.dot(v)
+    }
+
+    fn norm2_squared(&self) -> F {
+        self.0.norm2_squared() + self.1.norm2_squared()
+    }
+
     fn dist2_squared(&self, Pair(ref u, ref v) : &Self) -> F {
         self.0.dist2_squared(u) + self.1.dist2_squared(v)
     }

--- a/src/euclidean.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/euclidean.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -4,38 +4,35 @@
 
 use std::ops::{Mul, MulAssign, Div, DivAssign, Add, Sub, AddAssign, SubAssign, Neg};
 use crate::types::*;
-use crate::mapping::Space;
-
-/// Space (type) with a defined dot product.
-///
-/// `U` is the space of the multiplier, and `F` the space of scalars.
-/// Since `U` ≠ `Self`, this trait can also implement dual products.
-pub trait Dot<U, F> {
-    fn dot(&self, other : &U) -> F;
-}
+use crate::instance::Instance;
+use crate::norms::{HasDual, Reflexive};
 
 /// Space (type) with Euclidean and vector space structure
 ///
 /// The type should implement vector space operations (addition, subtraction, scalar
 /// multiplication and scalar division) along with their assignment versions, as well
-/// as the [`Dot`] product with respect to `Self`.
-pub trait Euclidean<F : Float> : Space + Dot<Self,F>
-        + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F>
-        + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F>
-        + Add<Self, Output=<Self as Euclidean<F>>::Output>
-        + Sub<Self, Output=<Self as Euclidean<F>>::Output>
-        + for<'b> Add<&'b Self, Output=<Self as Euclidean<F>>::Output>
-        + for<'b> Sub<&'b Self, Output=<Self as Euclidean<F>>::Output>
-        + AddAssign<Self> + for<'b> AddAssign<&'b Self>
-        + SubAssign<Self> + for<'b> SubAssign<&'b Self>
-        + Neg<Output=<Self as Euclidean<F>>::Output> {
+/// as an inner product.
+pub trait Euclidean<F : Float> : HasDual<F, DualSpace=Self> + Reflexive<F>
+    + Mul<F, Output=<Self as Euclidean<F>>::Output> + MulAssign<F>
+    + Div<F, Output=<Self as Euclidean<F>>::Output> + DivAssign<F>
+    + Add<Self, Output=<Self as Euclidean<F>>::Output>
+    + Sub<Self, Output=<Self as Euclidean<F>>::Output>
+    + for<'b> Add<&'b Self, Output=<Self as Euclidean<F>>::Output>
+    + for<'b> Sub<&'b Self, Output=<Self as Euclidean<F>>::Output>
+    + AddAssign<Self> + for<'b> AddAssign<&'b Self>
+    + SubAssign<Self> + for<'b> SubAssign<&'b Self>
+    + Neg<Output=<Self as Euclidean<F>>::Output>
+{
     type Output : Euclidean<F>;
 
+    // Inner product
+    fn dot<I : Instance<Self>>(&self, other : I) -> F;
+
     /// Calculate the square of the 2-norm, $\frac{1}{2}\\|x\\|_2^2$, where `self` is $x$.
-    #[inline]
-    fn norm2_squared(&self) -> F {
-        self.dot(self)
-    }
+    ///
+    /// This is not automatically implemented to avoid imposing
+    /// `for <'a> &'a Self : Instance<Self>` trait bound bloat.
+    fn norm2_squared(&self) -> F;
 
     /// Calculate the square of the 2-norm divided by 2, $\frac{1}{2}\\|x\\|_2^2$,
     /// where `self` is $x$.

--- a/src/fe_model/p2_local_model.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/fe_model/p2_local_model.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -6,7 +6,8 @@
 use crate::loc::Loc;
 use crate::sets::{Set,NPolygon,SpannedHalfspace};
 use crate::linsolve::*;
-use crate::euclidean::Dot;
+use crate::euclidean::Euclidean;
+use crate::instance::Instance;
 use super::base::{LocalModel,RealLocalModel};
 use crate::sets::Cube;
 use numeric_literals::replace_float_literals;
@@ -30,7 +31,8 @@
 
 impl<'a, F : Float> Set<Loc<F,1>> for RealInterval<F> {
     #[inline]
-    fn contains(&self, &Loc([x]) : &Loc<F,1>) -> bool {
+    fn contains<I : Instance<Loc<F, 1>>>(&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)
     }
@@ -38,11 +40,11 @@
 
 impl<'a, F : Float> Set<Loc<F,2>> for PlanarSimplex<F> {
     #[inline]
-    fn contains(&self, x : &Loc<F,2>) -> bool {
+    fn contains<I : Instance<Loc<F, 2>>>(&self, z : I) -> bool {
         let &[x0, x1, x2] = &self.0;
         NPolygon([[x0, x1].spanned_halfspace(),
                   [x1, x2].spanned_halfspace(),
-                  [x2, x0].spanned_halfspace()]).contains(x)
+                  [x2, x0].spanned_halfspace()]).contains(z)
     }
 }
 

--- a/src/loc.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/loc.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -429,19 +429,17 @@
 domination!(Linfinity, L2);
 domination!(L2, L1);
 
-impl<F : Num,const N : usize> Dot<Loc<F, N>,F> for Loc<F, N> {
+impl<F : Float,const N : usize> Euclidean<F> for Loc<F, N> {
+    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 : &Loc<F, N>) -> F {
+    fn dot<I : Instance<Self>>(&self, other : I) -> F {
         self.0.iter()
-              .zip(other.0.iter())
+              .zip(other.ref_instance().0.iter())
               .fold(F::ZERO, |m, (&v, &w)| m + v * w)
     }
-}
-
-impl<F : Float,const N : usize> Euclidean<F> for Loc<F, N> {
-    type Output = Self;
 
     /// This implementation is not stabilised as it's meant to be used for very small vectors.
     /// Use [`nalgebra`] for larger vectors.
@@ -696,7 +694,7 @@
     type Decomp = BasicDecomposition;
 }
 
-impl<F : Num, const N : usize> Mapping<Loc<F, N>> for Loc<F, N> {
+impl<F : Float, const N : usize> Mapping<Loc<F, N>> for Loc<F, N> {
     type Codomain = F;
 
     fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
@@ -704,9 +702,9 @@
     }
 }
 
-impl<F : Num, const N : usize> Linear<Loc<F, N>> for Loc<F, N> { }
+impl<F : Float, const N : usize> Linear<Loc<F, N>> for Loc<F, N> { }
 
-impl<F : Num, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> {
+impl<F : Float, const N : usize> AXPY<F, Loc<F, N>> for Loc<F, N> {
     type Owned = Self;
 
     #[inline]

--- a/src/nalgebra_support.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/nalgebra_support.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -1,7 +1,7 @@
 /*!
 Integration with nalgebra.
 
-This module mainly implements [`Euclidean`], [`Norm`], [`Dot`], [`Linear`], etc. for [`nalgebra`]
+This module mainly implements [`Euclidean`], [`Norm`], [`Linear`], etc. for [`nalgebra`]
 matrices and vectors.
 It also provides [`ToNalgebraRealField`] as a vomit-inducingly ugly workaround to nalgebra
 force-feeding its own versions of the same basic mathematical methods on `f32` and `f64` as
@@ -153,20 +153,6 @@
     }
 }
 
-impl<E,M,S,Si> Dot<Vector<E,M,Si>,E>
-for Vector<E,M,S>
-where M : Dim,
-      E : Float + Scalar + Zero + One,
-      S : Storage<E,M>,
-      Si : Storage<E,M>,
-      DefaultAllocator : Allocator<M> {
-
-    #[inline]
-    fn dot(&self, other : &Vector<E,M,Si>) -> E {
-        Vector::<E,M,S>::dot(self, other)
-    }
-}
-
 /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`.
 #[inline]
 fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>(
@@ -200,7 +186,12 @@
       DefaultAllocator : Allocator<M> {
 
     type Output = OVector<E, M>;
-    
+
+    #[inline]
+    fn dot<I : Instance<Self>>(&self, other : I) -> E {
+        Vector::<E,M,S>::dot(self, other.ref_instance())
+    }
+
     #[inline]
     fn norm2_squared(&self) -> E {
         Vector::<E,M,S>::norm_squared(self)

--- a/src/norms.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/norms.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -221,11 +221,16 @@
     type DualSpace : Normed<F>;
 }
 
+/// Automatically implemented trait for reflexive spaces
 pub trait Reflexive<F : Num = f64> : HasDual<F>
 where
     Self::DualSpace : HasDual<F, DualSpace = Self>
 { }
 
+impl<F : Num, X : HasDual<F>> Reflexive<F> for X
+where
+    X::DualSpace : HasDual<F, DualSpace = X>
+{ }
 
 pub trait HasDualExponent : NormExponent {
     type DualExp : NormExponent;

--- a/src/sets.rs	Sat Dec 21 23:32:20 2024 -0500
+++ b/src/sets.rs	Sun Dec 22 14:54:46 2024 -0500
@@ -3,19 +3,19 @@
 */
 
 use std::ops::{RangeFull,RangeFrom,Range,RangeInclusive,RangeTo,RangeToInclusive};
-use std::marker::PhantomData;
 use crate::types::*;
 use crate::loc::Loc;
-use crate::euclidean::Dot;
+use crate::euclidean::Euclidean;
+use crate::instance::{Space, Instance};
 use serde::Serialize;
 
 pub mod cube;
 pub use cube::Cube;
 
 /// Trait for arbitrary sets. The parameter `U` is the element type.
-pub trait Set<U> where U : ?Sized {
+pub trait Set<U> where U : Space {
     /// Check for element containment
-    fn contains(&self, item : &U) -> bool;
+    fn contains<I : Instance<U>>(&self, item : I) -> bool;
 }
 
 /// Additional ordering (besides [`PartialOrd`]) of a subfamily of sets:
@@ -31,22 +31,27 @@
 impl<U, const N : usize> Set<Loc<U, N>>
 for Cube<U,N>
 where U : Num + PartialOrd + Sized {
-    fn contains(&self, item : &Loc<U, N>) -> bool {
-        self.0.iter().zip(item.iter()).all(|(s, x)| s.contains(x))
+    fn contains<I : Instance<Loc<U, N>>>(&self, item : I) -> bool {
+        self.0.iter().zip(item.ref_instance().iter()).all(|(s, x)| s.contains(x))
     }
 }
 
-impl<U> Set<U> for RangeFull {
-    fn contains(&self, _item : &U) -> bool { true }
+impl<U : Space> Set<U> for RangeFull {
+    fn contains<I : Instance<U>>(&self, _item : I) -> bool { true }
 }
 
 macro_rules! impl_ranges {
     ($($range:ident),*) => { $(
-        impl<U,Idx> Set<U>
-        for $range<Idx>
-        where Idx : PartialOrd<U>, U : PartialOrd<Idx> + ?Sized, Idx : PartialOrd {
+        impl<U,Idx> Set<U> for $range<Idx>
+        where
+            Idx : PartialOrd<U>,
+            U : PartialOrd<Idx> + Space,
+            Idx : PartialOrd
+        {
             #[inline]
-            fn contains(&self, item : &U) -> bool { $range::contains(&self, item) }
+            fn contains<I : Instance<U>>(&self, item : I) -> bool {
+                item.eval(|x| $range::contains(&self, x))
+            }
         }
     )* }
 }
@@ -61,40 +66,39 @@
 /// `U` is the element type, `F` the floating point number type, and `A` the type of the
 /// orthogonal (dual) vectors. They need implement [`Dot<U, F>`].
 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
-pub struct Halfspace<A, F, U> where A : Dot<U, F>, F : Float {
+pub struct Halfspace<A, F> where A : Euclidean<F>, F : Float {
     pub orthogonal : A,
     pub offset : F,
-    _phantom : PhantomData<U>,
 }
 
-impl<A,F,U> Halfspace<A,F,U> where A : Dot<U,F>, F : Float {
+impl<A,F> Halfspace<A,F> where A : Euclidean<F>, F : Float {
     #[inline]
     pub fn new(orthogonal : A, offset : F) -> Self {
-        Halfspace{ orthogonal : orthogonal, offset : offset, _phantom : PhantomData }
+        Halfspace{ orthogonal : orthogonal, offset : offset }
     }
 }
 
 /// Trait for generating a halfspace spanned by another set `Self` of elements of type `U`.
-pub trait SpannedHalfspace<F, U> where  F : Float {
+pub trait SpannedHalfspace<F> where  F : Float {
     /// Type of the orthogonal vector describing the halfspace.
-    type A : Dot<U, F>;
+    type A : Euclidean<F>;
     /// Returns the halfspace spanned by this set.
-    fn spanned_halfspace(&self) -> Halfspace<Self::A, F, U>;
+    fn spanned_halfspace(&self) -> Halfspace<Self::A, F>;
 }
 
 // TODO: Gram-Schmidt for higher N.
-impl<F : Float> SpannedHalfspace<F,Loc<F, 1>> for [Loc<F, 1>; 2] {
+impl<F : Float> SpannedHalfspace<F> for [Loc<F, 1>; 2] {
     type A = Loc<F, 1>;
-    fn spanned_halfspace(&self) -> Halfspace<Self::A, F, Loc<F, 1>> {
+    fn spanned_halfspace(&self) -> Halfspace<Self::A, F> {
         let (x0, x1) = (self[0], self[1]);
         Halfspace::new(x1-x0, x0[0])
     }
 }
 
 // TODO: Gram-Schmidt for higher N.
-impl<F : Float> SpannedHalfspace<F,Loc<F, 2>> for [Loc<F, 2>; 2] {
+impl<F : Float> SpannedHalfspace<F> for [Loc<F, 2>; 2] {
     type A = Loc<F, 2>;
-    fn spanned_halfspace(&self) -> Halfspace<Self::A, F, Loc<F, 2>> {
+    fn spanned_halfspace(&self) -> Halfspace<Self::A, F> {
         let (x0, x1) = (&self[0], &self[1]);
         let d = x1 - x0;
         let orthog = loc![d[1], -d[0]]; // We don't normalise for efficiency
@@ -103,19 +107,29 @@
     }
 }
 
-impl<A,U,F> Set<U> for Halfspace<A,F,U> where A : Dot<U,F>, F : Float {
+impl<A,F> Set<A> for Halfspace<A,F>
+where
+    A : Euclidean<F>,
+    F : Float,
+{
     #[inline]
-    fn contains(&self, item : &U) -> bool {
+    fn contains<I : Instance<A>>(&self, item : I) -> bool {
         self.orthogonal.dot(item) >= self.offset
     }
 }
 
 /// Polygons defined by `N` `Halfspace`s.
 #[derive(Clone,Copy,Debug,Eq,PartialEq)]
-pub struct NPolygon<A, F, U, const N : usize>(pub [Halfspace<A,F,U>; N]) where A : Dot<U,F>, F : Float;
+pub struct NPolygon<A, F, const N : usize>(pub [Halfspace<A,F>; N])
+where A : Euclidean<F>, F : Float;
 
-impl<A,U,F,const N : usize> Set<U> for NPolygon<A,F,U,N> where A : Dot<U,F>, F : Float {
-    fn contains(&self, item : &U) -> bool {
-        self.0.iter().all(|halfspace| halfspace.contains(item))
+impl<A,F,const N : usize> Set<A> for NPolygon<A,F,N>
+where
+    A : Euclidean<F>,
+    F : Float,
+{
+    fn contains<I : Instance<A>>(&self, item : I) -> bool {
+        let r = item.ref_instance();
+        self.0.iter().all(|halfspace| halfspace.contains(r))
     }
 }