/* 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_BISTELLAR_COMPLEX_H
#define _POLYMAKE_BISTELLAR_COMPLEX_H "$Project: polymake $$Id: BistellarComplex.h 7579 2007-01-22 11:34:46Z gawrilow $"

#include <complex_tools.h>
#include <FacetList.h>
#include <RandomSubset.h>
#include <IncidenceMatrix.h>

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

namespace polymake { namespace topaz {
   
class BistellarComplex {
   
protected:
   typedef std::pair< Set<int>,Set<int> > option;
   
   class OptionsList {
      
   protected:
      int the_size;
      std_ext::hash_map< Set<int>,int > index_map;
      Array<option> the_options;

   public:
      OptionsList() : the_size(0) {};

      int size()
      {
	 return the_size;
      }

      void insert(const option& opt)
      {
	 if (the_options.size()==0)
	    the_options.resize(1);
	 if (the_size>=the_options.size())
	    the_options.resize(2*the_options.size());
	 the_options[the_size]=opt;
	 index_map[opt.first]=the_size;
	 ++the_size;
      }
      
      void insert(const Set<int>& f, const Set<int>& V)
      {
	 insert(option(f,V));
      }
      
      void remove(const Set<int>& f)
      {
	 std_ext::hash_map< Set<int>,int >::const_iterator find=index_map.find(f);
	 if (find!=index_map.end()) {
	    const int index = find->second;
	    the_options[index] = the_options[the_size-1];
	    index_map[the_options[the_size-1].first] = index;
	    index_map.erase(f);
	    --the_size;
	 }
      }
      
      Array<option> options() {
	 return Array<option>(the_size,the_options.begin());
      }

      /*
	void print_options()
	{
	for (int i=0; i<the_size; ++i)
	cout << "(" << the_options[i].first << "," << the_options[i].second << ") ";
	cout << endl;
	}
      */
   };
   
   FacetList the_facets;
   Array<OptionsList> raw_options;
   int dim;
   unsigned long seed;
   int verts;
   option next_move;
   Set<int> rev_move;
   bool allow_rev_move;
   Array<int> the_flip_vector; 
   bool verbose;
   int apex;
   bool closed;

public:
   BistellarComplex() : dim(0), seed(0), verts(0), allow_rev_move(false), verbose(false), apex(0), closed(true) {};

   BistellarComplex(const HasseDiagram<>& HD, const unsigned long in_seed,
		    const bool verb=false, const bool in_closed=false, const bool in_allow_rev_move=false) {
      verbose = verb;
      dim = HD.dim()-1;
      seed = in_seed;
      verts = 0;
      apex = 0;
      next_move = option();
      rev_move = Set<int>();
      allow_rev_move = in_allow_rev_move;
      the_flip_vector.resize((dim+1)/2);
      closed = in_closed;
      raw_options.resize(dim+1);

      // test if the complex is closed.
      if (!closed) {
	 const polymake::topaz::Boundary_of_PseudoManifold B = polymake::topaz::boundary_of_pseudo_manifold(HD);
	 bool closed = B.empty();

	 if (!closed) {

	    // compute C + cone(bound(C))
	    std::list< Set<int> > S;
	    
	    for (Entire<sequence>::const_iterator f=entire(HD.nodes_of_dim(HD.dim()-1));
		 !f.at_end(); ++f) {
	       S.push_back(HD.node(*f));
	       const int w = HD.node(*f).back();
	       if (w >= verts)
		  verts = w+1;
	    }
	    apex = verts;	    
	    ++verts;

	    for (Entire<polymake::topaz::Boundary_of_PseudoManifold>::const_iterator b=entire(B);
		 !b.at_end(); ++b)
	       S.push_back(*b+apex);
	    const HasseDiagram<> new_HD = polymake::topaz::pure_hasse_diagram(S);

	    // compute raw options
	    for (int d=0; d<=dim; ++d)
	       for (Entire<sequence>::const_iterator n=entire(new_HD.nodes_of_dim(d));
		    !n.at_end(); ++n) {
		  const Set<int> face = new_HD.node(*n);
		  
		  if (d==0 && face.front() == apex)  // the apex is not an option
		     continue;

		  if (d==dim) {  // each facet is an option
		     the_facets.insert(face);
		     raw_options[d].insert(face, Set<int>());
		     
		  } else {
		     const FacetList link = polymake::topaz::link_in_HD(new_HD,*n);
		     const Set<int> V = accumulate(link, operations::add());
		     
		     if (V.size()+face.size()==dim+2)  // face is raw option
			raw_options[d].insert(face,V);
		  }
	       }
	 }
      }
      
      if (closed)
	 for (int d=0; d<=dim; ++d)
	    for (Entire<sequence>::const_iterator n=entire(HD.nodes_of_dim(d));
		 !n.at_end(); ++n) {
	       const Set<int> face=HD.node(*n);
	       
	       if (d==0) {    // face is a vertex
		  const int v = face.front();
		  if (v >= verts)
		     verts = v+1;
	       }
	       
	       if (d==dim) {  // each facet is an option
		  the_facets.insert(face);
		  raw_options[d].insert(face, Set<int>());
		  
	       } else {
		  const FacetList link = polymake::topaz::link_in_HD(HD,*n);
		  const Set<int> V = accumulate(link, operations::add());
		  
		  if (V.size()+face.size()==dim+2)   // face is raw option
		     raw_options[d].insert(face,V);
	       }
	    }
   }
   
   int n_raw_options_of_dim(const int d) {
      return raw_options[d].options().size();
   }

   // Finds a reversed move of minimal dimension >= dim_min
   // and returns the dimension of the move. The move is != next_move.second
   // The move is stored in next_move.
   int find_move(const int dim_min) {
      for (int d=dim_min; d<=dim; ++d) {
	 const RandomPermutation< Array<option> > P(raw_options[d].options(),seed++);

	 for (Entire< RandomPermutation< Array<option> > >::const_iterator
		 opt=(entire(P));
	      !opt.at_end(); ++opt)
	    if ( (allow_rev_move || incl(opt->first,rev_move)!=0) &&
		 (d==dim || the_facets.findMax(opt->second).at_end()) ) {
	       next_move = *opt;
	       return opt->first.size() - 1;
	 }
      }
      
      throw std::runtime_error("BistellarComplex: No move found.");
   }
   
   int find_move() {
      return find_move(0);
   }

