/* Copyright (c) 1997-2006
   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_FACE_LATTICE_TOOLS_H
#define _POLYMAKE_FACE_LATTICE_TOOLS_H "$Project: polymake $$Project: polymake $$Id: face_lattice_tools.h 7315 2006-04-02 21:37:53Z gawrilow $"

#include <list>
#include <PowerSet.h>
#include <IncidenceMatrix.h>
#include <FaceMap.h>

namespace polymake { namespace polytope { namespace face_lattice {

using namespace pm;

/// \c first - all facets containing the given face
/// \c second - intersection of these facets = the enclosing face
typedef std::pair< Set<int>, Set<int> > face_closure;

/** Compute the closure of a face \c S
    @retval first all facets containing \c S
    @retval second intersection of these facets = the enclosing face
*/
template <typename _Set, typename _Matrix>
face_closure
closure(const GenericSet<_Set,int>& S, const GenericIncidenceMatrix<_Matrix>& M) {
   // intersection of all facets incident to vertices from S
   const Set<int> facets=accumulate(cols(M.minor(All,S)), operations::mul());

   // intersection of all these facets
   return face_closure(facets, accumulate(rows(M.minor(facets,All)), operations::mul()));
}


/// Compute a base for a given face H
template <typename _Set, typename _Matrix>
Set<int> c(const GenericSet<_Set,int>& H, const GenericIncidenceMatrix<_Matrix>& M) {
   if (H.top().empty())
      return Set<int>();

   typename Entire<_Set>::const_iterator H_i=entire(H);
   Set<int> C=scalar2set(*H_i);     // c(H) to be returned at the end
   Set<int> facets=M.col(*H_i);
   ++H_i;
   while (!H_i.at_end()) {
      int size = facets.size();
      facets *= M.col(*H_i);
      if (size > facets.size())
	 C.push_back(*H_i);
      ++H_i;
   }
   return C;
}


/// iterate over the (k+1)-faces containing a given k-face H
template <typename _Set, typename _Matrix>
class faces_one_above_iterator {
public:
   typedef forward_iterator_tag iterator_category;
   typedef face_closure value_type;
   typedef const value_type& reference;
   typedef const value_type* pointer;
   typedef ptrdiff_t difference_type;

   faces_one_above_iterator() { }

   faces_one_above_iterator(const GenericSet<_Set,int>& H_arg, const GenericIncidenceMatrix<_Matrix>& M_arg)
      : H(&H_arg), M(&M_arg), candidates(sequence(0,M->cols()) - *H), done(false)
   {
      find_next();
   }

   reference operator* () const { return result; }
   pointer operator-> () const { return &result; }

   faces_one_above_iterator& operator++ () { find_next(); return *this; }

   const faces_one_above_iterator operator++ (int) { faces_one_above_iterator copy=*this; operator++(); return copy; }

   bool operator== (const faces_one_above_iterator& it) const { return candidates==it.candidates; }
   bool operator!= (const faces_one_above_iterator& it) const { return !operator==(it); }
   bool at_end() const { return done; }

   void rewind()
   {
      candidates=sequence(0,M->cols()) - *H;
      minimal.clear();
      done=false;
      find_next();
   }
protected:
   void find_next()
   {
      while (!candidates.empty()) {
	 int v=candidates.front();  candidates.pop_front();
	 result=closure(*H+v, *M);
	 if (result.first.empty()) continue;	// the closure would be the whole polytope - absolutely dull
	 if ((result.second * candidates).empty() && (result.second * minimal).empty()) {
	    minimal.push_back(v);
	    return;
	 }
      }
      done=true;
   }

   const GenericSet<_Set> *H;
   const GenericIncidenceMatrix<_Matrix> *M;
   Set<int> candidates, minimal;
   face_closure result;
   bool done;
};

template <typename _Set, typename _Matrix> inline
faces_one_above_iterator<_Set,_Matrix>
all_faces_one_above(const GenericSet<_Set,int>& H, const GenericIncidenceMatrix<_Matrix>& M)
{
   return faces_one_above_iterator<_Set,_Matrix>(H,M);
}

/// Compute the face lattice
template <typename _Matrix, typename DiagrammFiller>
void compute(const GenericIncidenceMatrix<_Matrix>& VIF, DiagrammFiller HD)
{
   std::list< Set<int> > Q;    // queue of faces, which have been seen but who's faces above have not been computed yet.
   FaceMap<> Faces;

   // The bottom node: empty set (or dually: whole polytope)
   const int R=VIF.rows(), C=VIF.cols();
   HD.add_node(sequence(0,0), sequence(0,R));
   HD.increase_dim();
   int n, end_this_dim=0;

   if (C>1) {
      // The first level: vertices.
      std::copy(entire(all_subsets_of_1(sequence(0,C))), back_inserter(Q));
      n=HD.add_nodes(C, all_subsets_of_1(sequence(0,C)).begin(), cols(VIF).begin());
      end_this_dim=n+C;
      int end_next_dim=end_this_dim;
      HD.increase_dim();
      for (int i=n; i<end_this_dim; ++i)
	 HD.add_edge(0,i);

      if (C>2) {
	 bool facets_reached=false;
	 while (1) {
	    Set<int> H = Q.front(); Q.pop_front();
	    faces_one_above_iterator< Set<int>, _Matrix > faces(H, VIF);
	    while (!faces.at_end()) {
	       if (faces->first.size() == 1) {	// we have reached the facet level
		  if (!facets_reached) {
		     HD.add_nodes(R, rows(VIF).begin(), all_subsets_of_1(sequence(0,R)).begin());
		     facets_reached=true;
		  }
		  HD.add_edge(n, end_this_dim+faces->first.front());
	       } else {
		  int &node_ref = Faces[c(faces->second, VIF)];
		  if (node_ref==-1) {
		     node_ref=HD.add_node(faces->second, faces->first);
		     Q.push_back(faces->second);
		     ++end_next_dim;
		  }
		  HD.add_edge(n,node_ref);
	       }
	       ++faces;
	    }
	    if (++n == end_this_dim) {
	       HD.increase_dim();
	       if (Q.empty()) break;
	       end_this_dim = end_next_dim;
	    }
	 }
      } else {
	 // dim==1
	 HD.add_nodes(R, rows(VIF).begin(), all_subsets_of_1(sequence(0,R)).begin());
	 HD.add_edge(n, end_this_dim+VIF.col(0).front());
	 HD.add_edge(++n, end_this_dim+VIF.col(1).front());
	 HD.increase_dim();
      }
   }

   // The top node: whole polytope (or dually: empty set)
   n=HD.add_node(sequence(0,C), sequence(0,0));
   for (int i=end_this_dim; i<n; ++i)
      HD.add_edge(i,n);
}

} } } // end namespace face_lattice

namespace pm {

template <typename _Set, typename _Matrix>
struct check_iterator_feature<polymake::polytope::face_lattice::faces_one_above_iterator<_Set,_Matrix>, end_sensitive>
   : True { };
template <typename _Set, typename _Matrix>
struct check_iterator_feature<polymake::polytope::face_lattice::faces_one_above_iterator<_Set,_Matrix>, rewindable>
   : True { };

} // end namespace pm

#endif // _POLYMAKE_FACE_LATTICE_TOOLS_H

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