Merge branch 'main' into remove-39-retworkx

Author: IvanIsCoding

docs/source/api/algorithm_functions/geometry.rst

@@ -0,0 +1,10 @@
+.. _geometry:
+
+Geometry
+========
+
+.. autosummary::
+   :toctree: ../../apiref
+
+   rustworkx.hyperbolic_greedy_routing
+   rustworkx.hyperbolic_greedy_success_rate

docs/source/api/algorithm_functions/index.rst

@@ -11,6 +11,7 @@ Algorithm Functions
    connectivity_and_cycles
    dag_algorithms
    dominance
+   geometry
    graph_operations
    isomorphism
    link_analysis

releasenotes/notes/greedy-routing-2e539d3bd34c948b.yaml

@@ -0,0 +1,10 @@
+features:
+  - |
+    Adds the new ``geometry`` module to ``rustworkx_core``.
+    This module mainly defines the greedy routing algorithm which is useful when each node has a position in a metric space.
+    The ``greedy_routing`` function returns, if it exists, the greedy path between two nodes and its total length in the metric space.
+    The ``greedy_routing_success_rate`` returns the proportion of pairs of nodes for which the greedy path reaches the destination.
+    The greedy routing algorithm may be used with any distance.
+    The ``geometry`` module implements many distances: ``angular_distance``, ``euclidean_distance``, ``hyperboloid_hyperbolic_distance``, ``lp_distance``, ``maximum_distance`` and ``polar_hyperbolic_distance``.
+  - |
+    Adds the function :func:`.hyperbolic_greedy_routing` that performs the greedy routing algorithm in the hyperbolic space, and the function :func:`.hyperbolic_greedy_success_rate` that computes the proportion of pairs of nodes for which the greedy path reaches the destination.

rustworkx-core/src/generators/random_graph.rs

@@ -30,6 +30,7 @@ use rand_pcg::Pcg64;
 
 use super::star_graph;
 use super::InvalidInputError;
+use crate::geometry::hyperboloid_hyperbolic_distance;
 
 /// Generates a random regular graph
 ///
