/* -*- C++ -*-

 *

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

 *

 * Copyright (C) 2003-2008

 * 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/iterable_maps.h>

#include <lemon/edge_set.h>

#include <lemon/bucket_heap.h>

#include <lemon/ugraph_adaptor.h>

#include <lemon/color.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 planar

  ///

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

  ///

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

  /// 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.  

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

    bool run() {

      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, 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 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 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 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 The map for store of embedding

    typedef typename UGraph::template EdgeMap<Edge> EmbeddingMap;

    /// \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 the calculated embedding map

    ///

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

    /// graph.

    const EmbeddingMap& embedding() const {

      return _embedding;

    }

    /// \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;

	    }