isomorphism-impl.w

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

As discussed in \S\ref{sec:vertex-invariants}, we can quickly rule out
the possibility of any isomorphism between two graphs by checking to
see if the vertex invariants can match up. We sort both vectors of vertex
invariants, and then check to see if they are equal.

@d Quick return if the graph's invariants do not match
@{
{ // check if the graph's invariants do not match
  invar_vec1_t invar1_tmp(invar1_vec);
  invar_vec2_t invar2_tmp(invar2_vec);
  std::sort(invar1_tmp.begin(), invar1_tmp.end());
  std::sort(invar2_tmp.begin(), invar2_tmp.end());
  if (! std::equal(invar1_tmp.begin(), invar1_tmp.end(), 
                   invar2_tmp.begin()))
    return false;
}
@}

Next we compute the invariant multiplicity, the number of vertices
with the same invariant number. The \code{invar\_mult} vector is
indexed by invariant number. We loop through all the vertices in the
graph to record the multiplicity.

@d Compute invariant multiplicity
@{
std::vector<std::size_t> invar_mult(invariant1.max(), 0);
for (tie(i1, i1_end) = vertices(g1); i1 != i1_end; ++i1)      
  ++invar_mult[invar1[*i1]];
@}

\noindent We then order the vertices by their invariant multiplicity.
This will allow us to search the more constrained vertices first.
Since we will need to know the permutation from the original order to
the new order, we do not sort the vertices directly. Instead we sort
the vertex indices, creating the \code{perm} array.  Once sorted, this
array provides a mapping from the new index to the old index.
We then use the \code{permute} function to sort the vertices of
the graph, which we store in the \code{g1\_vertices} vector.

@d Sort vertices by invariant multiplicity
@{
std::vector<size_type> perm;
integer_range<size_type> range(0, num_vertices(g1));
std::copy(range.begin(), range.end(), std::back_inserter(perm));
std::sort(perm.begin(), perm.end(),
          detail::compare_invariant_multiplicity(invar1_vec.begin(),
                                                 invar_mult.begin()));

std::vector<vertex1_t> g1_vertices;
for (tie(i1, i1_end) = vertices(g1); i1 != i1_end; ++i1)
  g1_vertices.push_back(*i1);
permute(g1_vertices.begin(), g1_vertices.end(), perm.begin());
@}

\noindent The definition of the \code{compare\_multiplicity} predicate
is shown below. This predicate provides the glue that binds
\code{std::sort} to our current purpose.

@d Compare multiplicity predicate
@{
namespace detail {
  template <typename InvarMap, typename MultMap>
  struct compare_invariant_multiplicity_predicate
  {
    compare_invariant_multiplicity_predicate(InvarMap i, MultMap m)
      : m_invar(i), m_mult(m) { }

    template <typename Vertex>
    bool operator()(const Vertex& x, const Vertex& y) const
      { return m_mult[m_invar[x]] < m_mult[m_invar[y]]; }

    InvarMap m_invar;
    MultMap m_mult;
  };
  template <typename InvarMap, typename MultMap>
  compare_invariant_multiplicity_predicate<InvarMap, MultMap>
  compare_invariant_multiplicity(InvarMap i, MultMap m) {
    return compare_invariant_multiplicity_predicate<InvarMap, MultMap>(i,m);
  }
} // namespace detail
@}


\subsection{Ordering by DFS Discover Time}

To implement the ``visit adjacent vertices first'' heuristic, we order
the vertices according to DFS discover time.  We replace the ordering
in \code{perm} with the new DFS ordering. Again, we use \code{permute}
to sort the vertices of graph \code{g1}.

@d Order the vertices by DFS discover time
@{
{
  perm.clear();
  @<Compute DFS discover times@>
  g1_vertices.clear();
  for (tie(i1, i1_end) = vertices(g1); i1 != i1_end; ++i1)
    g1_vertices.push_back(*i1);
  permute(g1_vertices.begin(), g1_vertices.end(), perm.begin());
}
@}

