netoptim-rs 0.1.3

//! Network optimization algorithms in Rust.

/// Dijkstra's shortest path algorithm implementation.
pub mod dijkstra;

/// Error types for network optimization.
pub mod error;

/// Negative cycle detection using Howard's algorithm.
pub mod neg_cycle;

/// Maximum parametric optimization.
pub mod parametric;

/// Oracle for parametric network problems (cutting-plane / feasibility).
pub mod network_oracle;

/// Oracle for optimal matrix scaling.
pub mod optscaling_oracle;

/// Minimum cost-to-time cycle ratio solver.
pub mod min_cycle_ratio;

/// Graph utility functions.
pub mod utils;

pub use error::NetOptimError;
pub use utils::*;

#[cfg(test)]
mod integration_tests;

// Logging module - only available with std feature
#[cfg(feature = "std")]
pub mod logging;

use petgraph::prelude::*;

use petgraph::algo::{FloatMeasure, NegativeCycle};
use petgraph::visit::{
    IntoEdges, IntoNodeIdentifiers, NodeCount, NodeIndexable, VisitMap, Visitable,
};

/// Result of shortest path algorithms.
///
/// Contains the distances from the source node to all other nodes,
/// and the predecessor of each node along the shortest path.
#[derive(Debug, Clone)]
pub struct Paths<NodeId, EdgeWeight> {
    pub distances: Vec<EdgeWeight>,
    pub predecessors: Vec<Option<NodeId>>,
}

/// \[Generic\] Compute shortest paths from node `source` to all other.
///
/// $$ \text{dist}^{(k+1)}\[v\] = \min\left(\text{dist}^{(k)}\[v\], \min_{(u,v) \in E} (\text{dist}^{(k)}\[u\] + w(u,v))\right) $$
///
/// Using the [Bellman–Ford algorithm][bf]; negative edge costs are
/// permitted, but the graph must not have a cycle of negative weights
/// (in that case it will return an error).
///
/// On success, return one vec with path costs, and another one which points
/// out the predecessor of a node along a shortest path. The vectors
/// are indexed by the graph's node indices.
///
/// [bf]: https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm
///
/// # Example
/// ```rust
/// use petgraph::Graph;
/// use petgraph::algo::bellman_ford;
/// use petgraph::prelude::*;
///
/// let mut g = Graph::new();
/// let a = g.add_node(()); // node with no weight
/// let b = g.add_node(());
/// let c = g.add_node(());
/// let d = g.add_node(());
/// let edge = g.add_node(());
/// let f = g.add_node(());
/// g.extend_with_edges(&[
///     (0, 1, 2.0),
///     (0, 3, 4.0),
///     (1, 2, 1.0),
///     (1, 5, 7.0),
///     (2, 4, 5.0),
///     (4, 5, 1.0),
///     (3, 4, 1.0),
/// ]);
///
/// // Graph represented with the weight of each edge
/// //
/// //     2       1
/// // a ----- b ----- c
/// // | 4     | 7     |
/// // d       f       | 5
/// // | 1     | 1     |
/// // \------ edge ------/
///
/// let path = bellman_ford(&g, a);
/// assert!(path.is_ok());
/// let path = path.unwrap();
/// assert_eq!(path.distances, vec![    0.0,     2.0,    3.0,    4.0,     5.0,     6.0]);
/// assert_eq!(path.predecessors, vec![None, Some(a),Some(b),Some(a), Some(d), Some(edge)]);
///
/// // Node f (indice 5) can be reach from a with a path costing 6.
/// // Predecessor of f is Some(edge) which predecessor is Some(d) which predecessor is Some(a).
/// // Thus the path from a to f is a <-> d <-> edge <-> f
///
/// let graph_with_neg_cycle = Graph::<(), f32, Undirected>::from_edges(&[
///         (0, 1, -2.0),
///         (0, 3, -4.0),
///         (1, 2, -1.0),
///         (1, 5, -25.0),
///         (2, 4, -5.0),
///         (4, 5, -25.0),
///         (3, 4, -1.0),
/// ]);
///
/// assert!(bellman_ford(&graph_with_neg_cycle, NodeIndex::new(0)).is_err());
/// ```
pub fn bellman_ford<G>(
    g: G,
    source: G::NodeId,
) -> Result<Paths<G::NodeId, G::EdgeWeight>, NegativeCycle>
where
    G: NodeCount + IntoNodeIdentifiers + IntoEdges + NodeIndexable,
    G::EdgeWeight: FloatMeasure,
{
    let ix = |i| g.to_index(i);

    // Step 1 and Step 2: initialize and relax
    let (distances, predecessors) = bellman_ford_initialize_relax(g, source);

    // Step 3: check for negative weight cycle
    for i in g.node_identifiers() {
        for edge in g.edges(i) {
            let j = edge.target();
            let w = *edge.weight();
            if distances[ix(i)] + w < distances[ix(j)] {
                return Err(NegativeCycle(()));
            }
        }
    }

    Ok(Paths {
        distances,
        predecessors,
    })
}