protected:
   // Determines if a face is an option.
   bool is_option(const Set<int>& f, Set<int>& V) {
      if (!closed && f.size()==1 && f.front()==apex)  // apex is not an option
	 return false;

      for (FacetList::iteratorMax star = the_facets.findMax(f);
	   !star.at_end(); ++star)
	 V += *star;
      V -= f;

      return V.size()+f.size()==dim+2;
   }
   
public:
   // 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()
   {
      const Set<int> face = next_move.first;
      const int face_dim = face.size()-1;
      if (face_dim==dim)  // allocate an index for the new vertex
	 next_move.second = scalar2set(verts++);
      const Set<int> co_face= next_move.second;
      if (!allow_rev_move)  rev_move = co_face;      

      if (verbose)
	 cout << "BistellarComplex: executing move of dim "
	      << face.size()-1 << ": (" << face << "," << co_face << ")\n";

      // update the_flip_vector
      if (dim%2==1 || face_dim!=dim/2)
	 if (face_dim < (dim+1)/2)
	    --the_flip_vector[face_dim];
	 else
	    ++the_flip_vector[dim-face_dim];

      // remove star(face) from raw_options and the_facets
      std::list< Set<int> > star;
      the_facets.eraseMin(face, std::back_inserter(star));

      HasseDiagram<> star_HD = polymake::topaz::pure_hasse_diagram(star);
      for (int d=0; d<=dim; ++d)
	 for (Entire<sequence>::const_iterator n=entire(star_HD.nodes_of_dim(d));
	      !n.at_end(); ++n)
	    raw_options[d].remove(star_HD.node(*n));
      
      // add 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;
	 
	 the_facets.insert(f);
	 new_facets.push_back(f);
      }
      
      // find new raw_options
      HasseDiagram<> lokal_HD = polymake::topaz::pure_hasse_diagram(new_facets);
      for (int d=0; d<=dim; ++d)
	 for (Entire<sequence>::const_iterator n=entire(lokal_HD.nodes_of_dim(d));
	      !n.at_end(); ++n) {
	    const Set<int> f = lokal_HD.node(*n);
	    
	    if (d==dim)
	       raw_options[d].insert(f, Set<int>());
	    
	    else {
	       Set<int> V;
	       if( is_option(f,V) )
		  raw_options[d].insert(f,V);
	    }
	 }
   }
   
   // 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);
   }
   
   // The facets. Vertices are NOT numbered according to the topaz standard.
   FacetList facets() {
      if (closed)
	 return the_facets;

      // remove star(apex)
      FacetList real_facets = the_facets;
      real_facets.eraseMin(scalar2set(apex));

      return real_facets;
   }
  
  IncidenceMatrix<> as_incid_matrix() {
      FacetList F=facets();
      F.squeeze();
      return IncidenceMatrix<>(F.size(), F.cols(), F.begin());
   }
  
   int n_facets() {
      if (closed)
	 return the_facets.size();

      return facets().size();
   }
   
   const Array<int>& flip_vector() {
      return the_flip_vector;
   }
};
   
} } // end namespace polymake::topaz

#endif // _POLYMAKE_BISTELLAR_COMPLEX_H


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