Fix is_subgraph_isomorphic with parallel edges.

Closes Qiskit#429.
The solution is build an adjacency matrix with entries equal
to the number of edges connecting a pair of nodes instead of
0/1 values

releasenotes/notes/bugfix-vf2-parallel-edges-f8be125e255753f7.yaml

@@ -0,0 +1,5 @@
+---
+fixes:
+  - |
+    Fixes the output of :func:`~retworkx.is_subgraph_isomorphic`
+    if the input graphs have parallel edges.

src/isomorphism/vf2.rs

@@ -15,7 +15,6 @@
 // to handle PyDiGraph inputs instead of petgraph's generic Graph. However it has
 // since diverged significantly from the original petgraph implementation.
 
-use fixedbitset::FixedBitSet;
 use std::cmp::{Ordering, Reverse};
 use std::iter::Iterator;
 use std::marker;
@@ -29,9 +28,7 @@ use pyo3::PyTraverseError;
 
 use petgraph::stable_graph::NodeIndex;
 use petgraph::stable_graph::StableGraph;
-use petgraph::visit::{
-    EdgeRef, GetAdjacencyMatrix, IntoEdgeReferences, NodeIndexable,
-};
+use petgraph::visit::{EdgeRef, IntoEdgeReferences, NodeIndexable};
 use petgraph::EdgeType;
 use petgraph::{Directed, Incoming, Outgoing, Undirected};
 
@@ -41,6 +38,48 @@ use crate::iterators::NodeMap;
 
 type StablePyGraph<Ty> = StableGraph<PyObject, PyObject, Ty>;
 
+/// Returns the adjacency_matrix of a graph with `(i, j)` entry equal
+/// to number of edges from node `i` to node `j`.
+fn adjacency_matrix<Ty: EdgeType>(graph: &StablePyGraph<Ty>) -> Vec<usize> {
+    let n = graph.node_bound();
+    let mut matrix = vec![0; n * n];
+    for edge in graph.edge_references() {
+        let i = edge.source().index() * n + edge.target().index();
+        matrix[i] += 1;
+        if !graph.is_directed() {
+            let j = edge.source().index() + n * edge.target().index();
+            matrix[j] += 1;
+        }
+    }
+    matrix
+}
+
+/// Returns the number of edges from node `a` to node `b`.
+fn edge_multiplicity<Ty: EdgeType>(
+    graph: &StablePyGraph<Ty>,
+    matrix: &[usize],
+    a: NodeIndex,
+    b: NodeIndex,
+) -> usize {
+    let n = graph.node_bound();
+    let index = n * a.index() + b.index();
+    matrix[index]
+}
+
+/// Nodes `a`, `b` are adjacent if the number of edges
+/// from node `a` to node `b` are greater than `val`.
+fn is_adjacent<Ty: EdgeType>(
+    graph: &StablePyGraph<Ty>,
+    matrix: &[usize],
+    a: NodeIndex,
+    b: NodeIndex,
+    val: usize,
+) -> bool {
+    let n = graph.node_bound();
+    let index = n * a.index() + b.index();
+    matrix[index] >= val
+}
+
 trait NodeSorter<Ty>
 where
     Ty: EdgeType,
@@ -54,8 +93,12 @@ where
     ) -> (StablePyGraph<Ty>, HashMap<usize, usize>) {
         let order = self.sort(graph);
 
-        let mut new_graph = StablePyGraph::<Ty>::default();
-        let mut id_map: HashMap<NodeIndex, NodeIndex> = HashMap::new();
+        let mut new_graph = StablePyGraph::<Ty>::with_capacity(
+            graph.node_count(),
+            graph.edge_count(),
+        );
+        let mut id_map: HashMap<NodeIndex, NodeIndex> =
+            HashMap::with_capacity(graph.node_count());
         for node_index in order {
             let node_data = graph.node_weight(node_index).unwrap();
             let new_index = new_graph.add_node(node_data.clone_ref(py));
@@ -74,7 +117,7 @@ where
     }
 }
 
-// Sort nodes based on node ids.
+/// Sort nodes based on node ids.
 struct DefaultIdSorter;
 
 impl<Ty> NodeSorter<Ty> for DefaultIdSorter
@@ -86,7 +129,7 @@ where
     }
 }
 
-// Sort nodes based on  VF2++ heuristic.
+/// Sort nodes based on  VF2++ heuristic.
 struct Vf2ppSorter;
 
 impl<Ty> NodeSorter<Ty> for Vf2ppSorter
@@ -224,7 +267,7 @@ where
     ins: Vec<usize>,
     out_size: usize,
     ins_size: usize,
-    adjacency_matrix: FixedBitSet,
+    adjacency_matrix: Vec<usize>,
     generation: usize,
     _etype: marker::PhantomData<Directed>,
 }
