initial commit

Includes an initial implementation, some tests, and some documentation.

Signed-off-by: Stephen Waits <steve@waits.net>

src/height_vector.rs

@@ -0,0 +1,163 @@
+//! Implementation of a Vivaldi height vector or coordinate.
+//!
+//! A Vivaldi height vector consists of an N-dimensional Euclidean/Cartesian space vector and a
+//! greater-than-zero height. This vector is intented to represent latency; hence the distance
+//! between two Vivaldi coordinates should be close to the actual round-trip-time (i.e. ping time).
+//!
+//! The height component is included to help solve for the triangle inequality; the fact that three
+//! points in the real world won't necessarily follow the triangle equality principle. In other
+//! words, ping times aren't perfectly correlated to geographic distance. Thus, the height can be
+//! considered a representation of a sort of "stem time" which is the time from a node to the core
+//! internet; and the vectors can be considered locations on the core internet backbone where ping
+//! times are more consistently correlated with geographic distance.
+//!
+//! Quoting the paper:
+//!
+//! > A height vector consists of a Euclidean coordinate augmented
+//! > with a height. The Euclidean portion models a high-speed Internet
+//! > core with latencies proportional to geographic distance, while the
+//! > height models the time it takes packets to travel the access link
+//! > from the node to the core. The cause of the access link latency may
+//! > be queuing delay (as in the case of an oversubscribed cable line),
+//! > low bandwidth (as in the case of DSL, cable modems, or telephone
+//! > modems), or even the sheer length of the link (as in the case of
+//! > long-distance fiber-optic cables).
+//! > A packet sent from one node to another must travel the source
+//! > node’s height, then travel in the Euclidean space, then travel the
+//! > destination node’s height. Even if the two nodes have the same
+//! > height, the distance between them is their Euclidean distance plus
+//! > the two heights. This is the fundamental difference between height
+//! > vectors and adding a dimension to the Euclidean space. Intuitively,
+//! > packet transmission can only be done in the core, not above it.
+//! > The height vector model is implemented by redefining the usual
+//! > vector operations (note the + on the right hand side of the subtrac-
+//! > tion equation):
+//!
+//! ```text
+//! [x, xh] − [y, y] = [(x − y), xh + yh]
+//! ∥∥∥[x, xh]∥∥∥ = ∥∥∥x∥∥∥ + xh
+//! α × [x, xh] = [αx, αxh ]
+//! ```
+//!
+//! > Each node has a positive height element in its coordinates, so that
+//! > its height can always be scaled up or down.
+
+use std::ops::{Add, Mul, Sub};
+
+use nanorand::{Rng, WyRand};
+
+use crate::vector::Vector;
+
+//
+// **** Constants ****
+//
+
+//
+// **** Structs ****
+//
+
+/// The `HeightVector` is a tuple containing an `N`-dimensional Cartesian vector and a `> 0` height.
+///
+/// ## Generic Parameters
+///
+/// - `N`: the dimensionality of the vector portion (i.e. non-height) of the Vivaldi height vector
+#[derive(Copy, Clone, Debug)]
+pub(crate) struct HeightVector<const N: usize>(Vector<f32, N>, f32);
+
+//
+// **** Implementations ****
+//
+
+impl<const N: usize> HeightVector<N> {
+    /// A new height vector is a random unit vector plus a default height.
+    pub(crate) fn new() -> Self {
+        let mut rng = WyRand::new();
+        let mut vec = [0.0; N];
+        for i in vec.iter_mut().take(N) {
+            *i = rng.generate::<f32>();
+        }
+        let height = rng.generate::<f32>();
+        Self(Vector::<f32, N>::from(vec), height).normalized()
+    }
+
+    /// The magnitude of a Vivaldi height vector is defined as the magnitude of the vector plus the
+    /// height value.
+    pub(crate) fn len(&self) -> f32 {
+        self.0.len() + self.1
+    }
+
+    /// A normalized Vivaldi height vector is just like a normalized vector: the vector scaled
+    /// by the inveerse of its length.
+    pub(crate) fn normalized(&self) -> Self {
+        let len = self.len();
+        if len < f32::EPSILON {
+            // if we have a bad vector, generate a new random vector
+            Self::new()
+        } else {
+            // otherwise scale it by the inverse of its magnitude like any normal vector
+            Self(self.0 / len, self.1 / len)
+        }
+    }
+}
+
+//
+// **** Trait Implementations ****
+//
+
+impl<const N: usize> Default for HeightVector<N> {
+    /// Default value for a Vivaldi height vector is just the defaults of its children types.
+    fn default() -> Self {
+        Self(Default::default(), Default::default())
+    }
+}
+
+impl<const N: usize> From<([f32; N], f32)> for HeightVector<N> {
+    /// Convert from a `([f32; N], f32)` (vector, height) to a Vivaldi height vector type.
+    fn from(value: ([f32; N], f32)) -> Self {
+        Self(Vector::<f32, N>::from(value.0), value.1)
+    }
+}
+
+impl<const N: usize> Add for HeightVector<N> {
+    type Output = Self;
+
+    /// Add two Vivaldi height vectors.
+    fn add(self, rhs: Self) -> Self::Output {
+        Self(self.0 + rhs.0, self.1 + rhs.1)
+    }
+}
+
+impl<const N: usize> Sub for HeightVector<N> {
+    type Output = Self;
+
+    /// Subtract two Vivaldi height vectors. Note that this is the difference in the vectors
+    /// and the summation of the heights, as defined by Vivaldi's author.
+    fn sub(self, rhs: Self) -> Self::Output {
+        Self(self.0 - rhs.0, self.1 + rhs.1)
+    }
+}
+
+impl<const N: usize> Mul<f32> for HeightVector<N> {
+    type Output = Self;
+
+    /// Multiply a Vivaldi height vector by a scalar. Works the same as normal vector scaling.
+    fn mul(self, rhs: f32) -> Self::Output {
+        Self(self.0 * rhs, self.1 * rhs)
+    }
+}
+
+//
+// **** Tests ****
+//
+#[cfg(test)]
+mod tests {
+    use assert_approx_eq::assert_approx_eq;
+
+    use super::*;
+
+    #[test]
+    fn test_len() {
+        let a = HeightVector::<3>::from(([1.0, 2.0, 3.0], 4.0));
+        assert_approx_eq!(a.len(), 7.741_657, 0.00001);
+    }
+}

