Include support for vertex invariants in isomorphism()

This also properly handles self-loops and thus fixes #180.

src/graph/topology/graph_isomorphism.cc

@@ -25,14 +25,46 @@ using namespace boost;
 
 struct check_iso
 {
-    template <class Graph1, class Graph2, class IsoMap, class VertexIndexMap>
-    void operator()(Graph1& g1, Graph2* g2, IsoMap map, VertexIndexMap index1,
+
+    template <class Graph1, class Graph2, class IsoMap, class InvMap,
+              class VertexIndexMap>
+    void operator()(Graph1& g1, Graph2* g2, InvMap cinv_map1, InvMap cinv_map2,
+                    int64_t max_inv, IsoMap map, VertexIndexMap index1,
                     VertexIndexMap index2, bool& result) const
     {
-        result = isomorphism(g1, *g2, isomorphism_map(map).
+        auto inv_map1 = cinv_map1.get_unchecked(num_vertices(g1));
+        auto inv_map2 = cinv_map2.get_unchecked(num_vertices(*g2));
+
+        vinv_t<decltype(inv_map1)> vinv1(inv_map1, max_inv);
+        vinv_t<decltype(inv_map2)> vinv2(inv_map2, max_inv);
+
+        result = isomorphism(g1, *g2,
+                             isomorphism_map(map.get_unchecked(num_vertices(g1))).
+                             vertex_invariant1(vinv1).
+                             vertex_invariant2(vinv2).
                              vertex_index1_map(index1).
                              vertex_index2_map(index2));
     }
+
+    template <class Prop>
+    struct vinv_t
+    {
+        vinv_t(Prop& prop, int64_t max)
+            : _prop(prop), _max(max) {}
+        Prop& _prop;
+        int64_t _max;
+
+        template <class Vertex>
+        int64_t operator()(Vertex v) const
+        {
+            return _prop[v];
+        };
+
+        int64_t max() const { return _max; }
+
+        typedef int64_t result_type;
+        typedef size_t argument_type;
+    };
 };
 
 struct directed_graph_view_pointers:
@@ -49,10 +81,22 @@ typedef property_map_types::apply<integer_types,
     vertex_props_t;
 
 bool check_isomorphism(GraphInterface& gi1, GraphInterface& gi2,
-                       boost::any iso_map)
+                       boost::any ainv_map1, boost::any ainv_map2,
+                       int64_t max_inv, boost::any aiso_map)
 {
     bool result;
 
+    typedef property_map_type::apply<int32_t,
+                                     GraphInterface::vertex_index_map_t>::type
+        iso_map_t;
+    auto iso_map = any_cast<iso_map_t>(aiso_map);
+
+    typedef property_map_type::apply<int64_t,
+                                     GraphInterface::vertex_index_map_t>::type
+        inv_map_t;
+    auto inv_map1 = any_cast<inv_map_t>(ainv_map1);
+    auto inv_map2 = any_cast<inv_map_t>(ainv_map2);
+
     if (gi1.GetDirected() != gi2.GetDirected())
         return false;
     if (gi1.GetDirected())
@@ -60,20 +104,22 @@ bool check_isomorphism(GraphInterface& gi1, GraphInterface& gi2,
         run_action<graph_tool::detail::always_directed>()
             (gi1, std::bind(check_iso(),
                             placeholders::_1, placeholders::_2,
-                            placeholders::_3, gi1.GetVertexIndex(),
+                            inv_map1, inv_map2, max_inv, iso_map,
+                            gi1.GetVertexIndex(),
                             gi2.GetVertexIndex(), std::ref(result)),
-             directed_graph_view_pointers(), vertex_props_t())
-            (gi2.GetGraphView(), iso_map);
+             directed_graph_view_pointers())
+            (gi2.GetGraphView());
     }
     else
     {
         run_action<graph_tool::detail::never_directed>()
             (gi1, std::bind(check_iso(),
-                            placeholders::_1, placeholders::_2, placeholders::_3,
+                            placeholders::_1, placeholders::_2,
+                            inv_map1, inv_map2, max_inv, iso_map,
                             gi1.GetVertexIndex(),
                             gi2.GetVertexIndex(), std::ref(result)),
-             undirected_graph_view_pointers(), vertex_props_t())
-            (gi2.GetGraphView(), iso_map);
+             undirected_graph_view_pointers())
+            (gi2.GetGraphView());
     }
 
     return result;

