From 26a7b243dc079b10b36a47cbc582df578cce63ad Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Tue, 7 Feb 2023 20:25:10 -0700 Subject: [PATCH] fix: add numeric stability MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- proptest-regressions/vector.txt | 1 + src/height_vector.rs | 53 ++++++++++++++++++++++++++---- src/network_coordinate.rs | 30 +++++++++++------ src/vector.rs | 58 +++++++++++++++++---------------- 4 files changed, 98 insertions(+), 44 deletions(-) diff --git a/proptest-regressions/vector.txt b/proptest-regressions/vector.txt index 525328f..fc68e1d 100644 --- a/proptest-regressions/vector.txt +++ b/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 diff --git a/src/height_vector.rs b/src/height_vector.rs index d42f005..4b95dcb 100644 --- a/src/height_vector.rs +++ b/src/height_vector.rs @@ -70,14 +70,19 @@ pub(crate) struct HeightVector(Vector, f32); // impl HeightVector { - /// 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::(); } - let height = rng.generate::(); + let height = rng.generate::().abs(); Self(Vector::::from(vec), height).normalized() } @@ -99,6 +104,22 @@ impl HeightVector { 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 Default for HeightVector { impl From<([f32; N], f32)> for HeightVector { /// Convert from a `([f32; N], f32)` (vector, height) to a Vivaldi height vector type. fn from(value: ([f32; N], f32)) -> Self { - Self(Vector::::from(value.0), value.1) + let ret = Self(Vector::::from(value.0), value.1); + if ret.is_valid() { + ret + } else { + Self::random() + } } } @@ -124,7 +150,12 @@ impl Add for HeightVector { /// 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 Sub for HeightVector { /// 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 Mul for HeightVector { /// 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() + } } } diff --git a/src/network_coordinate.rs b/src/network_coordinate.rs index 674c7a6..12798c3 100644 --- a/src/network_coordinate.rs +++ b/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 { position: HeightVector, error: f32, @@ -168,9 +168,21 @@ impl NetworkCoordinate { // 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 NetworkCoordinate { // impl Default for NetworkCoordinate { + /// A default `NetworkCoordinate` has a random position and DEFAULT_ERROR fn default() -> Self { Self { - position: HeightVector::::new(), + position: HeightVector::::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] diff --git a/src/vector.rs b/src/vector.rs index 777a0f8..9adc8e4 100644 --- a/src/vector.rs +++ b/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::::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::::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()); } }