// -*- C++ -*-
/* Copyright (c) 1997-2005
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
   http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
*/

#ident "$Project: polymake $$Id: morse_matching_tools.h 6604 2005-12-08 19:15:48Z pfetsch $"

/**@file   morse_matching_tools.h
 * @brief  tools for computations with Morse matchings
 * @author Marc Pfetsch
 *
 * @todo Check whether overflow can occur in @a base (see checkAcyclicDFS() ).
 */

#ifndef MORSE_MATCHING_TOOLS_H
#define MORSE_MATCHING_TOOLS_H

#include <connected.h>
#include <Array.h>
#include <Map.h>
#include <PowerSet.h>

#include <vector>
#include <cassert>
#include <deque>

using namespace polymake;




//-------------------------------------------------------------------------------------------------------------
//------------------------------------- Check Morse properties ------------------------------------------------
//-------------------------------------------------------------------------------------------------------------


/**@brief Recursive function to detect cycles arising from reversing the arcs in @a M.
 * @param M         collection of arcs in Hasse diagram
 * @param marked    marker for DFS-search (node has been visited)
 * @param v         node to be searched
 * @param lower     whether we are on the lower part
 * @param base      base value for marker (see below)
 * @returns @c true if the layer is acyclic
 *
 * We use a depth-first-search (DFS) to detect (directed) cycles in the Hasse diagram
 * where arcs whose value is @c true are reversed. Such cycles cannot span several
 * levels. We therefore search each level individually. The parameter @c lower stores
 * whether the current node @a v is on the lower or the upper part of the level.
 *
 * To avoid reinitialization of @c marked we use the following convention for the markers:
 * - if the marker is less than @c base, the node has not been visited
 * - if the marker is equal to @c base, the node is in the current active subtree
 * - if the marker is equal to @c base + 1, the node and all its children have been processed.
 *
 * If we find an arc from @a v to some node marked with @c base we found a cycle.
 * If its mark is less than @c base we recurse on this node. If the mark is @c base + 1
 * we found a node in a different part of the tree, which does not yield a cycle.
 *
 * This function is called, e.g. from @p checkAcyclic.
 */
bool checkAcyclicDFS(const HasseDiagram<Set<int>,bool>& M, Array<int>& marked, int v, bool lower, int base)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::in_edge_list>::const_iterator HasseDiagramInConstIterator;
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   // mark node v as touched
   marked[v] = base;

   // if @v is on the lower side of a level
   if (lower)
   {
      // check all out-edges (pointing up in the Hasse diagram)
      for (HasseDiagramOutConstIterator e = entire(M.out_edges(v)); !e.at_end(); ++e)
      {
	 // if the arc is in the arc set, i.e. reversed (pointing down)
	 if ( *e )
	 {
	    int source = e.to_node();
	    // if the source-node is in the active tree we found a cycle
	    if (marked[source] == base)
	       return false;
	    else
	    {
	       // if source-node is unmarked recurse
	       if (marked[source] < base)
	       {
		  if (! checkAcyclicDFS(M, marked, source, false, base) )
		     return false;
	       }
	    }
	 }
      }
   }
   else // otherwise we are on the upper size of a level
   {
      // check all in-edges (pointing down in the Hasse diagram)
      for (HasseDiagramInConstIterator e = entire(M.in_edges(v)); !e.at_end(); ++e)
      {
	 if (! *e )
	 {
	    int target = e.from_node();
	    // if the target-node is in the active tree we found a cycle
	    if (marked[target] == base)
	       return false;
	    else
	    {
	       // if target-node is unmarked recurse
	       if (marked[target] < base)
	       {
		  if (! checkAcyclicDFS(M, marked, target, true, base) )
		     return false;
	       }
	    }
	 }
      }
   }
   marked[v] = base + 1;
   return true;
}





/**@brief Check whether a given solution is acyclic
 * @param M         collection of arcs in Hasse diagram
 * @param marked    marker for DFS-search (node has been visited)
 * @param base      base value for marker (see @c checkAcyclicDFS())
 * @returns @c true if the solution is acyclic
 *
 * We assume that @c marked has been initialized some time ago.
 * @warning @a base is changed!
 */
