Include is_planar() in topology module

This is a wrapper around BGL's boyer_myrvold_planarity_test.

src/graph/topology/Makefile.am

@@ -23,7 +23,9 @@ libgraph_tool_topology_la_SOURCES = \
     graph_transitive_closure.cc \
     graph_components.cc \
     graph_distance.cc \
-    graph_all_distances.cc
+    graph_all_distances.cc \
+    graph_planar.cc
+
 
 libgraph_tool_topology_la_include_HEADERS = \
     graph_components.hh

src/graph/topology/graph_planar.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007  Tiago de Paula Peixoto <tiago@forked.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#include "graph_filtering.hh"
+#include "graph.hh"
+#include "graph_properties.hh"
+
+#include <boost/graph/boyer_myrvold_planar_test.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+struct get_planar_embedding
+{
+    template <class EdgeMap>
+    class edge_inserter
+    {
+    public:
+        edge_inserter(EdgeMap edge_map): _edge_map(edge_map) {}
+
+        edge_inserter& operator++() { return *this; }
+        edge_inserter& operator++(int) { return *this; }
+        edge_inserter& operator*() { return *this; }
+
+        template <class Key>
+        edge_inserter& operator=(const Key& e)
+        {
+            _edge_map[e] = 1;
+            return *this;
+        }
+
+    private:
+        EdgeMap _edge_map;
+    };
+
+    template <class Graph, class VertexIndex, class EdgeIndex, class EmbedMap,
+              class KurMap>
+    void operator()(Graph& g, VertexIndex vertex_index, EdgeIndex edge_index,
+                    EmbedMap embed_map, KurMap kur_map, bool& is_planar) const
+    {
+        edge_inserter<KurMap> kur_insert(kur_map);
+        unchecked_vector_property_map
+            <vector<typename graph_traits<Graph>::edge_descriptor>, VertexIndex>
+            embedding(vertex_index, num_vertices(g));
+        is_planar = boyer_myrvold_planarity_test
+            (boyer_myrvold_params::graph = g,
+             boyer_myrvold_params::edge_index_map = edge_index,
+             boyer_myrvold_params::embedding = embedding,
+             boyer_myrvold_params::kuratowski_subgraph = kur_insert);
+
+        int i, N = num_vertices(g);
+        #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+        for (i = 0; i < N; ++i)
+        {
+            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
+            if (v == graph_traits<Graph>::null_vertex())
+                continue;
+            embed_map[v].resize(embedding[v].size());
+            for (size_t j = 0; j < embedding[v].size(); ++j)
+                embed_map[v][j] = edge_index[embedding[v][j]];
+        }
+    }
+
+    template <class Graph, class VertexIndex, class EdgeIndex, class KurMap>
+    void operator()(Graph& g, VertexIndex vertex_index, EdgeIndex edge_index,
+                    dummy_property_map embed_map, KurMap kur_map,
+                    bool& is_planar) const
+    {
+        edge_inserter<KurMap> kur_insert(kur_map);
+        is_planar = boyer_myrvold_planarity_test
+            (boyer_myrvold_params::graph = g,
+             boyer_myrvold_params::edge_index_map = edge_index,
+             boyer_myrvold_params::kuratowski_subgraph = kur_insert);
+    }
+
+};
+
+bool is_planar(GraphInterface& gi, boost::any embed_map, boost::any kur_map)
+{
+    bool is_planar;
+
+    if (embed_map.empty())
+        embed_map = dummy_property_map();
+    if (kur_map.empty())
+        kur_map = dummy_property_map();
+
+    typedef mpl::push_back<edge_scalar_properties,
+                           dummy_property_map>::type edge_map_types;
+    typedef mpl::push_back<vertex_scalar_vector_properties,
+                           dummy_property_map>::type vertex_map_types;
+
+    run_action<graph_tool::detail::never_directed>()
+        (gi, bind<void>(get_planar_embedding(), _1, gi.GetVertexIndex(),
+                        gi.GetEdgeIndex(), _2, _3, ref(is_planar)),
+         vertex_map_types(), edge_map_types())
+        (embed_map, kur_map);
+    return is_planar;
+}

src/graph/topology/graph_topology.cc

