/* -*- mode: C++; indent-tabs-mode: nil; -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library.

  *

  * Copyright (C) 2003-2009

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * precise terms see the accompanying LICENSE file.

  *

  * This software is provided "AS IS" with no warranty of any kind,

  * express or implied, and with no claim as to its suitability for any

  * purpose.

  *

  */

 #ifndef LEMON_PLANARITY_H

 #define LEMON_PLANARITY_H

 /// \ingroup planar

 /// \file

 /// \brief Planarity checking, embedding, drawing and coloring

 #include <vector>

 #include <list>

 #include <lemon/dfs.h>

 #include <lemon/bfs.h>

 #include <lemon/radix_sort.h>

 #include <lemon/maps.h>

 #include <lemon/path.h>

 #include <lemon/bucket_heap.h>

 #include <lemon/adaptors.h>

 #include <lemon/edge_set.h>

 #include <lemon/color.h>

 #include <lemon/dim2.h>

 namespace lemon {

   namespace _planarity_bits {

     template <typename Graph>

     struct PlanarityVisitor : DfsVisitor<Graph> {

       TEMPLATE_GRAPH_TYPEDEFS(Graph);

       typedef typename Graph::template NodeMap<Arc> PredMap;

       typedef typename Graph::template EdgeMap<bool> TreeMap;

       typedef typename Graph::template NodeMap<int> OrderMap;

       typedef std::vector<Node> OrderList;

       typedef typename Graph::template NodeMap<int> LowMap;

       typedef typename Graph::template NodeMap<int> AncestorMap;

       PlanarityVisitor(const Graph& graph,

                        PredMap& pred_map, TreeMap& tree_map,

                        OrderMap& order_map, OrderList& order_list,

                        AncestorMap& ancestor_map, LowMap& low_map)

         : _graph(graph), _pred_map(pred_map), _tree_map(tree_map),

           _order_map(order_map), _order_list(order_list),

           _ancestor_map(ancestor_map), _low_map(low_map) {}

       void reach(const Node& node) {

         _order_map[node] = _order_list.size();

         _low_map[node] = _order_list.size();

         _ancestor_map[node] = _order_list.size();

         _order_list.push_back(node);

       }

       void discover(const Arc& arc) {

         Node source = _graph.source(arc);

         Node target = _graph.target(arc);

         _tree_map[arc] = true;

         _pred_map[target] = arc;

       }

       void examine(const Arc& arc) {

         Node source = _graph.source(arc);

         Node target = _graph.target(arc);

         if (_order_map[target] < _order_map[source] && !_tree_map[arc]) {

           if (_low_map[source] > _order_map[target]) {

             _low_map[source] = _order_map[target];

           }

           if (_ancestor_map[source] > _order_map[target]) {

             _ancestor_map[source] = _order_map[target];

           }

         }

       }

       void backtrack(const Arc& arc) {

         Node source = _graph.source(arc);

         Node target = _graph.target(arc);

         if (_low_map[source] > _low_map[target]) {

           _low_map[source] = _low_map[target];

         }

       }

       const Graph& _graph;

       PredMap& _pred_map;

       TreeMap& _tree_map;

       OrderMap& _order_map;

       OrderList& _order_list;

       AncestorMap& _ancestor_map;

       LowMap& _low_map;

     };

     template <typename Graph, bool embedding = true>

     struct NodeDataNode {

       int prev, next;

       int visited;

       typename Graph::Arc first;

       bool inverted;

     };

     template <typename Graph>

     struct NodeDataNode<Graph, false> {

       int prev, next;

       int visited;

     };

     template <typename Graph>

     struct ChildListNode {

       typedef typename Graph::Node Node;

       Node first;

       Node prev, next;

     };

     template <typename Graph>

     struct ArcListNode {

       typename Graph::Arc prev, next;

     };

     template <typename Graph>

     class PlanarityChecking {

     private:

       TEMPLATE_GRAPH_TYPEDEFS(Graph);

       const Graph& _graph;

     private:

       typedef typename Graph::template NodeMap<Arc> PredMap;

       typedef typename Graph::template EdgeMap<bool> TreeMap;

       typedef typename Graph::template NodeMap<int> OrderMap;

       typedef std::vector<Node> OrderList;

       typedef typename Graph::template NodeMap<int> LowMap;

       typedef typename Graph::template NodeMap<int> AncestorMap;

       typedef _planarity_bits::NodeDataNode<Graph> NodeDataNode;

       typedef std::vector<NodeDataNode> NodeData;

       typedef _planarity_bits::ChildListNode<Graph> ChildListNode;

       typedef typename Graph::template NodeMap<ChildListNode> ChildLists;

       typedef typename Graph::template NodeMap<std::list<int> > MergeRoots;

       typedef typename Graph::template NodeMap<bool> EmbedArc;

     public:

       PlanarityChecking(const Graph& graph) : _graph(graph) {}

       bool run() {

         typedef _planarity_bits::PlanarityVisitor<Graph> Visitor;

         PredMap pred_map(_graph, INVALID);

         TreeMap tree_map(_graph, false);

         OrderMap order_map(_graph, -1);

         OrderList order_list;

         AncestorMap ancestor_map(_graph, -1);

         LowMap low_map(_graph, -1);

         Visitor visitor(_graph, pred_map, tree_map,

                         order_map, order_list, ancestor_map, low_map);

         DfsVisit<Graph, Visitor> visit(_graph, visitor);

         visit.run();

         ChildLists child_lists(_graph);

         createChildLists(tree_map, order_map, low_map, child_lists);

         NodeData node_data(2 * order_list.size());

         EmbedArc embed_arc(_graph, false);

         MergeRoots merge_roots(_graph);

         for (int i = order_list.size() - 1; i >= 0; --i) {

           Node node = order_list[i];

           Node source = node;

           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

             Node target = _graph.target(e);

             if (order_map[source] < order_map[target] && tree_map[e]) {

               initFace(target, node_data, order_map, order_list);

             }

           }

           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

             Node target = _graph.target(e);

             if (order_map[source] < order_map[target] && !tree_map[e]) {

               embed_arc[target] = true;

               walkUp(target, source, i, pred_map, low_map,

                      order_map, order_list, node_data, merge_roots);

             }

           }

           for (typename MergeRoots::Value::iterator it =

                  merge_roots[node].begin(); 

                it != merge_roots[node].end(); ++it) {

             int rn = *it;

             walkDown(rn, i, node_data, order_list, child_lists,

                      ancestor_map, low_map, embed_arc, merge_roots);

           }

           merge_roots[node].clear();

           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

             Node target = _graph.target(e);

             if (order_map[source] < order_map[target] && !tree_map[e]) {

               if (embed_arc[target]) {

                 return false;

               }

             }

           }

         }

         return true;

       }

     private:

       void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,

                             const LowMap& low_map, ChildLists& child_lists) {

         for (NodeIt n(_graph); n != INVALID; ++n) {

           Node source = n;

           std::vector<Node> targets;

           for (OutArcIt e(_graph, n); e != INVALID; ++e) {

             Node target = _graph.target(e);

             if (order_map[source] < order_map[target] && tree_map[e]) {

               targets.push_back(target);

             }

           }

           if (targets.size() == 0) {

             child_lists[source].first = INVALID;

           } else if (targets.size() == 1) {

             child_lists[source].first = targets[0];

             child_lists[targets[0]].prev = INVALID;

             child_lists[targets[0]].next = INVALID;

           } else {

             radixSort(targets.begin(), targets.end(), mapToFunctor(low_map));

             for (int i = 1; i < int(targets.size()); ++i) {

               child_lists[targets[i]].prev = targets[i - 1];

               child_lists[targets[i - 1]].next = targets[i];

             }

             child_lists[targets.back()].next = INVALID;

             child_lists[targets.front()].prev = INVALID;

             child_lists[source].first = targets.front();

           }

         }

       }

       void walkUp(const Node& node, Node root, int rorder,

                   const PredMap& pred_map, const LowMap& low_map,

                   const OrderMap& order_map, const OrderList& order_list,

                   NodeData& node_data, MergeRoots& merge_roots) {

         int na, nb;

         bool da, db;

         na = nb = order_map[node];

         da = true; db = false;

         while (true) {

           if (node_data[na].visited == rorder) break;

           if (node_data[nb].visited == rorder) break;

           node_data[na].visited = rorder;

           node_data[nb].visited = rorder;

           int rn = -1;

           if (na >= int(order_list.size())) {

             rn = na;

           } else if (nb >= int(order_list.size())) {

             rn = nb;

           }

           if (rn == -1) {

             int nn;

             nn = da ? node_data[na].prev : node_data[na].next;

             da = node_data[nn].prev != na;

             na = nn;

             nn = db ? node_data[nb].prev : node_data[nb].next;

             db = node_data[nn].prev != nb;

             nb = nn;

           } else {

             Node rep = order_list[rn - order_list.size()];

             Node parent = _graph.source(pred_map[rep]);

             if (low_map[rep] < rorder) {

               merge_roots[parent].push_back(rn);

             } else {

               merge_roots[parent].push_front(rn);

             }

             if (parent != root) {

               na = nb = order_map[parent];

               da = true; db = false;

             } else {

               break;

             }

           }

         }

       }

       void walkDown(int rn, int rorder, NodeData& node_data,

                     OrderList& order_list, ChildLists& child_lists,

                     AncestorMap& ancestor_map, LowMap& low_map,

                     EmbedArc& embed_arc, MergeRoots& merge_roots) {

         std::vector<std::pair<int, bool> > merge_stack;

         for (int di = 0; di < 2; ++di) {

           bool rd = di == 0;

           int pn = rn;

           int n = rd ? node_data[rn].next : node_data[rn].prev;

           while (n != rn) {

             Node node = order_list[n];

             if (embed_arc[node]) {

               // Merging components on the critical path

               while (!merge_stack.empty()) {

                 // Component root

                 int cn = merge_stack.back().first;

                 bool cd = merge_stack.back().second;

                 merge_stack.pop_back();

                 // Parent of component

                 int dn = merge_stack.back().first;

                 bool dd = merge_stack.back().second;

                 merge_stack.pop_back();

                 Node parent = order_list[dn];

                 // Erasing from merge_roots

                 merge_roots[parent].pop_front();

                 Node child = order_list[cn - order_list.size()];

                 // Erasing from child_lists

                 if (child_lists[child].prev != INVALID) {

                   child_lists[child_lists[child].prev].next =

                     child_lists[child].next;

                 } else {

                   child_lists[parent].first = child_lists[child].next;

                 }

                 if (child_lists[child].next != INVALID) {

                   child_lists[child_lists[child].next].prev =

                     child_lists[child].prev;

                 }

                 // Merging external faces

                 {

                   int en = cn;

                   cn = cd ? node_data[cn].prev : node_data[cn].next;

                   cd = node_data[cn].next == en;

                 }

                 if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;

                 if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;

               }

               bool d = pn == node_data[n].prev;

               if (node_data[n].prev == node_data[n].next &&

                   node_data[n].inverted) {

                 d = !d;

               }

               // Embedding arc into external face

               if (rd) node_data[rn].next = n; else node_data[rn].prev = n;

               if (d) node_data[n].prev = rn; else node_data[n].next = rn;

               pn = rn;

               embed_arc[order_list[n]] = false;

             }

             if (!merge_roots[node].empty()) {

               bool d = pn == node_data[n].prev;

               merge_stack.push_back(std::make_pair(n, d));

               int rn = merge_roots[node].front();

               int xn = node_data[rn].next;

               Node xnode = order_list[xn];

               int yn = node_data[rn].prev;

               Node ynode = order_list[yn];

               bool rd;

               if (!external(xnode, rorder, child_lists, 

                             ancestor_map, low_map)) {

                 rd = true;

               } else if (!external(ynode, rorder, child_lists,

                                    ancestor_map, low_map)) {

                 rd = false;

               } else if (pertinent(xnode, embed_arc, merge_roots)) {

                 rd = true;

               } else {

                 rd = false;

               }

               merge_stack.push_back(std::make_pair(rn, rd));

               pn = rn;

               n = rd ? xn : yn;

             } else if (!external(node, rorder, child_lists,

                                  ancestor_map, low_map)) {

               int nn = (node_data[n].next != pn ?

                         node_data[n].next : node_data[n].prev);

               bool nd = n == node_data[nn].prev;

               if (nd) node_data[nn].prev = pn;

               else node_data[nn].next = pn;

               if (n == node_data[pn].prev) node_data[pn].prev = nn;

               else node_data[pn].next = nn;

               node_data[nn].inverted =

                 (node_data[nn].prev == node_data[nn].next && nd != rd);

               n = nn;

             }

             else break;

           }

           if (!merge_stack.empty() || n == rn) {

             break;

           }

         }

       }

       void initFace(const Node& node, NodeData& node_data,

                     const OrderMap& order_map, const OrderList& order_list) {

         int n = order_map[node];

         int rn = n + order_list.size();

         node_data[n].next = node_data[n].prev = rn;

         node_data[rn].next = node_data[rn].prev = n;

         node_data[n].visited = order_list.size();

         node_data[rn].visited = order_list.size();

       }

       bool external(const Node& node, int rorder,

                     ChildLists& child_lists, AncestorMap& ancestor_map,

                     LowMap& low_map) {

         Node child = child_lists[node].first;

         if (child != INVALID) {

           if (low_map[child] < rorder) return true;

         }

         if (ancestor_map[node] < rorder) return true;

         return false;

       }

       bool pertinent(const Node& node, const EmbedArc& embed_arc,

                      const MergeRoots& merge_roots) {

         return !merge_roots[node].empty() || embed_arc[node];

       }

     };

   }

   /// \ingroup planar

   ///

   /// \brief Planarity checking of an undirected simple graph

   ///

   /// This function implements the Boyer-Myrvold algorithm for

   /// planarity checking of an undirected graph. It is a simplified

   /// version of the PlanarEmbedding algorithm class because neither

   /// the embedding nor the kuratowski subdivisons are not computed.

   template <typename GR>

   bool checkPlanarity(const GR& graph) {

     _planarity_bits::PlanarityChecking<GR> pc(graph);

     return pc.run();

   }

   /// \ingroup planar

   ///

   /// \brief Planar embedding of an undirected simple graph

   ///

   /// This class implements the Boyer-Myrvold algorithm for planar

   /// embedding of an undirected graph. The planar embedding is an

   /// ordering of the outgoing edges of the nodes, which is a possible

   /// configuration to draw the graph in the plane. If there is not

   /// such ordering then the graph contains a \f$ K_5 \f$ (full graph

   /// with 5 nodes) or a \f$ K_{3,3} \f$ (complete bipartite graph on

   /// 3 ANode and 3 BNode) subdivision.

   ///

   /// The current implementation calculates either an embedding or a

   /// Kuratowski subdivision. The running time of the algorithm is 

   /// \f$ O(n) \f$.

   template <typename Graph>

   class PlanarEmbedding {

   private:

     TEMPLATE_GRAPH_TYPEDEFS(Graph);

     const Graph& _graph;

     typename Graph::template ArcMap<Arc> _embedding;

     typename Graph::template EdgeMap<bool> _kuratowski;

   private:

     typedef typename Graph::template NodeMap<Arc> PredMap;

     typedef typename Graph::template EdgeMap<bool> TreeMap;

     typedef typename Graph::template NodeMap<int> OrderMap;

     typedef std::vector<Node> OrderList;

     typedef typename Graph::template NodeMap<int> LowMap;

     typedef typename Graph::template NodeMap<int> AncestorMap;

     typedef _planarity_bits::NodeDataNode<Graph> NodeDataNode;

     typedef std::vector<NodeDataNode> NodeData;

     typedef _planarity_bits::ChildListNode<Graph> ChildListNode;

     typedef typename Graph::template NodeMap<ChildListNode> ChildLists;

     typedef typename Graph::template NodeMap<std::list<int> > MergeRoots;

     typedef typename Graph::template NodeMap<Arc> EmbedArc;

     typedef _planarity_bits::ArcListNode<Graph> ArcListNode;

     typedef typename Graph::template ArcMap<ArcListNode> ArcLists;

     typedef typename Graph::template NodeMap<bool> FlipMap;

     typedef typename Graph::template NodeMap<int> TypeMap;

     enum IsolatorNodeType {

       HIGHX = 6, LOWX = 7,

       HIGHY = 8, LOWY = 9,

       ROOT = 10, PERTINENT = 11,

       INTERNAL = 12

     };

   public:

     /// \brief The map for store of embedding

     typedef typename Graph::template ArcMap<Arc> EmbeddingMap;

     /// \brief Constructor

     ///

     /// \note The graph should be simple, i.e. parallel and loop arc

     /// free.

     PlanarEmbedding(const Graph& graph)

       : _graph(graph), _embedding(_graph), _kuratowski(graph, false) {}

     /// \brief Runs the algorithm.

     ///

     /// Runs the algorithm.

     /// \param kuratowski If the parameter is false, then the

     /// algorithm does not compute a Kuratowski subdivision.

     ///\return %True when the graph is planar.

     bool run(bool kuratowski = true) {

       typedef _planarity_bits::PlanarityVisitor<Graph> Visitor;

       PredMap pred_map(_graph, INVALID);

       TreeMap tree_map(_graph, false);

       OrderMap order_map(_graph, -1);

       OrderList order_list;

       AncestorMap ancestor_map(_graph, -1);

       LowMap low_map(_graph, -1);

       Visitor visitor(_graph, pred_map, tree_map,

                       order_map, order_list, ancestor_map, low_map);

       DfsVisit<Graph, Visitor> visit(_graph, visitor);

       visit.run();

       ChildLists child_lists(_graph);

       createChildLists(tree_map, order_map, low_map, child_lists);

       NodeData node_data(2 * order_list.size());

       EmbedArc embed_arc(_graph, INVALID);

       MergeRoots merge_roots(_graph);

       ArcLists arc_lists(_graph);

       FlipMap flip_map(_graph, false);

       for (int i = order_list.size() - 1; i >= 0; --i) {

         Node node = order_list[i];

         node_data[i].first = INVALID;

         Node source = node;

         for (OutArcIt e(_graph, node); e != INVALID; ++e) {

           Node target = _graph.target(e);

           if (order_map[source] < order_map[target] && tree_map[e]) {

             initFace(target, arc_lists, node_data,

                      pred_map, order_map, order_list);

           }

         }

         for (OutArcIt e(_graph, node); e != INVALID; ++e) {

           Node target = _graph.target(e);

           if (order_map[source] < order_map[target] && !tree_map[e]) {

             embed_arc[target] = e;

             walkUp(target, source, i, pred_map, low_map,

                    order_map, order_list, node_data, merge_roots);

           }

         }

         for (typename MergeRoots::Value::iterator it =

                merge_roots[node].begin(); it != merge_roots[node].end(); ++it) {

           int rn = *it;

           walkDown(rn, i, node_data, arc_lists, flip_map, order_list,

                    child_lists, ancestor_map, low_map, embed_arc, merge_roots);

         }

         merge_roots[node].clear();

         for (OutArcIt e(_graph, node); e != INVALID; ++e) {

           Node target = _graph.target(e);

           if (order_map[source] < order_map[target] && !tree_map[e]) {

             if (embed_arc[target] != INVALID) {

               if (kuratowski) {

                 isolateKuratowski(e, node_data, arc_lists, flip_map,

                                   order_map, order_list, pred_map, child_lists,

                                   ancestor_map, low_map,

                                   embed_arc, merge_roots);

               }

               return false;

             }

           }

         }

       }

       for (int i = 0; i < int(order_list.size()); ++i) {

         mergeRemainingFaces(order_list[i], node_data, order_list, order_map,

                             child_lists, arc_lists);

         storeEmbedding(order_list[i], node_data, order_map, pred_map,

                        arc_lists, flip_map);

       }

       return true;

     }

     /// \brief Gives back the successor of an arc

     ///

     /// Gives back the successor of an arc. This function makes

     /// possible to query the cyclic order of the outgoing arcs from

     /// a node.

     Arc next(const Arc& arc) const {

       return _embedding[arc];

     }

     /// \brief Gives back the calculated embedding map

     ///

     /// The returned map contains the successor of each arc in the

     /// graph.

     const EmbeddingMap& embeddingMap() const {

       return _embedding;

     }

     /// \brief Gives back true if the undirected arc is in the

     /// kuratowski subdivision

     ///

     /// Gives back true if the undirected arc is in the kuratowski

     /// subdivision

     /// \note The \c run() had to be called with true value.

     bool kuratowski(const Edge& edge) {

       return _kuratowski[edge];

     }

   private:

     void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,

                           const LowMap& low_map, ChildLists& child_lists) {

       for (NodeIt n(_graph); n != INVALID; ++n) {

         Node source = n;

         std::vector<Node> targets;

         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

           Node target = _graph.target(e);

           if (order_map[source] < order_map[target] && tree_map[e]) {

             targets.push_back(target);

           }

         }

         if (targets.size() == 0) {

           child_lists[source].first = INVALID;

         } else if (targets.size() == 1) {

           child_lists[source].first = targets[0];

           child_lists[targets[0]].prev = INVALID;

           child_lists[targets[0]].next = INVALID;

         } else {

           radixSort(targets.begin(), targets.end(), mapToFunctor(low_map));

           for (int i = 1; i < int(targets.size()); ++i) {

             child_lists[targets[i]].prev = targets[i - 1];

             child_lists[targets[i - 1]].next = targets[i];

           }

           child_lists[targets.back()].next = INVALID;

           child_lists[targets.front()].prev = INVALID;

           child_lists[source].first = targets.front();

         }

       }

     }

     void walkUp(const Node& node, Node root, int rorder,

                 const PredMap& pred_map, const LowMap& low_map,

                 const OrderMap& order_map, const OrderList& order_list,

                 NodeData& node_data, MergeRoots& merge_roots) {

       int na, nb;

       bool da, db;

       na = nb = order_map[node];

       da = true; db = false;

       while (true) {

         if (node_data[na].visited == rorder) break;

         if (node_data[nb].visited == rorder) break;

         node_data[na].visited = rorder;

         node_data[nb].visited = rorder;

         int rn = -1;

         if (na >= int(order_list.size())) {

           rn = na;

         } else if (nb >= int(order_list.size())) {

           rn = nb;

         }

         if (rn == -1) {

           int nn;

           nn = da ? node_data[na].prev : node_data[na].next;

           da = node_data[nn].prev != na;

           na = nn;

           nn = db ? node_data[nb].prev : node_data[nb].next;

           db = node_data[nn].prev != nb;

           nb = nn;

         } else {

           Node rep = order_list[rn - order_list.size()];

           Node parent = _graph.source(pred_map[rep]);

           if (low_map[rep] < rorder) {

             merge_roots[parent].push_back(rn);

           } else {

             merge_roots[parent].push_front(rn);

           }

           if (parent != root) {

             na = nb = order_map[parent];

             da = true; db = false;

           } else {

             break;

           }

         }

       }

     }

     void walkDown(int rn, int rorder, NodeData& node_data,

                   ArcLists& arc_lists, FlipMap& flip_map,

                   OrderList& order_list, ChildLists& child_lists,

                   AncestorMap& ancestor_map, LowMap& low_map,

                   EmbedArc& embed_arc, MergeRoots& merge_roots) {

       std::vector<std::pair<int, bool> > merge_stack;

       for (int di = 0; di < 2; ++di) {

         bool rd = di == 0;

         int pn = rn;

         int n = rd ? node_data[rn].next : node_data[rn].prev;

         while (n != rn) {

           Node node = order_list[n];

           if (embed_arc[node] != INVALID) {

             // Merging components on the critical path

             while (!merge_stack.empty()) {

               // Component root

               int cn = merge_stack.back().first;

               bool cd = merge_stack.back().second;

               merge_stack.pop_back();

               // Parent of component

               int dn = merge_stack.back().first;

               bool dd = merge_stack.back().second;

               merge_stack.pop_back();

               Node parent = order_list[dn];

               // Erasing from merge_roots

               merge_roots[parent].pop_front();

               Node child = order_list[cn - order_list.size()];

               // Erasing from child_lists

               if (child_lists[child].prev != INVALID) {

                 child_lists[child_lists[child].prev].next =

                   child_lists[child].next;

               } else {

                 child_lists[parent].first = child_lists[child].next;

               }

               if (child_lists[child].next != INVALID) {

                 child_lists[child_lists[child].next].prev =

                   child_lists[child].prev;

               }

               // Merging arcs + flipping

               Arc de = node_data[dn].first;

               Arc ce = node_data[cn].first;

               flip_map[order_list[cn - order_list.size()]] = cd != dd;

               if (cd != dd) {

                 std::swap(arc_lists[ce].prev, arc_lists[ce].next);

                 ce = arc_lists[ce].prev;

                 std::swap(arc_lists[ce].prev, arc_lists[ce].next);

               }

               {

                 Arc dne = arc_lists[de].next;

                 Arc cne = arc_lists[ce].next;

                 arc_lists[de].next = cne;

                 arc_lists[ce].next = dne;

                 arc_lists[dne].prev = ce;

                 arc_lists[cne].prev = de;

               }

               if (dd) {

                 node_data[dn].first = ce;

               }

               // Merging external faces

               {

                 int en = cn;

                 cn = cd ? node_data[cn].prev : node_data[cn].next;

                 cd = node_data[cn].next == en;

                  if (node_data[cn].prev == node_data[cn].next &&

                     node_data[cn].inverted) {

                    cd = !cd;

                  }

               }

               if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;

               if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;

             }

             bool d = pn == node_data[n].prev;

             if (node_data[n].prev == node_data[n].next &&

                 node_data[n].inverted) {

               d = !d;

             }

             // Add new arc

             {

               Arc arc = embed_arc[node];

               Arc re = node_data[rn].first;

               arc_lists[arc_lists[re].next].prev = arc;

               arc_lists[arc].next = arc_lists[re].next;

               arc_lists[arc].prev = re;

               arc_lists[re].next = arc;

               if (!rd) {

                 node_data[rn].first = arc;

               }

               Arc rev = _graph.oppositeArc(arc);

               Arc e = node_data[n].first;

               arc_lists[arc_lists[e].next].prev = rev;

               arc_lists[rev].next = arc_lists[e].next;

               arc_lists[rev].prev = e;

               arc_lists[e].next = rev;

               if (d) {

                 node_data[n].first = rev;

               }

             }

             // Embedding arc into external face

             if (rd) node_data[rn].next = n; else node_data[rn].prev = n;

             if (d) node_data[n].prev = rn; else node_data[n].next = rn;

             pn = rn;

             embed_arc[order_list[n]] = INVALID;

           }

           if (!merge_roots[node].empty()) {

             bool d = pn == node_data[n].prev;

             if (node_data[n].prev == node_data[n].next &&

                 node_data[n].inverted) {

               d = !d;

             }

             merge_stack.push_back(std::make_pair(n, d));

             int rn = merge_roots[node].front();

             int xn = node_data[rn].next;

             Node xnode = order_list[xn];

             int yn = node_data[rn].prev;

             Node ynode = order_list[yn];

             bool rd;

             if (!external(xnode, rorder, child_lists, ancestor_map, low_map)) {

               rd = true;

             } else if (!external(ynode, rorder, child_lists,

                                  ancestor_map, low_map)) {

               rd = false;

             } else if (pertinent(xnode, embed_arc, merge_roots)) {

               rd = true;

             } else {

               rd = false;

             }

             merge_stack.push_back(std::make_pair(rn, rd));

             pn = rn;

             n = rd ? xn : yn;

           } else if (!external(node, rorder, child_lists,

                                ancestor_map, low_map)) {

             int nn = (node_data[n].next != pn ?

                       node_data[n].next : node_data[n].prev);

             bool nd = n == node_data[nn].prev;

             if (nd) node_data[nn].prev = pn;

             else node_data[nn].next = pn;

             if (n == node_data[pn].prev) node_data[pn].prev = nn;

             else node_data[pn].next = nn;

             node_data[nn].inverted =

               (node_data[nn].prev == node_data[nn].next && nd != rd);

             n = nn;

           }

           else break;

         }

         if (!merge_stack.empty() || n == rn) {

           break;

         }

       }

     }

     void initFace(const Node& node, ArcLists& arc_lists,

                   NodeData& node_data, const PredMap& pred_map,

                   const OrderMap& order_map, const OrderList& order_list) {

       int n = order_map[node];

       int rn = n + order_list.size();

       node_data[n].next = node_data[n].prev = rn;

       node_data[rn].next = node_data[rn].prev = n;

       node_data[n].visited = order_list.size();

       node_data[rn].visited = order_list.size();

       node_data[n].inverted = false;

       node_data[rn].inverted = false;

       Arc arc = pred_map[node];

       Arc rev = _graph.oppositeArc(arc);

       node_data[rn].first = arc;

       node_data[n].first = rev;

       arc_lists[arc].prev = arc;

       arc_lists[arc].next = arc;

       arc_lists[rev].prev = rev;

       arc_lists[rev].next = rev;

     }

     void mergeRemainingFaces(const Node& node, NodeData& node_data,

                              OrderList& order_list, OrderMap& order_map,

                              ChildLists& child_lists, ArcLists& arc_lists) {

       while (child_lists[node].first != INVALID) {

         int dd = order_map[node];

         Node child = child_lists[node].first;

         int cd = order_map[child] + order_list.size();

         child_lists[node].first = child_lists[child].next;

         Arc de = node_data[dd].first;

         Arc ce = node_data[cd].first;

         if (de != INVALID) {

           Arc dne = arc_lists[de].next;

           Arc cne = arc_lists[ce].next;

           arc_lists[de].next = cne;

           arc_lists[ce].next = dne;

           arc_lists[dne].prev = ce;

           arc_lists[cne].prev = de;

         }

         node_data[dd].first = ce;

       }

     }

     void storeEmbedding(const Node& node, NodeData& node_data,

                         OrderMap& order_map, PredMap& pred_map,

                         ArcLists& arc_lists, FlipMap& flip_map) {

       if (node_data[order_map[node]].first == INVALID) return;

       if (pred_map[node] != INVALID) {

         Node source = _graph.source(pred_map[node]);

         flip_map[node] = flip_map[node] != flip_map[source];

       }

       Arc first = node_data[order_map[node]].first;

       Arc prev = first;

       Arc arc = flip_map[node] ?

         arc_lists[prev].prev : arc_lists[prev].next;

       _embedding[prev] = arc;

       while (arc != first) {

         Arc next = arc_lists[arc].prev == prev ?

           arc_lists[arc].next : arc_lists[arc].prev;

         prev = arc; arc = next;

         _embedding[prev] = arc;

       }

     }

     bool external(const Node& node, int rorder,

                   ChildLists& child_lists, AncestorMap& ancestor_map,

                   LowMap& low_map) {

       Node child = child_lists[node].first;

       if (child != INVALID) {

         if (low_map[child] < rorder) return true;

       }

       if (ancestor_map[node] < rorder) return true;

       return false;

     }

     bool pertinent(const Node& node, const EmbedArc& embed_arc,

                    const MergeRoots& merge_roots) {

       return !merge_roots[node].empty() || embed_arc[node] != INVALID;

     }

     int lowPoint(const Node& node, OrderMap& order_map, ChildLists& child_lists,

                  AncestorMap& ancestor_map, LowMap& low_map) {

       int low_point;

       Node child = child_lists[node].first;

       if (child != INVALID) {

         low_point = low_map[child];

       } else {

         low_point = order_map[node];

       }

       if (low_point > ancestor_map[node]) {

         low_point = ancestor_map[node];

       }

       return low_point;

     }

     int findComponentRoot(Node root, Node node, ChildLists& child_lists,

                           OrderMap& order_map, OrderList& order_list) {

       int order = order_map[root];

       int norder = order_map[node];

       Node child = child_lists[root].first;

       while (child != INVALID) {

         int corder = order_map[child];

         if (corder > order && corder < norder) {

           order = corder;

         }

         child = child_lists[child].next;

       }

       return order + order_list.size();

     }

     Node findPertinent(Node node, OrderMap& order_map, NodeData& node_data,

                        EmbedArc& embed_arc, MergeRoots& merge_roots) {

       Node wnode =_graph.target(node_data[order_map[node]].first);

       while (!pertinent(wnode, embed_arc, merge_roots)) {

         wnode = _graph.target(node_data[order_map[wnode]].first);

       }

       return wnode;

     }

     Node findExternal(Node node, int rorder, OrderMap& order_map,

                       ChildLists& child_lists, AncestorMap& ancestor_map,

                       LowMap& low_map, NodeData& node_data) {

       Node wnode =_graph.target(node_data[order_map[node]].first);

       while (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {

         wnode = _graph.target(node_data[order_map[wnode]].first);

       }

       return wnode;

     }

     void markCommonPath(Node node, int rorder, Node& wnode, Node& znode,

                         OrderList& order_list, OrderMap& order_map,

                         NodeData& node_data, ArcLists& arc_lists,

                         EmbedArc& embed_arc, MergeRoots& merge_roots,

                         ChildLists& child_lists, AncestorMap& ancestor_map,

                         LowMap& low_map) {

       Node cnode = node;

       Node pred = INVALID;

       while (true) {

         bool pert = pertinent(cnode, embed_arc, merge_roots);

         bool ext = external(cnode, rorder, child_lists, ancestor_map, low_map);

         if (pert && ext) {

           if (!merge_roots[cnode].empty()) {

             int cn = merge_roots[cnode].back();

             if (low_map[order_list[cn - order_list.size()]] < rorder) {

               Arc arc = node_data[cn].first;

               _kuratowski.set(arc, true);

               pred = cnode;

               cnode = _graph.target(arc);

               continue;

             }

           }

           wnode = znode = cnode;

           return;

         } else if (pert) {

           wnode = cnode;

           while (!external(cnode, rorder, child_lists, ancestor_map, low_map)) {

             Arc arc = node_data[order_map[cnode]].first;

             if (_graph.target(arc) == pred) {

               arc = arc_lists[arc].next;

             }

             _kuratowski.set(arc, true);

             Node next = _graph.target(arc);

             pred = cnode; cnode = next;

           }

           znode = cnode;

           return;

         } else if (ext) {

           znode = cnode;

           while (!pertinent(cnode, embed_arc, merge_roots)) {

             Arc arc = node_data[order_map[cnode]].first;

             if (_graph.target(arc) == pred) {

               arc = arc_lists[arc].next;

             }

             _kuratowski.set(arc, true);

             Node next = _graph.target(arc);

             pred = cnode; cnode = next;

           }

           wnode = cnode;

           return;

         } else {

           Arc arc = node_data[order_map[cnode]].first;

           if (_graph.target(arc) == pred) {

             arc = arc_lists[arc].next;

           }

           _kuratowski.set(arc, true);

           Node next = _graph.target(arc);

           pred = cnode; cnode = next;

         }

       }

     }

     void orientComponent(Node root, int rn, OrderMap& order_map,

                          PredMap& pred_map, NodeData& node_data,

                          ArcLists& arc_lists, FlipMap& flip_map,

                          TypeMap& type_map) {

       node_data[order_map[root]].first = node_data[rn].first;

       type_map[root] = 1;

       std::vector<Node> st, qu;

       st.push_back(root);

       while (!st.empty()) {

         Node node = st.back();

         st.pop_back();

         qu.push_back(node);

         Arc arc = node_data[order_map[node]].first;

         if (type_map[_graph.target(arc)] == 0) {

           st.push_back(_graph.target(arc));

           type_map[_graph.target(arc)] = 1;

         }

         Arc last = arc, pred = arc;

         arc = arc_lists[arc].next;

         while (arc != last) {

           if (type_map[_graph.target(arc)] == 0) {

             st.push_back(_graph.target(arc));

             type_map[_graph.target(arc)] = 1;

           }

           Arc next = arc_lists[arc].next != pred ?

             arc_lists[arc].next : arc_lists[arc].prev;

           pred = arc; arc = next;

         }

       }

       type_map[root] = 2;

       flip_map[root] = false;

       for (int i = 1; i < int(qu.size()); ++i) {

         Node node = qu[i];

         while (type_map[node] != 2) {

           st.push_back(node);

           type_map[node] = 2;

           node = _graph.source(pred_map[node]);

         }

         bool flip = flip_map[node];

         while (!st.empty()) {

           node = st.back();

           st.pop_back();

           flip_map[node] = flip != flip_map[node];

           flip = flip_map[node];

           if (flip) {

             Arc arc = node_data[order_map[node]].first;

             std::swap(arc_lists[arc].prev, arc_lists[arc].next);

             arc = arc_lists[arc].prev;

             std::swap(arc_lists[arc].prev, arc_lists[arc].next);

             node_data[order_map[node]].first = arc;

           }

         }

       }

       for (int i = 0; i < int(qu.size()); ++i) {

         Arc arc = node_data[order_map[qu[i]]].first;

         Arc last = arc, pred = arc;

         arc = arc_lists[arc].next;

         while (arc != last) {

           if (arc_lists[arc].next == pred) {

             std::swap(arc_lists[arc].next, arc_lists[arc].prev);

           }

           pred = arc; arc = arc_lists[arc].next;

         }

       }

     }

     void setFaceFlags(Node root, Node wnode, Node ynode, Node xnode,

                       OrderMap& order_map, NodeData& node_data,

                       TypeMap& type_map) {

       Node node = _graph.target(node_data[order_map[root]].first);

       while (node != ynode) {

         type_map[node] = HIGHY;

         node = _graph.target(node_data[order_map[node]].first);

       }

       while (node != wnode) {

         type_map[node] = LOWY;

         node = _graph.target(node_data[order_map[node]].first);

       }

       node = _graph.target(node_data[order_map[wnode]].first);

       while (node != xnode) {

         type_map[node] = LOWX;

         node = _graph.target(node_data[order_map[node]].first);

       }

       type_map[node] = LOWX;

       node = _graph.target(node_data[order_map[xnode]].first);

       while (node != root) {

         type_map[node] = HIGHX;

         node = _graph.target(node_data[order_map[node]].first);

       }

       type_map[wnode] = PERTINENT;

       type_map[root] = ROOT;

     }

     void findInternalPath(std::vector<Arc>& ipath,

                           Node wnode, Node root, TypeMap& type_map,

                           OrderMap& order_map, NodeData& node_data,

                           ArcLists& arc_lists) {

       std::vector<Arc> st;

       Node node = wnode;

       while (node != root) {

         Arc arc = arc_lists[node_data[order_map[node]].first].next;

         st.push_back(arc);

         node = _graph.target(arc);

       }

       while (true) {

         Arc arc = st.back();

         if (type_map[_graph.target(arc)] == LOWX ||

             type_map[_graph.target(arc)] == HIGHX) {

           break;

         }

         if (type_map[_graph.target(arc)] == 2) {

           type_map[_graph.target(arc)] = 3;

           arc = arc_lists[_graph.oppositeArc(arc)].next;

           st.push_back(arc);

         } else {

           st.pop_back();

           arc = arc_lists[arc].next;

           while (_graph.oppositeArc(arc) == st.back()) {

             arc = st.back();

             st.pop_back();

             arc = arc_lists[arc].next;

           }

           st.push_back(arc);

         }

       }

       for (int i = 0; i < int(st.size()); ++i) {

         if (type_map[_graph.target(st[i])] != LOWY &&

             type_map[_graph.target(st[i])] != HIGHY) {

           for (; i < int(st.size()); ++i) {

             ipath.push_back(st[i]);

           }

         }

       }

     }

     void setInternalFlags(std::vector<Arc>& ipath, TypeMap& type_map) {

       for (int i = 1; i < int(ipath.size()); ++i) {

         type_map[_graph.source(ipath[i])] = INTERNAL;

       }

     }

     void findPilePath(std::vector<Arc>& ppath,

                       Node root, TypeMap& type_map, OrderMap& order_map,

                       NodeData& node_data, ArcLists& arc_lists) {

       std::vector<Arc> st;

       st.push_back(_graph.oppositeArc(node_data[order_map[root]].first));

       st.push_back(node_data[order_map[root]].first);

       while (st.size() > 1) {

         Arc arc = st.back();

         if (type_map[_graph.target(arc)] == INTERNAL) {

           break;

         }

         if (type_map[_graph.target(arc)] == 3) {

           type_map[_graph.target(arc)] = 4;

           arc = arc_lists[_graph.oppositeArc(arc)].next;

           st.push_back(arc);

         } else {

           st.pop_back();

           arc = arc_lists[arc].next;

           while (!st.empty() && _graph.oppositeArc(arc) == st.back()) {

             arc = st.back();

             st.pop_back();

             arc = arc_lists[arc].next;

           }

           st.push_back(arc);

         }

       }

       for (int i = 1; i < int(st.size()); ++i) {

         ppath.push_back(st[i]);

       }

     }

     int markExternalPath(Node node, OrderMap& order_map,

                          ChildLists& child_lists, PredMap& pred_map,

                          AncestorMap& ancestor_map, LowMap& low_map) {

       int lp = lowPoint(node, order_map, child_lists,

                         ancestor_map, low_map);

       if (ancestor_map[node] != lp) {

         node = child_lists[node].first;

         _kuratowski[pred_map[node]] = true;

         while (ancestor_map[node] != lp) {

           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

             Node tnode = _graph.target(e);

             if (order_map[tnode] > order_map[node] && low_map[tnode] == lp) {

               node = tnode;

               _kuratowski[e] = true;

               break;

             }

           }

         }

       }

       for (OutArcIt e(_graph, node); e != INVALID; ++e) {

         if (order_map[_graph.target(e)] == lp) {

           _kuratowski[e] = true;

           break;

         }

       }

       return lp;

     }

     void markPertinentPath(Node node, OrderMap& order_map,

                            NodeData& node_data, ArcLists& arc_lists,

                            EmbedArc& embed_arc, MergeRoots& merge_roots) {

       while (embed_arc[node] == INVALID) {

         int n = merge_roots[node].front();

         Arc arc = node_data[n].first;

         _kuratowski.set(arc, true);

         Node pred = node;

         node = _graph.target(arc);

         while (!pertinent(node, embed_arc, merge_roots)) {

           arc = node_data[order_map[node]].first;

           if (_graph.target(arc) == pred) {

             arc = arc_lists[arc].next;

           }

           _kuratowski.set(arc, true);

           pred = node;

           node = _graph.target(arc);

         }

       }

       _kuratowski.set(embed_arc[node], true);

     }

     void markPredPath(Node node, Node snode, PredMap& pred_map) {

       while (node != snode) {

         _kuratowski.set(pred_map[node], true);

         node = _graph.source(pred_map[node]);

       }

     }

     void markFacePath(Node ynode, Node xnode,

                       OrderMap& order_map, NodeData& node_data) {

       Arc arc = node_data[order_map[ynode]].first;

       Node node = _graph.target(arc);

       _kuratowski.set(arc, true);

       while (node != xnode) {

         arc = node_data[order_map[node]].first;

         _kuratowski.set(arc, true);

         node = _graph.target(arc);

       }

     }

     void markInternalPath(std::vector<Arc>& path) {

       for (int i = 0; i < int(path.size()); ++i) {

         _kuratowski.set(path[i], true);

       }

     }

     void markPilePath(std::vector<Arc>& path) {

       for (int i = 0; i < int(path.size()); ++i) {

         _kuratowski.set(path[i], true);

       }

     }

     void isolateKuratowski(Arc arc, NodeData& node_data,

                            ArcLists& arc_lists, FlipMap& flip_map,

                            OrderMap& order_map, OrderList& order_list,

                            PredMap& pred_map, ChildLists& child_lists,

                            AncestorMap& ancestor_map, LowMap& low_map,

                            EmbedArc& embed_arc, MergeRoots& merge_roots) {

       Node root = _graph.source(arc);

       Node enode = _graph.target(arc);

       int rorder = order_map[root];

       TypeMap type_map(_graph, 0);

       int rn = findComponentRoot(root, enode, child_lists,

                                  order_map, order_list);

       Node xnode = order_list[node_data[rn].next];

       Node ynode = order_list[node_data[rn].prev];

       // Minor-A

       {

         while (!merge_roots[xnode].empty() || !merge_roots[ynode].empty()) {

           if (!merge_roots[xnode].empty()) {

             root = xnode;

             rn = merge_roots[xnode].front();

           } else {

             root = ynode;

             rn = merge_roots[ynode].front();

           }

           xnode = order_list[node_data[rn].next];

           ynode = order_list[node_data[rn].prev];

         }

         if (root != _graph.source(arc)) {

           orientComponent(root, rn, order_map, pred_map,

                           node_data, arc_lists, flip_map, type_map);

           markFacePath(root, root, order_map, node_data);

           int xlp = markExternalPath(xnode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           int ylp = markExternalPath(ynode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

           Node lwnode = findPertinent(ynode, order_map, node_data,

                                       embed_arc, merge_roots);

           markPertinentPath(lwnode, order_map, node_data, arc_lists,

                             embed_arc, merge_roots);

           return;

         }

       }

       orientComponent(root, rn, order_map, pred_map,

                       node_data, arc_lists, flip_map, type_map);

       Node wnode = findPertinent(ynode, order_map, node_data,

                                  embed_arc, merge_roots);

       setFaceFlags(root, wnode, ynode, xnode, order_map, node_data, type_map);

       //Minor-B

       if (!merge_roots[wnode].empty()) {

         int cn = merge_roots[wnode].back();

         Node rep = order_list[cn - order_list.size()];

         if (low_map[rep] < rorder) {

           markFacePath(root, root, order_map, node_data);

           int xlp = markExternalPath(xnode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           int ylp = markExternalPath(ynode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           Node lwnode, lznode;

           markCommonPath(wnode, rorder, lwnode, lznode, order_list,

                          order_map, node_data, arc_lists, embed_arc,

                          merge_roots, child_lists, ancestor_map, low_map);

           markPertinentPath(lwnode, order_map, node_data, arc_lists,

                             embed_arc, merge_roots);

           int zlp = markExternalPath(lznode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           int minlp = xlp < ylp ? xlp : ylp;

           if (zlp < minlp) minlp = zlp;

           int maxlp = xlp > ylp ? xlp : ylp;

           if (zlp > maxlp) maxlp = zlp;

           markPredPath(order_list[maxlp], order_list[minlp], pred_map);

           return;

         }

       }

       Node pxnode, pynode;

       std::vector<Arc> ipath;

       findInternalPath(ipath, wnode, root, type_map, order_map,

                        node_data, arc_lists);

       setInternalFlags(ipath, type_map);

       pynode = _graph.source(ipath.front());

       pxnode = _graph.target(ipath.back());

       wnode = findPertinent(pynode, order_map, node_data,

                             embed_arc, merge_roots);

       // Minor-C

       {

         if (type_map[_graph.source(ipath.front())] == HIGHY) {

           if (type_map[_graph.target(ipath.back())] == HIGHX) {

             markFacePath(xnode, pxnode, order_map, node_data);

           }

           markFacePath(root, xnode, order_map, node_data);

           markPertinentPath(wnode, order_map, node_data, arc_lists,

                             embed_arc, merge_roots);

           markInternalPath(ipath);

           int xlp = markExternalPath(xnode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           int ylp = markExternalPath(ynode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

           return;

         }

         if (type_map[_graph.target(ipath.back())] == HIGHX) {

           markFacePath(ynode, root, order_map, node_data);

           markPertinentPath(wnode, order_map, node_data, arc_lists,

                             embed_arc, merge_roots);

           markInternalPath(ipath);

           int xlp = markExternalPath(xnode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           int ylp = markExternalPath(ynode, order_map, child_lists,

                                      pred_map, ancestor_map, low_map);

           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

           return;

         }

       }

       std::vector<Arc> ppath;

       findPilePath(ppath, root, type_map, order_map, node_data, arc_lists);

       // Minor-D

       if (!ppath.empty()) {

         markFacePath(ynode, xnode, order_map, node_data);

         markPertinentPath(wnode, order_map, node_data, arc_lists,

                           embed_arc, merge_roots);

         markPilePath(ppath);

         markInternalPath(ipath);

         int xlp = markExternalPath(xnode, order_map, child_lists,

                                    pred_map, ancestor_map, low_map);

         int ylp = markExternalPath(ynode, order_map, child_lists,

                                    pred_map, ancestor_map, low_map);

         markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

         return;

       }

       // Minor-E*

       {

         if (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {

           Node znode = findExternal(pynode, rorder, order_map,

                                     child_lists, ancestor_map,

                                     low_map, node_data);

           if (type_map[znode] == LOWY) {

             markFacePath(root, xnode, order_map, node_data);

             markPertinentPath(wnode, order_map, node_data, arc_lists,

                               embed_arc, merge_roots);

             markInternalPath(ipath);

             int xlp = markExternalPath(xnode, order_map, child_lists,

                                        pred_map, ancestor_map, low_map);

             int zlp = markExternalPath(znode, order_map, child_lists,

                                        pred_map, ancestor_map, low_map);

             markPredPath(root, order_list[xlp < zlp ? xlp : zlp], pred_map);

           } else {

             markFacePath(ynode, root, order_map, node_data);

             markPertinentPath(wnode, order_map, node_data, arc_lists,

                               embed_arc, merge_roots);

             markInternalPath(ipath);

             int ylp = markExternalPath(ynode, order_map, child_lists,

                                        pred_map, ancestor_map, low_map);

             int zlp = markExternalPath(znode, order_map, child_lists,

                                        pred_map, ancestor_map, low_map);

             markPredPath(root, order_list[ylp < zlp ? ylp : zlp], pred_map);

           }

           return;

         }

         int xlp = markExternalPath(xnode, order_map, child_lists,

                                    pred_map, ancestor_map, low_map);

         int ylp = markExternalPath(ynode, order_map, child_lists,

                                    pred_map, ancestor_map, low_map);

         int wlp = markExternalPath(wnode, order_map, child_lists,

                                    pred_map, ancestor_map, low_map);

         if (wlp > xlp && wlp > ylp) {

           markFacePath(root, root, order_map, node_data);

           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

           return;

         }

         markInternalPath(ipath);

         markPertinentPath(wnode, order_map, node_data, arc_lists,

                           embed_arc, merge_roots);

         if (xlp > ylp && xlp > wlp) {

           markFacePath(root, pynode, order_map, node_data);

           markFacePath(wnode, xnode, order_map, node_data);

           markPredPath(root, order_list[ylp < wlp ? ylp : wlp], pred_map);

           return;

         }

         if (ylp > xlp && ylp > wlp) {

           markFacePath(pxnode, root, order_map, node_data);

           markFacePath(ynode, wnode, order_map, node_data);

           markPredPath(root, order_list[xlp < wlp ? xlp : wlp], pred_map);

           return;

         }

         if (pynode != ynode) {

           markFacePath(pxnode, wnode, order_map, node_data);

           int minlp = xlp < ylp ? xlp : ylp;

           if (wlp < minlp) minlp = wlp;

           int maxlp = xlp > ylp ? xlp : ylp;

           if (wlp > maxlp) maxlp = wlp;

           markPredPath(order_list[maxlp], order_list[minlp], pred_map);

           return;

         }

         if (pxnode != xnode) {

           markFacePath(wnode, pynode, order_map, node_data);

           int minlp = xlp < ylp ? xlp : ylp;

           if (wlp < minlp) minlp = wlp;

           int maxlp = xlp > ylp ? xlp : ylp;

           if (wlp > maxlp) maxlp = wlp;

           markPredPath(order_list[maxlp], order_list[minlp], pred_map);

           return;

         }

         markFacePath(root, root, order_map, node_data);

         int minlp = xlp < ylp ? xlp : ylp;

         if (wlp < minlp) minlp = wlp;

         markPredPath(root, order_list[minlp], pred_map);

         return;

       }

     }

   };

   namespace _planarity_bits {

     template <typename Graph, typename EmbeddingMap>

     void makeConnected(Graph& graph, EmbeddingMap& embedding) {

       DfsVisitor<Graph> null_visitor;

       DfsVisit<Graph, DfsVisitor<Graph> > dfs(graph, null_visitor);

       dfs.init();

       typename Graph::Node u = INVALID;

       for (typename Graph::NodeIt n(graph); n != INVALID; ++n) {

         if (!dfs.reached(n)) {

           dfs.addSource(n);

           dfs.start();

           if (u == INVALID) {

             u = n;

           } else {

             typename Graph::Node v = n;

             typename Graph::Arc ue = typename Graph::OutArcIt(graph, u);

             typename Graph::Arc ve = typename Graph::OutArcIt(graph, v);

             typename Graph::Arc e = graph.direct(graph.addEdge(u, v), true);

             if (ue != INVALID) {

               embedding[e] = embedding[ue];

               embedding[ue] = e;

             } else {

               embedding[e] = e;

             }

             if (ve != INVALID) {

               embedding[graph.oppositeArc(e)] = embedding[ve];

               embedding[ve] = graph.oppositeArc(e);

             } else {

               embedding[graph.oppositeArc(e)] = graph.oppositeArc(e);

             }

           }

         }

       }

     }

     template <typename Graph, typename EmbeddingMap>

     void makeBiNodeConnected(Graph& graph, EmbeddingMap& embedding) {

       typename Graph::template ArcMap<bool> processed(graph);

       std::vector<typename Graph::Arc> arcs;

       for (typename Graph::ArcIt e(graph); e != INVALID; ++e) {

         arcs.push_back(e);

       }

       IterableBoolMap<Graph, typename Graph::Node> visited(graph, false);

       for (int i = 0; i < int(arcs.size()); ++i) {

         typename Graph::Arc pp = arcs[i];

         if (processed[pp]) continue;

         typename Graph::Arc e = embedding[graph.oppositeArc(pp)];

         processed[e] = true;

         visited.set(graph.source(e), true);

         typename Graph::Arc p = e, l = e;

         e = embedding[graph.oppositeArc(e)];

         while (e != l) {

           processed[e] = true;

           if (visited[graph.source(e)]) {

             typename Graph::Arc n =

               graph.direct(graph.addEdge(graph.source(p),

                                            graph.target(e)), true);

             embedding[n] = p;

             embedding[graph.oppositeArc(pp)] = n;

             embedding[graph.oppositeArc(n)] =

               embedding[graph.oppositeArc(e)];

             embedding[graph.oppositeArc(e)] =

               graph.oppositeArc(n);

             p = n;

             e = embedding[graph.oppositeArc(n)];

           } else {

             visited.set(graph.source(e), true);

             pp = p;

             p = e;

             e = embedding[graph.oppositeArc(e)];

           }

         }

         visited.setAll(false);

       }

     }

     template <typename Graph, typename EmbeddingMap>

     void makeMaxPlanar(Graph& graph, EmbeddingMap& embedding) {

       typename Graph::template NodeMap<int> degree(graph);

       for (typename Graph::NodeIt n(graph); n != INVALID; ++n) {

         degree[n] = countIncEdges(graph, n);

       }

       typename Graph::template ArcMap<bool> processed(graph);

       IterableBoolMap<Graph, typename Graph::Node> visited(graph, false);

       std::vector<typename Graph::Arc> arcs;

       for (typename Graph::ArcIt e(graph); e != INVALID; ++e) {

         arcs.push_back(e);

       }

       for (int i = 0; i < int(arcs.size()); ++i) {

         typename Graph::Arc e = arcs[i];

         if (processed[e]) continue;

         processed[e] = true;

         typename Graph::Arc mine = e;

         int mind = degree[graph.source(e)];

         int face_size = 1;

         typename Graph::Arc l = e;

         e = embedding[graph.oppositeArc(e)];

         while (l != e) {

           processed[e] = true;

           ++face_size;

           if (degree[graph.source(e)] < mind) {

             mine = e;

             mind = degree[graph.source(e)];

           }

           e = embedding[graph.oppositeArc(e)];

         }

         if (face_size < 4) {

           continue;

         }

         typename Graph::Node s = graph.source(mine);

         for (typename Graph::OutArcIt e(graph, s); e != INVALID; ++e) {

           visited.set(graph.target(e), true);

         }

         typename Graph::Arc oppe = INVALID;

         e = embedding[graph.oppositeArc(mine)];

         e = embedding[graph.oppositeArc(e)];

         while (graph.target(e) != s) {

           if (visited[graph.source(e)]) {

             oppe = e;

             break;

           }

           e = embedding[graph.oppositeArc(e)];

         }

         visited.setAll(false);

         if (oppe == INVALID) {

           e = embedding[graph.oppositeArc(mine)];

           typename Graph::Arc pn = mine, p = e;

           e = embedding[graph.oppositeArc(e)];

           while (graph.target(e) != s) {

             typename Graph::Arc n =

               graph.direct(graph.addEdge(s, graph.source(e)), true);

             embedding[n] = pn;

             embedding[graph.oppositeArc(n)] = e;

             embedding[graph.oppositeArc(p)] = graph.oppositeArc(n);

             pn = n;

             p = e;

             e = embedding[graph.oppositeArc(e)];

           }

           embedding[graph.oppositeArc(e)] = pn;

         } else {

           mine = embedding[graph.oppositeArc(mine)];

           s = graph.source(mine);

           oppe = embedding[graph.oppositeArc(oppe)];

           typename Graph::Node t = graph.source(oppe);

           typename Graph::Arc ce = graph.direct(graph.addEdge(s, t), true);

           embedding[ce] = mine;

           embedding[graph.oppositeArc(ce)] = oppe;

           typename Graph::Arc pn = ce, p = oppe;

           e = embedding[graph.oppositeArc(oppe)];

           while (graph.target(e) != s) {

             typename Graph::Arc n =

               graph.direct(graph.addEdge(s, graph.source(e)), true);

             embedding[n] = pn;

             embedding[graph.oppositeArc(n)] = e;

             embedding[graph.oppositeArc(p)] = graph.oppositeArc(n);

             pn = n;

             p = e;

             e = embedding[graph.oppositeArc(e)];

           }

           embedding[graph.oppositeArc(e)] = pn;

           pn = graph.oppositeArc(ce), p = mine;

           e = embedding[graph.oppositeArc(mine)];

           while (graph.target(e) != t) {

             typename Graph::Arc n =

               graph.direct(graph.addEdge(t, graph.source(e)), true);

             embedding[n] = pn;

             embedding[graph.oppositeArc(n)] = e;

             embedding[graph.oppositeArc(p)] = graph.oppositeArc(n);

             pn = n;

             p = e;

             e = embedding[graph.oppositeArc(e)];

           }

           embedding[graph.oppositeArc(e)] = pn;

         }

       }

     }

   }

   /// \ingroup planar

   ///

   /// \brief Schnyder's planar drawing algorithm

   ///

   /// The planar drawing algorithm calculates positions for the nodes

   /// in the plane which coordinates satisfy that if the arcs are

   /// represented with straight lines then they will not intersect

   /// each other.

   ///

   /// Scnyder's algorithm embeds the graph on \c (n-2,n-2) size grid,

   /// i.e. each node will be located in the \c [0,n-2]x[0,n-2] square.

   /// The time complexity of the algorithm is O(n).

   template <typename Graph>

   class PlanarDrawing {

   public:

     TEMPLATE_GRAPH_TYPEDEFS(Graph);

     /// \brief The point type for store coordinates

     typedef dim2::Point<int> Point;

     /// \brief The map type for store coordinates

     typedef typename Graph::template NodeMap<Point> PointMap;

     /// \brief Constructor

     ///

     /// Constructor

     /// \pre The graph should be simple, i.e. loop and parallel arc free.

     PlanarDrawing(const Graph& graph)

       : _graph(graph), _point_map(graph) {}

   private:

     template <typename AuxGraph, typename AuxEmbeddingMap>

     void drawing(const AuxGraph& graph,

                  const AuxEmbeddingMap& next,

                  PointMap& point_map) {

       TEMPLATE_GRAPH_TYPEDEFS(AuxGraph);

       typename AuxGraph::template ArcMap<Arc> prev(graph);

       for (NodeIt n(graph); n != INVALID; ++n) {

         Arc e = OutArcIt(graph, n);

         Arc p = e, l = e;

         e = next[e];

         while (e != l) {

           prev[e] = p;

           p = e;

           e = next[e];

         }

         prev[e] = p;

       }

       Node anode, bnode, cnode;

       {

         Arc e = ArcIt(graph);

         anode = graph.source(e);

         bnode = graph.target(e);

         cnode = graph.target(next[graph.oppositeArc(e)]);

       }

       IterableBoolMap<AuxGraph, Node> proper(graph, false);

       typename AuxGraph::template NodeMap<int> conn(graph, -1);

       conn[anode] = conn[bnode] = -2;

       {

         for (OutArcIt e(graph, anode); e != INVALID; ++e) {

           Node m = graph.target(e);

           if (conn[m] == -1) {

             conn[m] = 1;

           }

         }

         conn[cnode] = 2;

         for (OutArcIt e(graph, bnode); e != INVALID; ++e) {

           Node m = graph.target(e);

           if (conn[m] == -1) {

             conn[m] = 1;

           } else if (conn[m] != -2) {

             conn[m] += 1;

             Arc pe = graph.oppositeArc(e);

             if (conn[graph.target(next[pe])] == -2) {

               conn[m] -= 1;

             }

             if (conn[graph.target(prev[pe])] == -2) {

               conn[m] -= 1;

             }

             proper.set(m, conn[m] == 1);

           }

         }

       }

       typename AuxGraph::template ArcMap<int> angle(graph, -1);

       while (proper.trueNum() != 0) {

         Node n = typename IterableBoolMap<AuxGraph, Node>::TrueIt(proper);

         proper.set(n, false);

         conn[n] = -2;

         for (OutArcIt e(graph, n); e != INVALID; ++e) {

           Node m = graph.target(e);

           if (conn[m] == -1) {

             conn[m] = 1;

           } else if (conn[m] != -2) {

             conn[m] += 1;

             Arc pe = graph.oppositeArc(e);

             if (conn[graph.target(next[pe])] == -2) {

               conn[m] -= 1;

             }

             if (conn[graph.target(prev[pe])] == -2) {

               conn[m] -= 1;

             }

             proper.set(m, conn[m] == 1);

           }

         }

         {

           Arc e = OutArcIt(graph, n);

           Arc p = e, l = e;

           e = next[e];

           while (e != l) {

             if (conn[graph.target(e)] == -2 && conn[graph.target(p)] == -2) {

               Arc f = e;

               angle[f] = 0;

               f = next[graph.oppositeArc(f)];

               angle[f] = 1;

               f = next[graph.oppositeArc(f)];

               angle[f] = 2;

             }

             p = e;

             e = next[e];

           }

           if (conn[graph.target(e)] == -2 && conn[graph.target(p)] == -2) {

             Arc f = e;

             angle[f] = 0;

             f = next[graph.oppositeArc(f)];

             angle[f] = 1;

             f = next[graph.oppositeArc(f)];

             angle[f] = 2;

           }

         }

       }

       typename AuxGraph::template NodeMap<Node> apred(graph, INVALID);

       typename AuxGraph::template NodeMap<Node> bpred(graph, INVALID);

       typename AuxGraph::template NodeMap<Node> cpred(graph, INVALID);

       typename AuxGraph::template NodeMap<int> apredid(graph, -1);

       typename AuxGraph::template NodeMap<int> bpredid(graph, -1);

       typename AuxGraph::template NodeMap<int> cpredid(graph, -1);

       for (ArcIt e(graph); e != INVALID; ++e) {

         if (angle[e] == angle[next[e]]) {

           switch (angle[e]) {

           case 2:

             apred[graph.target(e)] = graph.source(e);

             apredid[graph.target(e)] = graph.id(graph.source(e));

             break;

           case 1:

             bpred[graph.target(e)] = graph.source(e);

             bpredid[graph.target(e)] = graph.id(graph.source(e));

             break;

           case 0:

             cpred[graph.target(e)] = graph.source(e);

             cpredid[graph.target(e)] = graph.id(graph.source(e));

             break;

           }

         }

       }

       cpred[anode] = INVALID;

       cpred[bnode] = INVALID;

       std::vector<Node> aorder, border, corder;

       {

         typename AuxGraph::template NodeMap<bool> processed(graph, false);

         std::vector<Node> st;

         for (NodeIt n(graph); n != INVALID; ++n) {

           if (!processed[n] && n != bnode && n != cnode) {

             st.push_back(n);

             processed[n] = true;

             Node m = apred[n];

             while (m != INVALID && !processed[m]) {

               st.push_back(m);

               processed[m] = true;

               m = apred[m];

             }

             while (!st.empty()) {

               aorder.push_back(st.back());

               st.pop_back();

             }

           }

         }

       }

       {

         typename AuxGraph::template NodeMap<bool> processed(graph, false);

         std::vector<Node> st;

         for (NodeIt n(graph); n != INVALID; ++n) {

           if (!processed[n] && n != cnode && n != anode) {

             st.push_back(n);

             processed[n] = true;

             Node m = bpred[n];

             while (m != INVALID && !processed[m]) {

               st.push_back(m);

               processed[m] = true;

               m = bpred[m];

             }

             while (!st.empty()) {

               border.push_back(st.back());

               st.pop_back();

             }

           }

         }

       }

       {

         typename AuxGraph::template NodeMap<bool> processed(graph, false);

         std::vector<Node> st;

         for (NodeIt n(graph); n != INVALID; ++n) {

           if (!processed[n] && n != anode && n != bnode) {

             st.push_back(n);

             processed[n] = true;

             Node m = cpred[n];

             while (m != INVALID && !processed[m]) {

               st.push_back(m);

               processed[m] = true;

               m = cpred[m];

             }

             while (!st.empty()) {

               corder.push_back(st.back());

               st.pop_back();

             }

           }

         }

       }

       typename AuxGraph::template NodeMap<int> atree(graph, 0);

       for (int i = aorder.size() - 1; i >= 0; --i) {

         Node n = aorder[i];

         atree[n] = 1;

         for (OutArcIt e(graph, n); e != INVALID; ++e) {

           if (apred[graph.target(e)] == n) {

             atree[n] += atree[graph.target(e)];

           }

         }

       }

       typename AuxGraph::template NodeMap<int> btree(graph, 0);

       for (int i = border.size() - 1; i >= 0; --i) {

         Node n = border[i];

         btree[n] = 1;

         for (OutArcIt e(graph, n); e != INVALID; ++e) {

           if (bpred[graph.target(e)] == n) {

             btree[n] += btree[graph.target(e)];

           }

         }

       }

       typename AuxGraph::template NodeMap<int> apath(graph, 0);

       apath[bnode] = apath[cnode] = 1;

       typename AuxGraph::template NodeMap<int> apath_btree(graph, 0);

       apath_btree[bnode] = btree[bnode];

       for (int i = 1; i < int(aorder.size()); ++i) {

         Node n = aorder[i];

         apath[n] = apath[apred[n]] + 1;

         apath_btree[n] = btree[n] + apath_btree[apred[n]];

       }

       typename AuxGraph::template NodeMap<int> bpath_atree(graph, 0);

       bpath_atree[anode] = atree[anode];

       for (int i = 1; i < int(border.size()); ++i) {

         Node n = border[i];

         bpath_atree[n] = atree[n] + bpath_atree[bpred[n]];

       }

       typename AuxGraph::template NodeMap<int> cpath(graph, 0);

       cpath[anode] = cpath[bnode] = 1;

       typename AuxGraph::template NodeMap<int> cpath_atree(graph, 0);

       cpath_atree[anode] = atree[anode];

       typename AuxGraph::template NodeMap<int> cpath_btree(graph, 0);

       cpath_btree[bnode] = btree[bnode];

       for (int i = 1; i < int(corder.size()); ++i) {

         Node n = corder[i];

         cpath[n] = cpath[cpred[n]] + 1;

         cpath_atree[n] = atree[n] + cpath_atree[cpred[n]];

         cpath_btree[n] = btree[n] + cpath_btree[cpred[n]];

       }

       typename AuxGraph::template NodeMap<int> third(graph);

       for (NodeIt n(graph); n != INVALID; ++n) {

         point_map[n].x =

           bpath_atree[n] + cpath_atree[n] - atree[n] - cpath[n] + 1;

         point_map[n].y =

           cpath_btree[n] + apath_btree[n] - btree[n] - apath[n] + 1;

       }

     }

   public:

     /// \brief Calculates the node positions

     ///

     /// This function calculates the node positions.

     /// \return %True if the graph is planar.

     bool run() {

       PlanarEmbedding<Graph> pe(_graph);

       if (!pe.run()) return false;

       run(pe);

       return true;

     }

     /// \brief Calculates the node positions according to a

     /// combinatorical embedding

     ///

     /// This function calculates the node locations. The \c embedding

     /// parameter should contain a valid combinatorical embedding, i.e.

     /// a valid cyclic order of the arcs.

     template <typename EmbeddingMap>

     void run(const EmbeddingMap& embedding) {

       typedef SmartEdgeSet<Graph> AuxGraph;

       if (3 * countNodes(_graph) - 6 == countEdges(_graph)) {

         drawing(_graph, embedding, _point_map);

         return;

       }

       AuxGraph aux_graph(_graph);

       typename AuxGraph::template ArcMap<typename AuxGraph::Arc>

         aux_embedding(aux_graph);

       {

         typename Graph::template EdgeMap<typename AuxGraph::Edge>

           ref(_graph);

         for (EdgeIt e(_graph); e != INVALID; ++e) {

           ref[e] = aux_graph.addEdge(_graph.u(e), _graph.v(e));

         }

         for (EdgeIt e(_graph); e != INVALID; ++e) {

           Arc ee = embedding[_graph.direct(e, true)];

           aux_embedding[aux_graph.direct(ref[e], true)] =

             aux_graph.direct(ref[ee], _graph.direction(ee));

           ee = embedding[_graph.direct(e, false)];

           aux_embedding[aux_graph.direct(ref[e], false)] =

             aux_graph.direct(ref[ee], _graph.direction(ee));

         }

       }

       _planarity_bits::makeConnected(aux_graph, aux_embedding);

       _planarity_bits::makeBiNodeConnected(aux_graph, aux_embedding);

       _planarity_bits::makeMaxPlanar(aux_graph, aux_embedding);

       drawing(aux_graph, aux_embedding, _point_map);

     }

     /// \brief The coordinate of the given node

     ///

     /// The coordinate of the given node.

     Point operator[](const Node& node) const {

       return _point_map[node];

     }

     /// \brief Returns the grid embedding in a \e NodeMap.

     ///

     /// Returns the grid embedding in a \e NodeMap of \c dim2::Point<int> .

     const PointMap& coords() const {

       return _point_map;

     }

   private:

     const Graph& _graph;

     PointMap _point_map;

   };

   namespace _planarity_bits {

     template <typename ColorMap>

     class KempeFilter {

     public:

       typedef typename ColorMap::Key Key;

       typedef bool Value;

       KempeFilter(const ColorMap& color_map,

                   const typename ColorMap::Value& first,

                   const typename ColorMap::Value& second)

         : _color_map(color_map), _first(first), _second(second) {}

       Value operator[](const Key& key) const {

         return _color_map[key] == _first || _color_map[key] == _second;

       }

     private:

       const ColorMap& _color_map;

       typename ColorMap::Value _first, _second;

     };

   }

   /// \ingroup planar

   ///

   /// \brief Coloring planar graphs

   ///

   /// The graph coloring problem is the coloring of the graph nodes

   /// that there are not adjacent nodes with the same color. The

   /// planar graphs can be always colored with four colors, it is

   /// proved by Appel and Haken and their proofs provide a quadratic

   /// time algorithm for four coloring, but it could not be used to

   /// implement efficient algorithm. The five and six coloring can be

   /// made in linear time, but in this class the five coloring has

   /// quadratic worst case time complexity. The two coloring (if

   /// possible) is solvable with a graph search algorithm and it is

   /// implemented in \ref bipartitePartitions() function in LEMON. To

   /// decide whether the planar graph is three colorable is

   /// NP-complete.

   ///

   /// This class contains member functions for calculate colorings

   /// with five and six colors. The six coloring algorithm is a simple

   /// greedy coloring on the backward minimum outgoing order of nodes.

   /// This order can be computed as in each phase the node with least

   /// outgoing arcs to unprocessed nodes is chosen. This order

   /// guarantees that when a node is chosen for coloring it has at

   /// most five already colored adjacents. The five coloring algorithm

   /// use the same method, but if the greedy approach fails to color

   /// with five colors, i.e. the node has five already different

   /// colored neighbours, it swaps the colors in one of the connected

   /// two colored sets with the Kempe recoloring method.

   template <typename Graph>

   class PlanarColoring {

   public:

     TEMPLATE_GRAPH_TYPEDEFS(Graph);

     /// \brief The map type for store color indexes

     typedef typename Graph::template NodeMap<int> IndexMap;

     /// \brief The map type for store colors

     typedef ComposeMap<Palette, IndexMap> ColorMap;

     /// \brief Constructor

     ///

     /// Constructor

     /// \pre The graph should be simple, i.e. loop and parallel arc free.

     PlanarColoring(const Graph& graph)

       : _graph(graph), _color_map(graph), _palette(0) {

       _palette.add(Color(1,0,0));

       _palette.add(Color(0,1,0));

       _palette.add(Color(0,0,1));

       _palette.add(Color(1,1,0));

       _palette.add(Color(1,0,1));

       _palette.add(Color(0,1,1));

     }

     /// \brief Returns the \e NodeMap of color indexes

     ///

     /// Returns the \e NodeMap of color indexes. The values are in the

     /// range \c [0..4] or \c [0..5] according to the coloring method.

     IndexMap colorIndexMap() const {

       return _color_map;

     }

     /// \brief Returns the \e NodeMap of colors

     ///

     /// Returns the \e NodeMap of colors. The values are five or six

     /// distinct \ref lemon::Color "colors".

     ColorMap colorMap() const {

       return composeMap(_palette, _color_map);

     }

     /// \brief Returns the color index of the node

     ///

     /// Returns the color index of the node. The values are in the

     /// range \c [0..4] or \c [0..5] according to the coloring method.

     int colorIndex(const Node& node) const {

       return _color_map[node];

     }

     /// \brief Returns the color of the node

     ///

     /// Returns the color of the node. The values are five or six

     /// distinct \ref lemon::Color "colors".

     Color color(const Node& node) const {

       return _palette[_color_map[node]];

     }

     /// \brief Calculates a coloring with at most six colors

     ///

     /// This function calculates a coloring with at most six colors. The time

     /// complexity of this variant is linear in the size of the graph.

     /// \return %True when the algorithm could color the graph with six color.

     /// If the algorithm fails, then the graph could not be planar.

     /// \note This function can return true if the graph is not

     /// planar but it can be colored with 6 colors.

     bool runSixColoring() {

       typename Graph::template NodeMap<int> heap_index(_graph, -1);

       BucketHeap<typename Graph::template NodeMap<int> > heap(heap_index);

       for (NodeIt n(_graph); n != INVALID; ++n) {

         _color_map[n] = -2;

         heap.push(n, countOutArcs(_graph, n));

       }

       std::vector<Node> order;

       while (!heap.empty()) {

         Node n = heap.top();

         heap.pop();

         _color_map[n] = -1;

         order.push_back(n);

         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

           Node t = _graph.runningNode(e);

           if (_color_map[t] == -2) {

             heap.decrease(t, heap[t] - 1);

           }

         }

       }

       for (int i = order.size() - 1; i >= 0; --i) {

         std::vector<bool> forbidden(6, false);

         for (OutArcIt e(_graph, order[i]); e != INVALID; ++e) {

           Node t = _graph.runningNode(e);

           if (_color_map[t] != -1) {

             forbidden[_color_map[t]] = true;

           }

         }

                for (int k = 0; k < 6; ++k) {

           if (!forbidden[k]) {

             _color_map[order[i]] = k;

             break;

           }

         }

         if (_color_map[order[i]] == -1) {

           return false;

         }

       }

       return true;

     }

   private:

     bool recolor(const Node& u, const Node& v) {

       int ucolor = _color_map[u];

       int vcolor = _color_map[v];

       typedef _planarity_bits::KempeFilter<IndexMap> KempeFilter;

       KempeFilter filter(_color_map, ucolor, vcolor);

       typedef FilterNodes<const Graph, const KempeFilter> KempeGraph;

       KempeGraph kempe_graph(_graph, filter);

       std::vector<Node> comp;

       Bfs<KempeGraph> bfs(kempe_graph);

       bfs.init();

       bfs.addSource(u);

       while (!bfs.emptyQueue()) {

         Node n = bfs.nextNode();

         if (n == v) return false;

         comp.push_back(n);

         bfs.processNextNode();

       }

       int scolor = ucolor + vcolor;

       for (int i = 0; i < static_cast<int>(comp.size()); ++i) {

         _color_map[comp[i]] = scolor - _color_map[comp[i]];

       }

       return true;

     }

     template <typename EmbeddingMap>

     void kempeRecoloring(const Node& node, const EmbeddingMap& embedding) {

       std::vector<Node> nodes;

       nodes.reserve(4);

       for (Arc e = OutArcIt(_graph, node); e != INVALID; e = embedding[e]) {

         Node t = _graph.target(e);

         if (_color_map[t] != -1) {

           nodes.push_back(t);

           if (nodes.size() == 4) break;

         }

       }