bool checkAcyclic(const HasseDiagram<Set<int>,bool>& M, Array<int>& marked, int& base)
{
   int d = M.dim() - 1;      // do not count empty face and top face

   // check each level in turn
   for (int k = 0; k < d; ++k)
   {
      // increase base
      base += 2;

      // check each face of dim k
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 int v = *f;
	 assert( v > 0 && v <= (M.nodes()-2) );

	 // for each unmarked face
	 if ( marked[v] < base )
	 {
	    bool acyclic = checkAcyclicDFS(M, marked, v, true, base);
	    if (! acyclic)
	       return false;
	 }
      }
   }
   return true;
}




/**@brief Check whether a given solution is acyclic
 * @param M         collection of arcs in Hasse diagram
 * @param k	    dimension of the level to be checked
 * @param marked    marker for DFS-search (node has been visited)
 * @param base      base value for marker (see @c checkAcyclicDFS())
 * @returns @c true if the solution is acyclic
 *
 * We assume that @c marked has been initialized some time ago.
 */
bool checkAcyclic(const HasseDiagram<Set<int>,bool>& M, int k, Array<int>& marked, int base)
{
   // check each face of dim k
   for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
   {
      int v = *f;
      assert( v > 0 && v <= (M.nodes()-2) );

      // for each unmarked face
      if ( marked[v] < base )
      {
	 bool acyclic = checkAcyclicDFS(M, marked, v, true, base);
	 if (! acyclic)
	    return false;
      }
   }
   return true;
}



/**@brief Check whether a given solution is acyclic
 * @param M  collection of arcs in Hasse diagram
 * @returns @c true if the solution is acyclic
 */
bool checkAcyclic(const HasseDiagram<Set<int>,bool>& M)
{
   int d = M.dim() - 1;      // do not count empty face
   int n = M.nodes() - 2;    // and top face

   Array<int> marked(n+1);
   // initialize markers
   for (int i = 0; i <= n; ++i)
      marked[i] = 0;

   // for each level in turn
   int base = -1;
   for (int k = 0; k < d; ++k)
   {
      // increase base
      base += 2;

      // check each face of dim k
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 int v = *f;
	 assert( (v > 0) && (v <= n) );

	 // for each unmarked face
	 if ( marked[v] < base )
	 {
	    bool acyclic = checkAcyclicDFS(M, marked, v, true, base);
	    if (! acyclic)
	       return false;
	 }
      }
   }
   return true;
}



/**@brief Check whether the given collection of arcs satisfies the matching constraints
 * @param M  collection of arcs in Hasse diagram
 * @returns @c true if the given arc set is a matching
 */
bool checkMatching(const HasseDiagram<Set<int>,bool>& M)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::in_edge_list>::const_iterator HasseDiagramInConstIterator;
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   int d = M.dim() - 1;

   // for each level in turn
   for (int k = 0; k <= d; ++k)
   {
      // check each face of dim k
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 int inc = 0;  // count incident arcs
	 if ( k > 0 )
	 {
	    for (HasseDiagramInConstIterator e = entire(M.in_edges(*f)); !e.at_end(); ++e)
	    {
	       // if arc is in collection, i.e. reverse we have an incident arc
	       if ( *e )
		  ++inc;
	       // if a node has more than two incident arcs the matching property is violated
	       if (inc >= 2)
		  return false;
	    }
	 }
	 if ( k < d )
	 {
	    for (HasseDiagramOutConstIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	    {
	       // if arc is in collection, i.e. reverse we have an incident arc
	       if ( *e )
		  ++inc;
	       // if a node has more than two incident arcs the matching property is violated
	       if (inc >= 2)
		  return false;
	    }
	 }
      }
   }
   return true;
}







//-------------------------------------------------------------------------------------------------------------
//------------------------------------- Other functions -------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------


/**@brief Find the critical faces w.r.t the Morse matching in @a M
 * @param M	    Morse matching in Hasse diagram
 * @returns bitset indicating the critical faces (nodes)
 *
 * On output an element (node) is set to true in the bitset if and only if the face is critical.
 */
Bitset collectCriticalFaces(const HasseDiagram<Set<int>,bool>& M)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::in_edge_list>::const_iterator HasseDiagramInConstIterator;
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   int d = M.dim() - 1;      // do not count empty face
   int n = M.nodes() - 2;    // and top face

   // ensure space
   Bitset critical(n+1);

   // loop over all levels
   for (int k = 0; k <= d; ++k)
   {
      // pass through all faces of dimension k
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 bool isCritical = true;
	 // if the dimension is larger than 0, we can look at in-arcs
	 if (k > 0)
	 {
	    // pass through all in-arcs
	    for (HasseDiagramInConstIterator e = entire(M.in_edges(*f)); !e.at_end(); ++e)
	    {
	       // if the arc is in the matching (i.e., reversed), f is not critical
	       if ( *e )
	       {
		  isCritical = false;
		  break;
	       }
	    }
	 }
	 // if the dimension is smaller than d, we can look at out-arcs
	 if ( isCritical && k < d )
	 {
	    // pass through all out-arcs
	    for (HasseDiagramOutConstIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	    {
	       // if the arc is in the matching (i.e., reversed), f is not critical
	       if ( *e )
	       {
		  isCritical = false;
		  break;
	       }
	    }
	 }
	 if ( isCritical )
	    critical += *f;
      }
   }
   return critical;
}