We implement the outer-loop of the DFS here, instead of calling the
\code{depth\_first\_search} function, because we want the roots of the
DFS tree's to be ordered by invariant multiplicity. We call
\code{depth\_\-first\_\-visit} to implement the recursive portion of
the DFS. The \code{record\_dfs\_order} adapts the DFS to record
the order in which DFS discovers the vertices.

@d Compute DFS discover times
@{
std::vector<default_color_type> color_vec(num_vertices(g1));
for (typename std::vector<vertex1_t>::iterator ui = g1_vertices.begin();
     ui != g1_vertices.end(); ++ui) {
  if (color_vec[get(index_map1, *ui)] 
      == color_traits<default_color_type>::white()) {
    depth_first_visit
      (g1, *ui, detail::record_dfs_order<Graph1, IndexMap1>(perm, 
                                                       index_map1), 
       make_iterator_property_map(&color_vec[0], index_map1, 
                                  color_vec[0]));
  }
}
@}

\noindent The definition of the \code{record\_dfs\_order} visitor
class is as follows. The index of each vertex is recorded in the
\code{dfs\_order} vector (which is the \code{perm} vector) in the
\code{discover\_vertex} event point.

@d Record DFS ordering visitor
@{
namespace detail {
  template <typename Graph1, typename IndexMap1>
  struct record_dfs_order : public default_dfs_visitor {
    typedef typename graph_traits<Graph1>::vertices_size_type size_type;
    typedef typename graph_traits<Graph1>::vertex_descriptor vertex;

    record_dfs_order(std::vector<size_type>& dfs_order, IndexMap1 index) 
      : dfs_order(dfs_order), index(index) { }

    void discover_vertex(vertex v, const Graph1& g) const {
      dfs_order.push_back(get(index, v));
    }
    std::vector<size_type>& dfs_order; 
    IndexMap1 index;
  };
} // namespace detail
@}


In the MATCH operation, we need to examine all the edges in the set
$E_1[k] - E_1[k-1]$. That is, we need to loop through all the edges of
the form $(k,j)$ or $(j,k)$ where $j \leq k$. To do this efficiently,
we create an array of all the edges in $G_1$ that has been sorted so
that $E_1[k] - E_1[k-1]$ forms a contiguous range.  To each edge
$e=(u,v)$ we assign the number $\max(u,v)$, and then sort the edges by
this number. All the edges $(u,v) \in E_1[k] - E_1[k-1]$ can then be
identified because $\max(u,v) = k$. The following code creates an
array of edges and then sorts them. The \code{edge\_\-ordering\_\-fun}
function object is described next.

@d Order the edges by DFS discover time
@{
typedef typename graph_traits<Graph1>::edge_descriptor edge1_t;
std::vector<edge1_t> edge_set;
std::copy(edges(g1).first, edges(g1).second, std::back_inserter(edge_set));

std::sort(edge_set.begin(), edge_set.end(), 
          detail::edge_ordering
          (make_iterator_property_map(perm.begin(), index_map1, perm[0]), g1));
@}

\noindent The \code{edge\_order} function computes the ordering number
for an edge, which for edge $e=(u,v)$ is $\max(u,v)$. The
\code{edge\_\-ordering\_\-fun} function object simply returns
comparison of two edge's ordering numbers.

@d Isomorph edge ordering predicate
@{
namespace detail {

  template <typename VertexIndexMap, typename Graph>
  std::size_t edge_order(const typename graph_traits<Graph>::edge_descriptor e,
                         VertexIndexMap index_map, const Graph& g) {
    return std::max(get(index_map, source(e, g)), get(index_map, target(e, g)));    
  }

  template <typename VertexIndexMap, typename Graph>
  class edge_ordering_fun {
  public:
    edge_ordering_fun(VertexIndexMap vip, const Graph& g)
      : m_index_map(vip), m_g(g) { }
    template <typename Edge>
    bool operator()(const Edge& e1, const Edge& e2) const {
      return edge_order(e1, m_index_map, m_g) < edge_order(e2, m_index_map, m_g);
    }
    VertexIndexMap m_index_map;
    const Graph& m_g;
  };
  template <class VertexIndexMap, class G>
  inline edge_ordering_fun<VertexIndexMap,G>
  edge_ordering(VertexIndexMap vip, const G& g)
  {
    return edge_ordering_fun<VertexIndexMap,G>(vip, g);
  }
} // namespace detail
@}


