retworkx_core/

max_weight_matching.rs

// 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.

// Needed to pass shared state between functions
// closures don't work because of recurssion
#![allow(clippy::too_many_arguments)]
// Allow single character names to match naming convention from
// paper
#![allow(clippy::many_single_char_names)]

use std::cmp::max;
use std::mem;

use hashbrown::{HashMap, HashSet};

use petgraph::graph::NodeIndex;
use petgraph::visit::{
    EdgeCount, EdgeRef, GraphBase, GraphProp, IntoEdges, IntoNodeIdentifiers,
    NodeCount,
};
use petgraph::Undirected;

/// Return 2 * slack of edge k (does not work inside blossoms).
fn slack(
    edge_index: usize,
    dual_var: &[i128],
    edges: &[(usize, usize, i128)],
) -> i128 {
    let (source_index, target_index, weight) = edges[edge_index];
    dual_var[source_index] + dual_var[target_index] - 2 * weight
}

/// Generate the leaf vertices of a blossom.
fn blossom_leaves<E>(
    blossom: usize,
    num_nodes: usize,
    blossom_children: &[Vec<usize>],
) -> Result<Vec<usize>, E> {
    let mut out_vec: Vec<usize> = Vec::new();
    if blossom < num_nodes {
        out_vec.push(blossom);
    } else {
        let child_blossom = &blossom_children[blossom];
        for c in child_blossom {
            let child = *c;
            if child < num_nodes {
                out_vec.push(child);
            } else {
                for v in blossom_leaves(child, num_nodes, blossom_children)? {
                    out_vec.push(v);
                }
            }
        }
    }
    Ok(out_vec)
}

/// Assign label t to the top-level blossom containing vertex w
/// and record the fact that w was reached through the edge with
/// remote endpoint p.
fn assign_label<E>(
    w: usize,
    t: usize,
    p: Option<usize>,
    num_nodes: usize,
    in_blossoms: &[usize],
    labels: &mut Vec<Option<usize>>,
    label_ends: &mut Vec<Option<usize>>,
    best_edge: &mut Vec<Option<usize>>,
    queue: &mut Vec<usize>,
    blossom_children: &[Vec<usize>],
    blossom_base: &[Option<usize>],
    endpoints: &[usize],
    mate: &HashMap<usize, usize>,
) -> Result<(), E> {
    let b = in_blossoms[w];
    assert!(labels[w] == Some(0) && labels[b] == Some(0));
    labels[w] = Some(t);
    labels[b] = Some(t);
    label_ends[b] = p;
    label_ends[w] = p;
    best_edge[w] = None;
    best_edge[b] = None;
    if t == 1 {
        // b became an S-vertex/blossom; add it(s verticies) to the queue
        queue.append(&mut blossom_leaves(b, num_nodes, blossom_children)?);
    } else if t == 2 {
        // b became a T-vertex/blossom; assign label S to its mate.
        // (If b is a non-trivial blossom, its base is the only vertex
        // with an external mate.)
        let blossom_index: usize = b;
        let base: usize = blossom_base[blossom_index].unwrap();
        assert!(mate.get(&base).is_some());
        assign_label(
            endpoints[mate[&base]],
            1,
            mate.get(&base).map(|p| (p ^ 1)),
            num_nodes,
            in_blossoms,
            labels,
            label_ends,
            best_edge,
            queue,
            blossom_children,
            blossom_base,
            endpoints,
            mate,
        )?;
    }
    Ok(())
}

/// Trace back from vertices v and w to discover either a new blossom
/// or an augmenting path. Return the base vertex of the new blossom or None.
fn scan_blossom(
    node_a: usize,
    node_b: usize,
    in_blossoms: &[usize],
    blossom_base: &[Option<usize>],
    endpoints: &[usize],
    labels: &mut Vec<Option<usize>>,
    label_ends: &[Option<usize>],
    mate: &HashMap<usize, usize>,
) -> Option<usize> {
    let mut v: Option<usize> = Some(node_a);
    let mut w: Option<usize> = Some(node_b);
    // Trace back from v and w, placing breadcrumbs as we go
    let mut path: Vec<usize> = Vec::new();
    let mut base: Option<usize> = None;
    while v.is_some() || w.is_some() {
        // Look for a breadcrumb in v's blossom or put a new breadcrump
        let mut blossom = in_blossoms[v.unwrap()];
        if labels[blossom].is_none() || labels[blossom].unwrap() & 4 > 0 {
            base = blossom_base[blossom];
            break;
        }
        assert!(labels[blossom] == Some(1));
        path.push(blossom);
        labels[blossom] = Some(5);
        // Trace one step bacl.
        assert!(
            label_ends[blossom]
                == mate.get(&blossom_base[blossom].unwrap()).copied()
        );
        if label_ends[blossom].is_none() {
            // The base of blossom is single; stop tracing this path
            v = None;
        } else {
            let tmp = endpoints[label_ends[blossom].unwrap()];
            blossom = in_blossoms[tmp];
            assert!(labels[blossom] == Some(2));
            // blossom is a T-blossom; trace one more step back.
            assert!(label_ends[blossom].is_some());
            v = Some(endpoints[label_ends[blossom].unwrap()]);
        }
        // Swap v and w so that we alternate between both paths.
        if w.is_some() {
            mem::swap(&mut v, &mut w);
        }
    }
    // Remvoe breadcrumbs.
    for blossom in path {
        labels[blossom] = Some(1);
    }
    // Return base vertex, if we found one.
    base
}