/**@brief Find the critical faces w.r.t the Morse matching in @a M
 * @param M	    Morse matching in Hasse diagram
 * @returns set containing critical faces
 */
PowerSet<int> findCriticalFaces(const HasseDiagram<Set<int>,bool>& M)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::in_edge_list>::const_iterator HasseDiagramInConstIterator;
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   int d = M.dim() - 1;      // do not count empty face

   // ensure space
   PowerSet<int> critical;

   // loop over all levels
   for (int k = 0; k <= d; ++k)
   {
      // pass through all faces of dimension k
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 bool isCritical = true;
	 // if the dimension is larger than 0, we can look at in-arcs
	 if (k > 0)
	 {
	    // pass through all in-arcs
	    for (HasseDiagramInConstIterator e = entire(M.in_edges(*f)); !e.at_end(); ++e)
	    {
	       // if the arc is in the matching (i.e., reversed), f is not critical
	       if ( *e )
	       {
		  isCritical = false;
		  break;
	       }
	    }
	 }
	 // if the dimension is smaller than d, we can look at out-arcs
	 if ( isCritical && k < d )
	 {
	    // pass through all out-arcs
	    for (HasseDiagramOutConstIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	    {
	       // if the arc is in the matching (i.e., reversed), f is not critical
	       if ( *e )
	       {
		  isCritical = false;
		  break;
	       }
	    }
	 }
	 if ( isCritical )
	    critical += M.node(*f);
      }
   }
   return critical;
}





/**@brief Compute the size of the subgraph in @a M
 * @param M  subset of the arcs
 * @returns the size of the arcs marked with @c true
 */
int subgraphSize(const HasseDiagram<Set<int>,bool>& M)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   int d = M.dim() - 1;      // do not count empty face
   int size = 0;

   // loop over all levels
   for (int k = 0; k < d; ++k)
   {
      // pass through all faces of dimension k
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 // pass through all out-arcs
	 for (HasseDiagramOutConstIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	 {
	    // if the arc is in the subset
	    if ( *e )
	       ++size;
	 }
      }
   }
   return size;
}



/**@class CompareByProperty
 * @brief Compares two objects by their value given by a property
 */
template <class Object, class Property>
class CompareByProperty
{
public:
   CompareByProperty(const Property& P) : P_(P)
   {}

   bool operator() (const Object &p, const Object &q) const
   {
      if ( P_[p] < P_[q] )
	 return true;
      return false;
   }

private:
   const Property& P_;
};





/**@brief Compute lexicographic ordering of the edges
 * @param M            subset of the arcs
 * @param order        output: order of arcs
 * @param bottomLevel  lowest level
 * @param topLevel     hightest level
 *
 * - We loop through each level, beginning at the top
 * - We sort the faces in a level lexicographically.
 * - For each face @a F in the lexicographic order,
 *   we sort its facets lexicographically and insert
 *   the corresponding arc numbers into the order.
 *
 * @a order is a std::vector, because we want to change its size dynamically.
 * Might be replaced by @c Array and precomputation of the size.
 */