src/lib.rs

@@ -0,0 +1,4 @@
+mod height_vector;
+mod vector;
+
+pub mod network_coordinate;

src/network_coordinate.rs

@@ -0,0 +1,279 @@
+//! A Vivaldi [`NetworkCoordinate`]
+//!
+//! This is an implementation of Vivaldi NCs per the original paper. It implements the following
+//! alogirthm (quoting from paper):
+//!
+//! ```text
+//! // Incorporate new information: node j has been
+//! // measured to be rtt ms away, has coordinates x j,
+//! // and an error estimate of e j .
+//! //
+//! // Our own coordinates and error estimate are xi and ei.
+//! //
+//! // The constants ce and cc are tuning parameters.
+//!
+//! vivaldi(rtt, xj, ej)
+//!     // Sample weight balances local and remote error. (1)
+//!     w = ei /(ei + ej)
+//!     // Compute relative error of this sample. (2)
+//!     es = ∣∣∣‖xi − xj‖ − rtt∣∣∣/rtt
+//!     // Update weighted moving average of local error. (3)
+//!     ei = es × ce × w + ei × (1 − ce × w)
+//!     // Update local coordinates. (4)
+//!     δ = cc × w
+//!     xi = xi + δ × (rtt − ‖xi − xj ‖) × u(xi − xj)
+//! ```
+//!
+
+use std::time::Duration;
+
+use crate::height_vector::HeightVector;
+
+//
+// **** Constants ****
+//
+
+// Vivaldi tuning parameters
+const C_ERROR: f32 = 0.25;
+const C_DELTA: f32 = 0.25;
+
+// initial error value
+const DEFAULT_ERROR: f32 = 1000.0;
+
+// error should always be greater than zero
+const MIN_ERROR: f32 = f32::EPSILON;
+
+//
+// **** Structs ****
+//
+
+#[derive(Debug)]
+pub struct NetworkCoordinate<const N: usize> {
+    position: HeightVector<N>,
+    error: f32,
+}
+
+//
+// **** Implementations ****
+//
+
+impl<const N: usize> NetworkCoordinate<N> {
+    /// Creates a new random [`NetworkCoordinate`]
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use vivaldi_nc::network_coordinate::NetworkCoordinate;
+    ///
+    /// // create a new 3-dimensional random NC
+    /// let a: NetworkCoordinate<3> = NetworkCoordinate::new();
+    ///
+    /// // print the NC
+    /// println!("Our new NC is: {:#?}", a);
+    /// ```
+    pub fn new() -> Self {
+        Self::default()
+    }
+
+    /// Given another Vivaldi [`NetworkCoordinate`], estimate the round trip time (ie ping) between them.
+    ///
+    /// This is done by computing the height vector distance between between the two coordinates.
+    /// Vivaldi uses this distance as a representation of estimated round trip time.
+    ///
+    /// # Parameters
+    ///
+    /// - `rhs`: the other coordinate
+    ///
+    /// # Returns
+    ///
+    /// - the estimated round trip time as a `Duration`
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use vivaldi_nc::network_coordinate::NetworkCoordinate;
+    ///
+    /// // create some 2-dimensional NCs for the sake of this example. These will just be random
+    /// // NCs. In a real usecase these would have meaningful values.
+    /// let a: NetworkCoordinate<2> = NetworkCoordinate::new();
+    /// let b: NetworkCoordinate<2> = NetworkCoordinate::new();
+    ///
+    /// // get the estimated RTT, convert to milliseconds, and print
+    /// println!("Estimated RTT: {}", a.estimated_rtt(&b).as_millis());
+    /// ```
+    ///
+    pub fn estimated_rtt(&self, rhs: &Self) -> Duration {
+        // estimated rss is euclidean distance between the two plus the sum of the heights
+        Duration::from_secs_f32((self.position - rhs.position).len() / 1000.0)
+    }
+
+    /// Given another Vivaldi [`NetworkCoordinate`], adjust our coordinateto better represent the actual round
+    /// trip time (aka distance) between us.
+    ///
+    /// # Parameters
+    ///
+    /// - `rhs`: the other coordinate
+    /// - `rtt`: the measured round trip time between `self` and `rhs`
+    ///
+    /// # Returns
+    ///
+    /// - a reference to `self`
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use std::time::Duration;
+    /// use vivaldi_nc::network_coordinate::NetworkCoordinate;
+    ///
+    /// // We always have our own NC:
+    /// let mut local: NetworkCoordinate<2> = NetworkCoordinate::new();
+    ///
+    /// // Assume we received a NC from a remote node:
+    /// let remote: NetworkCoordinate<2> = NetworkCoordinate::new();
+    ///
+    /// // And we measured the RTT between us and the remote node:
+    /// let rtt = Duration::from_millis(100);
+    ///
+    /// // Now we can update our NC to adjust our position relative to the remote node:
+    /// local.update(&remote, rtt);
+    /// ```
+    pub fn update(&mut self, rhs: &Self, rtt: Duration) -> &Self {
+        // rtt needs to be positive
+        if rtt.as_secs_f32() < 0.0 {
+            return self;
+        }
+
+        // convert Durations into f32s as fractional milliseconds for convenience
+        let rtt_ms = rtt.as_secs_f32() * 1000.0;
+        let rtt_estimated_ms = self.estimated_rtt(rhs).as_secs_f32() * 1000.0;
+
+        // Sample weight balances local and remote error. (1)
+        // w = ei /(ei + ej )
+        let w = self.error / (self.error + rhs.error);
+
+        // Compute relative error of this sample. (2)
+        // es = ∣∣∣‖xi − xj‖ − rtt∣∣∣/rtt
+        let error = rtt_ms - rtt_estimated_ms;
+        let es = error.abs() / rtt_ms;
+
+        // Update weighted moving average of local error. (3)
+        // ei = es × ce × w + ei × (1 − ce × w)
+        self.error = (es * C_ERROR * w + self.error * (1.0 - C_ERROR * w)).max(MIN_ERROR);
+
+        // Update local coordinates. (4)
+        // δ = cc × w
+        let delta = C_DELTA * w;
+        // xi = xi + δ × (rtt − ‖xi − xj ‖) × u(xi − xj)
+        self.position = self.position + (self.position - rhs.position).normalized() * delta * error;
+
+        // return reference to updated self
+        self
+    }
+}
+
+//
+// **** Trait Implementations ****
+//
+
+impl<const N: usize> Default for NetworkCoordinate<N> {
+    fn default() -> Self {
+        Self {
+            position: HeightVector::<N>::new(),
+            error: DEFAULT_ERROR,
+        }
+    }
+}
+
+//
+// **** Tests ****
+//
+#[cfg(test)]
+mod tests {
+    use assert_approx_eq::assert_approx_eq;
+
+    use super::*;
+
+    #[test]
+    fn test_convergence() {
+        let mut a = NetworkCoordinate::<3>::new();
+        let mut b = NetworkCoordinate::<3>::new();
+        let t = Duration::from_millis(250);
+        (0..20).for_each(|_| {
+            a.update(&b, t);
+            b.update(&a, t);
+        });
+        let rtt = a.estimated_rtt(&b);
+        assert_approx_eq!(rtt.as_secs_f32() * 1000.0, 250.0, 1.0);
+    }
+
+    #[test]
+    fn test_mini_network() {
+        // define a little network with these nodes:
+        //
+        // slc has 80ms stem time to core entry Seattle
+        // nyc has 30ms stem time to core entry Virginia
+        // lax has 15ms stem time to core entry Los Angeles
+        // mad has 60ms stem time to core entry London
+        //
+        // we'll assume that traffic in the core moves at 50% the speed of light
+        //
+        // that gives us this grid of RTTs (in ms) for the core:
+        //
+        // |             | Seattle | Virgina | Los Angeles | London |
+        // |-------------|---------|---------|-------------|--------|
+        // | Seattle     |       - |      52 |          20 |    102 |
+        // | Virginia    |         |       - |          50 |     78 |
+        // | Los Angeles |         |         |           - |    116 |
+        // | London      |         |         |             |      - |
+        //
+        // Which gives us these routes (plus their reverse) and times (ms):
+        //
+        // SLC -> Seattle -> Virginia -> NYC = 80 + 52 + 30 = 162
+        // SLC -> Seattle -> Los Angeles -> LAX = 80 + 20 + 15 = 115
+        // SLC -> Seattle -> Londong -> MAD = 80 + 102 + 60 = 242
+        // NYC -> Virginia -> Los Angeles -> LAX = 30 + 50 + 15 = 95
+        // NYC -> Virginia -> London -> MAD = 30 + 78 + 60 = 168
+        // LAX -> Los Angeles -> London -> MAD = 15 + 116 + 60 = 192
+
+        // create the NCs for each endpoint
+        let mut slc = NetworkCoordinate::<2>::new();
+        let mut nyc = NetworkCoordinate::<2>::new();
+        let mut lax = NetworkCoordinate::<2>::new();
+        let mut mad = NetworkCoordinate::<2>::new();
+
+        // verify the initial error
+        let error =
+            (slc.error.powf(2.0) + nyc.error.powf(2.0) + lax.error.powf(2.0) + mad.error.powf(2.0))
+                .sqrt();
+        println!("error = {error}");
+        assert_eq!(error, 2_000.0);
+
+        // iterate plenty of times to converge and minimize error
+        (0..50).for_each(|_| {
+            slc.update(&nyc, Duration::from_millis(162));
+            nyc.update(&slc, Duration::from_millis(162));
+
+            slc.update(&lax, Duration::from_millis(115));
+            lax.update(&slc, Duration::from_millis(115));
+
+            slc.update(&mad, Duration::from_millis(242));
+            mad.update(&slc, Duration::from_millis(242));
+
+            nyc.update(&lax, Duration::from_millis(95));
+            lax.update(&nyc, Duration::from_millis(95));
+
+            nyc.update(&mad, Duration::from_millis(168));
+            mad.update(&nyc, Duration::from_millis(168));
+
+            lax.update(&mad, Duration::from_millis(192));
+            mad.update(&lax, Duration::from_millis(192));
+        });
+
+        // compute and test the root mean squared error
+        let error =
+            (slc.error.powf(2.0) + nyc.error.powf(2.0) + lax.error.powf(2.0) + mad.error.powf(2.0))
+                .sqrt();
+        println!("error = {error}");
+        assert!(error < 0.07);
+    }
+}

src/vector.rs

@@ -0,0 +1,186 @@
+//! A minimal Euclidean/Cartesian vector implementation.
+//!
+//! This is the simplest possible implementation of a vector to support the needs of a
+//! Vivaldi `HeightVector`.
+
+use num_traits::Float;
+
+use std::ops::{Add, Div, Index, Mul, Sub};
+
+//
+// **** Structs ****
+//
+
+/// A minimal Euclidean/Cartesian vector implementation
+///
+/// ## Generic parameters:
+/// - `T`: the type to use, must satisfy the `Float` trait (i.e. `f32` or `f64`)
+/// - `N`: the number of dimensions (i.e. 2 for 2D vectors, 3 for 3D, etc)
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+pub(crate) struct Vector<T: Float, const N: usize>([T; N]);
+
+//
+// **** Implementations ****
+//
+
+impl<T: Float, const N: usize> Vector<T, N> {
+    /// Create a new Zero vector.
+    pub(crate) fn new() -> Self {
+        Self::default()
+    }
+
+    /// Compute magnitude of vector.
+    pub(crate) fn len(&self) -> T {
+        (0..N)
+            .map(|i| self.0[i] * self.0[i])
+            .fold(T::zero(), |acc, x| acc + x)
+            .sqrt()
+    }
+}
+
+//
+// **** Trait Implementations ****
+//
+
+impl<T: Float, const N: usize> Default for Vector<T, N> {
+    /// A default vector is all zeroes
+    fn default() -> Self {
+        Self([T::zero(); N])
+    }
+}
+
+impl<T: Float, const N: usize> From<[T; N]> for Vector<T, N> {
+    /// Create a new Vector from an array of type T, N dimensions.
+    fn from(value: [T; N]) -> Self {
+        Self(value)
+    }
+}
+
+impl<T: Float, const N: usize> Add for Vector<T, N> {
+    type Output = Self;
+
+    /// Vector+Vector addition.
+    fn add(self, rhs: Self) -> Self::Output {
+        let mut ret = [T::zero(); N]; //Self::Output::new();
+        for i in 0..N {
+            ret[i] = self[i] + rhs[i];
+        }
+        Self::Output::from(ret)
+    }
+}
+
+impl<T: Float, const N: usize> Sub for Vector<T, N> {
+    type Output = Self;
+
+    /// Vector-Vector subtraction.
+    fn sub(self, rhs: Self) -> Self::Output {
+        let mut ret = [T::zero(); N]; //Self::Output::new();
+        for i in 0..N {
+            ret[i] = self[i] - rhs[i];
+        }
+        Self::Output::from(ret)
+    }
+}
+
+impl<T: Float, const N: usize> Mul<T> for Vector<T, N> {
+    type Output = Self;
+
+    /// Vector*scalar multiplication.
+    fn mul(self, rhs: T) -> Self::Output {
+        let mut ret = Self::Output::new();
+        for i in 0..N {
+            ret.0[i] = self.0[i] * rhs
+        }
+        ret
+    }
+}
+
+impl<T: Float, const N: usize> Div<T> for Vector<T, N> {
+    type Output = Self;
+
+    /// Vector/scalar division.
+    fn div(self, rhs: T) -> Self::Output {
+        let mut ret = Self::Output::new();
+        for i in 0..N {
+            ret.0[i] = self.0[i] / rhs
+        }
+        ret
+    }
+}
+
+impl<T: Float, const N: usize> Index<usize> for Vector<T, N> {
+    type Output = T;
+
+    /// Nicer indexing [] for read only references.
+    fn index(&self, index: usize) -> &Self::Output {
+        &self.0[index]
+    }
+}
+
+//
+// **** Tests ****
+//
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    use assert_approx_eq::assert_approx_eq;
+
+    #[test]
+    fn test_default_new() {
+        let a = Vector::<f32, 3>::default();
+        let b = Vector::<f32, 3>::new();
+        let c = Vector::<f32, 3>::from([0.0; 3]);
+        assert_eq!(a, b);
+        assert_eq!(b, c);
+    }
+
+    #[test]
+    fn test_add() {
+        let a = Vector::<f32, 3>::from([1.0, 2.0, 3.0]);
+        let b = Vector::<f32, 3>::from([3.0, 2.0, 1.0]);
+        let c = a + b;
+        assert_eq!(c, Vector::<f32, 3>::from([4.0, 4.0, 4.0]));
+    }
+
+    #[test]
+    fn test_sub() {
+        let a = Vector::<f32, 3>::from([1.0, 2.0, 3.0]);
+        let b = Vector::<f32, 3>::from([3.0, 2.0, 1.0]);
+        let c = a - b;
+        assert_eq!(c, Vector::<f32, 3>::from([-2.0, 0.0, 2.0]));
+    }
+
+    #[test]
+    fn test_mul_scalar() {
+        let a = Vector::<f32, 3>::from([1.0, 2.0, 3.0]);
+        let b = 10.0;
+        let c = a * b;
+        assert_eq!(c, Vector::<f32, 3>::from([10.0, 20.0, 30.0]));
+    }
+
+    #[test]
+    fn test_div_scalar() {
+        let a = Vector::<f32, 3>::from([1.0, 2.0, 3.0]);
+        let b = 10.0;
+        let c = a / b;
+        assert_eq!(c, Vector::<f32, 3>::from([0.1, 0.2, 0.3]));
+    }
+
+    #[test]
+    fn test_index() {
+        let a = Vector::<f32, 3>::from([1.0, 2.0, 3.0]);
+        assert_eq!(a[0], 1.0);
+        assert_eq!(a[1], 2.0);
+        assert_eq!(a[2], 3.0);
+    }
+    #[test]
+    fn test_length() {
+        let a = Vector::<f32, 3>::from([1.0, 2.0, 3.0]);
+        assert_approx_eq!(a.len(), 3.741_65, 0.0001);
+
+        let b = a * 10.0;
+        assert_approx_eq!(b.len(), 37.416_57, 0.0001);
+    }
+}