bump crate to 0.6.0

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

.build.yml

@@ -1,23 +1,48 @@
 image: archlinux
 packages:
   - rustup
+  - cargo-tarpaulin
 sources:
   - https://git.sr.ht/~swaits/vivaldi-nc
+artifacts:
+  - vivaldi-nc/tarpaulin-report.html
 tasks:
   - setup: |
       rustup toolchain install nightly stable
       cd vivaldi-nc/
       rustup run stable cargo fetch
+  - format: |
+      rustup default stable
+      cd vivaldi-nc/
+      cargo fmt --verbose --check
+  - lint: |
+      rustup default stable
+      cd vivaldi-nc/
+      cargo clippy --verbose
+  - coverage: |
+      rustup default stable
+      cd vivaldi-nc/
+      cargo tarpaulin --verbose --all-features --out Html
   - stable: |
       rustup default stable
       cd vivaldi-nc/
       cargo build --verbose
-      cargo test --verbose
+      cargo test --verbose --all-features
   - nightly: |
       rustup default nightly
       cd vivaldi-nc/
-      cargo build --verbose ||:
+      cargo clippy --verbose ||:
-      cargo test --verbose  ||:
+      cargo build --verbose  ||:
+      cargo test --verbose --all-features ||:
+  - bench: |
+      rustup default stable
+      cd vivaldi-nc/
+      cargo bench --verbose --all-features
+  - examples: |
+      rustup default stable
+      cd vivaldi-nc/
+      cargo run --verbose --all-features --example planetlab
+      cargo run --verbose --all-features --example planetlab --release
   - docs: |
       cd vivaldi-nc/
       rustup run stable cargo doc --no-deps

.gitignore

@@ -1,2 +1,4 @@
+/flamegraph.svg
 /target
 Cargo.lock
+/tarpaulin-report.html

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "vivaldi-nc"
-version = "0.5.2"
+version = "0.6.0"
 edition = "2021"
 authors = ["Stephen Waits <steve@waits.net>"]
 description = "Vivaldi network coordinates for fast, distributed latency estimates in multinode networks, with a clean and simple interface."
@@ -16,17 +16,27 @@ exclude = [
 ]
 
 [features]
-default = ["f32"]
+default = []
 f32 = []
-f64 = []
+
+[profile.bench]
+debug = true
+
+[profile.release]
+debug = true
 
 [dev-dependencies]
 assert_approx_eq = "1.1.0"
+criterion = "0.4.0"
 proptest = "1.0.0"
 serde_json = "1.0.91"
 
+[[bench]]
+name = "million-updates"
+harness = false
+
 [dependencies]
-nanorand = {version = "0.7.0", features = ["wyrand", "getrandom"]}
 num-traits = "0.2.15"
+rand = "0.8.5"
 serde = {version = "1.0.152", features = ["serde_derive"]}
 serde_with = "2.2.0"

README.md

@@ -112,11 +112,13 @@ 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:
+By default, the internal data structures and operations all use `f64`, which is
+slightly faster on modern architectures. If you want to use `f32` instead, you
+can enable it as a cargo feature:
 
 ```toml
 [dependencies]
-vivaldi-nc = { version = "0.5.2", features = ["f64"] }
+vivaldi-nc = { version = "(version)", features = ["f32"] }
 ```
 
 ### Examples
@@ -172,10 +174,11 @@ One design goal of this crate is to minimize dependencies. When dependencies
 are required, I try to be very selective about them (re: bloat and licensing).
 This crate depends on:
 