/// Construct a new blossom with given base, containing edge k which
/// connects a pair of S vertices. Label the new blossom as S; set its dual
/// variable to zero; relabel its T-vertices to S and add them to the queue.
fn add_blossom<E>(
    base: usize,
    edge: usize,
    blossom_children: &mut Vec<Vec<usize>>,
    num_nodes: usize,
    edges: &[(usize, usize, i128)],
    in_blossoms: &mut Vec<usize>,
    dual_var: &mut Vec<i128>,
    labels: &mut Vec<Option<usize>>,
    label_ends: &mut Vec<Option<usize>>,
    best_edge: &mut Vec<Option<usize>>,
    queue: &mut Vec<usize>,
    blossom_base: &mut Vec<Option<usize>>,
    endpoints: &[usize],
    blossom_endpoints: &mut Vec<Vec<usize>>,
    unused_blossoms: &mut Vec<usize>,
    blossom_best_edges: &mut Vec<Vec<usize>>,
    blossom_parents: &mut Vec<Option<usize>>,
    neighbor_endpoints: &[Vec<usize>],
    mate: &HashMap<usize, usize>,
) -> Result<(), E> {
    let (mut v, mut w, _weight) = edges[edge];
    let blossom_b = in_blossoms[base];
    let mut blossom_v = in_blossoms[v];
    let mut blossom_w = in_blossoms[w];
    // Create blossom
    let blossom = unused_blossoms.pop().unwrap();
    blossom_base[blossom] = Some(base);
    blossom_parents[blossom] = None;
    blossom_parents[blossom_b] = Some(blossom);
    // Make list of sub-blossoms and their interconnecting edge endpoints.
    blossom_children[blossom].clear();
    blossom_endpoints[blossom].clear();
    // Trace back from blossom_v to base.
    while blossom_v != blossom_b {
        // Add blossom_v to the new blossom
        blossom_parents[blossom_v] = Some(blossom);
        blossom_children[blossom].push(blossom_v);
        let blossom_v_endpoint_label = label_ends[blossom_v].unwrap();
        blossom_endpoints[blossom].push(blossom_v_endpoint_label);
        assert!(
            labels[blossom_v] == Some(2)
                || (labels[blossom_v] == Some(1)
                    && label_ends[blossom_v]
                        == mate
                            .get(&blossom_base[blossom_v].unwrap())
                            .copied())
        );
        // Trace one step back.
        assert!(label_ends[blossom_v].is_some());
        v = endpoints[blossom_v_endpoint_label];
        blossom_v = in_blossoms[v];
    }
    // Reverse lists, add endpoint that connects the pair of S vertices.
    blossom_children[blossom].push(blossom_b);
    blossom_children[blossom].reverse();
    blossom_endpoints[blossom].reverse();
    blossom_endpoints[blossom].push(2 * edge);
    // Trace back from w to base.
    while blossom_w != blossom_b {
        // Add blossom_w to the new blossom
        blossom_parents[blossom_w] = Some(blossom);
        blossom_children[blossom].push(blossom_w);
        let blossom_w_endpoint_label = label_ends[blossom_w].unwrap();
        blossom_endpoints[blossom].push(blossom_w_endpoint_label ^ 1);
        assert!(
            labels[blossom_w] == Some(2)
                || (labels[blossom_w] == Some(1)
                    && label_ends[blossom_w]
                        == mate
                            .get(&blossom_base[blossom_w].unwrap())
                            .copied())
        );
        // Trace one step back
        assert!(label_ends[blossom_w].is_some());
        w = endpoints[blossom_w_endpoint_label];
        blossom_w = in_blossoms[w];
    }
    // Set label to S.
    assert!(labels[blossom_b] == Some(1));
    labels[blossom] = Some(1);
    label_ends[blossom] = label_ends[blossom_b];
    // Set dual variable to 0
    dual_var[blossom] = 0;
    // Relabel vertices
    for node in blossom_leaves(blossom, num_nodes, blossom_children)? {
        if labels[in_blossoms[node]] == Some(2) {
            // This T-vertex now turns into an S-vertex because it becomes
            // part of an S-blossom; add it to the queue
            queue.push(node);
        }
        in_blossoms[node] = blossom;
    }
    // Compute blossom_best_edges[blossom]
    let mut best_edge_to: HashMap<usize, usize> =
        HashMap::with_capacity(2 * num_nodes);
    for bv_ref in &blossom_children[blossom] {
        let bv = *bv_ref;
        // This sub-blossom does not have a list of least-slack edges;
        // get the information from the vertices.
        let nblists: Vec<Vec<usize>> = if blossom_best_edges[bv].is_empty() {
            let mut tmp: Vec<Vec<usize>> = Vec::new();
            for node in blossom_leaves(bv, num_nodes, blossom_children)? {
                tmp.push(
                    neighbor_endpoints[node].iter().map(|p| p / 2).collect(),
                );
            }
            tmp
        } else {
            // Walk this sub-blossom's least-slack edges.
            vec![blossom_best_edges[bv].clone()]
        };
        for nblist in nblists {
            for edge_index in nblist {
                let (mut i, mut j, _edge_weight) = edges[edge_index];
                if in_blossoms[j] == blossom {
                    mem::swap(&mut i, &mut j);
                }
                let blossom_j = in_blossoms[j];
                if blossom_j != blossom
                    && labels[blossom_j] == Some(1)
                    && (best_edge_to.get(&blossom_j).is_none()
                        || slack(edge_index, dual_var, edges)
                            < slack(best_edge_to[&blossom_j], dual_var, edges))
                {
                    best_edge_to.insert(blossom_j, edge_index);
                }
            }
        }
        // Forget about least-slack edges of the sub-blossom.
        blossom_best_edges[bv].clear();
        best_edge[bv] = None;
    }
    blossom_best_edges[blossom] = best_edge_to.values().copied().collect();
    //select best_edge[blossom]
    best_edge[blossom] = None;
    for edge_index in &blossom_best_edges[blossom] {
        if best_edge[blossom].is_none()
            || slack(*edge_index, dual_var, edges)
                < slack(best_edge[blossom].unwrap(), dual_var, edges)
        {
            best_edge[blossom] = Some(*edge_index);
        }
    }
    Ok(())
}