void orderEdgesLex(const HasseDiagram<Set<int>,bool>& M, std::vector<int>& order, int bottomLevel, int topLevel)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::in_edge_list>::const_iterator HasseDiagramInConstIterator;
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   int d = M.dim() - 1;      // do not count empty face

   order.clear();

   // create edge-number map
   Map<std::pair<int,int>, int> E;   // map to store numbers for edges

   int m = 0;
   for (int k = 0; k < d; ++k)
   {
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 for (HasseDiagramOutConstIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	 {
	    int source = e.from_node();
	    int target = e.to_node();
	    E[std::make_pair(source, target)] = m++;
	 }
      }
   }

   // loop over levels, starting from the top most
   for (int k = topLevel; k > bottomLevel; --k)
   {
      // collect all nodes of dimension k (use std::vector, because we don't know the size in advance)
      std::vector<Set<int> > FS;    // vertex sets of faces
      std::vector<int> FN;          // face numbers
      std::vector<int> index;       // indices to sort
      int i = 0;
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 Set<int> F = M.node(*f);  // get vertex set corr. to face
	 FS.push_back(F);          // store set and
	 FN.push_back(*f);         // face number
	 index.push_back(i++);
      }

      // sort
      CompareByProperty<int,std::vector<Set<int> > > cmp(FS);
      std::sort(index.begin(), index.end(), cmp);

      // pass through faces of dim k in lexicographic order
      std::vector<int>::const_iterator fit, fend = index.end();
      for (fit = index.begin(); fit != fend; ++fit)
      {
	 Set<int> F = FS[*fit];   // get current face

	 // collect neighbors of F
	 std::vector<Set<int> > N;
	 std::vector<int> EN;
	 std::vector<int> index2;
	 int l = 0;
	 for (HasseDiagramInConstIterator e = entire(M.in_edges(FN[*fit])); !e.at_end(); ++e)
	 {
	    Set<int> S = M.node(e.from_node());   // get vertex set corr. to neighbor
	    N.push_back(S);                       // store it
	    EN.push_back(E[std::make_pair(e.from_node(), e.to_node())]);  // store edge number
	    index2.push_back(l++);
	 }

	 // sort
	 CompareByProperty<int,std::vector<Set<int> > > Cmp(N);
	 std::sort(index2.begin(), index2.end(), Cmp);

	 /// pass through neighbors in lexicographic order
	 std::vector<int>::const_iterator it, iend = index2.end();
	 for (it = index2.begin(); it != iend; ++it)
	    order.push_back(EN[*it]);
      }
   }
}






/**@brief Set Morse matching to the empty set
 * @param M         Morse matching
 */
inline
void clearMorseMatching(HasseDiagram<Set<int>,bool>& M)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::iterator HasseDiagramOutIterator;

   int d = M.dim() - 1;
   for (int k = 0; k < d; ++k)
   {
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 for (HasseDiagramOutIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	    *e = false;
      }
   }
}




//-------------------------------------------------------------------------------------------------------------
//------------------------------------- Greedy Algorithm ------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------



/**@brief Compute a Morse matching by a simple greedy heuristic
 * @param M         Hasse diagram
 * @param varLevel  level of arcs
 * @param varOrder  order of arcs
 *
 * @note @a varOrder can contain a subset of the arc indices, but
 * @a varLevel stores the level of @b all arcs.
 */
template <class LevelMap, class Iterator>
int greedyHeuristic(HasseDiagram<Set<int>,bool>& M, const LevelMap& varLevel, Iterator orderIt, Iterator orderEnd)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::iterator HasseDiagramOutIterator;

   int d = M.dim() - 1;
   int n = M.nodes() - 2;
   int size = 0;
   int m = varLevel.size();

   Array<bool> matched(n+1);
   Array<int> marked(n+1);
   Array<HasseDiagramOutIterator> V(m);

   // init solution to the empty set and collect edge iterators
   int i = 0;
   for (int k = 0; k < d; ++k)
   {
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k)); !f.at_end(); ++f)
      {
	 for (HasseDiagramOutIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	 {
	    *e = false;
	    V[i] = e;
	    ++i;
	 }
      }
   }
   assert(i <= m);

   // init to no matched or marked faces
   for (int i = 0; i <= n; ++i)
   {
      matched[i] = false;
      marked[i] = 0;
   }

   int base = 1;
   for (; orderIt != orderEnd; ++orderIt)
   {
      int var = *orderIt;
      HasseDiagramOutIterator e = V[var];
      int source = e.from_node();
      int target = e.to_node();

      if ((! matched[source]) && (! matched[target]) )
      {
	 // tentatively include var in solution and check ...
	 *e = true;

	 assert( M.dim_of_node(source) < M.dim_of_node(target) );

	 if ( checkAcyclicDFS(M, marked, source, true, base) )
	 {
	    ++size;
	    matched[source] = true;
	    matched[target] = true;
	 }
	 else
	    *e = false;
	 // increase base
	 base += 2;
      }
   }
   return size;
}







//-------------------------------------------------------------------------------------------------------------
//------------------------------------- Find alternating paths ------------------------------------------------
//-------------------------------------------------------------------------------------------------------------



