fix: add numeric stability

Made `HeightVector` numerically safe; while `Vector` remains numerically unsafe (i.e. it can divide by 0, NaN, Inf, etc).

When HeightVector runs into a problem, it just re-randomizes itself. Since the algorithm converges relatively quickly, this seems to work fine.

Also fixed a bug so that now height will always start as a positive number (and stay that way, else get re-randomized).

Added safety to `NetworkCoordinate` which now checks that its error didn’t pull a `HeightVector` into a negative height; in which case it also re-randomizes.

Clarified the `random` construction, tweaked the DEFAULT_ERROR constant, and updated tests accordingly.

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

proptest-regressions/vector.txt

@@ -7,3 +7,4 @@
 cc 9746321b74151b1897633df1bb1acba31bd29b780939385d9bece281322379fb # shrinks to x = 0.0, y = 0.0, z = 0.0
 cc 92f14e65f9256a1213cabcecf6f45152c96f39193a1bdedb7a5cdfcaf74b8550 # shrinks to x = 1.345588e30, y = 0.0, z = 0.0, m = 8.635428e-39
 cc 6bd8ea879a7c332aea9f429cf4d6d9bd61c96661edfe8875c733aea750e0c5d8 # shrinks to x = 0.0, y = 0.0, z = 9.319843e-7, d = 0.0
+cc 0c5776cf71a7170fda31bd586efd8d07d8fba2c24431a58f5b722bd9a7f5621a # shrinks to x = 0.0, y = 0.0, z = 0.0, d = -0.0

src/height_vector.rs

@@ -70,14 +70,19 @@ pub(crate) struct HeightVector<const N: usize>(Vector<f32, N>, f32);
 //
 
 impl<const N: usize> HeightVector<N> {
-    /// A new height vector is a random unit vector plus a default height.
+    /// A new height vector is a random unit vector
     pub(crate) fn new() -> Self {
+        Self::random()
+    }
+
+    /// A new height vector is a random unit vector
+    pub(crate) fn random() -> 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>();
+        let height = rng.generate::<f32>().abs();
         Self(Vector::<f32, N>::from(vec), height).normalized()
     }
 
@@ -99,6 +104,22 @@ impl<const N: usize> HeightVector<N> {
             Self(self.0 / len, self.1 / len)
         }
     }
+
+    /// Checks whether the `HeightVector` is valid.
+    ///
+    /// In this case, valid means the height is positive, and none of the components are NaN or
+    /// Inf.
+    pub(crate) fn is_valid(&self) -> bool {
+        !self.is_invalid()
+    }
+
+    /// Checks whether the `HeightVector` is invalid.
+    ///
+    /// 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.1.is_sign_negative() || self.1.is_nan() || self.1.is_infinite() || self.0.is_invalid()
+    }
 }
 
 //
@@ -115,7 +136,12 @@ impl<const N: usize> Default for HeightVector<N> {
 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)
+        let ret = Self(Vector::<f32, N>::from(value.0), value.1);
+        if ret.is_valid() {
+            ret
+        } else {
+            Self::random()
+        }
     }
 }
 
@@ -124,7 +150,12 @@ impl<const N: usize> Add for HeightVector<N> {
 
     /// Add two Vivaldi height vectors.
     fn add(self, rhs: Self) -> Self::Output {
-        Self(self.0 + rhs.0, self.1 + rhs.1)
+        let ret = Self(self.0 + rhs.0, self.1 + rhs.1);
+        if ret.is_valid() {
+            ret
+        } else {
+            Self::random()
+        }
     }
 }
 
@@ -134,7 +165,12 @@ 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 {
-        Self(self.0 - rhs.0, self.1 + rhs.1)
+        let ret = Self(self.0 - rhs.0, self.1 + rhs.1);
+        if ret.is_valid() {
+            ret
+        } else {
+            Self::random()
+        }
     }
 }
 
@@ -143,7 +179,12 @@ impl<const N: usize> Mul<f32> for HeightVector<N> {
 
     /// 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)
+        let ret = Self(self.0 * rhs, self.1 * rhs);
+        if ret.is_valid() {
+            ret
+        } else {
+            Self::random()
+        }
     }
 }
 

src/network_coordinate.rs

@@ -40,7 +40,7 @@ const C_ERROR: f32 = 0.25;
 const C_DELTA: f32 = 0.25;
 
 // initial error value
-const DEFAULT_ERROR: f32 = 1000.0;
+const DEFAULT_ERROR: f32 = 200.0;
 
 // error should always be greater than zero
 const MIN_ERROR: f32 = f32::EPSILON;
@@ -49,7 +49,7 @@ const MIN_ERROR: f32 = f32::EPSILON;
 // **** Structs ****
 //
 