/// Expand the given top level blossom
fn expand_blossom<E>(
    blossom: usize,
    end_stage: bool,
    num_nodes: usize,
    blossom_children: &mut Vec<Vec<usize>>,
    blossom_parents: &mut Vec<Option<usize>>,
    in_blossoms: &mut Vec<usize>,
    dual_var: &[i128],
    labels: &mut Vec<Option<usize>>,
    label_ends: &mut Vec<Option<usize>>,
    best_edge: &mut Vec<Option<usize>>,
    queue: &mut Vec<usize>,
    blossom_base: &mut Vec<Option<usize>>,
    endpoints: &[usize],
    mate: &HashMap<usize, usize>,
    blossom_endpoints: &mut Vec<Vec<usize>>,
    allowed_edge: &mut Vec<bool>,
    unused_blossoms: &mut Vec<usize>,
) -> Result<(), E> {
    // Convert sub-blossoms into top-level blossoms.
    for s in blossom_children[blossom].clone() {
        blossom_parents[s] = None;
        if s < num_nodes {
            in_blossoms[s] = s
        } else if end_stage && dual_var[s] == 0 {
            // Recursively expand this sub-blossom
            expand_blossom(
                s,
                end_stage,
                num_nodes,
                blossom_children,
                blossom_parents,
                in_blossoms,
                dual_var,
                labels,
                label_ends,
                best_edge,
                queue,
                blossom_base,
                endpoints,
                mate,
                blossom_endpoints,
                allowed_edge,
                unused_blossoms,
            )?;
        } else {
            for v in blossom_leaves(s, num_nodes, blossom_children)? {
                in_blossoms[v] = s;
            }
        }
    }
    // if we expand a T-blossom during a stage, its a sub-blossoms must be
    // relabeled
    if !end_stage && labels[blossom] == Some(2) {
        // start at the sub-blossom through which the expanding blossom
        // obtained its label, and relabel sub-blossoms until we reach the
        // base.
        assert!(label_ends[blossom].is_some());
        let entry_child =
            in_blossoms[endpoints[label_ends[blossom].unwrap() ^ 1]];
        // Decied in which direction we will go around the blossom.
        let i = blossom_children[blossom]
            .iter()
            .position(|x| *x == entry_child)
            .unwrap();
        let mut j = i as i128;
        let j_step: i128;
        let endpoint_trick: usize = if i & 1 != 0 {
            // Start index is odd; go forward and wrap
            j -= blossom_children[blossom].len() as i128;
            j_step = 1;
            0
        } else {
            // start index is even; go backward
            j_step = -1;
            1
        };
        // Move along the blossom until we get to the base
        let mut p = label_ends[blossom].unwrap();
        while j != 0 {
            // Relabel the T-sub-blossom.
            labels[endpoints[p ^ 1]] = Some(0);
            if j < 0 {
                let length = blossom_endpoints[blossom].len();
                let index = length - j.abs() as usize;
                labels[endpoints[blossom_endpoints[blossom]
                    [index - endpoint_trick]
                    ^ endpoint_trick
                    ^ 1]] = Some(0);
            } else {
                labels[endpoints[blossom_endpoints[blossom]
                    [j as usize - endpoint_trick]
                    ^ endpoint_trick
                    ^ 1]] = Some(0);
            }
            assign_label(
                endpoints[p ^ 1],
                2,
                Some(p),
                num_nodes,
                in_blossoms,
                labels,
                label_ends,
                best_edge,
                queue,
                blossom_children,
                blossom_base,
                endpoints,
                mate,
            )?;
            // Step to the next S-sub-blossom and note it's forwward endpoint.
            let endpoint_index = if j < 0 {
                let tmp = j - endpoint_trick as i128;
                let length = blossom_endpoints[blossom].len();
                let index = length - tmp.abs() as usize;
                blossom_endpoints[blossom][index]
            } else {
                blossom_endpoints[blossom][j as usize - endpoint_trick]
            };
            allowed_edge[endpoint_index / 2] = true;
            j += j_step;
            p = if j < 0 {
                let tmp = j - endpoint_trick as i128;
                let length = blossom_endpoints[blossom].len();
                let index = length - tmp.abs() as usize;
                blossom_endpoints[blossom][index] ^ endpoint_trick
            } else {
                blossom_endpoints[blossom][j as usize - endpoint_trick]
                    ^ endpoint_trick
            };
            // Step to the next T-sub-blossom.
            allowed_edge[p / 2] = true;
            j += j_step;
        }
        // Relabel the base T-sub-blossom WITHOUT stepping through to
        // its mate (so don't call assign_label())
        let blossom_v = if j < 0 {
            let length = blossom_children[blossom].len();
            let index = length - j.abs() as usize;
            blossom_children[blossom][index]
        } else {
            blossom_children[blossom][j as usize]
        };
        labels[endpoints[p ^ 1]] = Some(2);
        labels[blossom_v] = Some(2);
        label_ends[endpoints[p ^ 1]] = Some(p);
        label_ends[blossom_v] = Some(p);
        best_edge[blossom_v] = None;
        // Continue along the blossom until we get back to entry_child
        j += j_step;
        let mut j_index = if j < 0 {
            let length = blossom_children[blossom].len();
            length - j.abs() as usize
        } else {
            j as usize
        };
        while blossom_children[blossom][j_index] != entry_child {
            // Examine the vertices of the sub-blossom to see whether
            // it is reachable from a neighboring S-vertex outside the
            // expanding blososm.
            let bv = blossom_children[blossom][j_index];
            if labels[bv] == Some(1) {
                // This sub-blossom just got label S through one of its
                // neighbors; leave it
                j += j_step;
                if j < 0 {
                    let length = blossom_children[blossom].len();
                    j_index = length - j.abs() as usize;
                } else {
                    j_index = j as usize;
                }
                continue;
            }
            let mut v: usize = 0;
            for tmp in blossom_leaves(bv, num_nodes, blossom_children)? {
                v = tmp;
                if labels[v].unwrap() != 0 {
                    break;
                }
            }
            // If the sub-blossom contains a reachable vertex, assign label T
            // to the sub-blossom
            if labels[v] != Some(0) {
                assert!(labels[v] == Some(2));
                assert!(in_blossoms[v] == bv);
                labels[v] = Some(0);
                labels[endpoints[mate[&blossom_base[bv].unwrap()]]] = Some(0);
                assign_label(
                    v,
                    2,
                    label_ends[v],
                    num_nodes,
                    in_blossoms,
                    labels,
                    label_ends,
                    best_edge,
                    queue,
                    blossom_children,
                    blossom_base,
                    endpoints,
                    mate,
                )?;
            }
            j += j_step;
            if j < 0 {
                let length = blossom_children[blossom].len();
                j_index = length - j.abs() as usize;
            } else {
                j_index = j as usize;
            }
        }
    }
    // Recycle the blossom number.
    labels[blossom] = None;
    label_ends[blossom] = None;
    blossom_children[blossom].clear();
    blossom_endpoints[blossom].clear();
    blossom_base[blossom] = None;
    best_edge[blossom] = None;
    unused_blossoms.push(blossom);
    Ok(())
}

