/* Copyright (c) 1997-2004
   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.
*/

#ifndef _POLYMAKE_PS_BISTELLAR_COMPLEX_H
#define _POLYMAKE_PS_BISTELLAR_COMPLEX_H "$Project: polymake $$Id: psBistellarComplex.h,v 1.11 2004/12/08 09:34:02 gawrilow Exp $"

#include <complex_tools.h>
#include <RandomSubset.h>
#include <PseudoComplex.h>
#include <Bitset.h>

#include <list>
#include <ext/hash_map>
#include <ext/hash_set>

namespace polymake { namespace topaz {
   
class psBistellarComplex {
   
protected:

  class OptionsList {
    
  protected:
    int the_size;
    Array<int> the_options;
    std_ext::hash_map< int,int > index_map;

  public:
     OptionsList() : the_size(0) {};
     
     void clear()
     {
	the_size = 0;
	index_map.clear();
     }
	
     int size()
     {
	return the_size;
     }
     
     void insert(const int opt)
     {
	if (the_size >= the_options.size())
	   the_options.resize(2*the_size+1);
	the_options[the_size]=opt;
	index_map[opt]=the_size;
	++the_size;
     }
     
     void remove(const int opt)
     {
	std_ext::hash_map< int,int >::const_iterator find=index_map.find(opt);
	if (find!=index_map.end()) {
	   const int index = find->second;
	   the_options[index] = the_options[the_size-1];
	   index_map[ the_options[index] ] = index;
	   index_map.erase(opt);
	   --the_size;
	}
     }
     
     Array<int> options() {
	return Array<int>(the_size,the_options.begin());
     }
  };  // end class OptionsList
   
   PseudoComplex C;
   Array<OptionsList> options;
   int dim;
   long seed;
   std::list<int> dead_nodes;
   int next_move;
   Set<int> rev_move;
   bool allow_rev_move;
   Array<int> the_f_vector;
   bool verbose;
   
public:
   psBistellarComplex() : dim(0), seed(0), next_move(-1), allow_rev_move(false), verbose(false) {};
   
protected:
   template <typename NodesType>
   Bitset search_up(const NodesType& nodes, Set<int>& facets)
   {
      std::list<int> queue;
      Bitset visited(C.nodes());
      for (typename Entire<NodesType>::const_iterator n=entire(nodes);
	   !n.at_end(); ++n) {
	 queue.push_back(*n);
	 visited += *n;
      }

      while (!queue.empty()) {
	 const int n=queue.front();
	 queue.pop_front();
	 
	 for (Entire<PseudoComplex::out_edge_list>::const_iterator e=entire(C.out_edges(n));
	      !e.at_end(); ++e) {
	    const int nn=e.to_node();

	    if (nn==C.top_node())  // n is facet
	       facets += n;
	    else
	       if (!visited.contains(nn)) {
		  visited += nn;
		  queue.push_back(nn);
	       }
	 }
      }
      
      return visited;
   }
   
   template <typename NodesType>
   Set<int> search_up(NodesType nodes)
   {
      Set<int> facets;
      search_up(nodes,facets);
      return facets;
   }

   Bitset search_up(const int n, Set<int>& facets)
   {
      std::list<int> nodes;
      nodes.push_back(n);
      return search_up(nodes,facets);
   }

   Set<int> search_up(const int n)
   {
      std::list<int> nodes;
      nodes.push_back(n);
      return search_up(nodes);
   }
   
   template <typename NodesType>
   Bitset search_down(const NodesType& nodes, Set<int>& vertices)
   {
      std::list<int> queue;
      Bitset visited(C.nodes());
      for (typename Entire<NodesType>::const_iterator n=entire(nodes);
	   !n.at_end(); ++n) {
	 queue.push_back(*n);
	 visited += *n;
      }

      while (!queue.empty()) {
	 const int n=queue.front();
	 queue.pop_front();
	 
	 for (Entire<PseudoComplex::in_edge_list>::const_iterator e=entire(C.in_edges(n));
	      !e.at_end(); ++e) {
	    const int nn=e.from_node();
	    
	   if (nn==C.bottom_node())  // n is vertex
	      vertices += n;
	   else
	      if (!visited.contains(nn)) {
		 visited += nn;
		 queue.push_back(nn);
	      }
	 }
       }

      return visited;
   }
   
   template <typename NodesType>
   Set<int> search_down(NodesType nodes)
   {
      Set<int> vertices;
      search_down(nodes,vertices);
      return vertices;
   }

   Bitset search_down(const int n, Set<int>& vertices)
   {
      std::list<int> nodes;
      nodes.push_back(n);
      return search_down(nodes,vertices);
   }
   
   Set<int> search_down(const int n)
   {
      std::list<int> nodes;
      nodes.push_back(n);
      return search_down(nodes);
   }
   
