hypergraphx 0.0.5

//! Bellman-Ford algorithms.

use super::measures::BoundedMeasure;
use crate::prelude::*;
use core::marker::Copy;

#[derive(Debug, Clone)]
pub struct Paths<K: BoundedMeasure> {
    pub distances: Vec<K>,
    pub predecessors: Vec<Option<usize>>,
}

/// \[Generic\] Compute shortest paths from node `source` to all other.
///
/// 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
///
pub fn bellman_ford<'a, G: ?Sized, E, F, K, T: GraphType>(
    g: &'a G,
    source: usize,
    edge_cost: F,
) -> Result<Paths<K>, HypergraphErrors>
where
    G: CommonProperties<'a, T, NodeIndex = usize, EdgeIndex = usize>,
    F: Fn(usize) -> K,
    K: BoundedMeasure + Copy,
{
    // Step 1 and Step 2: initialize and relax
    let (distances, predecessors) = bellman_ford_initialize_relax(g, source, &edge_cost);

    // Step 3: check for negative weight cycle
    for i in 0..g.node_count() {
        for (k, nodes) in g.costed_neighbours(i, &edge_cost) {
            for j in nodes {
                if distances[i] + k < distances[j] {
                    return Err(HypergraphErrors::NegativeCycle);
                }
            }
        }
    }

    Ok(Paths {
        distances,
        predecessors,
    })
}

/// \[Generic\] Find the path of a negative cycle reachable from node `source`.
///
/// 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 hypergraphx::algo::bellman_ford::find_negative_cycle;
/// use hypergraphx::prelude::*;
/// hypergraph! {
///     let mut graph: DirectedHypergraph<_, i32> {
///         nodes: {
///             let v0 = ();
///             let v1 = ();
///             let v2 = ();
///             let v3 = ();
///         };
///         edges: {
///             ([v0] -> [v1]) = 1;
///             ([v0] -> [v2]) = 1;
///             ([v0] -> [v3]) = 1;
///             ([v1] -> [v3]) = 1;
///             ([v2] -> [v1]) = 1;
///             ([v3] -> [v2]) = -3;
///         };
///     }
/// }
/// let path = find_negative_cycle(&graph, 0, &|e| graph.edge_weight_copied_unchecked(e));
/// assert_eq!(
///     path,
///     Some([1, 3, 2].to_vec())
/// );
/// ```
pub fn find_negative_cycle<'a, G, F, T: GraphType, K: BoundedMeasure + Copy>(
    g: &'a G,
    source: usize,
    edge_cost: &F,
) -> Option<Vec<usize>>
where
    G: CommonProperties<'a, T, NodeIndex = usize, EdgeIndex = usize>,
    // If you want contents of the edge, capture them in the closure
    // In general, unless you have a reference to the graph, you can only deal in indices, and not refs.
    F: Fn(usize) -> K,
{
    // let ix = |i| g.to_index(i);
    let mut path = Vec::<usize>::new();

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

    // Step 2: Check for negative weight cycle
    'outer: for i in 0..g.node_count() {
        // for edge in g.incident_edges(i).unwrap() {
        for (w, nodes) in g.costed_neighbours(i, edge_cost) {
            // let nb = g.contained_nodes(edge).unwrap();
            // let w = edge_cost(edge, g.edge(edge)?);
            for j in nodes {
                if i != j && distance[i] + w < distance[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[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.contains(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.insert(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
    }
}

/// \[Generic\] Find the path of a negative cycle reachable from node `source`.
///
/// 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())
/// );
/// ```
pub fn find_negative_directed_cycle<'a, G, F>(
    g: &'a G,
    source: usize,
    edge_cost: &F,
) -> Option<Vec<usize>>
where
    G: DiGraphProperties<'a, NodeIndex = usize, EdgeIndex = usize>,
    F: Fn(usize, <G as GraphBasics<'_>>::EdgeRef) -> usize,
{
    // let ix = |i| g.to_index(i);
    let mut path = Vec::<usize>::new();

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

    // Step 2: Check for negative weight cycle
    'outer: for i in 0..g.node_count() {
        for edge in g.out_edges(i).unwrap() {
            let nb = g.target(edge).unwrap();
            let w = edge_cost(edge, g.edge(edge)?);
            for j in nb {
                if i != j && distance[i] + w < distance[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[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.contains(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.insert(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<'a, G: ?Sized, F, K, T: GraphType>(
    g: &'a G,
    source: usize,
    edge_cost: F,
) -> (Vec<K>, Vec<Option<usize>>)
where
    G: CommonProperties<'a, T, NodeIndex = usize, EdgeIndex = usize>,
    F: Fn(usize) -> K,
    K: BoundedMeasure + Copy,
{
    // Step 1: initialize graph
    let mut predecessor = vec![None; g.node_count()];
    let mut distance = vec![K::max(); g.node_count()];
    // let ix = |i| g.to_index(i);
    distance[source] = K::default();

    // Step 2: relax edges repeatedly
    for _ in 1..g.node_count() {
        let mut did_update = false;
        for i in 0..g.node_count() {
            // for edge_index in g.incident_edges(i).unwrap() {
            //     // let j = edge_ref.target();
            //     let w = edge_cost(edge_index, g.edge(edge_index).unwrap());
            //     for j in g.contained_nodes(edge_index).unwrap() {
            for (w, nodes) in g.costed_neighbours(i, &edge_cost) {
                for j in nodes {
                    if distance[i] + w < distance[j] {
                        distance[j] = distance[i] + w;
                        predecessor[j] = Some(i);
                        did_update = true;
                    }
                }
            }
        }
        if !did_update {
            break;
        }
    }
    (distance, predecessor)
}

// Perform Step 1 and Step 2 of the Bellman-Ford algorithm.
#[inline(always)]
fn bellman_ford_initialize_relax_directed<'a, G, F, K>(
    g: &'a G,
    source: usize,
    mut edge_cost: F,
) -> (Vec<K>, Vec<Option<usize>>)
where
    G: DiGraphProperties<'a, NodeIndex = usize, EdgeIndex = usize>,
    F: FnMut(usize, <G as GraphBasics<'_>>::EdgeRef) -> K,
    K: BoundedMeasure + Copy,
{
    // Step 1: initialize graph
    let mut predecessor = vec![None; g.node_count()];
    let mut distance = vec![K::max(); g.node_count()];
    // let ix = |i| g.to_index(i);
    distance[source] = K::default();

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