/// \[Generic\] Find the path of a negative cycle reachable from node `source`.
///
/// $$ \text{A cycle } C \text{ is negative iff } \sum_{(u,v) \in C} w(u,v) < 0 $$
///
/// Using the [find_negative_cycle][nc]; will search the Graph for negative cycles using
/// [Bellman–Ford algorithm][bf]. If no negative cycle is found the function will return `None`.
///
/// If a negative cycle is found from source, return one vec with a path of `NodeId`s.
///
/// The time complexity of this algorithm should be the same as the Bellman-Ford (O(|V|·|E|)).
///
/// [nc]: https://blogs.asarkar.com/assets/docs/algorithms-curated/Negative-Weight%20Cycle%20Algorithms%20-%20Huang.pdf
/// [bf]: https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm
///
/// # Example
/// ```rust
/// use petgraph::Graph;
/// use petgraph::algo::find_negative_cycle;
/// use petgraph::prelude::*;
///
/// let graph_with_neg_cycle = Graph::<(), f32, Directed>::from_edges(&[
///         (0, 1, 1.),
///         (0, 2, 1.),
///         (0, 3, 1.),
///         (1, 3, 1.),
///         (2, 1, 1.),
///         (3, 2, -3.),
/// ]);
///
/// let path = find_negative_cycle(&graph_with_neg_cycle, NodeIndex::new(0));
/// assert_eq!(
///     path,
///     Some([NodeIndex::new(1), NodeIndex::new(3), NodeIndex::new(2)].to_vec())
/// );
/// ```
/// # Example: Graph with no negative cycle
/// ```rust
/// use petgraph::Graph;
/// use petgraph::algo::find_negative_cycle;
/// use petgraph::prelude::*;
///
/// let graph = Graph::<(), f32, Directed>::from_edges(&[
///     (0, 1, 1.0),
///     (1, 2, 1.0),
///     (2, 3, 1.0),
/// ]);
/// let result = find_negative_cycle(&graph, NodeIndex::new(0));
/// assert!(result.is_none());
/// ```
pub fn find_negative_cycle<G>(g: G, source: G::NodeId) -> Option<Vec<G::NodeId>>
where
    G: NodeCount + IntoNodeIdentifiers + IntoEdges + NodeIndexable + Visitable,
    G::EdgeWeight: FloatMeasure,
{
    let ix = |i| g.to_index(i);
    let mut path = Vec::<G::NodeId>::new();

    // Step 1: initialize and relax
    let (distance, predecessor) = bellman_ford_initialize_relax(g, source);

    // Step 2: Check for negative weight cycle
    'outer: for i in g.node_identifiers() {
        for edge in g.edges(i) {
            let j = edge.target();
            let w = *edge.weight();
            if distance[ix(i)] + w < distance[ix(j)] {
                // Step 3: negative cycle found
                let start = j;
                let mut node = start;
                let mut visited = g.visit_map();
                // Go backward in the predecessor chain
                loop {
                    let ancestor = match predecessor[ix(node)] {
                        Some(predecessor_node) => predecessor_node,
                        None => node, // no predecessor, self cycle
                    };
                    // We have only 2 ways to find the cycle and break the loop:
                    // 1. start is reached
                    if ancestor == start {
                        path.push(ancestor);
                        break;
                    }
                    // 2. some node was reached twice
                    else if visited.is_visited(&ancestor) {
                        // Drop any node in path that is before the first ancestor
                        let pos = path
                            .iter()
                            .position(|&p| p == ancestor)
                            .expect("we should always have a position");
                        path = path[pos..path.len()].to_vec();

                        break;
                    }

                    // None of the above, some middle path node
                    path.push(ancestor);
                    visited.visit(ancestor);
                    node = ancestor;
                }
                // We are done here
                break 'outer;
            }
        }
    }
    if !path.is_empty() {
        // Users will probably need to follow the path of the negative cycle
        // so it should be in the reverse order than it was found by the algorithm.
        path.reverse();
        Some(path)
    } else {
        None
    }
}