   Set<int> verts_of_link(const int f)
   {
      const Set<int> star_facets = search_up(f);
      return search_down(star_facets) - C.node(f);
   }

   bool is_option(const int f, const int d)
   {
      if (d==dim)
	 return true;

      // check #facets in star(f)
      const Set<int> star_facets = search_up(f);
      if (star_facets.size() != dim-d+1)
	 return false;

      // check if star is a simplicial complex
      Set<int> tmp;
      const Bitset star = search_down(star_facets, tmp);
      tmp.clear();

      std_ext::hash_set< Set<int> > faces;
      for (Entire<Bitset>::const_iterator s=entire(star); !s.at_end(); ++s) {
	 const Set<int>& face = C.node(*s);
	 std_ext::hash_set< Set<int> >::const_iterator find=faces.find(face);
	 if (find!=faces.end())  // face exists twice
	    return false;
	 else
	    faces.insert(face);
      }

      return true;
   }
   
public:
   psBistellarComplex(const PseudoComplex& C_in, const long in_seed,
		      const bool verb=false, const bool in_allow_rev_move=false) {
      C = C_in;
      C.restore_vertex_sets();  // make sure, vertices are indexed in C by themselves
      C.restore();
      verbose = verb;
      dim = C.dim();
      seed = in_seed;
      the_f_vector.resize(dim+1);
      next_move = -1;
      rev_move = Set<int>();
      allow_rev_move = in_allow_rev_move;
      options.resize(dim+1);
      
      // compute options
      for (int d=0; d<=dim; ++d) {
	 the_f_vector[d] = C.nodes_of_dim(d).size();
	 
	 for (Entire< std::list<int> >::const_iterator n=entire(C.nodes_of_dim(d));
	      !n.at_end(); ++n)    
	    if (is_option(*n,d))
	       options[d].insert(*n);
      }
   }
   
   int n_options_of_dim(const int d) {
      return options[d].options().size();
   }

   // Finds a reversed move of minimal dimension >= dim_min
   // and returns the dimension of the move.
   // The move is stored in next_move.
   int find_move(const int dim_min) {
      for (int d=dim_min; d<=dim; ++d) {
	 if (options[d].options().size()>1) {
	    const RandomPermutation< Array<int> > P(options[d].options(),seed++);
	    
	    for (Entire< RandomPermutation< Array<int> > >::const_iterator opt=(entire(P));
		 !opt.at_end(); ++opt)
	       if ( allow_rev_move || incl(C.node(*opt),rev_move)!=0 ) {
		  next_move = *opt;
		  return d;
	       }
	 }
      }
      
      throw std::runtime_error("psBistellarComplex: No move found.");
   } 
   
   /*
   // Finds the maximal down move
   int cool_move() {
   const int cool_dim = (dim-1)/2;  // maximal down move
   for (int d=cool_dim; d>=0; --d)
   if (options[d].options().size()>1) {
   const RandomPermutation< Array<int> > P(options[d].options(),seed++);
   const RandomPermutation< Array<int> >::const_iterator opt = P.begin();
   next_move = *opt;
   execute_move();
   return d;
   }

   const int d = find_move(cool_dim+1);
   execute_move();
   
   return d;
   }  
   */   

   int find_move() {
      return find_move(0);
   }
   