src/graph/topology/graph_topology.cc

@@ -24,7 +24,8 @@ using namespace boost::python;
 using namespace graph_tool;
 
 bool check_isomorphism(GraphInterface& gi1, GraphInterface& gi2,
-                       boost::any iso_map);
+                       boost::any ainv_map1, boost::any ainv_map2,
+                       int64_t max_inv, boost::any aiso_map);
 void get_kruskal_spanning_tree(GraphInterface& gi, boost::any weight_map,
                                boost::any tree_map);
 void get_prim_spanning_tree(GraphInterface& gi, size_t root,

src/graph_tool/topology/__init__.py

@@ -68,7 +68,9 @@ dl_import("from . import libgraph_tool_topology")
 from .. import _prop, Vector_int32_t, _check_prop_writable, \
      _check_prop_scalar, _check_prop_vector, Graph, PropertyMap, GraphView,\
      libcore, _get_rng, _degree, perfect_prop_hash
+from .. stats import label_self_loops
 import random, sys, numpy
+
 __all__ = ["isomorphism", "subgraph_isomorphism", "mark_subgraph",
            "max_cardinality_matching", "max_independent_vertex_set",
            "min_spanning_tree", "random_spanning_tree", "dominator_tree",
@@ -156,12 +158,34 @@ def similarity(g1, g2, label1=None, label2=None, norm=True):
     return s
 
 
-def isomorphism(g1, g2, isomap=False):
+def isomorphism(g1, g2, vertex_inv1=None, vertex_inv2=None, isomap=False):
     r"""Check whether two graphs are isomorphic.
 
-    If `isomap` is True, a vertex :class:`~graph_tool.PropertyMap` with the
+    Parameters
+    ----------
+    g1 : :class:`~graph_tool.Graph`
+        First graph.
+    g2 : :class:`~graph_tool.Graph`
+        Second graph.
+    vertex_inv1 : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
+        Vertex invariant of the first graph. Only vertices with with the same
+        invariants are considered in the isomorphism.
+    vertex_inv2 : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
+        Vertex invariant of the second graph. Only vertices with with the same
+        invariants are considered in the isomorphism.
+    isomap : ``bool`` (optional, default: ``False``)
+        If ``True``, a vertex :class:`~graph_tool.PropertyMap` with the
         isomorphism mapping is returned as well.
 
+    Returns
+    -------
+    is_isomorphism : ``bool``
+        ``True`` if both graphs are isomorphic, otherwise ``False``.
+    isomap : :class:`~graph_tool.PropertyMap`
+         Isomorphism mapping corresponding to a property map belonging to the
+         first graph which maps its vertices to their corresponding vertices of
+         the second graph.
+
     Examples
     --------
     .. testcode::
@@ -182,8 +206,34 @@ def isomorphism(g1, g2, isomap=False):
 
     """
     imap = g1.new_vertex_property("int32_t")
+    if vertex_inv1 is None:
+        vertex_inv1 = g1.degree_property_map("total").copy("int64_t")
+    else:
+        vertex_inv1 = vertex_inv1.copy("int64_t")
+        d = g1.degree_property_map("total")
+        vertex_inv1.fa += (vertex_inv1.fa.max() + 1) * d.a
+    if vertex_inv2 is None:
+        vertex_inv2 = g2.degree_property_map("total").copy("int64_t")
+    else:
+        vertex_inv2 = vertex_inv2.copy("int64_t")
+        d = g2.degree_property_map("total")
+        vertex_inv2.fa += (vertex_inv2.fa.max() + 1) * d.a
+
+    inv_max = max(vertex_inv1.fa.max(),vertex_inv2.fa.max()) + 1
+
+    l1 = label_self_loops(g1, mark_only=True)
+    if l1.fa.max() > 0:
+        g1 = GraphView(g1, efilt=1 - l1.fa)
+
+    l2 = label_self_loops(g2, mark_only=True)
+    if l2.fa.max() > 0:
+        g2 = GraphView(g2, efilt=1 - l2.fa)
+
     iso = libgraph_tool_topology.\
            check_isomorphism(g1._Graph__graph, g2._Graph__graph,
+                             _prop("v", g1, vertex_inv1),
+                             _prop("v", g2, vertex_inv2),
+                             inv_max,
                              _prop("v", g1, imap))
     if isomap:
         return iso, imap