/// Swap matched/unmatched edges over an alternating path through blossom b
/// between vertex v and the base vertex. Keep blossom bookkeeping consistent.
fn augment_blossom(
    blossom: usize,
    node: usize,
    num_nodes: usize,
    blossom_parents: &[Option<usize>],
    endpoints: &[usize],
    blossom_children: &mut Vec<Vec<usize>>,
    blossom_endpoints: &mut Vec<Vec<usize>>,
    blossom_base: &mut Vec<Option<usize>>,
    mate: &mut HashMap<usize, usize>,
) {
    // Bubble up through the blossom tree from vertex v to an immediate
    // sub-blossom of b.
    let mut tmp = node;
    while blossom_parents[tmp] != Some(blossom) {
        tmp = blossom_parents[tmp].unwrap();
    }
    // Recursively deal with the first sub-blossom.
    if tmp >= num_nodes {
        augment_blossom(
            tmp,
            node,
            num_nodes,
            blossom_parents,
            endpoints,
            blossom_children,
            blossom_endpoints,
            blossom_base,
            mate,
        );
    }
    // Decide in which direction we will go around the blossom.
    let i = blossom_children[blossom]
        .iter()
        .position(|x| *x == tmp)
        .unwrap();
    let mut j: i128 = i as i128;
    let j_step: i128;
    let endpoint_trick: usize = if i & 1 != 0 {
        // start index is odd; go forward and wrap.
        j -= blossom_children[blossom].len() as i128;
        j_step = 1;
        0
    } else {
        // Start index is even; go backward.
        j_step = -1;
        1
    };
    // Move along the blossom until we get to the base.
    while j != 0 {
        // Step to the next sub-blossom and augment it recursively.
        j += j_step;

        tmp = if j < 0 {
            let length = blossom_children[blossom].len();
            let index = length - j.abs() as usize;
            blossom_children[blossom][index]
        } else {
            blossom_children[blossom][j as usize]
        };
        let p = if j < 0 {
            let length = blossom_endpoints[blossom].len();
            let index = length - j.abs() as usize - endpoint_trick;
            blossom_endpoints[blossom][index] ^ endpoint_trick
        } else {
            blossom_endpoints[blossom][j as usize - endpoint_trick]
                ^ endpoint_trick
        };
        if tmp > num_nodes {
            augment_blossom(
                tmp,
                endpoints[p],
                num_nodes,
                blossom_parents,
                endpoints,
                blossom_children,
                blossom_endpoints,
                blossom_base,
                mate,
            );
        }
        j += j_step;
        if j < 0 {
            let length = blossom_children[blossom].len();
            let index = length - j.abs() as usize;
            tmp = blossom_children[blossom][index];
        } else {
            tmp = blossom_children[blossom][j as usize];
        }
        if tmp > num_nodes {
            augment_blossom(
                tmp,
                endpoints[p ^ 1],
                num_nodes,
                blossom_parents,
                endpoints,
                blossom_children,
                blossom_endpoints,
                blossom_base,
                mate,
            );
        }
        // Match the edge connecting those sub-blossoms.
        mate.insert(endpoints[p], p ^ 1);
        mate.insert(endpoints[p ^ 1], p);
    }
    // Rotate the list of sub-blossoms to put the new base at the front.
    let mut children: Vec<usize> = blossom_children[blossom][i..].to_vec();
    children.extend_from_slice(&blossom_children[blossom][..i]);
    blossom_children[blossom] = children;
    let mut endpoint: Vec<usize> = blossom_endpoints[blossom][i..].to_vec();
    endpoint.extend_from_slice(&blossom_endpoints[blossom][..i]);
    blossom_endpoints[blossom] = endpoint;
    blossom_base[blossom] = blossom_base[blossom_children[blossom][0]];
    assert!(blossom_base[blossom] == Some(node));
}

/// Swap matched/unmatched edges over an alternating path between two
/// single vertices. The augmenting path runs through edge k, which
/// connects a pair of S vertices.
fn augment_matching(
    edge: usize,
    num_nodes: usize,
    edges: &[(usize, usize, i128)],
    in_blossoms: &[usize],
    labels: &[Option<usize>],
    label_ends: &[Option<usize>],
    blossom_parents: &[Option<usize>],
    endpoints: &[usize],
    blossom_children: &mut Vec<Vec<usize>>,
    blossom_endpoints: &mut Vec<Vec<usize>>,
    blossom_base: &mut Vec<Option<usize>>,
    mate: &mut HashMap<usize, usize>,
) {
    let (v, w, _weight) = edges[edge];
    for (s_ref, p_ref) in [(v, 2 * edge + 1), (w, 2 * edge)].iter() {
        // Match vertex s to remote endpoint p. Then trace back from s
        // until we find a single vertex, swapping matched and unmatched
        // edges as we go.
        let mut s: usize = *s_ref;
        let mut p: usize = *p_ref;
        loop {
            let blossom_s = in_blossoms[s];
            assert!(labels[blossom_s] == Some(1));
            assert!(
                label_ends[blossom_s]
                    == mate.get(&blossom_base[blossom_s].unwrap()).copied()
            );
            // Augment through the S-blossom from s to base.
            if blossom_s >= num_nodes {
                augment_blossom(
                    blossom_s,
                    s,
                    num_nodes,
                    blossom_parents,
                    endpoints,
                    blossom_children,
                    blossom_endpoints,
                    blossom_base,
                    mate,
                );
            }
            // Update mate[s]
            mate.insert(s, p);
            // Trace one step back.
            if label_ends[blossom_s].is_none() {
                // Reached a single vertex; stop
                break;
            }
            let t = endpoints[label_ends[blossom_s].unwrap()];
            let blossom_t = in_blossoms[t];
            assert!(labels[blossom_t] == Some(2));
            // Trace one step back
            assert!(label_ends[blossom_t].is_some());
            s = endpoints[label_ends[blossom_t].unwrap()];
            let j = endpoints[label_ends[blossom_t].unwrap() ^ 1];
            // Augment through the T-blossom from j to base.
            assert!(blossom_base[blossom_t] == Some(t));
            if blossom_t >= num_nodes {
                augment_blossom(
                    blossom_t,
                    j,
                    num_nodes,
                    blossom_parents,
                    endpoints,
                    blossom_children,
                    blossom_endpoints,
                    blossom_base,
                    mate,
                );
            }
            // Update mate[j]
            mate.insert(j, label_ends[blossom_t].unwrap());
            // Keep the opposite endpoint;
            // it will be assigned to mate[s] in the next step.
            p = label_ends[blossom_t].unwrap() ^ 1;
        }
    }
}