// Perform Step 1 and Step 2 of the Bellman-Ford algorithm.
#[inline(always)]
fn bellman_ford_initialize_relax<G>(
    g: G,
    source: G::NodeId,
) -> (Vec<G::EdgeWeight>, Vec<Option<G::NodeId>>)
where
    G: NodeCount + IntoNodeIdentifiers + IntoEdges + NodeIndexable,
    G::EdgeWeight: FloatMeasure,
{
    // Step 1: initialize graph
    let mut predecessor = vec![None; g.node_bound()];
    let mut distance = vec![<_>::infinite(); g.node_bound()];
    let ix = |i| g.to_index(i);
    distance[ix(source)] = <_>::zero();

    // Step 2: relax edges repeatedly
    for _ in 1..g.node_count() {
        let mut did_update = false;
        for i in g.node_identifiers() {
            for edge in g.edges(i) {
                let j = edge.target();
                let w = *edge.weight();
                if distance[ix(i)] + w < distance[ix(j)] {
                    distance[ix(j)] = distance[ix(i)] + w;
                    predecessor[ix(j)] = Some(i);
                    did_update = true;
                }
            }
        }
        if !did_update {
            break;
        }
    }
    (distance, predecessor)
}

#[cfg(test)]
mod tests {
    use super::*;
    use petgraph::Graph;

    #[test]
    fn test_bellman_ford_negative_cycle() {
        let graph_with_neg_cycle =
            Graph::<(), f32, Directed>::from_edges([(0, 1, 1.0), (1, 2, 1.0), (2, 0, -3.0)]);
        let result = bellman_ford(&graph_with_neg_cycle, NodeIndex::new(0));
        assert!(result.is_err());
    }

    #[test]
    fn test_bellman_ford_no_edges() {
        let mut graph = Graph::<(), f32, Directed>::new();
        let n0 = graph.add_node(());
        let result = bellman_ford(&graph, n0);
        assert!(result.is_ok());
        let paths = result.unwrap();
        assert_eq!(paths.distances, vec![0.0]);
        assert_eq!(paths.predecessors, vec![None]);
    }

    #[test]
    fn test_bellman_ford_disconnected_components() {
        let mut graph = Graph::<(), f32, Directed>::new();
        let n0 = graph.add_node(());
        let n1 = graph.add_node(());
        let n2 = graph.add_node(());
        graph.add_edge(n0, n1, 1.0);

        let result = bellman_ford(&graph, n0);
        assert!(result.is_ok());
        let paths = result.unwrap();
        // Node 2 is unreachable, so its distance should be infinite
        assert_eq!(paths.distances.len(), 3);
        assert_eq!(paths.distances[n0.index()], 0.0);
        assert_eq!(paths.distances[n1.index()], 1.0);
        assert!(paths.distances[n2.index()].is_infinite());
        assert_eq!(paths.predecessors, vec![None, Some(n0), None]);
    }