@@ -1004,7 +1005,7 @@ where
     let between = Uniform::new(0.0, 1.0).unwrap();
     for (v, p1) in pos.iter().enumerate().take(num_nodes - 1) {
         for (w, p2) in pos.iter().enumerate().skip(v + 1) {
-            let dist = hyperbolic_distance(p1, p2);
+            let dist = hyperboloid_hyperbolic_distance(p1, p2).unwrap();
             let is_edge = match beta {
                 Some(b) => {
                     let prob_inverse = (b / 2. * (dist - r)).exp() + 1.;
@@ -1025,27 +1026,6 @@ where
     Ok(graph)
 }
 
-#[inline]
-fn hyperbolic_distance(x: &[f64], y: &[f64]) -> f64 {
-    let mut sum_squared_x = 0.;
-    let mut sum_squared_y = 0.;
-    let mut inner_product = 0.;
-    for (x_i, y_i) in x.iter().zip(y.iter()) {
-        if x_i.is_infinite() || y_i.is_infinite() || x_i.is_nan() || y_i.is_nan() {
-            return f64::NAN;
-        }
-        sum_squared_x = x_i.mul_add(*x_i, sum_squared_x);
-        sum_squared_y = y_i.mul_add(*y_i, sum_squared_y);
-        inner_product = x_i.mul_add(*y_i, inner_product);
-    }
-    let arg = (1. + sum_squared_x).sqrt() * (1. + sum_squared_y).sqrt() - inner_product;
-    if arg < 1. {
-        0.
-    } else {
-        arg.acosh()
-    }
-}
-
 #[cfg(test)]
 mod tests {
     use crate::generators::InvalidInputError;
@@ -1056,10 +1036,7 @@ mod tests {
     };
     use crate::petgraph;
 
-    use super::hyperbolic_distance;
-
     // Test gnp_random_graph
-
     #[test]
     fn test_gnp_random_graph_directed() {
         let g: petgraph::graph::DiGraph<(), ()> =
@@ -1546,18 +1523,6 @@ mod tests {
     // x = sinh(r)cos(theta)
     // y = sinh(r)sin(theta)
 
-    #[test]
-    fn test_hyperbolic_dist() {
-        assert_eq!(
-            hyperbolic_distance(&[3_f64.sinh(), 0.], &[-0.5_f64.sinh(), 0.]),
-            3.5
-        );
-    }
-    #[test]
-    fn test_hyperbolic_dist_inf() {
-        assert!(hyperbolic_distance(&[f64::INFINITY, 0.], &[0., 0.]).is_nan());
-    }
-
     #[test]
     fn test_hyperbolic_random_graph_seeded() {
         let g = hyperbolic_random_graph::<petgraph::graph::UnGraph<(), ()>, _, _, _, _>(

rustworkx-core/src/geometry/distances.rs

@@ -0,0 +1,279 @@
+// Licensed under the Apache License, Version 2.0 (the "License"); you may
+// not use this file except in compliance with the License. You may obtain
+// a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+// License for the specific language governing permissions and limitations
+// under the License.
+
+use num_traits::Float;
+use std::iter::Sum;
+
+/// Error returned when two points have a different dimension.
+#[derive(Debug, PartialEq, Eq)]
+pub struct IncompatiblePointsError;
+
+/// Computes the L^`p` distance between `x` and `y`.
+///
+/// Works for any `p`>0. An [`IncompatiblePointsError`] is returned when `x` and `y` have different
+/// lengths.
+pub fn lp_distance<T>(x: &[T], y: &[T], p: i32) -> Result<T, IncompatiblePointsError>
+where
+    T: Float + Sum + From<i32>,
+{
+    if x.len() != y.len() {
+        Err(IncompatiblePointsError {})
+    } else {
+        Ok(x.iter()
+            .zip(y.iter())
+            .map(|(&a, &b)| T::powi(T::abs(b - a), p))
+            .sum::<T>()
+            .powf(T::one() / p.into()))
+    }
+}
+
+/// Computes the Euclidean distance between `x` and `y`.
+///
+/// An [`IncompatiblePointsError`] is returned when `x` and `y` have different lengths.
+pub fn euclidean_distance<T>(x: &[T], y: &[T]) -> Result<T, IncompatiblePointsError>
+where
+    T: Float + Sum + From<i32>,
+{
+    lp_distance(x, y, 2)
+}
+
+/// Computes the maximum distance (Chebyshev distance or L^infinity distance) between `x` and `y`.
+///
+/// An [`IncompatiblePointsError`] is returned when `x` and `y` have different lengths.
+pub fn maximum_distance<T: Float>(x: &[T], y: &[T]) -> Result<T, IncompatiblePointsError> {
+    if x.len() != y.len() {
+        Err(IncompatiblePointsError {})
+    } else {
+        Ok(x.iter()
+            .zip(y.iter())
+            .map(|(&a, &b)| T::abs(b - a))
+            .reduce(Float::max)
+            .unwrap_or(T::zero()))
+    }
+}
+
+/// Computes the Euclidean dot product between points the unit n-sphere, where n is the length of
+/// `angles1` and `angles2`.
+///
+/// No check is done on the lengths of `angles1` and `angles2`. The Euclidean dot product is also
+/// the cosine of the angular distance between `angles1` and `angles2`.
+fn euclidean_dot_product<T: Float>(angles1: &[T], angles2: &[T]) -> T {
+    let mut total = T::zero();
+    let mut sin_prod = T::one();
+    let d = angles1.len();
+    for (i, (t1, t2)) in angles1.iter().zip(angles2.iter()).enumerate() {
+        total = T::mul_add(sin_prod, t1.cos() * t2.cos(), total);
+        sin_prod = sin_prod * t1.sin() * t2.sin();
+        if i == d - 1 {
+            total = total + sin_prod;
+        }
+    }
+    total
+}
+
+/// Computes the distance between the points `angles1` and `angles2` on the unit n-sphere, where n
+/// is the length of `angles1` and `angles2`.
+///
+/// The last element of `angles1` and `angles2` is assumed to be in [0, 2pi] or [-pi, pi] (and the
+/// other elements are in [0, pi]). An [`IncompatiblePointsError`] is returned when `angles1` and
+/// `angles2` have different lengths.
+pub fn angular_distance<T: Float>(
+    angles1: &[T],
+    angles2: &[T],
+) -> Result<T, IncompatiblePointsError> {
+    if angles1.len() != angles2.len() {
+        Err(IncompatiblePointsError {})
+    } else {
+        Ok(euclidean_dot_product(angles1, angles2).acos().abs())
+    }
+}
+
+/// Computes the hyperbolic distance between two points in polar coordinates.
+///
+/// `r1` and `r2` are the distances to the origin and the last element of `angles1` and `angles2`
+/// is assumed to be in [0, 2pi] or [-pi, pi] (and the other elements are in [0, pi]). An
+/// [`IncompatiblePointsError`] is returned when `angles1` and `angles2` have different lengths.
+pub fn polar_hyperbolic_distance<T: Float>(
+    r1: T,
+    angles1: &[T],
+    r2: T,
+    angles2: &[T],
+) -> Result<T, IncompatiblePointsError> {
+    if angles1.len() != angles2.len() {
+        Err(IncompatiblePointsError {})
+    } else {
+        let arg = (r1 - r2).cosh()
+            + (T::one() - euclidean_dot_product(angles1, angles2)) * r1.sinh() * r2.sinh();
+        Ok(if arg < T::one() {
+            T::zero()
+        } else {
+            arg.acosh()
+        })
+    }
+}
+
+/// Computes the hyperbolic distance between the points `x1` and `x2` in the hyperboloid model.
+///
+/// The "time" coordinate (opposite sign in the metric) is inferred from the others and should not
+/// be included in `x1` and `x2`. An [`IncompatiblePointsError`] is returned when `x` and `y` have
+/// different lengths.
+pub fn hyperboloid_hyperbolic_distance<T: Float>(
+    x: &[T],
+    y: &[T],
+) -> Result<T, IncompatiblePointsError> {
+    if x.len() != y.len() {
+        Err(IncompatiblePointsError {})
+    } else {
+        let mut sum_x_squared = T::zero();
+        let mut sum_y_squared = T::zero();
+        let mut sum_xy = T::zero();
+        for (x_i, y_i) in x.iter().zip(y.iter()) {
+            sum_x_squared = T::mul_add(*x_i, *x_i, sum_x_squared);
+            sum_y_squared = T::mul_add(*y_i, *y_i, sum_y_squared);
+            sum_xy = T::mul_add(*x_i, *y_i, sum_xy);
+        }
+        let arg = (T::one() + sum_x_squared).sqrt() * (T::one() + sum_y_squared).sqrt() - sum_xy;
+        Ok(if arg < T::one() {
+            T::zero()
+        } else {
+            arg.acosh()
+        })
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use std::f64::consts::PI;
+
+    use super::{
+        angular_distance, euclidean_distance, hyperboloid_hyperbolic_distance, lp_distance,
+        maximum_distance, polar_hyperbolic_distance, IncompatiblePointsError,
+    };
+
+    #[test]
+    fn test_l4_dist() {
+        assert!(
+            (lp_distance(&[1., 2., 3.], &[5., 3., 1.], 4).unwrap()
+                - (256_f64 + 1. + 16.).sqrt().sqrt())
+            .abs()
+                < 1e-15
+        );
+    }
+
+    #[test]
+    fn test_l1_dist() {
+        assert!(
+            (lp_distance(&[1., 2., 3.], &[5., 3., 1.], 1).unwrap() - (4. + 1. + 2_f64)).abs()
+                < 1e-15
+        );
+    }
+
+    #[test]
+    fn test_lp_dist_incompatible_error() {
+        assert_eq!(
+            lp_distance(&[0., 0.], &[0., 0., 0.], 1),
+            Err(IncompatiblePointsError {})
+        );
+        assert_eq!(
+            lp_distance(&[0., 0.], &[0., 0., 0.], 4),
+            Err(IncompatiblePointsError {})
+        );
+    }
+
+    #[test]
+    fn test_euclidean_dist() {
+        assert!(
+            (euclidean_distance(&[1., 2., 3.], &[5., 3., 1.]).unwrap() - (16_f64 + 1. + 4.).sqrt())
+                .abs()
+                < 1e-15
+        );
+    }
+
+    #[test]
+    fn test_euclidean_dist_incompatible_error() {
+        assert_eq!(
+            euclidean_distance(&[0., 0.], &[0., 0., 0.]),
+            Err(IncompatiblePointsError {})
+        );
+    }
+
+    #[test]
+    fn test_maximum_dist() {
+        assert!((maximum_distance(&[1., 2., 3.], &[5., 3., 1.]).unwrap() - 4_f64).abs() < 1e-15);
+    }
+
+    #[test]
+    fn test_maximum_dist_incompatible_error() {
+        assert_eq!(
+            maximum_distance(&[0., 0.], &[0., 0., 0.]),
+            Err(IncompatiblePointsError {})
+        );
+    }
+
+    #[test]
+    fn test_angular_dist() {
+        assert!((angular_distance(&[0.3], &[0.5]).unwrap() - 0.2_f64).abs() < 1e-15);
+        assert!((angular_distance(&[0.5 * PI, PI], &[0.5 * PI, 0.]).unwrap() - PI).abs() < 1e-15);
+        assert!((angular_distance(&[0.2, 0., 1.], &[0., 0., 1.]).unwrap() - 0.2_f64).abs() < 1e-15);
+        assert!(
+            (angular_distance(&[2. * PI, 0., 1.], &[0., 0., 1.]).unwrap() - 0_f64).abs() < 1e-15
+        );
+    }
+
+    #[test]
+    fn test_angular_dist_incompatible_error() {
+        assert_eq!(
+            angular_distance(&[0.], &[0., 0.]),
+            Err(IncompatiblePointsError {})
+        );
+    }
+
+    #[test]
+    fn test_polar_hyperbolic_dist() {
+        assert_eq!(
+            polar_hyperbolic_distance(3., &[0.], 0.5, &[PI]).unwrap(),
+            3.5
+        );
+    }
+
+    #[test]
+    fn test_polar_hyperbolic_dist_incompatible_error() {
+        assert_eq!(
+            polar_hyperbolic_distance(1., &[0.], 1., &[0., 0.]),
+            Err(IncompatiblePointsError {})
+        );
+    }
+
+    #[test]
+    fn test_hyperboloid_dist() {
+        assert_eq!(
+            hyperboloid_hyperbolic_distance(&[3_f64.sinh(), 0.], &[-0.5_f64.sinh(), 0.]).unwrap(),
+            3.5
+        );
+    }
+
+    #[test]
+    fn test_hyperboloid_dist_inf() {
+        assert!(
+            hyperboloid_hyperbolic_distance(&[f64::INFINITY, 0.], &[0., 0.])
+                .unwrap()
+                .is_nan()
+        );
+    }
+    #[test]
+    fn test_hyperboloid_dist_length_error() {
+        assert_eq!(
+            euclidean_distance(&[0., 0.], &[0., 0., 0.]),
+            Err(IncompatiblePointsError {})
+        );
+    }
+}