/**@brief  Depth first search to find alternating paths
 * @param M	    Morse matching in Hasse diagram
 * @param marked    marker for dfs-search (node has been visited)
 * @param p	    stores the DFS-tree edges
 * @param v	    start node
 * @param lower     whether we are on the lower part
 *
 * This is a modified DFS, see processAlternatingPaths for a description.
 */
void findAlternatingPathDFS(const HasseDiagram<Set<int>,bool>& M, Array<int>& marked, Array<int>& p, int v, bool lower)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::in_edge_list>::const_iterator HasseDiagramInConstIterator;
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   marked[v] = 1;
   if (lower)  // if we are one the lower side of a level
   {
      for (HasseDiagramOutConstIterator e = entire(M.out_edges(v)); !e.at_end(); ++e)
      {
	 if ( *e )
	 {
	    int u = e.to_node();
	    if (marked[u] == 0)
	    {
	       p[u] = v;
	       findAlternatingPathDFS(M, marked, p, u, false);
	    }
	    else
	       ++(marked[u]);
	 }
      }
   }
   else
   {
      for (HasseDiagramInConstIterator e = entire(M.in_edges(v)); ! e.at_end(); ++e)
      {
	 if (! *e )
	 {
	    int u = e.from_node();
	    if (marked[u] == 0)
	    {
	       p[u] = v;
	       findAlternatingPathDFS(M, marked, p, u, true);
	    }
	    else
	       ++(marked[u]);
	 }
      }
   }
}




/**@brief Exchange matching arcs on alternating paths
 * @param M	Morse matching in Hasse diagram
 * @param p	stores the DFS-tree edges
 * @param v	start node
 * @param v	end node
 * @param size	size of the matching
 */
void exchangePath(HasseDiagram<Set<int>,bool>& M, const Array<int>& p, int u, int v, int& size)
{
   int w = v;   // start node
   do
   {
      int ww = p[w];
      if ( M.edge_exists(w, ww) )
      {
	 bool r = M.edge(w,ww);
	 M.edge(w,ww) = !r;
	 if ( r )
	    --size;
	 else
	    ++size;
      }
      else
      {
	 bool r = M.edge(ww,w);
	 M.edge(ww,w) = !r;
	 if ( r )
	    --size;
	 else
	    ++size;
      }
      w = ww;
   }
   while (w != u);
}







/**@brief Improve a Morse matching by canceling critical faces along alternating paths.
 * @param M	       Morse matching in Hasse diagram
 * @param size	       size of the computed solution
 * @param bottomLevel  lowest level to consider
 * @param topLevel     highest level to consider
 *
 * An alternating path w.r.t. a given Morse matching @a M is a path in the modified
 * Hasse diagram @a H(M) starting at a critical face and ending at a critical face;
 * it is called alternating because normal arcs and flipped arcs alternate on this
 * path. Forman observed that if there exists a @a unique alternating path then
 * one can flip the direction of every arc in the path to again obtain a Morse
 * matching. The two critical faces at the end of the path are not critical w.r.t.
 * the new Morse matching.
 *
 * Alternating paths can be found by a modified DFS as follows. We work on the
 * acyclic graph @a H(M) and start at a critical face @a u and then recurse. Initially
 * all faces have been marked 0. Let @a v be the current face and @a a = (@a v,@a w)
 * be an out-arc of @a v. If @a w has not been visited (its mark is 0), then we mark
 * @a w with 1 and set the DFS tree arc to be @a a. If @a w has been visited (marked
 * with a positive number) then we increase its mark by one.
 *
 * After running this algorithm we can check each critical face @a v (not equal
 * to @a u). If its mark is 1, we check whether each face on the path from @a v to
 * @a u is marked with 1. If this is the case, then we found a unique path from @a u
 * to @a v. This can be seen as follows:
 *
 * If there exists a unique path, then every face @a w on the path to @a v is visited
 * exactly once (so it is marked with 1).
 *
 * Conversely, assume that the path is not unique, i.e. there exist at least two
 * distinct paths from @a u to @a v (if there exists no path, then @a v is not marked).
 * DFS follows each path until it reaches an already visited face @a w; then the mark
 * of @a w is increased by one. Hence, at least one face on every path from @a u to
 * @a v has mark at least 2 (twice visisted).
 *
 * The algorithm is then simple to describe. We first compute all critical faces
 * w.r.t. the current solution. We loop over all critical faces @a u and run the
 * modified DFS. If there exists a critical face such that each face on the path to
 * @a u is marked with 1, we flip all arcs and set the two endpoints as non-critical.
 *
 * Note that the result may depend on the order of the paths flipped.
 */