We are now ready to enter the main part of the algorithm, the
backtracking search implemented by the \code{isomorph} function (which
corresponds to the ISOMORPH algorithm).  The set $S$ is not
represented directly; instead we represent $V_2 - S$.  Initially $S =
\emptyset$ so $V_2 - S = V_2$.  We use the permuted indices for the
vertices of graph \code{g1}. We represent $V_2 - S$ with a bitset.  We
use \code{std::vector} instead of \code{boost::dyn\_bitset} for speed
instead of space.

@d Invoke recursive \code{isomorph} function
@{
std::vector<char> not_in_S_vec(num_vertices(g2), true);
iterator_property_map<char*, IndexMap2, char, char&>
  not_in_S(&not_in_S_vec[0], index_map2);

return detail::isomorph(g1_vertices.begin(), g1_vertices.end(), 
      edge_set.begin(), edge_set.end(), g1, g2,
      make_iterator_property_map(perm.begin(), index_map1, perm[0]),
      index_map2, f, invar1, invar2, not_in_S);
@}


\subsection{Implementation of ISOMORPH}

The ISOMORPH algorithm is implemented with the \code{isomorph}
function. The vertices of $G_1$ are searched in the order specified by
the iterator range \code{[k\_iter,last)}. The function returns true if
a isomorphism is found between the vertices of $G_1$ in
\code{[k\_iter,last)} and the vertices of $G_2$ in \code{not\_in\_S}.
The mapping is recorded in the parameter \code{f}.

@d Signature for the recursive isomorph function
@{
template <class VertexIter, class EdgeIter, class Graph1, class Graph2,
  class IndexMap1, class IndexMap2, class IndexMapping, 
  class Invar1, class Invar2, class Set>
bool isomorph(VertexIter k_iter, VertexIter last,
              EdgeIter edge_iter, EdgeIter edge_iter_end,
              const Graph1& g1, const Graph2& g2,
              IndexMap1 index_map1,
              IndexMap2 index_map2,
              IndexMapping f, Invar1 invar1, Invar2 invar2,
              const Set& not_in_S)
@}

\noindent The steps for this function are as follows.

@d Body of the isomorph function
@{
{
  @<Some typedefs and variable declarations@>
  @<Return true if matching is complete@>
  @<Create a copy of $f_{k-1}$ which will become $f_k$@>
  @<Compute $M$, the potential matches for $k$@>
  @<Invoke isomorph for each vertex in $M$@>
}
@}

\noindent Here we create short names for some often-used types
and declare some variables.

@d Some typedefs and variable declarations
@{
typedef typename graph_traits<Graph1>::vertex_descriptor vertex1_t;
typedef typename graph_traits<Graph2>::vertex_descriptor vertex2_t;
typedef typename graph_traits<Graph1>::vertices_size_type size_type;

vertex1_t k = *k_iter;
@}

\noindent We have completed creating an isomorphism if \code{k\_iter == last}.

@d Return true if matching is complete
@{
if (k_iter == last) 
  return true;
@}