-#[derive(Debug, Serialize, Deserialize)]
+#[derive(Clone, Debug, Serialize, Deserialize)]
 pub struct NetworkCoordinate<const N: usize> {
     position: HeightVector<N>,
     error: f32,
@@ -168,9 +168,21 @@ impl<const N: usize> NetworkCoordinate<N> {
         // xi = xi + δ × (rtt − ‖xi − xj ‖) × u(xi − xj)
         self.position = self.position + (self.position - rhs.position).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()
+        }
+
         // return reference to updated self
         self
     }
+
+    /// getter for error value - useful for consumers to understand the estimated accuracty of this
+    /// `NetworkCoordinate`
+    pub fn error(&self) -> f32 {
+        self.error
+    }
 }
 
 //
@@ -178,9 +190,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
     fn default() -> Self {
         Self {
-            position: HeightVector::<N>::new(),
+            position: HeightVector::<N>::random(),
             error: DEFAULT_ERROR,
         }
     }
@@ -247,11 +260,10 @@ mod tests {
         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);
+        assert_eq!(error, 400.0);
 
         // iterate plenty of times to converge and minimize error
-        (0..50).for_each(|_| {
+        (0..20).for_each(|_| {
             slc.update(&nyc, Duration::from_millis(162));
             nyc.update(&slc, Duration::from_millis(162));
 
@@ -272,11 +284,9 @@ mod tests {
         });
 
         // 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();
+        let error = slc.error + nyc.error + lax.error + mad.error;
         println!("error = {error}");
-        assert!(error < 0.2);
+        assert!(error < 5.0);
     }
 
     #[test]

src/vector.rs

@@ -48,6 +48,13 @@ where
             .fold(T::zero(), |sum_sq, x| sum_sq + x)
             .sqrt()
     }
+
+    /// Checks whether the `Vector` is invalid.
+    ///
+    /// In this case, valid means tnone of the components are NaN or Inf.
+    pub(crate) fn is_invalid(&self) -> bool {
+        self.inner.iter().any(|x| x.is_nan() || x.is_infinite())
+    }
 }
 
 //
@@ -119,9 +126,6 @@ where
         let mut ret = Self::Output::new();
         for i in 0..N {
             ret.inner[i] = self.inner[i] * rhs;
-            if !ret.inner[i].is_normal() {
-                ret.inner[i] = T::zero();
-            }
         }
         ret
     }
@@ -138,9 +142,6 @@ where
         let mut ret = Self::Output::new();
         for i in 0..N {
             ret.inner[i] = self.inner[i] / rhs;
-            if !ret.inner[i].is_normal() {
-                ret.inner[i] = T::zero();
-            }
         }
         ret
     }
@@ -167,7 +168,6 @@ mod tests {
     use super::*;
 
     use assert_approx_eq::assert_approx_eq;
-    use num_traits::Zero;
     use proptest::prelude::*;
 
     proptest! {
@@ -186,33 +186,35 @@ 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;
-            if m.is_zero() {
-                assert_eq!(w[0], 0.0);
-                assert_eq!(w[1], 0.0);
-                assert_eq!(w[2], 0.0);
-            }
-            assert!(v[0].is_finite());
-            assert!(v[1].is_finite());
-            assert!(v[1].is_finite());
-            assert!(w[0].is_finite());
-            assert!(w[1].is_finite());
-            assert!(w[1].is_finite());
+            assert_eq!(v[0] * m, w[0]);
+            assert_eq!(v[1] * m, w[1]);
+            assert_eq!(v[2] * m, w[2]);
         }
+
         #[test]
         fn proptest_div(x: f32, y: f32, z: f32, d: f32) {
             let v = Vector::<f32, 3>::from([x,y,z]);
             let w = v / d;
-            if d.is_zero() {
-                assert_eq!(w[0], 0.0);
-                assert_eq!(w[1], 0.0);
-                assert_eq!(w[2], 0.0);
+
+            let (x,y,z) = (v[0]/d, v[1]/d, v[2]/d);
+
+            if x.is_nan() {
+                assert!(w[0].is_nan());
+            } else {
+                assert_eq!(v[0] / d, x);
+            }
+
+            if y.is_nan() {
+                assert!(w[1].is_nan());
+            } else {
+                assert_eq!(v[1] / d, y);
+            }
+
+            if z.is_nan() {
+                assert!(w[2].is_nan());
+            } else {
+                assert_eq!(v[2] / d, z);
             }
-            assert!(v[0].is_finite());
-            assert!(v[1].is_finite());
-            assert!(v[1].is_finite());
-            assert!(w[0].is_finite());
-            assert!(w[1].is_finite());
-            assert!(w[1].is_finite());
         }
     }