    #[test]
    fn test_find_negative_cycle_exists() {
        let graph_with_neg_cycle =
            Graph::<(), f32, Directed>::from_edges([(0, 1, 1.0), (1, 2, 1.0), (2, 0, -3.0)]);
        let result = find_negative_cycle(&graph_with_neg_cycle, NodeIndex::new(0));
        assert!(result.is_some());
        let cycle = result.unwrap();
        assert_eq!(cycle.len(), 3);
        assert!(cycle.contains(&NodeIndex::new(0)));
        assert!(cycle.contains(&NodeIndex::new(1)));
        assert!(cycle.contains(&NodeIndex::new(2)));
    }

    #[test]
    fn test_find_negative_cycle_none() {
        let graph = Graph::<(), f32, Directed>::from_edges([(0, 1, 1.0), (1, 2, 1.0), (2, 3, 1.0)]);
        let result = find_negative_cycle(&graph, NodeIndex::new(0));
        assert!(result.is_none());
    }

    #[test]
    fn test_find_negative_cycle_unreachable() {
        let graph =
            Graph::<(), f32, Directed>::from_edges([(0, 1, 1.0), (2, 3, -1.0), (3, 2, -1.0)]);
        let result = find_negative_cycle(&graph, NodeIndex::new(0));
        assert!(result.is_none());
    }
    use crate::neg_cycle::NegCycleFinder;
    use crate::parametric::{MaxParametricSolver, ParametricAPI};
    use num::rational::Ratio;
    use petgraph::graph::{DiGraph, EdgeReference};

    #[test]
    fn test_neg_cycle_multiple_neg_cycles() {
        let digraph = DiGraph::<(), Ratio<i32>>::from_edges([
            (0, 1, Ratio::new(1, 1)),
            (1, 0, Ratio::new(-2, 1)), // Cycle 1: 0 -> 1 -> 0 (weight -1)
            (2, 3, Ratio::new(1, 1)),
            (3, 2, Ratio::new(-3, 1)), // Cycle 2: 2 -> 3 -> 2 (weight -2)
        ]);

        let mut ncf = NegCycleFinder::new(&digraph);
        let mut dist = [
            Ratio::new(0, 1),
            Ratio::new(0, 1),
            Ratio::new(0, 1),
            Ratio::new(0, 1),
        ];
        let result = ncf.howard(&mut dist, |e| *e.weight());
        assert!(result.is_some());
        let cycle = result.unwrap();
        let cycle_weight: Ratio<i32> = cycle.iter().map(|e| *e.weight()).sum();
        assert!(cycle_weight < Ratio::new(0, 1));
    }

    #[test]
    fn test_neg_cycle_not_reachable() {
        let digraph = DiGraph::<(), Ratio<i32>>::from_edges([
            (0, 1, Ratio::new(1, 1)),
            (2, 3, Ratio::new(-1, 1)),
            (3, 2, Ratio::new(-1, 1)),
        ]);

        let mut ncf = NegCycleFinder::new(&digraph);
        let mut dist = [
            Ratio::new(0, 1),
            Ratio::new(0, 1),
            Ratio::new(0, 1),
            Ratio::new(0, 1),
        ];
        let result = ncf.howard(&mut dist, |e| *e.weight());
        assert!(result.is_some());
    }

    struct TestParametricAPI;

    impl ParametricAPI<(), Ratio<i32>> for TestParametricAPI {
        fn distance(&self, ratio: &Ratio<i32>, edge: &EdgeReference<Ratio<i32>>) -> Ratio<i32> {
            *edge.weight() - *ratio
        }

        fn zero_cancel(&self, cycle: &[EdgeReference<Ratio<i32>>]) -> Ratio<i32> {
            let mut sum_a = Ratio::new(0, 1);
            let mut sum_b = Ratio::new(0, 1);
            for edge in cycle {
                sum_a += *edge.weight();
                sum_b += Ratio::new(1, 1);
            }
            sum_a / sum_b
        }
    }