@@ -31,6 +31,8 @@ void get_prim_spanning_tree(GraphInterface& gi, size_t root,
 void topological_sort(GraphInterface& gi, vector<int32_t>& sort);
 void dominator_tree(GraphInterface& gi, size_t entry, boost::any pred_map);
 void transitive_closure(GraphInterface& gi, GraphInterface& tcgi);
+bool is_planar(GraphInterface& gi, boost::any embed_map, boost::any kur_map);
+
 void export_components();
 void export_dists();
 void export_all_dists();
@@ -43,6 +45,7 @@ BOOST_PYTHON_MODULE(libgraph_tool_topology)
     def("topological_sort", &topological_sort);
     def("dominator_tree", &dominator_tree);
     def("transitive_closure", &transitive_closure);
+    def("is_planar", &is_planar);
     export_components();
     export_dists();
     export_all_dists();

src/graph_tool/topology/__init__.py

@@ -34,6 +34,7 @@ Summary
    transitive_closure
    label_components
    label_biconnected_components
+   is_planar
 
 Contents
 ++++++++
@@ -47,7 +48,7 @@ from .. core import _prop, Vector_int32_t, _check_prop_writable, \
 import random, sys, numpy
 __all__ = ["isomorphism", "min_spanning_tree", "dominator_tree",
            "topological_sort", "transitive_closure", "label_components",
-           "label_biconnected_components", "shortest_distance"]
+           "label_biconnected_components", "shortest_distance", "is_planar"]
 
 def isomorphism(g1, g2, isomap=False):
     """Check whether two graphs are isomorphisms. If `isomap` is True, a vertex
@@ -560,3 +561,100 @@ def shortest_distance(g, source=None, weights=None, max_dist=None,
         if directed != None:
             g.pop_filter(directed=True)
     return dist_map
+
+def is_planar(g, embedding=False, kuratowski=False):
+    """
+    Test if the graph is planar.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    embedding : bool (optional, default: False)
+        If true, return a mapping from vertices to the clockwise order of
+        out-edges in the planar embedding.
+    kuratowski : bool (optional, default: False)
+        If true, the minimal set of edges that form the obstructing Kuratowski
+        subgraph will be returned as a property map, if the graph is not planar.
+
+    Returns
+    -------
+    is_planar : bool
+        Whether or not the graph is planar.
+    embedding : :class:`~graph_tool.PropertyMap` (only if `embedding=True`)
+        A vertex property map with the out-edges indexes in clockwise order in
+        the planar embedding,
+    kuratowski : :class:`~graph_tool.PropertyMap` (only if `kuratowski=True`)
+        An edge property map with the minimal set of edges that form the
+        obstructing Kuratowski subgraph (if the value of kuratowski[e] is 1,
+        the edge belongs to the set)
+
+    Notes
+    -----
+
+    A graph is planar if it can be drawn in two-dimensional space without any of
+    its edges crossing. This algorithm performs the Boyer-Myrvold planarity
+    testing [boyer-myrvold]_. See [boost-planarity]_ for more details.
+
+    This algorithm runs in :math:`O(V)` time.
+
+    Examples
+    --------
+    >>> from numpy.random import seed, random
+    >>> seed(42)
+    >>> g = gt.triangulation(random((100,2)))[0]
+    >>> p, embed_order = gt.is_planar(g, embedding=True)
+    >>> print p
+    True
+    >>> print list(embed_order[g.vertex(0)])
+    [0, 3, 6, 17, 5, 13, 1, 20, 7, 23, 10, 22, 14, 2, 24, 8, 4, 15, 11, 12, 9, 18, 19, 21, 16]
+    >>> g = gt.random_graph(100, lambda: 4, directed=False)
+    >>> p, kur = gt.is_planar(g, kuratowski=True)
+    >>> print p
+    False
+    >>> g.set_edge_filter(kur, True)
+    >>> gt.graph_draw(g, layout="arf",  size=(6,6), output="kuratowski.png")
+    <...>
+
+    .. figure:: kuratowski.png
+        :align: center
+
+        Obstructing Kuratowski subgraph of a random graph.
+
+    References
+    ----------
+    .. [boyer-myrvold] John M. Boyer and Wendy J. Myrvold, "On the Cutting Edge:
+       Simplified O(n) Planarity by Edge Addition Journal of Graph Algorithms
+       and Applications", 8(2): 241-273, 2004.
+    .. [boost-planarity] http://www.boost.org/libs/graph/doc/boyer_myrvold.html
+    """
+
+    g.stash_filter(directed=True)
+    g.set_directed(False)
+
+    if embedding:
+        embed = g.new_vertex_property("vector<int>")
+    else:
+        embed = None
+
+    if kuratowski:
+        kur = g.new_edge_property("bool")
+    else:
+        kur = None
+
+    try:
+        is_planar = libgraph_tool_topology.is_planar(g._Graph__graph,
+                                                     _prop("v", g, embed),
+                                                     _prop("e", g, kur))
+    finally:
+        g.pop_filter(directed=True)
+
+    ret = [is_planar]
+    if embed != None:
+        ret.append(embed)
+    if kur != None:
+        ret.append(kur)
+    if len(ret) == 1:
+        return ret[0]
+    else:
+        return tuple(ret)