void processAlternatingPaths(HasseDiagram<Set<int>,bool>& M, int& size, int bottomLevel, int topLevel)
{
   int n = M.nodes() - 2;
   int cnt = 0; // number of alternating paths

   // compute critical faces
   Bitset critical = collectCriticalFaces(M);

   // find alternating paths
   Array<int> marked(n+1);
   Array<int> dfsTree(n+1);
   for (int k = bottomLevel; k < topLevel; ++k)
   {
      // check only nodes at the upper level
      for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(k+1)); !f.at_end(); ++f)
      {
	 int u = *f;
	 if ( critical.contains(u) )
	 {
	    for (int i = 0; i <= n; ++i)
	    {
	       dfsTree[i] = -1;
	       marked[i] = 0;
	    }

	    findAlternatingPathDFS(M, marked, dfsTree, u, false);

	    // check whether we found a path: lower part
	    for (Entire<sequence>::const_iterator f2 = entire(M.nodes_of_dim(k)); !f2.at_end(); ++f2)
	    {
	       int v = *f2;

	       //if ( (v != u) && critical.contains(v) && (marked[v] == 1) )
	       if ( critical.contains(v) && (marked[v] == 1) )
	       {
		  // check whether path contains only nodes visited once
		  int w = v;
		  do
		  {
		     w = dfsTree[w];
		     assert( w >= 0 );
		     if ( marked[w] != 1 )
			break;
		  }
		  while (w != u);

		  if (w == u)
		  {
		     // exchange path
		     exchangePath(M, dfsTree, u, v, size);
		     ++cnt;

		     critical -= u;  // remove u and v
		     critical -= v;
		     break;
		  }
	       }
	    }
	 }
      }
   }
#if POLY_DEBUG
   cout << "Changed alternating paths: " << cnt << endl;
#endif
}












//-------------------------------------------------------------------------------------------------------------
//------------------------------------- Compute Maximum Forests -----------------------------------------------
//-------------------------------------------------------------------------------------------------------------



/**@brief Modified Prim algorithm to find a maximum forest in a graph
 * @param G        graph
 * @param p        map to store edge to predecessor
 * @param visited  visted map
 *
 * We do not care about weights, we just want to produce an arbitrary spanning forest.
 *
 * - visited = -1 :  node deleted
 * - visited = 0  :  node not visited
 * - visited = 1  :  node visited and not root in the forest
 * - visited = 2  :  node visited and root in the forest
 *
 * The map @a p stores the predecessor of a node. Its type must be
 * equal to the type stored on each edge. The map should provide
 * operator[].
 *
 * @todo Use concept checks to guarantee that the type stored on each edge is
 *       the same as the one store in the predecessor map.
 *
 * @see completeToBottomLevel(), completToTopLevel()
 */
template <class Graph, class PMap>
void findMaximumForest(const Graph& G, PMap& p, Array<int>& visited)
{
   typedef typename Entire<typename Graph::out_edge_list>::const_iterator GraphOutConstIterator;

   // initialize visited:
   for (int i = 0; i < G.nodes(); ++i)
      visited[i] = 0;

   // search for start vertex and then start BFS iteration
   for (int start = 0; start < G.nodes(); ++start)
   {
      if ( visited[start] == 0 )   // if vertex has not be seen so far
      {
	 visited[start] = 2;   // mark vertex as root
	 std::deque<int> Q;

	 // initialize edges from start vertex
	 for (GraphOutConstIterator e = entire(G.out_edges(start)); !e.at_end(); ++e)
	 {
	    int w = e.to_node();
	    if ( visited[w] == 0 )
	    {
	       p[w] = *e;
	       Q.push_back(w);
	    }
	 }

	 // start BFS iterations
	 while (! Q.empty() )
	 {
	    int v = Q.front();
	    Q.pop_front();
	    visited[v] = 1;

	    for (GraphOutConstIterator e = entire(G.out_edges(v)); !e.at_end(); ++e)
	    {
	       int w = e.to_node();
	       if ( visited[w] == 0 )
	       {
		  p[w] = *e;
		  Q.push_back(w);
	       }
	    }
	 }
      }
   }
}