    #[test]
    fn test_max_parametric_solver_no_neg_cycle() {
        let digraph = DiGraph::<(), Ratio<i32>>::from_edges([
            (0, 1, Ratio::new(1, 1)),
            (1, 2, Ratio::new(1, 1)),
            (2, 0, Ratio::new(1, 1)),
        ]);

        let mut solver = MaxParametricSolver::new(&digraph, TestParametricAPI);
        let mut dist = [Ratio::new(0, 1), Ratio::new(0, 1), Ratio::new(0, 1)];
        let mut ratio = Ratio::new(0, 1);

        let cycle = solver.run(&mut dist, &mut ratio);
        assert!(cycle.is_empty());
    }

    #[test]
    fn test_max_parametric_solver_multiple_neg_cycles() {
        let digraph = DiGraph::<(), Ratio<i32>>::from_edges([
            (0, 1, Ratio::new(1, 1)),
            (1, 0, Ratio::new(-2, 1)), // Cycle 1: ratio -1/1
            (2, 3, Ratio::new(1, 1)),
            (3, 2, Ratio::new(-4, 1)), // Cycle 2: ratio -3/1
        ]);

        let mut solver = MaxParametricSolver::new(&digraph, TestParametricAPI);
        let mut dist = [
            Ratio::new(0, 1),
            Ratio::new(0, 1),
            Ratio::new(0, 1),
            Ratio::new(0, 1),
        ];
        let mut ratio = Ratio::new(0, 1);

        let cycle = solver.run(&mut dist, &mut ratio);
        assert!(!cycle.is_empty());
        assert_eq!(ratio, Ratio::new(-3, 2));
    }

    #[test]
    fn test_bellman_ford_neg_cycle() {
        let graph_with_neg_cycle =
            Graph::<(), f32, Directed>::from_edges([(0, 1, 1.0), (1, 2, 1.0), (2, 0, -3.0)]);
        let result = bellman_ford(&graph_with_neg_cycle, NodeIndex::new(0));
        assert!(result.is_err());
    }

    #[test]
    fn test_bellman_ford_no_edge() {
        let mut graph = Graph::<(), f32, Directed>::new();
        let n0 = graph.add_node(());
        let result = bellman_ford(&graph, n0);
        assert!(result.is_ok());
        let paths = result.unwrap();
        assert_eq!(paths.distances, vec![0.0]);
        assert_eq!(paths.predecessors, vec![None]);
    }

    #[test]
    fn test_bellman_ford_disconnected() {
        let mut graph = Graph::<(), f32, Directed>::new();
        let n0 = graph.add_node(());
        let n1 = graph.add_node(());
        let n2 = graph.add_node(());
        graph.add_edge(n0, n1, 1.0);

        let result = bellman_ford(&graph, n0);
        assert!(result.is_ok());
        let paths = result.unwrap();
        // Node 2 is unreachable, so its distance should be infinite
        assert_eq!(paths.distances.len(), 3);
        assert_eq!(paths.distances[n0.index()], 0.0);
        assert_eq!(paths.distances[n1.index()], 1.0);
        assert!(paths.distances[n2.index()].is_infinite());
        assert_eq!(paths.predecessors, vec![None, Some(n0), None]);
    }

    #[test]
    fn test_find_negative_cycle_multiple() {
        let graph_with_neg_cycle = Graph::<(), f32, Directed>::from_edges([
            (0, 1, 1.0),
            (1, 0, -2.0),
            (2, 3, 1.0),
            (3, 2, -3.0),
        ]);
        let result = find_negative_cycle(&graph_with_neg_cycle, NodeIndex::new(0));
        assert!(result.is_some());
    }

    #[test]
    fn test_find_negative_cycle_no_neg_cycle() {
        let graph = Graph::<(), f32, Directed>::from_edges([(0, 1, 1.0), (1, 2, 1.0), (2, 3, 1.0)]);
        let result = find_negative_cycle(&graph, NodeIndex::new(0));
        assert!(result.is_none());
    }

    #[test]
    fn test_find_negative_cycle_unreachable_neg_cycle() {
        let graph =
            Graph::<(), f32, Directed>::from_edges([(0, 1, 1.0), (2, 3, -1.0), (3, 2, -1.0)]);
        let result = find_negative_cycle(&graph, NodeIndex::new(0));
        assert!(result.is_none());
    }
}