-- [`nanorand`](https://crates.io/crates/nanorand): A fast PRNG based on WyRand.
+- [`rand`](https://crates.io/crates/rand): I started with `nanorand`; however,
-  It pulls in `getrandom` as a portable source of entropy. I chose this over
+  this is faster because of `rand::thread_rng()`; whereas with `nanorand` we
-  the more commonly used, but heavyweight `rand` because it's significantly
+  ended up needing to reimplement that same functionality through something
-  smaller and does just what I need, and not much more.
+  like `lazy_static!` plus `Mutex`, which is not worth it for a few percentage
+  points of performance.
 - [`num-traits`](https://crates.io/crates/num-traits): A brilliant crate which
   makes it easier to operate on numbers in generics; like using `Float` as a
   constraint on a generic type. Its convenience outweighs its cost.
@@ -190,7 +193,8 @@ This crate depends on:
 
 Several crates implemented Vivaldi NC before this one. So, why another?
 
-I had several design goals which the existing crates didn't satisfy:
+I had several design goals which the existing crates didn't satisfy, in order
+of priority:
 
 1. Provide the *simplest interface possible*. I just want to have some sort of
    Coordinate struct and be able to update it, and then use it to estimate
@@ -202,7 +206,10 @@ I had several design goals which the existing crates didn't satisfy:
 3. Serializable/Deserializable by default. These things are useful for sending
    across networks. The intent is to be able to do that without a bunch of
    rigamarole. Support `serde` traits by default.
-4. Well documented. Well tested. This varies by the crates.
+4. Well documented. Well tested. This varies in the existing crates.
+5. Performance. While it's a low priority, it should be reasonable. The library
+   currently does about 10,000 `NetworkCoordinate::update()` calls per
+   millisecond on my Macbook Pro (Intel i7-9750H (12) @ 2.60GHz).
 
 ### Other Vivaldi NC Implementations
 

benches/million-updates.rs

@@ -0,0 +1,52 @@
+use std::time::Duration;
+
+use criterion::{criterion_group, criterion_main, Criterion};
+use vivaldi_nc::{NetworkCoordinate2D, NetworkCoordinate3D};
+
+const NUM_NODES: usize = 1_000;
+
+fn million_vivaldi_2d_updates() {
+    // create NUM_NODES NCs
+    let mut nc: Vec<NetworkCoordinate2D> =
+        (0..NUM_NODES).map(|_| NetworkCoordinate2D::new()).collect();
+
+    // loop, update each item for every other item
+    for i in 0..NUM_NODES {
+        let nc_i = nc[i].clone();
+        for (j, nc_j) in nc.iter_mut().enumerate().take(NUM_NODES) {
+            if i == j {
+                continue;
+            }
+            let rtt = if i > j { i - j } else { j - i };
+            // let nc_j = nc[j].clone();
+            nc_j.update(&nc_i, Duration::from_millis(rtt.try_into().unwrap()));
+        }
+    }
+}
+
+fn million_vivaldi_3d_updates() {
+    // create NUM_NODES NCs
+    let mut nc: Vec<NetworkCoordinate3D> =
+        (0..NUM_NODES).map(|_| NetworkCoordinate3D::new()).collect();
+
+    // loop, update each item for every other item
+    for i in 0..NUM_NODES {
+        let nc_i = nc[i].clone();
+        for (j, nc_j) in nc.iter_mut().enumerate().take(NUM_NODES) {
+            if i == j {
+                continue;
+            }
+            let rtt = if i > j { i - j } else { j - i };
+            // let nc_j = nc[j].clone();
+            nc_j.update(&nc_i, Duration::from_millis(rtt.try_into().unwrap()));
+        }
+    }
+}
+
+fn run_benchmarks(c: &mut Criterion) {
+    c.bench_function("million 2D updates", |b| b.iter(million_vivaldi_2d_updates));
+    c.bench_function("million 3D updates", |b| b.iter(million_vivaldi_3d_updates));
+}
+
+criterion_group!(benches, run_benchmarks);
+criterion_main!(benches);

examples/planetlab.rs

@@ -2,13 +2,25 @@
 //
 // for info on the format, see `NetLatency-Data/README.md`
 
+//
+// **** 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(not(feature = "f32"))]
+type FloatType = f64;
+
 use std::{
     fs::File,
     io::{BufReader, Read},
     time::Duration,
 };
 
-use nanorand::{Rng, WyRand};
+use rand::Rng;
 use vivaldi_nc::NetworkCoordinate;
 
 // this will return a Vec<Vec<Duration>>
@@ -66,10 +78,10 @@ fn main() {
     // 4. adjust local NC
     let mut error = 0.0;
     for _ in 0..30_000 {
-        let mut rng = WyRand::new();
+        let mut rng = rand::thread_rng();
         // 1. choose a local and remote NC
-        let i_local = rng.generate::<usize>() % NUM_NODES;
+        let i_local = rng.gen_range(0..NUM_NODES);
-        let i_remote = rng.generate::<usize>() % NUM_NODES;
+        let i_remote = rng.gen_range(0..NUM_NODES);
         if i_local == i_remote {
             continue;
         }
@@ -85,7 +97,7 @@ fn main() {
         nc[i_local].update(&nc_remote, rtt_measured);
 
         // complete, now let's sum up the error
-        error = nc.iter().map(|n| n.error()).sum::<f32>() / NUM_NODES as f32;
+        error = nc.iter().map(|n| n.error()).sum::<FloatType>() / NUM_NODES as FloatType;
         assert!(error.is_finite());
         if error < 5.0 {
             break;

proptest-regressions/height_vector.txt

@@ -8,3 +8,5 @@ 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
+cc 7fa6f2f56f2398bae6804c5f0c0f3f47a7304e505c34f3addb676546cefd97f2 # shrinks to x = -0.0, y = 1.500389e38, h = 2.828499e38

src/height_vector.rs

@@ -42,20 +42,23 @@
 //! > 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 rand::prelude::*;
 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;
 
-#[cfg(feature = "f64")]
+/// `FloatType` is a type alias for either `f32` or `f64` depending on cargo features
+#[cfg(not(feature = "f32"))]
 type FloatType = f64;
 
 //
@@ -72,9 +75,15 @@ type FloatType = f64;
 ///
 /// - `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> {
+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,
 }
 
@@ -90,12 +99,12 @@ impl<const N: usize> HeightVector<N> {
 
     /// A new height vector is a random unit vector
     pub(crate) fn random() -> Self {
-        let mut rng = WyRand::new();
+        let mut rng = rand::thread_rng();
         let mut vec = [0.0; N];
         for i in vec.iter_mut().take(N) {
-            *i = rng.generate::<FloatType>();
+            *i = rng.gen::<FloatType>() - 0.5;
         }
-        let height = rng.generate::<FloatType>().abs();
+        let height = rng.gen::<FloatType>().abs();
         Self {
             position: Vector::<FloatType, N>::from(vec),
             height,
@@ -113,15 +122,16 @@ impl<const N: usize> HeightVector<N> {
     /// by the inveerse of its length.
     pub(crate) fn normalized(&self) -> Self {
         let len = self.len();
-        if len < FloatType::EPSILON {
+        // scale it by the inverse of its magnitude like any normal vector
+        let ret = Self {
+            position: self.position / len,
+            height: self.height / len,
+        };
+        if ret.is_valid() {
+            ret
+        } else {
             // 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 {
-                position: self.position / len,
-                height: self.height / len,
-            }
         }
     }
 
@@ -153,7 +163,7 @@ 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 {
-            position: Default::default(),
+            position: vector::Vector::default(),
             height: Default::default(),
         }
     }
@@ -238,14 +248,19 @@ mod tests {
     proptest! {
         #[test]
         fn proptest_len(x: FloatType, y: FloatType, h: FloatType) {
-            let len = ((x*x) + (y*y)).sqrt() + h.abs();
+            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);
+                if len.is_nan() || len.is_infinite() {
+                    assert!(a.len().is_nan() || a.len().is_infinite());
+                } else {
+                    assert_approx_eq!(a.len(), len);
+                }
             }
         }
 
@@ -258,9 +273,9 @@ mod tests {
             let a = HeightVector::<2>::from(([fx0 ,fy0 ],fh0 ));
             let b = HeightVector::<2>::from(([fx1 ,fy1 ],fh1 ));
             let c = a+b;
-            assert_eq!(c.position[0], fx0 + fx1);
+            assert_approx_eq!(c.position[0], fx0 + fx1);
-            assert_eq!(c.position[1], fy0 + fy1);
+            assert_approx_eq!(c.position[1], fy0 + fy1);
-            assert_eq!(c.height, fh0 + fh1);
+            assert_approx_eq!(c.height, fh0 + fh1);
         }
 
         #[test]
@@ -272,9 +287,9 @@ mod tests {
             let a = HeightVector::<2>::from(([fx0 ,fy0 ],fh0 ));
             let b = HeightVector::<2>::from(([fx1 ,fy1 ],fh1 ));
             let c = a-b;
-            assert_eq!(c.position[0], fx0 - fx1);
+            assert_approx_eq!(c.position[0], fx0 - fx1);
-            assert_eq!(c.position[1], fy0 - fy1);
+            assert_approx_eq!(c.position[1], fy0 - fy1);
-            assert_eq!(c.height, fh0 + fh1);
+            assert_approx_eq!(c.height, fh0 + fh1);
         }
 
         #[test]
@@ -285,11 +300,11 @@ mod tests {
             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.position[0], fx * fm);
+                assert_approx_eq!(b.position[0], fx * fm);
-                assert_eq!(b.position[1], fy * fm);
+                assert_approx_eq!(b.position[1], fy * fm);
-                assert_eq!(b.height, fh * fm);
+                assert_approx_eq!(b.height, fh * fm);
             }
         }
     }
@@ -317,9 +332,9 @@ mod tests {
     #[test]
     fn test_from() {
         let a = HeightVector::<2>::from(([1.0, 2.0], 3.0));
-        assert_eq!(a.position[0], 1.0);
+        assert_approx_eq!(a.position[0], 1.0);
-        assert_eq!(a.position[1], 2.0);
+        assert_approx_eq!(a.position[1], 2.0);
-        assert_eq!(a.height, 3.0);
+        assert_approx_eq!(a.height, 3.0);
     }
 
     #[test]
@@ -327,9 +342,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.position[0], 4.0);
+        assert_approx_eq!(c.position[0], 4.0);
-        assert_eq!(c.position[1], 4.0);
+        assert_approx_eq!(c.position[1], 4.0);
-        assert_eq!(c.height, 4.0);
+        assert_approx_eq!(c.height, 4.0);
     }
 
     #[test]
@@ -337,17 +352,17 @@ 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.position[0], -2.0);
+        assert_approx_eq!(c.position[0], -2.0);
-        assert_eq!(c.position[1], 0.0);
+        assert_approx_eq!(c.position[1], 0.0);
-        assert_eq!(c.height, 4.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.position[0], 10.0);
+        assert_approx_eq!(a.position[0], 10.0);
-        assert_eq!(a.position[1], 20.0);
+        assert_approx_eq!(a.position[1], 20.0);
-        assert_eq!(a.height, 30.0);
+        assert_approx_eq!(a.height, 30.0);
     }
 
     #[test]
@@ -363,4 +378,36 @@ mod tests {
         let t = serde_json::to_string(&a);
         assert_eq!(t.as_ref().unwrap(), s);
     }
+
+    #[test]
+    fn test_sub_invalid() {
+        let valid = HeightVector::<2>::from(([1.0, 2.0], 3.0));
+        let mut invalid = HeightVector::<2>::from(([1.0, 2.0], 3.0));
+        // force invalidity here - because `::from()` actually also catches this
+        invalid.height = 1.0 / 0.0;
+        let result = valid - invalid;
+        assert_approx_eq!(result.len(), 1.0);
+    }
+
+    #[test]
+    fn test_add_invalid() {
+        let valid = HeightVector::<2>::from(([1.0, 2.0], 3.0));
+        let mut invalid = HeightVector::<2>::from(([1.0 / 0.0, 2.0], 3.0));
+
+        // make sure `from()` caught and replaced the invalid vec
+        assert_approx_eq!(invalid.len(), 1.0);
+
+        // force invalidity here - because `::from()` actually also catches this
+        invalid.height = 1.0 / 0.0;
+        let result = valid + invalid;
+        assert_approx_eq!(result.len(), 1.0);
+    }
+
+    #[test]
+    fn test_zero_norm() {
+        let a = HeightVector::<2>::from(([0.0, 0.0], 0.0));
+        let b = a.normalized();
+        let len = b.len();
+        assert_approx_eq!(len, 1.0);
+    }
 }

src/lib.rs

@@ -93,13 +93,21 @@
 //! ```
 //!
 
+#![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;

src/network_coordinate.rs

@@ -1,29 +1,6 @@
-//! A Vivaldi [`NetworkCoordinate`]
+//! Main interface module for Vivaldi network coordinates.
-//!
-//! 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)
-//! ```
 //!
+//! For usage explanation and examples, please see the main [`crate`] documentation.
 
 use std::time::Duration;
 
@@ -37,7 +14,7 @@ use crate::height_vector::HeightVector;
 #[cfg(feature = "f32")]
 type FloatType = f32;
 
-#[cfg(feature = "f64")]
+#[cfg(not(feature = "f32"))]
 type FloatType = f64;
 
 //
@@ -58,6 +35,21 @@ 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> {
     #[serde(flatten)]
@@ -70,13 +62,13 @@ pub struct NetworkCoordinate<const N: usize> {
 /// 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>`.
+/// `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>`.
+/// `NetworkCoordinate<3>`. For more information, see [`NetworkCoordinate`].
 pub type NetworkCoordinate3D = NetworkCoordinate<3>;
 
 //
@@ -97,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()
     }
@@ -128,12 +121,13 @@ 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
         #[cfg(feature = "f32")]
         return Duration::from_secs_f32((self.heightvec - rhs.heightvec).len() / 1000.0);
 
-        #[cfg(feature = "f64")]
+        #[cfg(not(feature = "f32"))]
         return Duration::from_secs_f64((self.heightvec - rhs.heightvec).len() / 1000.0);
     }
 
@@ -167,6 +161,33 @@ 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")]
@@ -174,14 +195,17 @@ impl<const N: usize> NetworkCoordinate<N> {
         #[cfg(feature = "f32")]
         let rtt_estimated_ms = self.estimated_rtt(rhs).as_secs_f32() * 1000.0;
 
-        #[cfg(feature = "f64")]
+        #[cfg(not(feature = "f32"))]
         let rtt_ms = rtt.as_secs_f64() * 1000.0;
-        #[cfg(feature = "f64")]
+        #[cfg(not(feature = "f32"))]
         let rtt_estimated_ms = self.estimated_rtt(rhs).as_secs_f64() * 1000.0;
 
         // rtt needs to be positive
         if rtt_ms < 0.0 {
-            return self;
+            // Note: `rtt` is guaranteed to be positive because `Duration` enforces it.
+            //       If this panics, something changed where `Duration` now allows for negative
+            //       values.
+            unreachable!();
         }
 
         // Sample weight balances local and remote error. (1)
@@ -195,7 +219,11 @@ 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
@@ -207,7 +235,11 @@ impl<const N: usize> NetworkCoordinate<N> {
         // if we ended up with an invalid coordinate, return a new random coordinate with default
         // error
         if self.heightvec.is_invalid() {
-            *self = Self::new()
+            *self = Self::new();
+
+            // We should never get here because the call to `normalized()` above should catch any
+            // invalid `heightvec`
+            unreachable!();
         }
 
         // return reference to updated self
@@ -216,7 +248,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) -> FloatType {
+    #[must_use]
+    pub const fn error(&self) -> FloatType {
         self.error
     }
 }
@@ -226,7 +259,7 @@ 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 {
             heightvec: HeightVector::<N>::random(),
@@ -293,10 +326,8 @@ mod tests {
         let mut mad = NetworkCoordinate::<2>::new();
 
         // verify the initial error
-        let error =
+        let error = slc.error.hypot(nyc.error.hypot(lax.error.hypot(mad.error)));
-            (slc.error.powf(2.0) + nyc.error.powf(2.0) + lax.error.powf(2.0) + mad.error.powf(2.0))
+        assert_approx_eq!(error, 400.0);
-                .sqrt();
-        assert_eq!(error, 400.0);
 
         // iterate plenty of times to converge and minimize error
         (0..20).for_each(|_| {
@@ -333,7 +364,7 @@ mod tests {
 
         // make sure it's the right length and works like we expect a normal NC
         assert_approx_eq!(a.heightvec.len(), 2.649_509, 0.001);
-        assert_eq!(a.error, 1.0);
+        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
@@ -352,4 +383,11 @@ mod tests {
         let estimate = a.estimated_rtt(&b);
         assert_approx_eq!(estimate.as_secs_f32(), 0.080_099);
     }
+
+    #[test]
+    fn test_error_getter() {
+        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();
+        assert_approx_eq!(a.error(), 1.0);
+    }
 }

src/vector.rs

@@ -20,7 +20,7 @@ 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>,
 {
@@ -44,10 +44,8 @@ where
 
     /// Compute magnitude of vector.
     pub(crate) fn len(&self) -> T {
-        (0..N)
+        // use hypot() instead of sqrt() of sum of squares to avoid overflows
-            .map(|i| self.inner[i] * self.inner[i])
+        self.inner.iter().fold(T::zero(), |acc, x| acc.hypot(*x))
-            .fold(T::zero(), |sum_sq, x| sum_sq + x)
-            .sqrt()
     }
 
     /// Checks whether the `Vector` is invalid.
@@ -175,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_approx_eq!(v[0], x);
-            assert_eq!(v[1], y);
+            assert_approx_eq!(v[1], y);
-            assert_eq!(v[2], z);
+            assert_approx_eq!(v[2], z);
             assert!(v[0].is_finite());
             assert!(v[1].is_finite());
             assert!(v[1].is_finite());
@@ -187,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]);
+            let (x,y,z) = (v[0]*m, v[1]*m, v[2]*m);
-            assert_eq!(v[2] * m, w[2]);
+            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]
@@ -198,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() || x.is_infinite() {
-            if x.is_nan() {
+                assert!(w[0].is_nan() ||w[0].is_infinite());
-                assert!(w[0].is_nan());
             } else {
-                assert_eq!(v[0] / d, x);
+                assert_approx_eq!(w[0], x);
             }
 
-            if y.is_nan() {
+            if y.is_nan() || y.is_infinite() {
-                assert!(w[1].is_nan());
+                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() {
+            if z.is_nan() || z.is_infinite() {
-                assert!(w[2].is_nan());
+                assert!(w[2].is_nan() || w[2].is_infinite());
             } else {
-                assert_eq!(v[2] / d, z);
+                assert_approx_eq!(w[2], z);
             }
         }
     }
@@ -263,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_approx_eq!(a[0], 1.0);
-        assert_eq!(a[1], 2.0);
+        assert_approx_eq!(a[1], 2.0);
-        assert_eq!(a[2], 3.0);
+        assert_approx_eq!(a[2], 3.0);
     }
     #[test]
     fn test_length() {