/**@brief Modified Prim algorithm to find a maximum forest in a graph, first trying to included marked nodes
 * @param G        graph
 * @param marked   map to mark vertices
 * @param p        map to store edge to predecessor
 * @param visited  visted map
 *
 * This function is similar to findMaximumForest(), but first tries to include nodes that are
 * marked by @a marked.
 *
 * We do not care about weights, we just want to produce an arbitrary spanning forest.
 *
 * - visited = -1 :  node deleted
 * - visited = 0  :  node not visited
 * - visited = 1  :  node visited and not root in the forest
 * - visited = 2  :  node visited and root in the forest
 *
 * The map @a p stores the predecessor of a node. Its type must be
 * equal to the type stored on each edge. The map should provide
 * operator[].
 *
 * @todo Use concept checks to guarantee that the type stored on each edge is
 *       the same as the one store in the predecessor map.
 *
 * @see completeToBottomLevel(), completToTopLevel()
 */
template <class Graph, class PMap, class MarkedMap>
void findMaximumForestMarked(const Graph& G, MarkedMap& marked, PMap& p, Array<int>& visited)
{
   typedef typename Entire<typename Graph::out_edge_list>::const_iterator GraphOutConstIterator;

   // initialize visited:
   for (int i = 0; i < G.nodes(); ++i)
      visited[i] = 0;

   // first check for marked nodes
   for (int start = 0; start < G.nodes(); ++start)
   {
      if ( visited[start] == 0 && marked[start] )   // if loop is not visited
      {
	 visited[start] = 2;
	 std::deque<int> Q;

	 // initialize edges from start vertex
	 for (GraphOutConstIterator e = entire(G.out_edges(start)); !e.at_end(); ++e)
	 {
	    int w = e.to_node();
	    if ( visited[w] == 0 )
	    {
	       p[w] = *e;
	       Q.push_back(w);
	    }
	 }

	 // start BFS iterations
	 while (! Q.empty() )
	 {
	    int v = Q.front();
	    Q.pop_front();
	    visited[v] = 1;

	    for (GraphOutConstIterator e = entire(G.out_edges(v)); !e.at_end(); ++e)
	    {
	       int w = e.to_node();
	       if ( visited[w] == 0 )
	       {
		  p[w] = *e;
		  Q.push_back(w);
	       }
	    }
	 }
      }
   }

   // search for start vertex and then start BFS iteration
   for (int start = 0; start < G.nodes(); ++start)
   {
      if ( visited[start] == 0 )   // if vertex has not be seen so far
      {
	 visited[start] = 2;   // mark vertex as root
	 std::deque<int> Q;

	 // initialize edges from start vertex
	 for (GraphOutConstIterator e = entire(G.out_edges(start)); !e.at_end(); ++e)
	 {
	    int w = e.to_node();
	    if ( visited[w] == 0 )
	    {
	       p[w] = *e;
	       Q.push_back(w);
	    }
	 }

	 // start BFS iterations
	 while (! Q.empty() )
	 {
	    int v = Q.front();
	    Q.pop_front();
	    visited[v] = 1;

	    for (GraphOutConstIterator e = entire(G.out_edges(v)); !e.at_end(); ++e)
	    {
	       int w = e.to_node();
	       if ( visited[w] == 0 )
	       {
		  p[w] = *e;
		  Q.push_back(w);
	       }
	    }
	 }
      }
   }
}





//-------------------------------------------------------------------------------------------------------------
//------------------------------------- Complete Solution -----------------------------------------------------
//-------------------------------------------------------------------------------------------------------------



/**@brief Complete the Morse matching to the bottom level by a maximum forest computation
 * @param M   Morse matching
 *
 * The current Morse matching in @a M is extended to the bottom level as follows:
 *
 * -# Compute the currently critical faces.
 * -# Generate the graph of the complex with non-critical edges removed.
 * -# Compute a maximum forest on the graph.
 * -# The forest is given by edges oriented towards a root.
 * -# The matching is extended by the vertices matched with the corresponding
 *    edges leading in the direction of the root.
 *
 * Currently we assume that the complex is connected. It follows that the graph
 * is connected as well (see Joswig & Pfetsch [2005]).
 */