   // Executes what ever move is set by find_move() or find_move(const int).
   // You MUST set a move by using find_move() or find_move(const int) befor executing it. 
   void execute_move()
   {
#ifdef POLY_DEBUG
      Array<int> f_vec = f_vector();
      int e=0;
      for (int i=0; i<f_vec.size(); ++i)
	 e += f_vec[i]*(1-2*(i%2));
      
      if (!is_option(next_move,C.node(next_move).size()-1))
      	 throw std::runtime_error("next_move is not an option");
#endif

      const Set<int> face = C.node(next_move);
      const int face_dim = face.size()-1;
      
      Set<int> co_face;
      if (face_dim==dim)  // face is a facet
	 if (dead_nodes.empty()) {
	    co_face = scalar2set(C.nodes());
	    C.add_node(co_face);
	 } else {
	    co_face = scalar2set(dead_nodes.front());
	    C.node(dead_nodes.front()) = co_face;
	    dead_nodes.pop_front();
	 }
      else
	 co_face = verts_of_link(next_move);

      if (!allow_rev_move)  rev_move = co_face;
      
      if (verbose)
	 cout << "psBistellarComplex: executing move of dim "
	      << face.size()-1 << ": (" << face << "," << co_face << ")\n";

      // remove star(face) from options and C
      // and compute boundary(star(face))
      const Set<int> star_facets = search_up(next_move);
      std::list<int> boundary;
      Set<int> tmp;
      const Bitset star = search_down(star_facets, tmp);
      tmp.clear();

      for (Entire<Bitset>::const_iterator s=entire(star);
	   !s.at_end(); ++s) {
	 const Set<int> f = C.node(*s);
	 const int f_dim = f.size()-1;
	 options[f_dim].remove(*s);

	 if (incl(face,f)<=0) {  // face is contained in vertset of *s
	    C.in_edges(*s).clear();
	    C.out_edges(*s).clear();
	    dead_nodes.push_back(*s);
	    --the_f_vector[f_dim];
	 }
	 else
	    boundary.push_back(*s);
	 
      }

      // compute co_face * boundary(face)
      std::list< Set<int> > new_facets;
      for (Entire< Set<int> >::const_iterator w=entire(face);
	   !w.at_end(); ++w) {
	 Set<int> f=face;
	 f-=*w;
	 f+=co_face;
	 
	 new_facets.push_back(f);
      }
      
      // compute boundary_map and new_star_map
      HasseDiagram<> new_star = pure_hasse_diagram(new_facets);

      std_ext::hash_map< Set<int>, int > boundary_map;
      std_ext::hash_map< int, int > new_star_map;

      for (Entire< std::list<int> >::const_iterator b=entire(boundary);
	   !b.at_end(); ++b)
	 boundary_map[C.node(*b)] = *b;

      new_star_map[new_star.top_node()] = C.top_node();
      new_star_map[new_star.bottom_node()] = C.bottom_node();
      for (int d=0; d<=dim; ++d)
	 for (Entire<sequence>::const_iterator n=entire(new_star.nodes_of_dim(d));
	      !n.at_end(); ++n) {
	    const Set<int> new_face = new_star.node(*n);

	    std_ext::hash_map< Set<int>,int >::const_iterator find=boundary_map.find(new_face);
	    if (find==boundary_map.end()) {  // new_face is not in boundary
	       ++the_f_vector[d];
	       if (d==0) // a new vertex
		  new_star_map[*n] = new_face.front();
	       else
		  if (dead_nodes.empty()) {
		     new_star_map[*n] = C.nodes();
		     C.add_node(new_face);
		  } else {
		     new_star_map[*n] = dead_nodes.front();
		     C.node(dead_nodes.front()) = new_face;
		     dead_nodes.pop_front();
		  }
	    }		  
	    else  // new_face is in boundary
	       new_star_map[*n] = find->second;
	 }
      
      // merge new_star and C according to new_star_map
      for (Entire<HasseDiagram<>::edge_container>::const_iterator e=entire(edges(new_star));
	   !e.at_end(); ++e)
	 C.edge(new_star_map[ e.from_node() ],new_star_map[ e.to_node() ]);
      
      // find new options
      for (int d=0; d<=dim; ++d)
	 for (Entire<sequence>::const_iterator n=entire(new_star.nodes_of_dim(d));
	      !n.at_end(); ++n) {
	    const int f = new_star_map[*n];
	    if (is_option(f,d))
	       options[d].insert(f);
	 }

#ifdef POLY_DEBUG
      C.restore_dim_map();
      bool not_pmf = false;
      for (Entire< std::list<int> >::const_iterator n=entire(C.nodes_of_dim(dim-1));
	   !n.at_end(); ++n)
	 if (C.out_edges(*n).size()!=2) {
	    not_pmf = true;
	    cerr << *n << " = " << C.node(*n) << " is contained in " << C.out_edges(*n).size() << " facets.\n";
	 }
      
      f_vec = f_vector();
      int ee=0;
      for (int i=0; i<f_vec.size(); ++i)
         ee += f_vec[i]*(1-2*(i%2));
      if (ee!=e)
	 cerr << "euler has changed: " << f_vec << " and " << e << " -> " << ee << endl;  
      
      if (ee!=e || not_pmf)
	 throw std::runtime_error("ERROR");
#endif
   }
   
   // Finds minimal revers move >= dim_min and executes it. Return the dimension of the move.
   int min_rev_move(const int dim_min) {
      const int d = find_move(dim_min);
      execute_move();
      
      return d;
   }
   
   int min_rev_move() {
      return min_rev_move(0);
   }
   
   // A zero move.
   void zero_move() {
      min_rev_move(dim);
   }

   // C is unsqueezed and without compleat dim_map !!!
   PseudoComplex pseudo_complex()
     {
	return C;
     }
   
   int n_facets() {
     return C.in_edges(C.top_node()).size();
   }
   
   const Array<int>& f_vector() {
     return the_f_vector;
   }
};
   
} } // end namespace polymake::topaz

#endif // _POLYMAKE_PS_BISTELLAR_COMPLEX_H


// Local Variables:
// mode:C++
// c-basic-offset:3
// End:
