bump crate to 0.5.3

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

.build.yml

@@ -11,13 +11,16 @@ tasks:
   - stable: |
       rustup default stable
       cd vivaldi-nc/
+      cargo fmt --verbose --check
+      cargo clippy --verbose
       cargo build --verbose
       cargo test --verbose
   - nightly: |
       rustup default nightly
       cd vivaldi-nc/
-      cargo build --verbose ||:
-      cargo test --verbose  ||:
+      cargo clippy --verbose ||:
+      cargo build --verbose  ||:
+      cargo test --verbose   ||:
   - docs: |
       cd vivaldi-nc/
       rustup run stable cargo doc --no-deps

Cargo.toml

@@ -1,9 +1,9 @@
 [package]
 name = "vivaldi-nc"
-version = "0.4.1"
+version = "0.5.3"
 edition = "2021"
 authors = ["Stephen Waits <steve@waits.net>"]
-description = "Vivaldi network coordinates with a simple interface."
+description = "Vivaldi network coordinates for fast, distributed latency estimates in multinode networks, with a clean and simple interface."
 license = "MIT"
 keywords = ["vivaldi", "network", "coordinates", "distributed", "p2p"]
 categories = ["algorithms", "network-programming"]
@@ -15,6 +15,11 @@ exclude = [
 	"examples/*",
 ]
 
+[features]
+default = ["f32"]
+f32 = []
+f64 = []
+
 [dev-dependencies]
 assert_approx_eq = "1.1.0"
 proptest = "1.0.0"

README.md

@@ -110,28 +110,49 @@ and compact binary formats like [bincode](https://crates.io/crates/bincode) and
 See the [module documentation](https://docs.rs/vivaldi-nc) for more detailed
 usage examples.
 
+### Cargo Features
+
+By default, the internal data structures and operations all use `f32`. If you want to use `f64`, you can enable it as a cargo feature:
+
+```toml
+[dependencies]
+vivaldi-nc = { version = "0.5.2", features = ["f64"] }
+```
+
 ### Examples
 
 The repository includes an example which loads a n-to-n latency sample from
 PlanetLab and iterates on `NetworkCoordinate`s until some low-enough mean error
-is reached. The output is a JSON array of elements that look like this:
+is reached. The output is a JSON array of elements which contain the NC's
+position, its height (or stem latency estimate), and its estimation of error
+(lower is better) that look like this:
 
 ```json
-  {
-    "position": [
-      {
-        "inner": [     # Euclidean position (3D in this case):
-          9.249476,      # first dimension; "x"
-          23.747225,     # second dimension; "x"
-          -6.473008      # third dimension; "x"
-        ]
-      },
-      72.10599         # Height component (estimate of stem latency in milliseconds)
-    ],
-    "error": 2.9166403 # Estimated error (lower is better)
-  },
-
+{
+  "position": [
+    5.563593,
+    -2.332495,
+    7.3957834
+  ],
+  "height": 55.56651,
+  "error": 3.241348
+}
 ```
+Adding some detail to these three fields:
+
+- `position`: Estimated position of this node (or `NetworkCoordinate`) on the
+  network core. Think of this like the node's endpoint into the Internet
+  backbone. This is a cartesian coordinate in a multi-dimensional latency space
+  (in this case 3D), measured in milliseconds.
+- `height`: Estimated stem time, which is the time from the node itself to the
+  backbone. So, for example this might represent the latency between your
+  node's actual location, like your home, and the point at which it gets into
+  the highest-speed, core part of the Internet. Height helps adjust for
+  [triangle inequality violations](https://en.wikipedia.org/wiki/Triangle_inequality),
+  which are common on the Internet.
+- `error`: Estimated error of the current `position`, and `height`. Lower is
+  better.
+
 To run the example, first clone the repository locally:
 
 ```bash

proptest-regressions/height_vector.txt

@@ -8,3 +8,4 @@ cc 4a5cb6ca844d5306203ae00a99375c83a20ddf4ca427f0ae0b541697efa3d9b7 # shrinks to
 cc 38ed1e2c313317809f5469d856d2c1c61e4851c64a156cd1bdd0757fa7c77128 # shrinks to x = 0.0, y = -0.0, h = 0.0, m = -0.0
 cc b7d1b5e9fe3e33066b583c4532865373179972e9b5786ecc62438f23bc939e87 # shrinks to x = -0.0, y = 0.0, h = -3.448118e-33, m = -0.0
 cc 3c4fc4f21c82ffe5bcbc740ca9ab598e0894226acb86a3046a3a24cbbcdbf8e6 # shrinks to x = 0.0, y = 0.0, h = -0.0
+cc d1a551d7eb92bf69e3d0cd0b677baac5aaa5a95d6a016fccf67cf57af71c8b0c # shrinks to x = -5.098216e26, y = 0.0, h = 0.0

src/height_vector.rs

@@ -42,12 +42,24 @@
 //! > 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 core::ops::{Add, Mul, Sub};
 
 use nanorand::{Rng, WyRand};
 use serde::{Deserialize, Serialize};
 
-use crate::vector::Vector;
+use crate::vector::{self, Vector};
+
+//
+// **** Features ****
+//
+
+/// `FloatType` is a type alias for either `f32` or `f64` depending on cargo features
+#[cfg(feature = "f32")]
+type FloatType = f32;
+
+/// `FloatType` is a type alias for either `f32` or `f64` depending on cargo features
+#[cfg(feature = "f64")]
+type FloatType = f64;
 
 //
 // **** Constants ****
@@ -63,7 +75,17 @@ use crate::vector::Vector;
 ///
 /// - `N`: the dimensionality of the vector portion (i.e. non-height) of the Vivaldi height vector
 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
-pub(crate) struct HeightVector<const N: usize>(Vector<f32, N>, f32);
+pub struct HeightVector<const N: usize> {
+    /// `position` is the Euclidean coordinate part of the `HeightVector`, representing a position
+    /// in the network's latency space
+    #[serde(flatten)]
+    position: Vector<FloatType, N>,
+
+    /// `height` is a representation of a node's stemp time, or its latency to get into the network
+    /// core (or backbone). It's a way to account for triangle inequality violations.
+    /// in the network's latency space
+    height: FloatType,
+}
 
 //
 // **** Implementations ****
@@ -80,28 +102,35 @@ impl<const N: usize> HeightVector<N> {
         let mut rng = WyRand::new();
         let mut vec = [0.0; N];
         for i in vec.iter_mut().take(N) {
-            *i = rng.generate::<f32>();
+            *i = rng.generate::<FloatType>();
         }
-        let height = rng.generate::<f32>().abs();
-        Self(Vector::<f32, N>::from(vec), height).normalized()
+        let height = rng.generate::<FloatType>().abs();
+        Self {
+            position: Vector::<FloatType, 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
+    pub(crate) fn len(&self) -> FloatType {
+        self.position.len() + self.height
     }
 
     /// 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 len < FloatType::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)
+            Self {
+                position: self.position / len,
+                height: self.height / len,
+            }
         }
     }
 
@@ -118,7 +147,10 @@ impl<const N: usize> HeightVector<N> {
     /// In this case, valid means the height is positive, and none of the components are NaN or
     /// Inf.
     pub(crate) fn is_invalid(&self) -> bool {
-        self.0.is_invalid() || self.1.is_nan() || self.1.is_infinite() || self.1 < 0.0
+        self.position.is_invalid()
+            || self.height.is_nan()
+            || self.height.is_infinite()
+            || self.height < 0.0
     }
 }
 
@@ -129,14 +161,20 @@ impl<const N: usize> HeightVector<N> {
 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())
+        Self {
+            position: vector::Vector::default(),
+            height: 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 {
-        let ret = Self(Vector::<f32, N>::from(value.0), value.1);
+impl<const N: usize> From<([FloatType; N], FloatType)> for HeightVector<N> {
+    /// Convert from a `([FloatType; N], FloatType)` (vector, height) to a Vivaldi height vector type.
+    fn from(value: ([FloatType; N], FloatType)) -> Self {
+        let ret = Self {
+            position: Vector::<FloatType, N>::from(value.0),
+            height: value.1,
+        };
         if ret.is_valid() {
             ret
         } else {
@@ -150,7 +188,10 @@ impl<const N: usize> Add for HeightVector<N> {
 
     /// Add two Vivaldi height vectors.
     fn add(self, rhs: Self) -> Self::Output {
-        let ret = Self(self.0 + rhs.0, self.1 + rhs.1);
+        let ret = Self {
+            position: self.position + rhs.position,
+            height: self.height + rhs.height,
+        };
         if ret.is_valid() {
             ret
         } else {
@@ -165,7 +206,10 @@ impl<const N: usize> Sub for HeightVector<N> {
     /// 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 {
-        let ret = Self(self.0 - rhs.0, self.1 + rhs.1);
+        let ret = Self {
+            position: self.position - rhs.position,
+            height: self.height + rhs.height,
+        };
         if ret.is_valid() {
             ret
         } else {
@@ -174,12 +218,15 @@ impl<const N: usize> Sub for HeightVector<N> {
     }
 }
 
-impl<const N: usize> Mul<f32> for HeightVector<N> {
+impl<const N: usize> Mul<FloatType> 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 {
-        let ret = Self(self.0 * rhs, self.1 * rhs);
+    fn mul(self, rhs: FloatType) -> Self::Output {
+        let ret = Self {
+            position: self.position * rhs,
+            height: self.height * rhs,
+        };
         if ret.is_valid() {
             ret
         } else {
@@ -199,59 +246,60 @@ mod tests {
 
     proptest! {
         #[test]
-        fn proptest_len(x: f32, y: f32, h: f32) {
-            let len = ((x*x) + (y*y)).sqrt() + h.abs();
+        fn proptest_len(x: FloatType, y: FloatType, h: FloatType) {
+            let len = x.hypot(y) + h.abs();
             let a = HeightVector::<2>::from(([x,y],h));
+            println!("a = {a:#?}");
             if x.is_nan() || x.is_infinite() || y.is_nan() || y.is_infinite() || h.is_nan() || h.is_infinite() || h < 0.0 {
                 // we should've gottne a random univ vector here
                 assert_approx_eq!(a.len(), 1.0);
             } else {
                 // `HeightVector` we created from proptest values should be valid and have a length
-                assert_eq!(a.len(), len);
+                assert_approx_eq!(a.len(), len);
             }
         }
 
         #[test]
         fn proptest_add(x0 in -1_000_000_000..1_000_000_000i32, y0 in -1_000_000_000..1_000_000_000i32, h0 in 0..1_000_000_000i32, x1 in -1_000_000_000..1_000_000_000i32, y1 in -1_000_000_000..1_000_000_000i32, h1 in 0..1_000_000i32) {
-            // convert our integer range inputs to f32s
-            let (fx0,fy0,fh0) = (x0 as f32 / 1_000.0, y0 as f32 / 1_000.0, h0 as f32 / 1_000.0);
-            let (fx1,fy1,fh1) = (x1 as f32 / 1_000.0, y1 as f32 / 1_000.0, h1 as f32 / 1_000.0);
+            // convert our integer range inputs to FloatType
+            let (fx0,fy0,fh0) = (x0 as FloatType / 1_000.0, y0 as FloatType / 1_000.0, h0 as FloatType / 1_000.0);
+            let (fx1,fy1,fh1) = (x1 as FloatType / 1_000.0, y1 as FloatType / 1_000.0, h1 as FloatType / 1_000.0);
 
             let a = HeightVector::<2>::from(([fx0 ,fy0 ],fh0 ));
             let b = HeightVector::<2>::from(([fx1 ,fy1 ],fh1 ));
             let c = a+b;
-            assert_eq!(c.0[0], fx0 + fx1);
-            assert_eq!(c.0[1], fy0 + fy1);
-            assert_eq!(c.1, fh0 + fh1);
+            assert_approx_eq!(c.position[0], fx0 + fx1);
+            assert_approx_eq!(c.position[1], fy0 + fy1);
+            assert_approx_eq!(c.height, fh0 + fh1);
         }
 
         #[test]
         fn proptest_sub(x0 in -1_000_000_000..1_000_000_000i32, y0 in -1_000_000_000..1_000_000_000i32, h0 in 0..1_000_000_000i32, x1 in -1_000_000_000..1_000_000_000i32, y1 in -1_000_000_000..1_000_000_000i32, h1 in 0..1_000_000i32) {
-            // convert our integer range inputs to f32s
-            let (fx0,fy0,fh0) = (x0 as f32 / 1_000.0, y0 as f32 / 1_000.0, h0 as f32 / 1_000.0);
-            let (fx1,fy1,fh1) = (x1 as f32 / 1_000.0, y1 as f32 / 1_000.0, h1 as f32 / 1_000.0);
+            // convert our integer range inputs to FloatType
+            let (fx0,fy0,fh0) = (x0 as FloatType / 1_000.0, y0 as FloatType / 1_000.0, h0 as FloatType / 1_000.0);
+            let (fx1,fy1,fh1) = (x1 as FloatType / 1_000.0, y1 as FloatType / 1_000.0, h1 as FloatType / 1_000.0);
 
             let a = HeightVector::<2>::from(([fx0 ,fy0 ],fh0 ));
             let b = HeightVector::<2>::from(([fx1 ,fy1 ],fh1 ));
             let c = a-b;
-            assert_eq!(c.0[0], fx0 - fx1);
-            assert_eq!(c.0[1], fy0 - fy1);
-            assert_eq!(c.1, fh0 + fh1);
+            assert_approx_eq!(c.position[0], fx0 - fx1);
+            assert_approx_eq!(c.position[1], fy0 - fy1);
+            assert_approx_eq!(c.height, fh0 + fh1);
         }
 
         #[test]
         fn proptest_mul(x in -1_000_000_000..1_000_000_000i32, y in -1_000_000_000..1_000_000_000i32, h in 0..1_000_000_000i32, m in -1_000_000_000..1_000_000_000i32) {
-            // convert our integer range inputs to f32s
-            let (fx,fy,fh,fm) = (x as f32 / 1_000.0, y as f32 / 1_000.0, h as f32 / 1_000.0, m as f32 / 1_000.0);
+            // convert our integer range inputs to FloatType
+            let (fx,fy,fh,fm) = (x as FloatType / 1_000.0, y as FloatType / 1_000.0, h as FloatType / 1_000.0, m as FloatType / 1_000.0);
 
             let a = HeightVector::<2>::from(([fx, fy], fh));
             let b = a * fm;
             if fm < 0.0 {
-                assert_approx_eq!(b.len(), 1.0)
+                assert_approx_eq!(b.len(), 1.0);
             } else {
-                assert_eq!(b.0[0], fx * fm);
-                assert_eq!(b.0[1], fy * fm);
-                assert_eq!(b.1, fh * fm);
+                assert_approx_eq!(b.position[0], fx * fm);
+                assert_approx_eq!(b.position[1], fy * fm);
+                assert_approx_eq!(b.height, fh * fm);
             }
         }
     }
@@ -279,9 +327,9 @@ mod tests {
     #[test]
     fn test_from() {
         let a = HeightVector::<2>::from(([1.0, 2.0], 3.0));
-        assert_eq!(a.0[0], 1.0);
-        assert_eq!(a.0[1], 2.0);
-        assert_eq!(a.1, 3.0);
+        assert_approx_eq!(a.position[0], 1.0);
+        assert_approx_eq!(a.position[1], 2.0);
+        assert_approx_eq!(a.height, 3.0);
     }
 
     #[test]
@@ -289,9 +337,9 @@ mod tests {
         let a = HeightVector::<2>::from(([1.0, 2.0], 3.0));
         let b = HeightVector::<2>::from(([3.0, 2.0], 1.0));
         let c = a + b;
-        assert_eq!(c.0[0], 4.0);
-        assert_eq!(c.0[1], 4.0);
-        assert_eq!(c.1, 4.0);
+        assert_approx_eq!(c.position[0], 4.0);
+        assert_approx_eq!(c.position[1], 4.0);
+        assert_approx_eq!(c.height, 4.0);
     }
 
     #[test]
@@ -299,23 +347,23 @@ mod tests {
         let a = HeightVector::<2>::from(([1.0, 2.0], 3.0));
         let b = HeightVector::<2>::from(([3.0, 2.0], 1.0));
         let c = a - b;
-        assert_eq!(c.0[0], -2.0);
-        assert_eq!(c.0[1], 0.0);
-        assert_eq!(c.1, 4.0);
+        assert_approx_eq!(c.position[0], -2.0);
+        assert_approx_eq!(c.position[1], 0.0);
+        assert_approx_eq!(c.height, 4.0);
     }
 
     #[test]
     fn test_mul_scalar() {
         let a = HeightVector::<2>::from(([1.0, 2.0], 3.0)) * 10.0;
-        assert_eq!(a.0[0], 10.0);
-        assert_eq!(a.0[1], 20.0);
-        assert_eq!(a.1, 30.0);
+        assert_approx_eq!(a.position[0], 10.0);
+        assert_approx_eq!(a.position[1], 20.0);
+        assert_approx_eq!(a.height, 30.0);
     }
 
     #[test]
     fn test_serde() {
         // start with JSON, deserialize it
-        let s = "[{\"inner\":[1.0,2.0,3.0]},4.0]";
+        let s = "{\"position\":[1.0,2.0,3.0],\"height\":4.0}";
         let a: HeightVector<3> = serde_json::from_str(s).unwrap();
 
         // make sure it's the right length and works like we expect a normal NC

src/lib.rs

@@ -1,4 +1,5 @@
-//! Vivaldi Network Coordinates with a simple interface
+//! Vivaldi network coordinates for fast, distributed latency estimates in multinode networks, with
+//! a clean and simple interface.
 //!
 //! Network Coordinates (NC) are a way to represent a node's position in a network's
 //! latency space. Nodes with low latency (ping or round trip time) between each
@@ -19,7 +20,8 @@
 //!
 //! # Usage
 //!
-//! This crate exports a single struct, [`NetworkCoordinate`]. Typical use of Vivaldi NCs for a
+//! This crate exports a single struct, [`NetworkCoordinate`] and two type aliases
+//! ([`NetworkCoordinate2D`] and [`NetworkCoordinate3D`]). Typical use of Vivaldi NCs for a
 //! distributed network works like this:
 //!
 //! 1. Each node in the network has its own instance of [`NetworkCoordinate`]. *See "Note on
@@ -51,7 +53,7 @@
 //!
 //! // mock up a remote NC for the sake of this example (normally the node would actually receive
 //! // this from a remote node)
-//! let received_msg = "{\"position\":[{\"inner\":[1.5,0.5,2.0]},0.1],\"error\":1.0}";
+//! let received_msg = "{\"position\":[1.5,0.5,2.0],\"height\":0.1,\"error\":1.0}";
 //! let remote_nc: NetworkCoordinate<3> = serde_json::from_str(received_msg).unwrap();
 //!
 //! // mock up a measured RTT between us and that node (normally the nodes would coordinate
@@ -69,16 +71,19 @@
 //! forwarded to us by someone else. And in this case, we aren't measuring the actual RTT. Instead,
 //! we're using the NCs to estimate RTT.
 //!
+//! This example also uses the type alias [`NetworkCoordinate3D`] instead of its equivalent
+//! [`NetworkCoordinate<3>`].
+//!
 //! ```
-//! use vivaldi_nc::NetworkCoordinate;
+//! use vivaldi_nc::NetworkCoordinate3D;
 //!
 //! // create our local NC as a 3-dimensional Vivaldi coordinate
-//! let mut my_nc = NetworkCoordinate::<3>::new();
+//! let mut my_nc = NetworkCoordinate3D::new();
 //!
 //! // mock up a remote NC for the sake of this example (normally the node would actually receive
 //! // this from a remote node)
-//! let received_msg = "{\"position\":[{\"inner\":[1.5,0.5,2.0]},0.1],\"error\":1.0}";
-//! let remote_nc: NetworkCoordinate<3> = serde_json::from_str(received_msg).unwrap();
+//! let received_msg = "{\"position\":[1.5,0.5,2.0],\"height\":0.1,\"error\":1.0}";
+//! let remote_nc: NetworkCoordinate3D = serde_json::from_str(received_msg).unwrap();
 //!
 //! // estimate the RTT to this remote node without measuring it directly
 //! let rtt_estimated = my_nc.estimated_rtt(&remote_nc);
@@ -88,9 +93,27 @@
 //! ```
 //!
 
+#![deny(clippy::all)]
+#![deny(clippy::correctness)]
+#![deny(clippy::style)]
+#![deny(clippy::suspicious)]
+#![deny(clippy::complexity)]
+#![deny(clippy::perf)]
+#![deny(clippy::cargo)]
+#![deny(clippy::pedantic)]
+#![deny(clippy::nursery)]
+#![allow(clippy::type_repetition_in_bounds)]
+#![allow(clippy::module_name_repetitions)]
+
 mod height_vector;
 mod vector;
 
+// make sure we don't have conflicting features
+#[cfg(all(feature = "f32", feature = "f64"))]
+compile_error!("feature \"f32\" and feature \"f64\" cannot be enabled at the same time");
+
 // publish our interface
 pub mod network_coordinate;
 pub use network_coordinate::NetworkCoordinate;
+pub use network_coordinate::NetworkCoordinate2D;
+pub use network_coordinate::NetworkCoordinate3D;

src/network_coordinate.rs

@@ -1,29 +1,6 @@
-//! 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)
-//! ```
+//! Main interface module for Vivaldi network coordinates.
 //!
+//! For usage explanation and examples, please see the main [`crate`] documentation.
 
 use std::time::Duration;
 
@@ -31,30 +8,69 @@ use serde::{Deserialize, Serialize};
 
 use crate::height_vector::HeightVector;
 
+//
+// **** Features ****
+//
+#[cfg(feature = "f32")]
+type FloatType = f32;
+
+#[cfg(feature = "f64")]
+type FloatType = f64;
+
 //
 // **** Constants ****
 //
 
 // Vivaldi tuning parameters
-const C_ERROR: f32 = 0.25;
-const C_DELTA: f32 = 0.25;
+const C_ERROR: FloatType = 0.25;
+const C_DELTA: FloatType = 0.25;
 
 // initial error value
-const DEFAULT_ERROR: f32 = 200.0;
+const DEFAULT_ERROR: FloatType = 200.0;
 
 // error should always be greater than zero
-const MIN_ERROR: f32 = f32::EPSILON;
+const MIN_ERROR: FloatType = FloatType::EPSILON;
 
 //
 // **** Structs ****
 //
 
+/// A `NetworkCoordinate<N>` is the main interface to a Vivaldi network coordinate.
+///
+/// # Generic Parameters
+///
+/// - `N`: Const generic for number of dimensions. For example, `NetworkCoordinate<3>` is a
+/// 3-Dimentionsal Euclidean coordinate plus a height. Should be a positive number greater than
+/// zero.
+///
+/// **Note:** Dimensions other than 2D or 3D are usually not useful. If you want to use one of
+/// those dimensions, you can use type aliases ([`NetworkCoordinate2D`] or [`NetworkCoordinate3D`])
+/// which are a little more ergonomic than using the generic here.
+///
+/// # Examples
+///
+/// For an explanation and examples of usage, please see the main [`crate`] documentation.
 #[derive(Clone, Debug, Serialize, Deserialize)]
 pub struct NetworkCoordinate<const N: usize> {
-    position: HeightVector<N>,
-    error: f32,
+    #[serde(flatten)]
+    heightvec: HeightVector<N>,
+    error: FloatType,
 }
 
+// type aliases for convenience
+
+/// A 2D [`NetworkCoordinate`]. Includes a 2D Euclidean position and a height.
+///
+/// This type alias is just for convenience. It's functionally equivalent to
+/// `NetworkCoordinate<2>`. For more information, see [`NetworkCoordinate`].
+pub type NetworkCoordinate2D = NetworkCoordinate<2>;
+
+/// A 3D [`NetworkCoordinate`]. Includes a 3D Euclidean position and a height.
+///
+/// This type alias is just for convenience. It's functionally equivalent to
+/// `NetworkCoordinate<3>`. For more information, see [`NetworkCoordinate`].
+pub type NetworkCoordinate3D = NetworkCoordinate<3>;
+
 //
 // **** Implementations ****
 //
@@ -73,6 +89,7 @@ impl<const N: usize> NetworkCoordinate<N> {
     /// // print the NC
     /// println!("Our new NC is: {:#?}", a);
     /// ```
+    #[must_use]
     pub fn new() -> Self {
         Self::default()
     }
@@ -104,9 +121,14 @@ impl<const N: usize> NetworkCoordinate<N> {
     /// println!("Estimated RTT: {}", a.estimated_rtt(&b).as_millis());
     /// ```
     ///
+    #[must_use]
     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)
+        #[cfg(feature = "f32")]
+        return Duration::from_secs_f32((self.heightvec - rhs.heightvec).len() / 1000.0);
+
+        #[cfg(feature = "f64")]
+        return Duration::from_secs_f64((self.heightvec - rhs.heightvec).len() / 1000.0);
     }
 
     /// Given another Vivaldi [`NetworkCoordinate`], adjust our coordinateto better represent the actual round
@@ -139,16 +161,50 @@ impl<const N: usize> NetworkCoordinate<N> {
     /// // Now we can update our NC to adjust our position relative to the remote node:
     /// local.update(&remote, rtt);
     /// ```
+    ///
+    /// # Algorithm
+    ///
+    /// 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 xj,
+    /// // and an error estimate of ej .
+    /// //
+    /// // 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)
+    /// ```
+    ///
     pub fn update(&mut self, rhs: &Self, rtt: Duration) -> &Self {
+        // convert Durations into FloatType as fractional milliseconds for convenience
+        #[cfg(feature = "f32")]
+        let rtt_ms = rtt.as_secs_f32() * 1000.0;
+        #[cfg(feature = "f32")]
+        let rtt_estimated_ms = self.estimated_rtt(rhs).as_secs_f32() * 1000.0;
+
+        #[cfg(feature = "f64")]
+        let rtt_ms = rtt.as_secs_f64() * 1000.0;
+        #[cfg(feature = "f64")]
+        let rtt_estimated_ms = self.estimated_rtt(rhs).as_secs_f64() * 1000.0;
+
         // rtt needs to be positive
-        if rtt.as_secs_f32() < 0.0 {
+        if rtt_ms < 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);
@@ -160,18 +216,23 @@ impl<const N: usize> NetworkCoordinate<N> {
 
         // 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);
+        // self.error = (es * C_ERROR * w + self.error * (1.0 - C_ERROR * w)).max(MIN_ERROR);
+        // NOTE: using `mul_add()` which is a little safer (avoid overflows)
+        self.error = (es * C_ERROR)
+            .mul_add(w, self.error * C_ERROR.mul_add(-w, 1.0))
+            .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;
+        self.heightvec =
+            self.heightvec + (self.heightvec - rhs.heightvec).normalized() * delta * error;
 
         // if we ended up with an invalid coordinate, return a new random coordinate with default
         // error
-        if self.position.is_invalid() {
-            *self = Self::new()
+        if self.heightvec.is_invalid() {
+            *self = Self::new();
         }
 
         // return reference to updated self
@@ -180,7 +241,8 @@ impl<const N: usize> NetworkCoordinate<N> {
 
     /// getter for error value - useful for consumers to understand the estimated accuracty of this
     /// `NetworkCoordinate`
-    pub fn error(&self) -> f32 {
+    #[must_use]
+    pub const fn error(&self) -> FloatType {
         self.error
     }
 }
@@ -190,10 +252,10 @@ impl<const N: usize> NetworkCoordinate<N> {
 //
 
 impl<const N: usize> Default for NetworkCoordinate<N> {
-    /// A default `NetworkCoordinate` has a random position and DEFAULT_ERROR
+    /// A default `NetworkCoordinate` has a random position and `DEFAULT_ERROR`
     fn default() -> Self {
         Self {
-            position: HeightVector::<N>::random(),
+            heightvec: HeightVector::<N>::random(),
             error: DEFAULT_ERROR,
         }
     }
@@ -257,10 +319,8 @@ mod tests {
         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();
-        assert_eq!(error, 400.0);
+        let error = slc.error.hypot(nyc.error.hypot(lax.error.hypot(mad.error)));
+        assert_approx_eq!(error, 400.0);
 
         // iterate plenty of times to converge and minimize error
         (0..20).for_each(|_| {
@@ -292,12 +352,12 @@ mod tests {
     #[test]
     fn test_serde() {
         // start with JSON, deserialize it
-        let s = "{\"position\":[{\"inner\":[1.5,0.5,2.0]},0.1],\"error\":1.0}";
+        let s = "{\"position\":[1.5,0.5,2.0],\"height\":0.1,\"error\":1.0}";
         let a: NetworkCoordinate<3> = serde_json::from_str(s).unwrap();
 
         // make sure it's the right length and works like we expect a normal NC
-        assert_approx_eq!(a.position.len(), 2.649_509, 0.001);
-        assert_eq!(a.error, 1.0);
+        assert_approx_eq!(a.heightvec.len(), 2.649_509, 0.001);
+        assert_approx_eq!(a.error, 1.0);
         assert_eq!(a.estimated_rtt(&a).as_millis(), 0);
 
         // serialize it into a new JSON string and make sure it matches the original
@@ -308,9 +368,9 @@ mod tests {
     #[test]
     fn test_estimated_rtt() {
         // start with JSON, deserialize it
-        let s = "{\"position\":[{\"inner\":[1.5,0.5,2.0]},25.0],\"error\":1.0}";
+        let s = "{\"position\":[1.5,0.5,2.0],\"height\":25.0,\"error\":1.0}";
         let a: NetworkCoordinate<3> = serde_json::from_str(s).unwrap();
-        let s = "{\"position\":[{\"inner\":[-1.5,-0.5,-2.0]},50.0],\"error\":1.0}";
+        let s = "{\"position\":[-1.5,-0.5,-2.0],\"height\":50.0,\"error\":1.0}";
         let b: NetworkCoordinate<3> = serde_json::from_str(s).unwrap();
 
         let estimate = a.estimated_rtt(&b);

src/vector.rs

@@ -20,10 +20,11 @@ use std::ops::{Add, Div, Index, Mul, Sub};
 /// - `N`: the number of dimensions (i.e. 2 for 2D vectors, 3 for 3D, etc)
 #[serde_as]
 #[derive(Copy, Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
-pub(crate) struct Vector<T, const N: usize>
+pub struct Vector<T, const N: usize>
 where
     T: Float + Serialize + for<'d> Deserialize<'d>,
 {
+    #[serde(rename = "position")]
     #[serde_as(as = "[_; N]")]
     inner: [T; N],
 }
@@ -43,10 +44,8 @@ where
 
     /// Compute magnitude of vector.
     pub(crate) fn len(&self) -> T {
-        (0..N)
-            .map(|i| self.inner[i] * self.inner[i])
-            .fold(T::zero(), |sum_sq, x| sum_sq + x)
-            .sqrt()
+        // use hypot() instead of sqrt() of sum of squares to avoid overflows
+        self.inner.iter().fold(T::zero(), |acc, x| acc.hypot(*x))
     }
 
     /// Checks whether the `Vector` is invalid.
@@ -174,9 +173,9 @@ mod tests {
         #[test]
         fn proptest_new(x: f32, y: f32, z: f32){
             let v = Vector::<f32, 3>::from([x,y,z]);
-            assert_eq!(v[0], x);
-            assert_eq!(v[1], y);
-            assert_eq!(v[2], z);
+            assert_approx_eq!(v[0], x);
+            assert_approx_eq!(v[1], y);
+            assert_approx_eq!(v[2], z);
             assert!(v[0].is_finite());
             assert!(v[1].is_finite());
             assert!(v[1].is_finite());
@@ -186,9 +185,25 @@ mod tests {
         fn proptest_mul(x: f32, y: f32, z: f32, m: f32) {
             let v = Vector::<f32, 3>::from([x,y,z]);
             let w = v * m;
-            assert_eq!(v[0] * m, w[0]);
-            assert_eq!(v[1] * m, w[1]);
-            assert_eq!(v[2] * m, w[2]);
+
+            let (x,y,z) = (v[0]*m, v[1]*m, v[2]*m);
+            if x.is_nan() || x.is_infinite() {
+                assert!(w[0].is_nan() || w[0].is_infinite());
+            } else {
+                assert_approx_eq!(w[0], x);
+            }
+
+            if y.is_nan() || y.is_infinite() {
+                assert!(w[1].is_nan() || w[1].is_infinite());
+            } else {
+                assert_approx_eq!(w[1], y);
+            }
+
+            if z.is_nan() || z.is_infinite() {
+                assert!(w[2].is_nan() || w[2].is_infinite());
+            } else {
+                assert_approx_eq!(w[2], z);
+            }
         }
 
         #[test]
@@ -197,23 +212,22 @@ mod tests {
             let w = v / d;
 
             let (x,y,z) = (v[0]/d, v[1]/d, v[2]/d);
-
-            if x.is_nan() {
-                assert!(w[0].is_nan());
+            if x.is_nan() || x.is_infinite() {
+                assert!(w[0].is_nan() ||w[0].is_infinite());
             } else {
-                assert_eq!(v[0] / d, x);
+                assert_approx_eq!(w[0], x);
             }
 
-            if y.is_nan() {
-                assert!(w[1].is_nan());
+            if y.is_nan() || y.is_infinite() {
+                assert!(w[1].is_nan() || w[1].is_infinite());
             } else {
-                assert_eq!(v[1] / d, y);
+                assert_approx_eq!(w[1], y);
             }
 
-            if z.is_nan() {
-                assert!(w[2].is_nan());
+            if z.is_nan() || z.is_infinite() {
+                assert!(w[2].is_nan() || w[2].is_infinite());
             } else {
-                assert_eq!(v[2] / d, z);
+                assert_approx_eq!(w[2], z);
             }
         }
     }
@@ -262,9 +276,9 @@ mod tests {
     #[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);
+        assert_approx_eq!(a[0], 1.0);
+        assert_approx_eq!(a[1], 2.0);
+        assert_approx_eq!(a[2], 3.0);
     }
     #[test]
     fn test_length() {
@@ -279,7 +293,7 @@ mod tests {
     fn test_serde() {
         let a = Vector::<f32, 3>::from([1.0, 2.0, 3.0]);
         let s = serde_json::to_string(&a);
-        assert_eq!(s.as_ref().unwrap(), "{\"inner\":[1.0,2.0,3.0]}");
+        assert_eq!(s.as_ref().unwrap(), "{\"position\":[1.0,2.0,3.0]}");
         assert!(s.is_ok());
     }
 }