@@ -236,7 +279,7 @@ where
     pub fn new(graph: StablePyGraph<Ty>) -> Self {
         let c0 = graph.node_count();
         let is_directed = graph.is_directed();
-        let adjacency_matrix = graph.adjacency_matrix();
+        let adjacency_matrix = adjacency_matrix(&graph);
         Vf2State {
             graph,
             mapping: vec![NodeIndex::end(); c0],
@@ -593,10 +636,19 @@ where
                 if m_neigh == end {
                     continue;
                 }
-                let has_edge = st[1 - j].graph.is_adjacent(
+                let val = edge_multiplicity(
+                    &st[j].graph,
+                    &st[j].adjacency_matrix,
+                    nodes[j],
+                    n_neigh,
+                );
+
+                let has_edge = is_adjacent(
+                    &st[1 - j].graph,
                     &st[1 - j].adjacency_matrix,
                     nodes[1 - j],
                     m_neigh,
+                    val,
                 );
                 if !has_edge {
                     return Ok(false);
@@ -622,10 +674,19 @@ where
                     if m_neigh == end {
                         continue;
                     }
-                    let has_edge = st[1 - j].graph.is_adjacent(
+                    let val = edge_multiplicity(
+                        &st[j].graph,
+                        &st[j].adjacency_matrix,
+                        n_neigh,
+                        nodes[j],
+                    );
+
+                    let has_edge = is_adjacent(
+                        &st[1 - j].graph,
                         &st[1 - j].adjacency_matrix,
                         m_neigh,
                         nodes[1 - j],
+                        val,
                     );
                     if !has_edge {
                         return Ok(false);

tests/digraph/test_isomorphic.py

@@ -297,6 +297,17 @@ def test_digraph_non_isomorphic_rule_ins_incoming(self):
             retworkx.is_isomorphic(graph, second_graph, id_order=True)
         )
 
+    def test_isomorphic_parallel_edges(self):
+        first = retworkx.PyDiGraph()
+        first.extend_from_edge_list([
+            (0, 1), (0, 1), (1, 2), (2, 3)
+        ])
+        second = retworkx.PyDiGraph()
+        second.extend_from_edge_list([
+            (0, 1), (1, 2), (1, 2), (2, 3)
+        ])
+        self.assertFalse(retworkx.is_isomorphic(first, second))
+
     def test_digraph_isomorphic_insufficient_call_limit(self):
         graph = retworkx.generators.directed_path_graph(5)
         self.assertFalse(retworkx.is_isomorphic(graph, graph, call_limit=2))

tests/graph/test_isomorphic.py

@@ -264,6 +264,17 @@ def test_graph_isomorphic_self_loop(self):
         graph.add_edges_from_no_data([(0, 0), (0, 1)])
         self.assertTrue(retworkx.is_isomorphic(graph, graph))
 
+    def test_isomorphic_parallel_edges(self):
+        first = retworkx.PyGraph()
+        first.extend_from_edge_list([
+            (0, 1), (0, 1), (1, 2), (2, 3)
+        ])
+        second = retworkx.PyGraph()
+        second.extend_from_edge_list([
+            (0, 1), (1, 2), (1, 2), (2, 3)
+        ])
+        self.assertFalse(retworkx.is_isomorphic(first, second))
+
     def test_graph_isomorphic_insufficient_call_limit(self):
         graph = retworkx.generators.path_graph(5)
         self.assertFalse(retworkx.is_isomorphic(graph, graph, call_limit=2))

tests/graph/test_subgraph_isomorphic.py

@@ -219,6 +219,22 @@ def test_non_induced_subgraph_isomorphic(self):
                     )
                 )
 
+    def test_subgraph_isomorphic_parallel_edges(self):
+        first = retworkx.PyGraph()
+        first.extend_from_edge_list([
+            (0, 1), (1, 2), (2, 3)
+        ])
+        second = retworkx.PyGraph()
+        second.extend_from_edge_list([
+            (0, 1), (0, 1)
+        ])
+        self.assertFalse(
+            retworkx.is_subgraph_isomorphic(first, second, induced=True)
+        )
+        self.assertFalse(
+            retworkx.is_subgraph_isomorphic(first, second, induced=False)
+        )
+
     def test_non_induced_grid_subgraph_isomorphic(self):
         g_a = retworkx.generators.grid_graph(2, 2)
         g_b = retworkx.PyGraph()
@@ -233,6 +249,22 @@ def test_non_induced_grid_subgraph_isomorphic(self):
             retworkx.is_subgraph_isomorphic(g_a, g_b, induced=False)
         )
 
+    def test_non_induced_subgraph_isomorphic_parallel_edges(self):
+        first = retworkx.PyGraph()
+        first.extend_from_edge_list([
+            (0, 1), (0, 1), (1, 2), (1, 2)
+        ])
+        second = retworkx.PyGraph()
+        second.extend_from_edge_list([
+            (0, 1), (1, 2), (1, 2)
+        ])
+        self.assertFalse(
+            retworkx.is_subgraph_isomorphic(first, second, induced=True)
+        )
+        self.assertTrue(
+            retworkx.is_subgraph_isomorphic(first, second, induced=False)
+        )
+
     def test_subgraph_vf2_mapping(self):
         graph = retworkx.generators.grid_graph(10, 10)
         second_graph = retworkx.generators.grid_graph(2, 2)