/// Swap matched/unmatched edges over an alternating path between two
/// single vertices. The augmenting path runs through the edge, which
/// connects a pair of S vertices.
fn verify_optimum(
    max_cardinality: bool,
    num_nodes: usize,
    num_edges: usize,
    edges: &[(usize, usize, i128)],
    endpoints: &[usize],
    dual_var: &[i128],
    blossom_parents: &[Option<usize>],
    blossom_endpoints: &[Vec<usize>],
    blossom_base: &[Option<usize>],
    mate: &HashMap<usize, usize>,
) {
    let dual_var_node_min: i128 = *dual_var[..num_nodes].iter().min().unwrap();
    let node_dual_offset: i128 = if max_cardinality {
        // Vertices may have negative dual;
        // find a constant non-negative number to add to all vertex duals.
        max(0, -dual_var_node_min)
    } else {
        0
    };
    assert!(dual_var_node_min + node_dual_offset >= 0);
    assert!(*dual_var[num_nodes..].iter().min().unwrap() >= 0);
    // 0. all edges have non-negative slack and
    // 1. all matched edges have zero slack;
    for (edge, (i, j, weight)) in edges.iter().enumerate().take(num_edges) {
        let mut s = dual_var[*i] + dual_var[*j] - 2 * weight;
        let mut i_blossoms: Vec<usize> = vec![*i];
        let mut j_blossoms: Vec<usize> = vec![*j];
        while blossom_parents[*i_blossoms.last().unwrap()].is_some() {
            i_blossoms
                .push(blossom_parents[*i_blossoms.last().unwrap()].unwrap());
        }
        while blossom_parents[*j_blossoms.last().unwrap()].is_some() {
            j_blossoms
                .push(blossom_parents[*j_blossoms.last().unwrap()].unwrap());
        }
        i_blossoms.reverse();
        j_blossoms.reverse();
        for (blossom_i, blossom_j) in i_blossoms.iter().zip(j_blossoms.iter()) {
            if blossom_i != blossom_j {
                break;
            }
            s += 2 * dual_var[*blossom_i];
        }
        assert!(s >= 0);

        if (mate.get(i).is_some() && mate.get(i).unwrap() / 2 == edge)
            || (mate.get(j).is_some() && mate.get(j).unwrap() / 2 == edge)
        {
            assert!(mate[i] / 2 == edge && mate[j] / 2 == edge);
            assert!(s == 0);
        }
    }
    // 2. all single vertices have zero dual value;
    for (node, dual_var_node) in dual_var.iter().enumerate().take(num_nodes) {
        assert!(
            mate.get(&node).is_some() || dual_var_node + node_dual_offset == 0
        );
    }
    // 3. all blossoms with positive dual value are full.
    for blossom in num_nodes..2 * num_nodes {
        if blossom_base[blossom].is_some() && dual_var[blossom] > 0 {
            assert!(blossom_endpoints[blossom].len() % 2 == 1);
            for p in blossom_endpoints[blossom].iter().skip(1).step_by(2) {
                assert!(mate.get(&endpoints[*p]).copied() == Some(p ^ 1));
                assert!(mate.get(&endpoints[*p ^ 1]) == Some(p));
            }
        }
    }
}