void completeToBottomLevel(HasseDiagram<Set<int>,bool>& M)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::in_edge_list>::const_iterator HasseDiagramInConstIterator;

   // find critical faces of best_solution_
   Bitset critical = collectCriticalFaces(M);

   // build helper graph
   Graph<int,int> G;

   // to translate between face indices and nodes of the graph:
   Map<int, int> FTON;

   // create nodes of the graph ( = vertices (0-faces) of complex )
   for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(0)); !f.at_end(); ++f)
   {
      int v = G.add_node(*f);  // add 1-face
      FTON[*f] = v;
   }

   // create edges (arising from 1-faces in the complex)
   for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(1)); !f.at_end(); ++f)
   {
      if ( critical.contains(*f) )
      {
	 // find the two nodes incident to the edge
	 int v1 = -1;
	 int v2 = -1;
	 for (HasseDiagramInConstIterator e = entire(M.in_edges(*f)); !e.at_end(); ++e)
	 {
	    if (v1 == -1)
	       v1 = e.from_node();
	    else
	    {
	       assert( v2 == -1);
	       v2 = e.from_node();
	    }
	 }
	 assert(v2 != -1);

	 // create edge:
	 G.create_edge(FTON[v1], FTON[v2], *f);
      }
   }

   // prepare computation of spanning forest
   int numNodes = G.nodes();
   Array<int> visited(numNodes);
   Array<int> p(numNodes);

   // check whether graph is really connected (not strictly necessary!)
   assert( graph::is_connected(G) );

   findMaximumForest(G, p, visited);

   // now add corresponding arcs to the matching
   for (int i = 0; i < numNodes; ++i)
   {
      assert( visited[i] != 0);
      // if vertex is in the forest, but not the root
      if (visited[i] == 1)
      {
	 int fnum = p[i];       // index of 1-face in complex
	 int vnum = G.node(i);  // index of 0-face in complex
	 M.edge(vnum, fnum) = true;
      }
   }
}






/**@brief Complete the Morse matching to the top level by a maximum forest computation
 * @param M   Morse matching
 *
 * This is only possible if M arises from a pseudo-manifold.
 *
 * The current Morse matching in @a M is extended to the top level as in
 * @c completeToBottomLevel(), but with the following exceptions:
 *
 * - We work on the dual graph of the complex.
 * - The boundary faces produce a loop in the dual graph.
 * - When computing a spanning forest we first try to take vertices with a loop
 *   as the root, since this produces a component without critical face.
 * - When completing the matching, the arcs corresponding to loops are added.
 *
 */
void completeToTopLevel(HasseDiagram<Set<int>,bool>& M)
{
   typedef Entire<HasseDiagram<Set<int>,bool>::out_edge_list>::const_iterator HasseDiagramOutConstIterator;

   int d = M.dim() - 1;

   // find critical faces of best_solution_
   Bitset critical = collectCriticalFaces(M);

   // build helper graph
   Graph<int,int> G;

   // to translate between face indices and nodes of the graph:
   Map<int, int> FTON;

   // to store loops
   std::vector<bool> loop;

   // create nodes of the graph ( = vertices (d-faces) of complex )
   for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(d)); !f.at_end(); ++f)
   {
      int v = G.add_node(*f);  // add d-face
      FTON[*f] = v;
      loop.push_back(false);   // assume that node numbers a consecutive
   }

   // create edges (arising from 1-faces in the complex)
   for (Entire<sequence>::const_iterator f = entire(M.nodes_of_dim(d-1)); !f.at_end(); ++f)
   {
      if ( critical.contains(*f) )
      {
	 // find the two nodes incident to the edge
	 int v1 = -1;
	 int v2 = -1;
	 for (HasseDiagramOutConstIterator e = entire(M.out_edges(*f)); !e.at_end(); ++e)
	 {
	    if (v1 == -1)
	       v1 = e.to_node();
	    else
	    {
	       assert( v2 == -1);
	       v2 = e.to_node();
	    }
	 }
	 assert(v1 != -1);

	 // either produce an edge or a loop:
	 if (v2 != -1)
	    G.create_edge(FTON[v1], FTON[v2], *f);
	 else
	    loop[FTON[v1]] = true;   // store that the d-face corr. to a loop
      }
   }

   // prepare computation of spanning forest
   int numNodes = G.nodes();
   Array<int> visited(numNodes);
   Array<int> p(numNodes);

   findMaximumForestMarked(G, loop, p, visited);

   // now add corresponding arcs to the matching
   for (int i = 0; i < numNodes; ++i)
   {
      assert( visited[i] != 0);
      // if vertex is in the forest or the root and a loop
      if ( visited[i] == 1 || (visited[i] == 2 && loop[i]) )
      {
	 int rnum = p[i];       // index of (d-1)-face in complex
	 int fnum = G.node(i);  // index of d-face in complex
	 M.edge(rnum, fnum) = true;
      }
   }
}

#endif