In the pseudo-code for ISOMORPH, we iterate through each vertex in $v
\in V_2 - S$ and check if $k$ and $v$ can match.  A more efficient
approach is to directly iterate through the potential matches for $k$,
for this often is many fewer vertices than $V_2 - S$. Let $M$ be the
set of potential matches for $k$. $M$ consists of all the vertices $v
\in V_2 - S$ such that if $(k,j)$ or $(j,k) \in E_1[k] - E_1[k-1]$
then $(v,f(j)$ or $(f(j),v) \in E_2$ with $i(v) = i(k)$. Note that
this means if there are no edges in $E_1[k] - E_1[k-1]$ then $M = V_2
- S$. In the case where there are edges in $E_1[k] - E_1[k-1]$ we
break the computation of $M$ into two parts, computing $out$ sets
which are vertices that can match according to an out-edge of $k$, and
computing $in$ sets which are vertices that can match according to an
in-edge of $k$.

The implementation consists of a loop through the edges of $E_1[k] -
E_1[k-1]$. The straightforward implementation would initialize $M
\leftarrow V_2 - S$, and then intersect $M$ with the $out$ or $in$ set
for each edge. However, to reduce the cost of the intersection
operation, we start with $M \leftarrow \emptyset$, and on the first
iteration of the loop we do $M \leftarrow out$ or $M \leftarrow in$
instead of an intersection operation.

@d Compute $M$, the potential matches for $k$
@{
std::vector<vertex2_t> potential_matches;
bool some_edges = false;

for (; edge_iter != edge_iter_end; ++edge_iter) {
  if (get(index_map1, k) != edge_order(*edge_iter, index_map1, g1))
    break;      
  if (k == source(*edge_iter, g1)) { // (k,j)
    @<Compute the $out$ set@>
    if (some_edges == false) {
      @<Perform $M \leftarrow out$@>
    } else {
      @<Perform $M \leftarrow M \intersect out$@>
    }
    some_edges = true;
  } else { // (j,k)
    @<Compute the $in$ set@>
    if (some_edges == false) {
      @<Perform $M \leftarrow in$@>
    } else {
      @<Perform $M \leftarrow M \intersect in$@>
    }
    some_edges = true;
  }
  if (potential_matches.empty())
    break;
} // for edge_iter
if (some_edges == false) {
  @<Perform $M \leftarrow V_2 - S$@>
}
@}

To compute the $out$ set, we iterate through the out-edges $(k,j)$ of
$k$, and for each $j$ we iterate through the in-edges $(v,f(j))$ of
$f(j)$, putting all of the $v$'s in $out$ that have the same vertex
invariant as $k$, and which are in $V_2 - S$. Figure~\ref{fig:out}
depicts the computation of the $out$ set. The implementation is as
follows.

@d Compute the $out$ set
@{
vertex1_t j = target(*edge_iter, g1);
std::vector<vertex2_t> out;
typename graph_traits<Graph2>::in_edge_iterator ei, ei_end;
for (tie(ei, ei_end) = in_edges(get(f, j), g2); ei != ei_end; ++ei) {
  vertex2_t v = source(*ei, g2); // (v,f[j])
  if (invar1[k] == invar2[v] && not_in_S[v])
    out.push_back(v);
}
@}

\noindent Here initialize $M$ with the $out$ set. Since we are
representing sets with sorted vectors, we sort \code{out} before
copying to \code{potential\_matches}.

@d Perform $M \leftarrow out$
@{
indirect_cmp<IndexMap2,std::less<std::size_t> > cmp(index_map2);
std::sort(out.begin(), out.end(), cmp);
std::copy(out.begin(), out.end(), std::back_inserter(potential_matches));
@}

\noindent We use \code{std::set\_intersection} to implement $M
\leftarrow M \intersect out$. Since there is no version of
\code{std::set\_intersection} that works in-place, we create a
temporary for the result and then swap.

@d Perform $M \leftarrow M \intersect out$
@{
indirect_cmp<IndexMap2,std::less<std::size_t> > cmp(index_map2);
std::sort(out.begin(), out.end(), cmp);
std::vector<vertex2_t> tmp_matches;
std::set_intersection(out.begin(), out.end(),
                      potential_matches.begin(), potential_matches.end(),
                      std::back_inserter(tmp_matches), cmp);
std::swap(potential_matches, tmp_matches);
@}

% Shoot, there is some problem with f(j). Could have to do with the
% change from the edge set to just using out_edges and in_edges.
% Yes, have to visit edges in correct order to we don't hit
% part of f that is not yet defined.

\vizfig{out}{Computing the $out$ set.}

@c out.dot
@{
digraph G {
  node[shape=circle]
  size="4,2"