source: lemon-0.x/lemon/planarity.h @ 2480:eecaeab41472

/* -*- C++ -*-
 *
 * This file is a part of LEMON, a generic C++ optimization library
 *
 * Copyright (C) 2003-2007
 * 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 graph_prop
/// \file
/// \brief Planarity checking, embedding

#include <vector>
#include <list>

#include <lemon/dfs.h>
#include <lemon/radix_sort.h>
#include <lemon/maps.h>
#include <lemon/path.h>


namespace lemon {

  namespace _planarity_bits {

    template <typename UGraph>
    struct PlanarityVisitor : DfsVisitor<UGraph> {

      typedef typename UGraph::Node Node;
      typedef typename UGraph::Edge Edge;

      typedef typename UGraph::template NodeMap<Edge> PredMap;

      typedef typename UGraph::template UEdgeMap<bool> TreeMap;

      typedef typename UGraph::template NodeMap<int> OrderMap;
      typedef std::vector<Node> OrderList;

      typedef typename UGraph::template NodeMap<int> LowMap;
      typedef typename UGraph::template NodeMap<int> AncestorMap;

      PlanarityVisitor(const UGraph& ugraph,
                       PredMap& pred_map, TreeMap& tree_map,
                       OrderMap& order_map, OrderList& order_list,
                       AncestorMap& ancestor_map, LowMap& low_map)
        : _ugraph(ugraph), _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 Edge& edge) {
        Node source = _ugraph.source(edge);
        Node target = _ugraph.target(edge);

        _tree_map[edge] = true;
        _pred_map[target] = edge;
      }

      void examine(const Edge& edge) {
        Node source = _ugraph.source(edge);
        Node target = _ugraph.target(edge);
        
        if (_order_map[target] < _order_map[source] && !_tree_map[edge]) {
          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 Edge& edge) {
        Node source = _ugraph.source(edge);
        Node target = _ugraph.target(edge);

        if (_low_map[source] > _low_map[target]) {
          _low_map[source] = _low_map[target];
        }
      }

      const UGraph& _ugraph;
      PredMap& _pred_map;
      TreeMap& _tree_map;
      OrderMap& _order_map;
      OrderList& _order_list;
      AncestorMap& _ancestor_map;
      LowMap& _low_map;
    };

    template <typename UGraph, bool embedding = true>
    struct NodeDataNode {
      int prev, next;
      int visited;
      typename UGraph::Edge first;
      bool inverted;
    };

    template <typename UGraph>
    struct NodeDataNode<UGraph, false> {
      int prev, next;
      int visited;
    };

    template <typename UGraph>
    struct ChildListNode {
      typedef typename UGraph::Node Node;
      Node first;
      Node prev, next;
    };

    template <typename UGraph>
    struct EdgeListNode {
      typename UGraph::Edge prev, next;
    };

  }

  /// \ingroup  graph_prop
  ///
  /// \brief Planarity checking of an undirected simple graph
  ///
  /// This class implements the Boyer-Myrvold algorithm for planar
  /// checking of an undirected graph. This class is a simplified
  /// version of the PlanarEmbedding algorithm class, and it does
  /// provide neither embedding nor kuratowski subdivisons.
  template <typename UGraph>
  class PlanarityChecking {
  private:
    
    UGRAPH_TYPEDEFS(typename UGraph)
      
    const UGraph& _ugraph;

  private:

    typedef typename UGraph::template NodeMap<Edge> PredMap;
    
    typedef typename UGraph::template UEdgeMap<bool> TreeMap;

    typedef typename UGraph::template NodeMap<int> OrderMap;
    typedef std::vector<Node> OrderList;

    typedef typename UGraph::template NodeMap<int> LowMap;
    typedef typename UGraph::template NodeMap<int> AncestorMap;

    typedef _planarity_bits::NodeDataNode<UGraph> NodeDataNode;
    typedef std::vector<NodeDataNode> NodeData;
    
    typedef _planarity_bits::ChildListNode<UGraph> ChildListNode;
    typedef typename UGraph::template NodeMap<ChildListNode> ChildLists;

    typedef typename UGraph::template NodeMap<std::list<int> > MergeRoots;
 
    typedef typename UGraph::template NodeMap<bool> EmbedEdge;

  public:

    /// \brief Constructor
    ///
    /// \warining The graph should be simple, i.e. parallel and loop edge
    /// free.
    PlanarityChecking(const UGraph& ugraph) : _ugraph(ugraph) {}

    /// \brief Runs the algorithm.
    ///
    /// Runs the algorithm.  
    /// \param kuratowski If the parameter is false, then the
    /// algorithm does not calculate the isolate Kuratowski
    /// subdivisions.
    /// \return %True when the graph is planar.
    bool run(bool kuratowski = true) {
      typedef _planarity_bits::PlanarityVisitor<UGraph> Visitor;

      PredMap pred_map(_ugraph, INVALID);
      TreeMap tree_map(_ugraph, false);

      OrderMap order_map(_ugraph, -1);
      OrderList order_list;

      AncestorMap ancestor_map(_ugraph, -1);
      LowMap low_map(_ugraph, -1);

      Visitor visitor(_ugraph, pred_map, tree_map,
                      order_map, order_list, ancestor_map, low_map);
      DfsVisit<UGraph, Visitor> visit(_ugraph, visitor);
      visit.run();

      ChildLists child_lists(_ugraph);
      createChildLists(tree_map, order_map, low_map, child_lists);

      NodeData node_data(2 * order_list.size());
      
      EmbedEdge embed_edge(_ugraph, false);

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

        Node node = order_list[i];

        Node source = node;
        for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
          Node target = _ugraph.target(e);
          
          if (order_map[source] < order_map[target] && tree_map[e]) {
            initFace(target, node_data, pred_map, order_map, order_list);
          }
        }
        
        for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
          Node target = _ugraph.target(e);
          
          if (order_map[source] < order_map[target] && !tree_map[e]) {
            embed_edge[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_edge, merge_roots);
        }
        merge_roots[node].clear();
        
        for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
          Node target = _ugraph.target(e);
          
          if (order_map[source] < order_map[target] && !tree_map[e]) {
            if (embed_edge[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(_ugraph); n != INVALID; ++n) {
        Node source = n;
        
        std::vector<Node> targets;  
        for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
          Node target = _ugraph.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(), mapFunctor(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 = _ugraph.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,
                  EmbedEdge& embed_edge, 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_edge[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 edge 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_edge[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_edge, 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 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();
      
    }

    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 EmbedEdge& embed_edge,
                   const MergeRoots& merge_roots) {
      return !merge_roots[node].empty() || embed_edge[node];
    }

  };

  /// \ingroup graph_prop
  ///
  /// \brief Planar embedding of an undirected simple graph
  ///
  /// This class implements the Boyer-Myrvold algorithm for planar
  /// embedding of an undirected graph. The planar embeding is an
  /// ordering of the outgoing edges in each node, 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 an \f$ K_{3,3} \f$ (complete bipartite graph on
  /// 3 ANode and 3 BNode) subdivision.
  ///
  /// The current implementation calculates an embedding or an
  /// Kuratowski subdivision if the graph is not planar. The running
  /// time of the algorithm is \f$ O(n) \f$.
  template <typename UGraph>
  class PlanarEmbedding {
  private:
    
    UGRAPH_TYPEDEFS(typename UGraph)
      
    const UGraph& _ugraph;
    typename UGraph::template EdgeMap<Edge> _embedding;

    typename UGraph::template UEdgeMap<bool> _kuratowski;

  private:

    typedef typename UGraph::template NodeMap<Edge> PredMap;
    
    typedef typename UGraph::template UEdgeMap<bool> TreeMap;

    typedef typename UGraph::template NodeMap<int> OrderMap;
    typedef std::vector<Node> OrderList;

    typedef typename UGraph::template NodeMap<int> LowMap;
    typedef typename UGraph::template NodeMap<int> AncestorMap;

    typedef _planarity_bits::NodeDataNode<UGraph> NodeDataNode;
    typedef std::vector<NodeDataNode> NodeData;
    
    typedef _planarity_bits::ChildListNode<UGraph> ChildListNode;
    typedef typename UGraph::template NodeMap<ChildListNode> ChildLists;

    typedef typename UGraph::template NodeMap<std::list<int> > MergeRoots;
 
    typedef typename UGraph::template NodeMap<Edge> EmbedEdge;

    typedef _planarity_bits::EdgeListNode<UGraph> EdgeListNode;
    typedef typename UGraph::template EdgeMap<EdgeListNode> EdgeLists;

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

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

    enum IsolatorNodeType {
      HIGHX = 6, LOWX = 7,
      HIGHY = 8, LOWY = 9,
      ROOT = 10, PERTINENT = 11,
      INTERNAL = 12
    };

  public:

    /// \brief Constructor
    ///
    /// \warining The graph should be simple, i.e. parallel and loop edge
    /// free.
    PlanarEmbedding(const UGraph& ugraph)
      : _ugraph(ugraph), _embedding(_ugraph), _kuratowski(ugraph, false) {}

    /// \brief Runs the algorithm.
    ///
    /// Runs the algorithm.  
    /// \param kuratowski If the parameter is false, then the
    /// algorithm does not calculate the isolate Kuratowski
    /// subdivisions.
    ///\return %True when the graph is planar.
    bool run(bool kuratowski = true) {
      typedef _planarity_bits::PlanarityVisitor<UGraph> Visitor;

      PredMap pred_map(_ugraph, INVALID);
      TreeMap tree_map(_ugraph, false);

      OrderMap order_map(_ugraph, -1);
      OrderList order_list;

      AncestorMap ancestor_map(_ugraph, -1);
      LowMap low_map(_ugraph, -1);

      Visitor visitor(_ugraph, pred_map, tree_map,
                      order_map, order_list, ancestor_map, low_map);
      DfsVisit<UGraph, Visitor> visit(_ugraph, visitor);
      visit.run();

      ChildLists child_lists(_ugraph);
      createChildLists(tree_map, order_map, low_map, child_lists);

      NodeData node_data(2 * order_list.size());
      
      EmbedEdge embed_edge(_ugraph, INVALID);

      MergeRoots merge_roots(_ugraph);

      EdgeLists edge_lists(_ugraph);

      FlipMap flip_map(_ugraph, 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 (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
          Node target = _ugraph.target(e);
          
          if (order_map[source] < order_map[target] && tree_map[e]) {
            initFace(target, edge_lists, node_data,
                      pred_map, order_map, order_list);
          }
        }
        
        for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
          Node target = _ugraph.target(e);
          
          if (order_map[source] < order_map[target] && !tree_map[e]) {
            embed_edge[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, edge_lists, flip_map, order_list, 
                   child_lists, ancestor_map, low_map, embed_edge, merge_roots);
        }
        merge_roots[node].clear();
        
        for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
          Node target = _ugraph.target(e);
          
          if (order_map[source] < order_map[target] && !tree_map[e]) {
            if (embed_edge[target] != INVALID) {
              if (kuratowski) {
                isolateKuratowski(e, node_data, edge_lists, flip_map,
                                  order_map, order_list, pred_map, child_lists,
                                  ancestor_map, low_map, 
                                  embed_edge, 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, edge_lists);
        storeEmbedding(order_list[i], node_data, order_map, pred_map,
                       edge_lists, flip_map);
      }

      return true;
    }

    /// \brief Gives back the successor of an edge
    ///
    /// Gives back the successor of an edge. This function makes
    /// possible to query the cyclic order of the outgoing edges from
    /// a node.
    Edge next(const Edge& edge) const {
      return _embedding[edge];
    }

    /// \brief Gives back true when the undirected edge is in the
    /// kuratowski subdivision
    ///
    /// Gives back true when the undirected edge is in the kuratowski
    /// subdivision
    bool kuratowski(const UEdge& uedge) {
      return _kuratowski[uedge];
    }

  private:

    void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,
                          const LowMap& low_map, ChildLists& child_lists) {

      for (NodeIt n(_ugraph); n != INVALID; ++n) {
        Node source = n;
        
        std::vector<Node> targets;  
        for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
          Node target = _ugraph.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(), mapFunctor(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 = _ugraph.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,
                  EdgeLists& edge_lists, FlipMap& flip_map, 
                  OrderList& order_list, ChildLists& child_lists,
                  AncestorMap& ancestor_map, LowMap& low_map,
                  EmbedEdge& embed_edge, 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_edge[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 edges + flipping
              Edge de = node_data[dn].first;
              Edge ce = node_data[cn].first;

              flip_map[order_list[cn - order_list.size()]] = cd != dd;
              if (cd != dd) {
                std::swap(edge_lists[ce].prev, edge_lists[ce].next);
                ce = edge_lists[ce].prev;
                std::swap(edge_lists[ce].prev, edge_lists[ce].next);
              }

              {
                Edge dne = edge_lists[de].next; 
                Edge cne = edge_lists[ce].next; 

                edge_lists[de].next = cne;
                edge_lists[ce].next = dne;
              
                edge_lists[dne].prev = ce;
                edge_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 edge
            {
              Edge edge = embed_edge[node];
              Edge re = node_data[rn].first;

              edge_lists[edge_lists[re].next].prev = edge;
              edge_lists[edge].next = edge_lists[re].next;
              edge_lists[edge].prev = re;
              edge_lists[re].next = edge;

              if (!rd) {
                node_data[rn].first = edge;
              }

              Edge rev = _ugraph.oppositeEdge(edge);
              Edge e = node_data[n].first;

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

              if (d) {
                node_data[n].first = rev;
              }
              
            }

            // Embedding edge 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_edge[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_edge, 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, EdgeLists& edge_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;

      Edge edge = pred_map[node];
      Edge rev = _ugraph.oppositeEdge(edge);

      node_data[rn].first = edge;
      node_data[n].first = rev;

      edge_lists[edge].prev = edge;
      edge_lists[edge].next = edge;

      edge_lists[rev].prev = rev;
      edge_lists[rev].next = rev;

    }

    void mergeRemainingFaces(const Node& node, NodeData& node_data,
                             OrderList& order_list, OrderMap& order_map,
                             ChildLists& child_lists, EdgeLists& edge_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;

        Edge de = node_data[dd].first;
        Edge ce = node_data[cd].first;

        if (de != INVALID) {
          Edge dne = edge_lists[de].next; 
          Edge cne = edge_lists[ce].next; 
          
          edge_lists[de].next = cne;
          edge_lists[ce].next = dne;
          
          edge_lists[dne].prev = ce;
          edge_lists[cne].prev = de;
        }
        
        node_data[dd].first = ce;

      }
    }

    void storeEmbedding(const Node& node, NodeData& node_data,
                        OrderMap& order_map, PredMap& pred_map,
                        EdgeLists& edge_lists, FlipMap& flip_map) {

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

      if (pred_map[node] != INVALID) {
        Node source = _ugraph.source(pred_map[node]);
        flip_map[node] = flip_map[node] != flip_map[source];
      }
      
      Edge first = node_data[order_map[node]].first;
      Edge prev = first;

      Edge edge = flip_map[node] ?
        edge_lists[prev].prev : edge_lists[prev].next;

      _embedding[prev] = edge;
      
      while (edge != first) {
        Edge next = edge_lists[edge].prev == prev ?
          edge_lists[edge].next : edge_lists[edge].prev;
        prev = edge; edge = next;
        _embedding[prev] = edge;
      }
    }


    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 EmbedEdge& embed_edge,
                   const MergeRoots& merge_roots) {
      return !merge_roots[node].empty() || embed_edge[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,
                       EmbedEdge& embed_edge, MergeRoots& merge_roots) {
      Node wnode =_ugraph.target(node_data[order_map[node]].first);
      while (!pertinent(wnode, embed_edge, merge_roots)) {
        wnode = _ugraph.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 =_ugraph.target(node_data[order_map[node]].first);
      while (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {
        wnode = _ugraph.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, EdgeLists& edge_lists, 
                        EmbedEdge& embed_edge, 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_edge, 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) {
              Edge edge = node_data[cn].first;
              _kuratowski.set(edge, true);
              
              pred = cnode;
              cnode = _ugraph.target(edge);
            
              continue;
            }
          }
          wnode = znode = cnode;
          return;

        } else if (pert) {
          wnode = cnode;
          
          while (!external(cnode, rorder, child_lists, ancestor_map, low_map)) {
            Edge edge = node_data[order_map[cnode]].first;
          
            if (_ugraph.target(edge) == pred) {
              edge = edge_lists[edge].next;
            }
            _kuratowski.set(edge, true);
            
            Node next = _ugraph.target(edge);
            pred = cnode; cnode = next;
          }
          
          znode = cnode;
          return;

        } else if (ext) {
          znode = cnode;
          
          while (!pertinent(cnode, embed_edge, merge_roots)) {
            Edge edge = node_data[order_map[cnode]].first;
          
            if (_ugraph.target(edge) == pred) {
              edge = edge_lists[edge].next;
            }
            _kuratowski.set(edge, true);
            
            Node next = _ugraph.target(edge);
            pred = cnode; cnode = next;
          }
          
          wnode = cnode;
          return;
          
        } else {
          Edge edge = node_data[order_map[cnode]].first;
          
          if (_ugraph.target(edge) == pred) {
            edge = edge_lists[edge].next;
          }
          _kuratowski.set(edge, true);

          Node next = _ugraph.target(edge);
          pred = cnode; cnode = next;
        }
        
      }

    }

    void orientComponent(Node root, int rn, OrderMap& order_map,
                         PredMap& pred_map, NodeData& node_data, 
                         EdgeLists& edge_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);
        
        Edge edge = node_data[order_map[node]].first;
        
        if (type_map[_ugraph.target(edge)] == 0) {
          st.push_back(_ugraph.target(edge));
          type_map[_ugraph.target(edge)] = 1;
        } 
        
        Edge last = edge, pred = edge;
        edge = edge_lists[edge].next;
        while (edge != last) {

          if (type_map[_ugraph.target(edge)] == 0) {
            st.push_back(_ugraph.target(edge));
            type_map[_ugraph.target(edge)] = 1;
          } 
          
          Edge next = edge_lists[edge].next != pred ? 
            edge_lists[edge].next : edge_lists[edge].prev;
          pred = edge; edge = 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 = _ugraph.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) {
            Edge edge = node_data[order_map[node]].first;
            std::swap(edge_lists[edge].prev, edge_lists[edge].next);
            edge = edge_lists[edge].prev;
            std::swap(edge_lists[edge].prev, edge_lists[edge].next);
            node_data[order_map[node]].first = edge;
          }
        }
      }

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

        Edge edge = node_data[order_map[qu[i]]].first;
        Edge last = edge, pred = edge;

        edge = edge_lists[edge].next;
        while (edge != last) {

          if (edge_lists[edge].next == pred) {
            std::swap(edge_lists[edge].next, edge_lists[edge].prev);
          } 
          pred = edge; edge = edge_lists[edge].next;
        }
        
      }
    }

    void setFaceFlags(Node root, Node wnode, Node ynode, Node xnode,
                      OrderMap& order_map, NodeData& node_data, 
                      TypeMap& type_map) {
      Node node = _ugraph.target(node_data[order_map[root]].first);

      while (node != ynode) {
        type_map[node] = HIGHY;
        node = _ugraph.target(node_data[order_map[node]].first);
      }

      while (node != wnode) {
        type_map[node] = LOWY;
        node = _ugraph.target(node_data[order_map[node]].first);
      }
      
      node = _ugraph.target(node_data[order_map[wnode]].first);
      
      while (node != xnode) {
        type_map[node] = LOWX;
        node = _ugraph.target(node_data[order_map[node]].first);
      }
      type_map[node] = LOWX;

      node = _ugraph.target(node_data[order_map[xnode]].first);
      while (node != root) {
        type_map[node] = HIGHX;
        node = _ugraph.target(node_data[order_map[node]].first);
      }

      type_map[wnode] = PERTINENT;
      type_map[root] = ROOT;
    }

    void findInternalPath(std::vector<Edge>& ipath,
                          Node wnode, Node root, TypeMap& type_map, 
                          OrderMap& order_map, NodeData& node_data, 
                          EdgeLists& edge_lists) {
      std::vector<Edge> st;

      Node node = wnode;
      
      while (node != root) {
        Edge edge = edge_lists[node_data[order_map[node]].first].next;
        st.push_back(edge);
        node = _ugraph.target(edge);
      }
      
      while (true) {
        Edge edge = st.back();
        if (type_map[_ugraph.target(edge)] == LOWX ||
            type_map[_ugraph.target(edge)] == HIGHX) {
          break;
        }
        if (type_map[_ugraph.target(edge)] == 2) {
          type_map[_ugraph.target(edge)] = 3;
          
          edge = edge_lists[_ugraph.oppositeEdge(edge)].next;
          st.push_back(edge);
        } else {
          st.pop_back();
          edge = edge_lists[edge].next;
          
          while (_ugraph.oppositeEdge(edge) == st.back()) {
            edge = st.back();
            st.pop_back();
            edge = edge_lists[edge].next;
          }
          st.push_back(edge);
        }
      }
      
      for (int i = 0; i < int(st.size()); ++i) {
        if (type_map[_ugraph.target(st[i])] != LOWY &&
            type_map[_ugraph.target(st[i])] != HIGHY) {
          for (; i < int(st.size()); ++i) {
            ipath.push_back(st[i]);
          }
        }
      }
    }

    void setInternalFlags(std::vector<Edge>& ipath, TypeMap& type_map) {
      for (int i = 1; i < int(ipath.size()); ++i) {
        type_map[_ugraph.source(ipath[i])] = INTERNAL;
      }
    }

    void findPilePath(std::vector<Edge>& ppath,
                      Node root, TypeMap& type_map, OrderMap& order_map, 
                      NodeData& node_data, EdgeLists& edge_lists) {
      std::vector<Edge> st;

      st.push_back(_ugraph.oppositeEdge(node_data[order_map[root]].first));
      st.push_back(node_data[order_map[root]].first);
      
      while (st.size() > 1) {
        Edge edge = st.back();
        if (type_map[_ugraph.target(edge)] == INTERNAL) {
          break;
        }
        if (type_map[_ugraph.target(edge)] == 3) {
          type_map[_ugraph.target(edge)] = 4;
          
          edge = edge_lists[_ugraph.oppositeEdge(edge)].next;
          st.push_back(edge);
        } else {
          st.pop_back();
          edge = edge_lists[edge].next;
          
          while (!st.empty() && _ugraph.oppositeEdge(edge) == st.back()) {
            edge = st.back();
            st.pop_back();
            edge = edge_lists[edge].next;
          }
          st.push_back(edge);
        }
      }
      
      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 (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
            Node tnode = _ugraph.target(e); 
            if (order_map[tnode] > order_map[node] && low_map[tnode] == lp) {
              node = tnode;
              _kuratowski[e] = true;
              break;
            }
          }
        }
      }

      for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
        if (order_map[_ugraph.target(e)] == lp) {
          _kuratowski[e] = true;
          break;
        }
      }
      
      return lp;
    }

    void markPertinentPath(Node node, OrderMap& order_map, 
                           NodeData& node_data, EdgeLists& edge_lists,
                           EmbedEdge& embed_edge, MergeRoots& merge_roots) {
      while (embed_edge[node] == INVALID) {
        int n = merge_roots[node].front();
        Edge edge = node_data[n].first;

        _kuratowski.set(edge, true);
        
        Node pred = node;
        node = _ugraph.target(edge);
        while (!pertinent(node, embed_edge, merge_roots)) {
          edge = node_data[order_map[node]].first;
          if (_ugraph.target(edge) == pred) {
            edge = edge_lists[edge].next;
          }
          _kuratowski.set(edge, true);
          pred = node;
          node = _ugraph.target(edge);
        }
      }
      _kuratowski.set(embed_edge[node], true);
    } 

    void markPredPath(Node node, Node snode, PredMap& pred_map) {
      while (node != snode) {
        _kuratowski.set(pred_map[node], true);
        node = _ugraph.source(pred_map[node]);
      }
    }

    void markFacePath(Node ynode, Node xnode, 
                      OrderMap& order_map, NodeData& node_data) {
      Edge edge = node_data[order_map[ynode]].first;
      Node node = _ugraph.target(edge);
      _kuratowski.set(edge, true);
        
      while (node != xnode) {
        edge = node_data[order_map[node]].first;
        _kuratowski.set(edge, true);
        node = _ugraph.target(edge);
      }
    }

    void markInternalPath(std::vector<Edge>& path) {
      for (int i = 0; i < int(path.size()); ++i) {
        _kuratowski.set(path[i], true);
      }
    }

    void markPilePath(std::vector<Edge>& path) {
      for (int i = 0; i < int(path.size()); ++i) {
        _kuratowski.set(path[i], true);
      }
    }

    void isolateKuratowski(Edge edge, NodeData& node_data, 
                           EdgeLists& edge_lists, FlipMap& flip_map,
                           OrderMap& order_map, OrderList& order_list, 
                           PredMap& pred_map, ChildLists& child_lists,
                           AncestorMap& ancestor_map, LowMap& low_map, 
                           EmbedEdge& embed_edge, MergeRoots& merge_roots) {

      Node root = _ugraph.source(edge);
      Node enode = _ugraph.target(edge);

      int rorder = order_map[root];

      TypeMap type_map(_ugraph, 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 != _ugraph.source(edge)) {
          orientComponent(root, rn, order_map, pred_map, 
                          node_data, edge_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_edge, merge_roots);
          
          markPertinentPath(lwnode, order_map, node_data, edge_lists,
                            embed_edge, merge_roots);
          
          return;
        }
      }
      
      orientComponent(root, rn, order_map, pred_map, 
                      node_data, edge_lists, flip_map, type_map);

      Node wnode = findPertinent(ynode, order_map, node_data,
                                 embed_edge, 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, edge_lists, embed_edge, 
                         merge_roots, child_lists, ancestor_map, low_map);
                  
          markPertinentPath(lwnode, order_map, node_data, edge_lists,
                            embed_edge, 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<Edge> ipath;
      findInternalPath(ipath, wnode, root, type_map, order_map,
                       node_data, edge_lists);
      setInternalFlags(ipath, type_map);
      pynode = _ugraph.source(ipath.front());
      pxnode = _ugraph.target(ipath.back());

      wnode = findPertinent(pynode, order_map, node_data,
                            embed_edge, merge_roots);
      
      // Minor-C
      {
        if (type_map[_ugraph.source(ipath.front())] == HIGHY) {
          if (type_map[_ugraph.target(ipath.back())] == HIGHX) {
            markFacePath(xnode, pxnode, order_map, node_data);
          }
          markFacePath(root, xnode, order_map, node_data);
          markPertinentPath(wnode, order_map, node_data, edge_lists,
                            embed_edge, 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[_ugraph.target(ipath.back())] == HIGHX) {
          markFacePath(ynode, root, order_map, node_data);
          markPertinentPath(wnode, order_map, node_data, edge_lists,
                            embed_edge, 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<Edge> ppath;
      findPilePath(ppath, root, type_map, order_map, node_data, edge_lists);
      
      // Minor-D
      if (!ppath.empty()) {
        markFacePath(ynode, xnode, order_map, node_data);
        markPertinentPath(wnode, order_map, node_data, edge_lists,
                          embed_edge, 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, edge_lists,
                              embed_edge, 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, edge_lists,
                              embed_edge, 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, edge_lists,
                          embed_edge, 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;
      }
      
    }

  };

}

#endif