/// Compute a maximum-weighted matching in the general undirected weighted
/// graph given by "edges". If `max_cardinality` is true, only
/// maximum-cardinality matchings are considered as solutions.
///
/// The algorithm is based on "Efficient Algorithms for Finding Maximum
/// Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986.
///
/// Based on networkx implementation
/// <https://github.com/networkx/networkx/blob/3351206a3ce5b3a39bb2fc451e93ef545b96c95b/networkx/algorithms/matching.py>
///
/// With reference to the standalone protoype implementation from:
/// <http://jorisvr.nl/article/maximum-matching>
///
/// <http://jorisvr.nl/files/graphmatching/20130407/mwmatching.py>
///
/// The function takes time O(n**3)
///
/// Arguments:
///
/// * `graph` - The undirected graph to compute the maximum weight matching for
/// * `max_cardinality` - If set to true compute the maximum-cardinality matching
///     with maximum weight among all maximum-cardinality matchings
/// * `weight_fn` - A callback function that will be give a edge reference and
///     expected to return a `i128` representing the weight of the edge
/// * `verify_optimum_flag`: If true an prior to returning an additional routine
///     to verify the optimal solution was found will be run after computing
///     the maximum weight matching. If it's true and the found matching is not
///     an optimal solution this function will panic. This option should
///     normally be only set true during testing.
///
/// # Example
/// ```rust
/// use retworkx_core::petgraph;
/// use retworkx_core::max_weight_matching::max_weight_matching;
/// use retworkx_core::Result;
///
/// use hashbrown::HashSet;
///
/// // Create a path graph
/// let g = petgraph::graph::UnGraph::<i32, i128>::from_edges(&[
///     (1, 2, 5), (2, 3, 11), (3, 4, 5)
/// ]);
///
/// // Run max weight matching with max cardinality set to false
/// let res: Result<HashSet<(usize, usize)>> = max_weight_matching(
///     &g, false, |e| Ok(*e.weight()), true
/// );
/// // Run max weight matching with max cardinality set to true
/// let maxc_res: Result<HashSet<(usize, usize)>> = max_weight_matching(
///     &g, true, |e| Ok(*e.weight()), true
/// );
///
/// let matching = res.unwrap();
/// let maxc_matching = maxc_res.unwrap();
/// // Check output
/// assert_eq!(matching.len(), 1);
/// assert!(matching.contains(&(2, 3)) || matching.contains(&(3, 2)));
/// assert_eq!(maxc_matching.len(), 2);
/// assert!(maxc_matching.contains(&(1, 2)) || maxc_matching.contains(&(2, 1)));
/// assert!(maxc_matching.contains(&(3, 4)) || maxc_matching.contains(&(4, 3)));
/// ```
pub fn max_weight_matching<G, F, E>(
    graph: G,
    max_cardinality: bool,
    mut weight_fn: F,
    verify_optimum_flag: bool,
) -> Result<HashSet<(usize, usize)>, E>
where
    G: EdgeCount
        + NodeCount
        + IntoNodeIdentifiers
        + GraphProp<EdgeType = Undirected>
        + GraphBase<NodeId = NodeIndex>
        + IntoEdges,
    F: FnMut(G::EdgeRef) -> Result<i128, E>,
{
    let num_edges = graph.edge_count();
    let num_nodes = graph.node_count();
    let mut out_set: HashSet<(usize, usize)> =
        HashSet::with_capacity(num_nodes);
    // Exit fast for graph without edges
    if num_edges == 0 {
        return Ok(HashSet::new());
    }
    // Node indicies in the PyGraph may not be contiguous however the
    // algorithm operates on contiguous indices 0..num_nodes. node_map maps
    // the PyGraph's NodeIndex to the contingous usize used inside the
    // algorithm
    let node_map: HashMap<NodeIndex, usize> = graph
        .node_identifiers()
        .enumerate()
        .map(|(index, node_index)| (node_index, index))
        .collect();
    let mut edges: Vec<(usize, usize, i128)> = Vec::with_capacity(num_edges);
    let mut max_weight: i128 = 0;
    for edge in graph.edge_references() {
        let edge_weight: i128 = weight_fn(edge)?;
        if edge_weight > max_weight {
            max_weight = edge_weight;
        };
        edges.push((
            node_map[&edge.source()],
            node_map[&edge.target()],
            edge_weight,
        ));
    }
    // If p is an edge endpoint
    // endpoints[p] is the node index to which endpoint p is attached
    let endpoints: Vec<usize> = (0..2 * num_edges)
        .map(|endpoint| {
            let edge_tuple = edges[endpoint / 2];
            let out_value: usize = if endpoint % 2 == 0 {
                edge_tuple.0
            } else {
                edge_tuple.1
            };
            out_value
        })
        .collect();
    // If v is a node/vertex
    // neighbor_endpoints[v] is the list of remote endpoints of the edges
    // attached to v.
    // Not modified by algorithm (only mut to initially construct contents).
    let mut neighbor_endpoints: Vec<Vec<usize>> =
        (0..num_nodes).map(|_| Vec::new()).collect();
    for edge in 0..num_edges {
        neighbor_endpoints[edges[edge].0].push(2 * edge + 1);
        neighbor_endpoints[edges[edge].1].push(2 * edge);
    }
    // If v is a vertex,
    // mate[v] is the remote endpoint of its matched edge, or None if it is
    // single (i.e. endpoint[mate[v]] is v's partner vertex).
    // Initially all vertices are single; updated during augmentation.
    let mut mate: HashMap<usize, usize> = HashMap::with_capacity(num_nodes);
    // If b is a top-level blossom
    // label[b] is 0 if b is unlabeled (free);
    //             1 if b is a S-vertex/blossom;
    //             2 if b is a T-vertex/blossom;
    // The label of a vertex/node is found by looking at the label of its
    // top-level containing blossom
    // If v is a node/vertex inside a T-Blossom,
    // label[v] is 2 if and only if v is reachable from an S_Vertex outside
    // the blossom.
    let mut labels: Vec<Option<usize>>;
    // If b is a labeled top-level blossom,
    // label_ends[b] is the remote endpoint of the edge through which b
    // obtained its label, or None if b's base vertex is single.
    // If v is a vertex inside a T-blossom and label[v] == 2,
    // label_ends[v] is the remote endpoint of the edge through which v is
    // reachable from outside the blossom
    let mut label_ends: Vec<Option<usize>>;
    // If v is a vertex/node
    // in_blossoms[v] is the top-level blossom to which v belongs.
    // If v is a top-level vertex, v is itself a blossom (a trivial blossom)
    // and in_blossoms[v] == v.
    // Initially all nodes are top-level trivial blossoms.
    let mut in_blossoms: Vec<usize> = (0..num_nodes).collect();
    // if b is a sub-blossom
    // blossom_parents[b] is its immediate parent
    // If b is a top-level blossom, blossom_parents[b] is None
    let mut blossom_parents: Vec<Option<usize>> = vec![None; 2 * num_nodes];
    // If b is a non-trivial (sub-)blossom,
    // blossom_children[b] is an ordered list of its sub-blossoms, starting with
    // the base and going round the blossom.
    let mut blossom_children: Vec<Vec<usize>> =
        (0..2 * num_nodes).map(|_| Vec::new()).collect();
    // If b is a (sub-)blossom,
    // blossombase[b] is its base VERTEX (i.e. recursive sub-blossom).
    let mut blossom_base: Vec<Option<usize>> =
        (0..num_nodes).map(Some).collect();
    blossom_base.append(&mut vec![None; num_nodes]);
    // If b is a non-trivial (sub-)blossom,
    // blossom_endpoints[b] is a list of endpoints on its connecting edges,
    // such that blossom_endpoints[b][i] is the local endpoint of
    // blossom_children[b][i] on the edge that connects it to
    // blossom_children[b][wrap(i+1)].
    let mut blossom_endpoints: Vec<Vec<usize>> =
        (0..2 * num_nodes).map(|_| Vec::new()).collect();
    // If v is a free vertex (or an unreached vertex inside a T-blossom),
    // best_edge[v] is the edge to an S-vertex with least slack,
    // or None if there is no such edge. If b is a (possibly trivial)
    // top-level S-blossom,
    // best_edge[b] is the least-slack edge to a different S-blossom,
    // or None if there is no such edge.
    let mut best_edge: Vec<Option<usize>>;
    // If b is a non-trivial top-level S-blossom,
    // blossom_best_edges[b] is a list of least-slack edges to neighboring
    // S-blossoms, or None if no such list has been computed yet.
    // This is used for efficient computation of delta3.
    let mut blossom_best_edges: Vec<Vec<usize>> =
        (0..2 * num_nodes).map(|_| Vec::new()).collect();
    let mut unused_blossoms: Vec<usize> = (num_nodes..2 * num_nodes).collect();
    // If v is a vertex,
    // dual_var[v] = 2 * u(v) where u(v) is the v's variable in the dual
    // optimization problem (multiplication by two ensures integer values
    // throughout the algorithm if all edge weights are integers).
    // If b is a non-trivial blossom,
    // dual_var[b] = z(b) where z(b) is b's variable in the dual optimization
    // problem.
    // dual_var is for vertex v in 0..num_nodes and blossom b in
    // num_nodes..2* num nodes
    let mut dual_var: Vec<i128> = vec![max_weight; num_nodes];
    dual_var.append(&mut vec![0; num_nodes]);
    // If allowed_edge[k] is true, edge k has zero slack in the optimization
    // problem; if allowed_edge[k] is false, the edge's slack may or may not
    // be zero.
    let mut allowed_edge: Vec<bool>;
    // Queue of newly discovered S-vertices
    let mut queue: Vec<usize> = Vec::with_capacity(num_nodes);

    // Main loop: continue until no further improvement is possible
    for _ in 0..num_nodes {
        // Each iteration of this loop is a "stage".
        // A stage finds an augmenting path and uses that to improve
        // the matching.

        // Removal labels from top-level blossoms/vertices
        labels = vec![Some(0); 2 * num_nodes];
        label_ends = vec![None; 2 * num_nodes];
        // Forget all about least-slack edges.
        best_edge = vec![None; 2 * num_nodes];
        blossom_best_edges
            .splice(num_nodes.., (0..num_nodes).map(|_| Vec::new()));
        // Loss of labeling means that we can not be sure that currently
        // allowable edges remain allowable througout this stage.
        allowed_edge = vec![false; num_edges];
        // Make queue empty
        queue.clear();
        // Label single blossom/vertices with S and put them in queue.
        for v in 0..num_nodes {
            if mate.get(&v).is_none() && labels[in_blossoms[v]] == Some(0) {
                assign_label(
                    v,
                    1,
                    None,
                    num_nodes,
                    &in_blossoms,
                    &mut labels,
                    &mut label_ends,
                    &mut best_edge,
                    &mut queue,
                    &blossom_children,
                    &blossom_base,
                    &endpoints,
                    &mate,
                )?;
            }
        }
        // Loop until we succeed in augmenting the matching.
        let mut augmented = false;
        loop {
            // Each iteration of this loop is a "substage".
            // A substage tries to find an augmenting path;
            // if found, the path is used to improve the matching and
            // the stage ends. If there is no augmenting path, the
            // primal-dual method is used to find some slack out of
            // the dual variables.

            // Continue labeling until all vertices which are reachable
            // through an alternating path have got a label.
            while !queue.is_empty() && !augmented {
                // Take an S vertex from the queue
                let v = queue.pop().unwrap();
                assert!(labels[in_blossoms[v]] == Some(1));

                // Scan its neighbors
                for p in &neighbor_endpoints[v] {
                    let k = *p / 2;
                    let mut kslack = 0;
                    let w = endpoints[*p];
                    // w is a neighbor of v
                    if in_blossoms[v] == in_blossoms[w] {
                        // this edge is internal to a blossom; ignore it
                        continue;
                    }
                    if !allowed_edge[k] {
                        kslack = slack(k, &dual_var, &edges);
                        if kslack <= 0 {
                            // edge k has zero slack -> it is allowable
                            allowed_edge[k] = true;
                        }
                    }
                    if allowed_edge[k] {
                        if labels[in_blossoms[w]] == Some(0) {
                            // (C1) w is a free vertex;
                            // label w with T and label its mate with S (R12).
                            assign_label(
                                w,
                                2,
                                Some(*p ^ 1),
                                num_nodes,
                                &in_blossoms,
                                &mut labels,
                                &mut label_ends,
                                &mut best_edge,
                                &mut queue,
                                &blossom_children,
                                &blossom_base,
                                &endpoints,
                                &mate,
                            )?;
                        } else if labels[in_blossoms[w]] == Some(1) {
                            // (C2) w is an S-vertex (not in the same blossom);
                            // follow back-links to discover either an
                            // augmenting path or a new blossom.
                            let base = scan_blossom(
                                v,
                                w,
                                &in_blossoms,
                                &blossom_base,
                                &endpoints,
                                &mut labels,
                                &label_ends,
                                &mate,
                            );
                            match base {
                                // Found a new blossom; add it to the blossom
                                // bookkeeping and turn it into an S-blossom.
                                Some(base) => add_blossom(
                                    base,
                                    k,
                                    &mut blossom_children,
                                    num_nodes,
                                    &edges,
                                    &mut in_blossoms,
                                    &mut dual_var,
                                    &mut labels,
                                    &mut label_ends,
                                    &mut best_edge,
                                    &mut queue,
                                    &mut blossom_base,
                                    &endpoints,
                                    &mut blossom_endpoints,
                                    &mut unused_blossoms,
                                    &mut blossom_best_edges,
                                    &mut blossom_parents,
                                    &neighbor_endpoints,
                                    &mate,
                                )?,
                                // Found an augmenting path; augment the
                                // matching and end this stage.
                                None => {
                                    augment_matching(
                                        k,
                                        num_nodes,
                                        &edges,
                                        &in_blossoms,
                                        &labels,
                                        &label_ends,
                                        &blossom_parents,
                                        &endpoints,
                                        &mut blossom_children,
                                        &mut blossom_endpoints,
                                        &mut blossom_base,
                                        &mut mate,
                                    );
                                    augmented = true;
                                    break;
                                }
                            };
                        } else if labels[w] == Some(0) {
                            // w is inside a T-blossom, but w itself has not
                            // yet been reached from outside the blossom;
                            // mark it as reached (we need this to relabel
                            // during T-blossom expansion).
                            assert!(labels[in_blossoms[w]] == Some(2));
                            labels[w] = Some(2);
                            label_ends[w] = Some(*p ^ 1);
                        }
                    } else if labels[in_blossoms[w]] == Some(1) {
                        // Keep track of the least-slack non-allowable edge to
                        // a different S-blossom
                        let blossom = in_blossoms[v];
                        if best_edge[blossom].is_none()
                            || kslack
                                < slack(
                                    best_edge[blossom].unwrap(),
                                    &dual_var,
                                    &edges,
                                )
                        {
                            best_edge[blossom] = Some(k);
                        }
                    } else if labels[w] == Some(0) {
                        // w is a free vertex (or an unreached vertex inside
                        // a T-blossom) but we can not reach it yet;
                        // keep track of the least-slack edge that reaches w.
                        if best_edge[w].is_none()
                            || kslack
                                < slack(
                                    best_edge[w].unwrap(),
                                    &dual_var,
                                    &edges,
                                )
                        {
                            best_edge[w] = Some(k)
                        }
                    }
                }
            }
            if augmented {
                break;
            }
            // There is no augmenting path under these constraints;
            // compute delta and reduce slack in the optimization problem.
            // (Note that our vertex dual variables, edge slacks and delta's
            // are pre-multiplied by two.)
            let mut delta_type = -1;
            let mut delta: Option<i128> = None;
            let mut delta_edge: Option<usize> = None;
            let mut delta_blossom: Option<usize> = None;

            // Compute delta1: the minimum value of any vertex dual.
            if !max_cardinality {
                delta_type = 1;
                delta = Some(*dual_var[..num_nodes].iter().min().unwrap());
            }

            // Compute delta2: the minimum slack on any edge between
            // an S-vertex and a free vertex.
            for v in 0..num_nodes {
                if labels[in_blossoms[v]] == Some(0) && best_edge[v].is_some() {
                    let d = slack(best_edge[v].unwrap(), &dual_var, &edges);
                    if delta_type == -1 || Some(d) < delta {
                        delta = Some(d);
                        delta_type = 2;
                        delta_edge = best_edge[v];
                    }
                }
            }

            // Compute delta3: half the minimum slack on any edge between a
            // pair of S-blossoms
            for blossom in 0..2 * num_nodes {
                if blossom_parents[blossom].is_none()
                    && labels[blossom] == Some(1)
                    && best_edge[blossom].is_some()
                {
                    let kslack =
                        slack(best_edge[blossom].unwrap(), &dual_var, &edges);
                    assert!(kslack % 2 == 0);
                    let d = Some(kslack / 2);
                    if delta_type == -1 || d < delta {
                        delta = d;
                        delta_type = 3;
                        delta_edge = best_edge[blossom];
                    }
                }
            }

            // Compute delta4: minimum z variable of any T-blossom
            for blossom in num_nodes..2 * num_nodes {
                if blossom_base[blossom].is_some()
                    && blossom_parents[blossom].is_none()
                    && labels[blossom] == Some(2)
                    && (delta_type == -1 || dual_var[blossom] < delta.unwrap())
                {
                    delta = Some(dual_var[blossom]);
                    delta_type = 4;
                    delta_blossom = Some(blossom);
                }
            }
            if delta_type == -1 {
                // No further improvement possible; max-cardinality optimum
                // reached. Do a final delta update to make the optimum
                // verifyable
                assert!(max_cardinality);
                delta_type = 1;
                delta =
                    Some(max(0, *dual_var[..num_nodes].iter().min().unwrap()));
            }

            // Update dual variables according to delta.
            for v in 0..num_nodes {
                if labels[in_blossoms[v]] == Some(1) {
                    // S-vertex: 2*u = 2*u - 2*delta
                    dual_var[v] -= delta.unwrap();
                } else if labels[in_blossoms[v]] == Some(2) {
                    // T-vertex: 2*u = 2*u + 2*delta
                    dual_var[v] += delta.unwrap();
                }
            }
            for b in num_nodes..2 * num_nodes {
                if blossom_base[b].is_some() && blossom_parents[b].is_none() {
                    if labels[b] == Some(1) {
                        // top-level S-blossom: z = z + 2*delta
                        dual_var[b] += delta.unwrap();
                    } else if labels[b] == Some(2) {
                        // top-level T-blossom: z = z - 2*delta
                        dual_var[b] -= delta.unwrap();
                    }
                }
            }
            // Take action at the point where minimum delta occured.
            if delta_type == 1 {
                // No further improvement possible; optimum reached
                break;
            } else if delta_type == 2 {
                // Use the least-slack edge to continue the search.
                allowed_edge[delta_edge.unwrap()] = true;
                let (mut i, mut j, _weight) = edges[delta_edge.unwrap()];
                if labels[in_blossoms[i]] == Some(0) {
                    mem::swap(&mut i, &mut j);
                }
                assert!(labels[in_blossoms[i]] == Some(1));
                queue.push(i);
            } else if delta_type == 3 {
                // Use the least-slack edge to continue the search.
                allowed_edge[delta_edge.unwrap()] = true;
                let (i, _j, _weight) = edges[delta_edge.unwrap()];
                assert!(labels[in_blossoms[i]] == Some(1));
                queue.push(i);
            } else if delta_type == 4 {
                // Expand the least-z blossom
                expand_blossom(
                    delta_blossom.unwrap(),
                    false,
                    num_nodes,
                    &mut blossom_children,
                    &mut blossom_parents,
                    &mut in_blossoms,
                    &dual_var,
                    &mut labels,
                    &mut label_ends,
                    &mut best_edge,
                    &mut queue,
                    &mut blossom_base,
                    &endpoints,
                    &mate,
                    &mut blossom_endpoints,
                    &mut allowed_edge,
                    &mut unused_blossoms,
                )?;
            }
            // end of this substage
        }
        // Stop when no more augment paths can be found
        if !augmented {
            break;
        }

        // End of a stage; expand all S-blossoms which have a dual_var == 0
        for blossom in num_nodes..2 * num_nodes {
            if blossom_parents[blossom].is_none()
                && blossom_base[blossom].is_some()
                && labels[blossom] == Some(1)
                && dual_var[blossom] == 0
            {
                expand_blossom(
                    blossom,
                    true,
                    num_nodes,
                    &mut blossom_children,
                    &mut blossom_parents,
                    &mut in_blossoms,
                    &dual_var,
                    &mut labels,
                    &mut label_ends,
                    &mut best_edge,
                    &mut queue,
                    &mut blossom_base,
                    &endpoints,
                    &mate,
                    &mut blossom_endpoints,
                    &mut allowed_edge,
                    &mut unused_blossoms,
                )?;
            }
        }
    }
    if verify_optimum_flag {
        verify_optimum(
            max_cardinality,
            num_nodes,
            num_edges,
            &edges,
            &endpoints,
            &dual_var,
            &blossom_parents,
            &blossom_endpoints,
            &blossom_base,
            &mate,
        );
    }

    // Transform mate[] such that mate[v] is the vertex to which v is paired
    // Also handle holes in node indices from PyGraph node removals by mapping
    // linear index to node index.
    let mut seen: HashSet<(usize, usize)> =
        HashSet::with_capacity(2 * num_nodes);
    let node_list: Vec<NodeIndex> = graph.node_identifiers().collect();
    for (index, node) in mate.iter() {
        let tmp = (
            node_list[*index].index(),
            node_list[endpoints[*node]].index(),
        );
        let rev_tmp = (
            node_list[endpoints[*node]].index(),
            node_list[*index].index(),
        );
        if !seen.contains(&tmp) && !seen.contains(&rev_tmp) {
            out_set.insert(tmp);
            seen.insert(tmp);
            seen.insert(rev_tmp);
        }
    }
    Ok(out_set)
}