/* 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_BENEATH_BEYOND_TCC
#define _POLYMAKE_BENEATH_BEYOND_TCC "$Project: polymake $$Id: beneath_beyond.tcc 7315 2006-04-02 21:37:53Z gawrilow $"

namespace polymake { namespace polytope {

template <typename E> template <typename Iterator>
void beneath_beyond_algo<E>::compute(Iterator perm)
{
   if (perm.at_end()) return;

   // the first point
   int d2=points.cols()-2, p1=*perm;  ++perm;
   null_space(entire(item2container(points[p1])), black_hole<int>(), black_hole<int>(), AH);

   // look for the second point different from the first one
   while (true) {
      if (perm.at_end()) {
	 // the special case: a single point
	 triang_size=1;
	 triangulation.push_back(scalar2set(p1));
	 return;
      }
      int p2=*perm;  ++perm;
      null_space(entire(item2container(points[p2])), black_hole<int>(), black_hole<int>(), AH);
      if (AH.rows()==d2) {
	 // two different points found: initialize the polytope
	 start_with_points(p1,p2);
	 break;
      }
      if (!already_VERTICES) interior_points += p2;
   }

   // as long as the affine hull is not empty...
   while (AH.rows() && !perm.at_end()) {
      add_point_low_dim(*perm);
#ifdef POLY_DEBUG
      if (debug==do_check) {
	 check_p(*perm);
	 if (!pm::identical<Iterator,sequence::iterator>::value) points_so_far+=*perm;
      }
#endif
      ++perm;
   }

   // then take the shortcut
   while (!perm.at_end()) {
      add_point_full_dim(*perm);
#ifdef POLY_DEBUG
      if (debug==do_check) {
	 check_p(*perm);
	 if (!pm::identical<Iterator,sequence::iterator>::value) points_so_far+=*perm;
      }
#endif
      ++perm;
   }

   if (!facet_normals_valid) {
      facet_normals_low_dim();
      facet_normals_valid=true;
   }

   // sweep out the unneeded facets from the recycling list
   dual_graph.squeeze();
   if (dual_graph.nodes()==2)
      dual_graph.delete_edge(0,1);	// special case of dimension==1: the dual graph is not connected
   free_facet=-1;
#ifdef POLY_DEBUG
   active_facets=sequence(0,dual_graph.nodes());
   cout << "final "; dump();
#endif
}

template <typename E>
void beneath_beyond_algo<E>::start_with_points(int p1, int p2)
{
   new_facet();
   dual_graph.node(0).vertices=scalar2set(p1);
   new_facet();
   dual_graph.node(1).vertices=scalar2set(p2);
   dual_graph.edge(0,1);
   vertices_so_far=scalar2set(p1)+scalar2set(p2);
   triangulation.push_back(vertices_so_far);
   triang_size=1;
   dual_graph.node(0).simplices.push_back(incident_simplex(triangulation.front(),p2));
   dual_graph.node(1).simplices.push_back(incident_simplex(triangulation.front(),p1));
   valid_facet=0;
#ifdef POLY_DEBUG
   cout << "starting points: " << p1 << ' ' << p2 << "\nAH:\n" << AH << endl;
#endif
   if ((facet_normals_valid=!AH.rows())) {	// dimension==1, will need the facet normals immediately
      dual_graph.node(0).coord_full_dim(*this);
      dual_graph.node(1).coord_full_dim(*this);
   }
}

/** @param p the new point
    @param f the facet index to start from
    @retval index of the violated/incident facet or -1 if nothing found
*/
template <typename E>
int beneath_beyond_algo<E>::descend_to_violated_facet(int f, int p)
{

   visited_facets+=f;
   E fxp= dual_graph.node(f).normal * points[p];
   if ((dual_graph.node(f).orientation=sign(fxp))<=0) return f;

   // starting facet stays valid in this step: let's look for another one violated by p.
   // The search is performed in the dual graph, following the steepest descend of the
   // (square of the) distance between p and the facets

   if (!already_VERTICES) vertices_this_step += dual_graph.node(f).vertices;
   fxp=fxp*fxp/dual_graph.node(f).sqr_normal;	// square of the distance from p to the facet

   int nextf;
   do {
#ifdef POLY_DEBUG
      cout << " *" << f << '(' << fxp << ')';
#endif
      nextf=-1;
      for (typename Entire<typename DualGraph::neighbor_list>::iterator neighbor=entire(dual_graph.neighbors(f));
	   !neighbor.at_end(); ++neighbor) {
	 int f2=*neighbor;
	 if (visited_facets.contains(f2)) continue;

	 visited_facets+=f2;
	 E f2xp= dual_graph.node(f2).normal * points[p];
	 if ((dual_graph.node(f2).orientation=sign(f2xp))<=0) return f2;

	 if (!already_VERTICES) vertices_this_step += dual_graph.node(f2).vertices;
	 f2xp=f2xp*f2xp/dual_graph.node(f2).sqr_normal;
#ifdef POLY_DEBUG
	 cout << ' ' << f2 << '(' << f2xp << ')';
#endif
	 if (f2xp<=fxp) {
	    nextf=f2;
	    fxp=f2xp;
	 }
      }
   } while ((f=nextf)>=0);

   return f;	// -1 : local minimum of sqr(distance) reached
}

template <typename Set> static inline
int first_or_none(const Set& set)
{
   typename Entire<Set>::const_iterator s=entire(set);
   return s.at_end() ? -1 : *s;
}

template <typename E>
void beneath_beyond_algo<E>::add_point_full_dim(int p)
{
#ifdef POLY_DEBUG
   cout << "point " << p << "=[ " << points[p] << " ] : valid facets";
#endif
   // reset the working variables
   visited_facets.clear();
   if (!already_VERTICES) vertices_this_step.clear();

   // first try the facet added last in the previous step
   int f=valid_facet;
   do {
      if ((f=descend_to_violated_facet(f,p))>=0) {
	 update_facets(f,p);
	 return;
      }
      f=first_or_none(active_facets-visited_facets);
   } while (f>=0);

   // no violated facet found: p must be a redundant point

   if (!already_VERTICES) {
      interior_points += p;
#ifdef POLY_DEBUG
      if (debug>=do_dump)
	 cout << "\ninterior points: " << interior_points
	      << "\n=======================================" << endl;
#endif
   }
}

template <typename E>
void beneath_beyond_algo<E>::facet_normals_low_dim()
{
   // facets must be orthogonal to the affine hull
   facet_nullspace=unit_matrix<E>(points.cols());
   null_space(entire(rows( zero_vector<E>(AH.rows()) | AH.minor(All,range(1,points.cols()-1)) )),
	      black_hole<int>(), black_hole<int>(), facet_nullspace);

   for (typename Entire<ActiveFacets>::iterator f=entire(select(nodes(dual_graph), active_facets));  !f.at_end();  ++f) {
      f->coord_low_dim(*this);
#ifdef POLY_DEBUG
      if (debug>=do_dump) cout << f.index() << ": =[" << f->normal << " ]\n";
#endif
   }
}

template <typename E>
void beneath_beyond_algo<E>::add_point_low_dim(int p)
{
   // update the affine hull
   int codim=AH.rows();
   null_space(entire(item2container(points[p])), black_hole<int>(), black_hole<int>(), AH);

   if (AH.rows()<codim) {
      // point set dimension increased
      if (facet_nullspace.rows()) {
	 generic_position=false;	// the base facet is more than a simplex
	 facet_nullspace.clear();
      }

      // build a pyramid with the former polytope as a base and the point as an apex
      int nf_index=new_facet();
      facet_info& nf=dual_graph.node(nf_index);
      nf.vertices=vertices_so_far;
      vertices_so_far+=p;

      // triangulation simplices are 'pyramidized' too
      for (typename Entire<Triangulation>::iterator simplex=entire(triangulation); !simplex.at_end(); ++simplex) {
	 *simplex += p;
	 nf.simplices.push_back(incident_simplex(*simplex,p));
      }

      for (typename Entire< Edges<DualGraph> >::iterator e=entire(edges(dual_graph)); !e.at_end(); ++e)
	 *e += p;
      facet_normals_valid=!AH.rows();
      for (typename Entire<ActiveFacets>::iterator f=entire(select(nodes(dual_graph), active_facets));  !f.at_end();  ++f) {
	 // for all facets, except the new one
	 if (f.index() != nf_index) {
	    f.out_edges().insert(nf_index,f->vertices);
	    f->vertices+=p;
	 }
	 if (facet_normals_valid)
	    // the polytope became full-dimensional: will need the facet coordinates the whole rest of the time
	    f->coord_full_dim(*this);
      }
#ifdef POLY_DEBUG
      cout << "point " << p << ": dim increased\nAH:\n" << AH << endl;
#endif
   } else {
      // point set dimension not increased
      if (!facet_normals_valid) {
	 // the polytope was a simplex, the facet coordinates are still not computed;
	 // now we need them for the visibility region search.
	 facet_normals_low_dim();
	 facet_normals_valid=true;
      }
      add_point_full_dim(p);
   }
}

template <typename E> inline
void beneath_beyond_algo<E>::update_facets(int f, int p)
{
#ifdef POLY_DEBUG
   cout << "\nupdating:";
#endif
   // queue for BFS
   std::list<int> Q;
   Q.push_back(f);
   if (!already_VERTICES) interior_points_this_step.clear();
   std::list<int> incident_facets;
   if (dual_graph.node(f).orientation==0) {
      dual_graph.node(f).vertices += p;
      generic_position=false;
      incident_facets.push_back(f);
   }

   /* BFS in the visible hemisphere.
      We visit all facets violated by or incident with p.
      Incident facets are important since they can contain redundant points not discovered before this iteration.
   */
   while (!Q.empty()) {
      f=Q.front(); Q.pop_front();
      int f_orientation=dual_graph.node(f).orientation;
      // remember the position where the new simplices will end
      typename Triangulation::iterator new_simplex_end=triangulation.begin();

      if (f_orientation<0) {
	 // the facet is violated
#ifdef POLY_DEBUG
	 cout << " -" << f;
#endif
	 if (!already_VERTICES) interior_points_this_step += dual_graph.node(f).vertices;

	 // build new triangulation simplices using the triangulation of the facet
	 for (typename Entire<typename facet_info::simplex_list>::iterator is=entire(dual_graph.node(f).simplices);
	      !is.at_end(); ++is) {
	    triangulation.push_front(*is->simplex);
	    ++triang_size;
	    // just take the existing simplex and replace the vertex behind the facet by the new point
	    (triangulation.front() -= is->opposite_vertex) += p;
	 }
#ifdef POLY_DEBUG
      } else {
	 cout << " " << f;
#endif
      }

      // check the neighbor facets
      for (typename Entire<typename DualGraph::out_edge_list>::iterator e=entire(dual_graph.out_edges(f));
	   !e.at_end(); ++e) {
	 int f2=e.to_node();
	 facet_info& nbf=dual_graph.node(f2);
	 if (!visited_facets.contains(f2)) {
	    visited_facets+=f2;
	    nbf.orientation=sign(nbf.normal * points[p]);
	    if (nbf.orientation==0) {
	       // incident facet
	       nbf.vertices += p;
	       generic_position=false;
	       incident_facets.push_back(f2);
	    }
	    if (nbf.orientation<=0)
	       Q.push_back(f2);
	    else if (!already_VERTICES)
	       vertices_this_step += nbf.vertices;
	 }

	 if (f_orientation<0) {
	    if (nbf.orientation>0) {
	       // found a ridge on the visibility border: create a new facet
	       int nf_index=new_facet();
	       facet_info& nf=dual_graph.node(nf_index);
	       nf.vertices=*e + p;
	       if (AH.rows())
		  nf.coord_low_dim(*this);
	       else
		  nf.coord_full_dim(*this);
#ifdef POLY_DEBUG
	       if (debug==do_check) check_f(nf_index, p);
#endif
	       dual_graph.edge(nf_index,f2)=*e;
	       incident_facets.push_back(nf_index);
	       nf.add_incident_simplices(triangulation.begin(), new_simplex_end);
	    } else if (nbf.orientation==0) {
	       nbf.add_incident_simplices(triangulation.begin(), new_simplex_end);
	    }
	 } else if (nbf.orientation==0) {
	    *e += p;	// include the point into the edge, since it's incident to both facets
	 }
      }

      if (f_orientation<0) delete_facet(f);
   }

   if (!already_VERTICES) {
      if (interior_points_this_step.empty()) { // = no violated facets visited
	 interior_points += p;
#ifdef POLY_DEBUG
	 if (debug>=do_dump)
	    cout << "\ninterior points: " << interior_points
		 << "\n=======================================" << endl;
#endif
	 return;
      } else {
	 interior_points_this_step -= vertices_this_step;
	 interior_points += interior_points_this_step;
      }
   }

   valid_facet=incident_facets.front();
   update_dual_graph(select(nodes(dual_graph),incident_facets));
   if (AH.rows()) vertices_so_far+=p;
#ifdef POLY_DEBUG
   if (debug>=do_dump) {
      dump();
      cout << "\ninterior points: " << interior_points
	   << "\n=======================================" << endl;
   }
#endif
}

template <typename E> template <typename FacetSubset>
void beneath_beyond_algo<E>::update_dual_graph (FacetSubset facets)
{
   int min_ridge=points.cols()-AH.rows()-2;

   for (typename Entire<FacetSubset>::iterator f=entire(facets); !f.at_end(); ++f) {
      const bool vis=visited_facets.contains(f.index());

      typename Entire<FacetSubset>::iterator f2=f;
      for (++f2; !f2.at_end(); ++f2) {
	 // if both facets are incident to p, they could already have a connecting edge
	 if (vis && visited_facets.contains(f2.index()) && f.out_edges().exists(f2.index())) continue;

	 const Set<int> ridge= f->vertices * f2->vertices;
	 if (ridge.size()>=min_ridge) {
	    bool add=true;
	    typename Entire<typename DualGraph::out_edge_list>::iterator e=entire(f.out_edges());
	    while (!e.at_end()) {
	       int inc=incl(*e,ridge);
	       if (inc==2) {
		  ++e;
	       } else {
		  if (inc<=0) f.out_edges().erase(e++);
		  if (inc>=0) { add=false; break; }
	       }
	    }
	    if (add) f.out_edges().insert(f2.index(),ridge);
	 }
      }
   }
}

template <typename E> inline
void beneath_beyond_algo<E>::facet_info::coord_full_dim (const beneath_beyond_algo<E>& A)
{
   normal=rows(null_space(A.points.minor(vertices,All))).front();
   if (normal * A.points[(A.vertices_so_far - vertices).front()] < 0) inv_sign(normal);
   sqr_normal=sqr(normal);
}

template <typename E>
void beneath_beyond_algo<E>::facet_info::coord_low_dim (const beneath_beyond_algo<E>& A)
{
   ListMatrix< SparseVector<E> > Fn=A.facet_nullspace;
   null_space(entire(rows(A.points.minor(vertices,All))), black_hole<int>(), black_hole<int>(), Fn);
   normal=rows(Fn).front();
   if (normal * A.points[(A.vertices_so_far - vertices).front()] < 0) inv_sign(normal);
   sqr_normal=sqr(normal);
}

template <typename Top> inline
int single_or_nothing(const GenericSet<Top,int>& s)
{
   int x=-1;
   typename Entire<Top>::const_iterator e=entire(s);
   if (!e.at_end()) {
      x=*e; ++e;
      if (!e.at_end()) x=-1;
   }
   return x;
}

template <typename E> template <typename Iterator> inline
void beneath_beyond_algo<E>::facet_info::add_incident_simplices(Iterator s, Iterator s_end)
{
   for (; s != s_end; ++s) {
      int opv=single_or_nothing(*s-vertices);
      if (opv>=0) simplices.push_back(incident_simplex(*s,opv));
   }
}

template <typename E>
int beneath_beyond_algo<E>::new_facet()
{
   int f=free_facet;
   if (f>=0) {
      free_facet=dual_graph.node(f).orientation;
   } else {
      f=dual_graph.nodes();
      dual_graph.resize(f+1);
   }
   active_facets+=f;
   return f;
}

template <typename E>
void beneath_beyond_algo<E>::delete_facet(int f)
{
   active_facets-=f;
   facet_info& F=dual_graph.node(f);
   dual_graph.out_edges(f).clear();
   F.normal.clear(); F.vertices.clear(); F.simplices.clear();
   F.orientation=free_facet;
   free_facet=f;
}

#ifdef POLY_DEBUG
template <typename E>
void beneath_beyond_algo<E>::dump() const
{
   cout << "dual_graph:\n";
   const bool show_normals= debug==do_full_dump && (!AH.rows() || facet_nullspace.rows());
   for (typename Entire<ActiveFacets>::const_iterator f=entire(select(nodes(dual_graph), active_facets)); !f.at_end(); ++f) {
      cout << f.index() << ": " << f->vertices;
      if (show_normals) cout << "=[ " << f->normal << " ]";
      if (debug==do_full_dump) {
	 cout << ' ' << f.out_edges();
	 cout << " <<";
	 for (typename Entire<typename facet_info::simplex_list>::const_iterator s=entire(f->simplices); !s.at_end(); ++s)
	 cout << ' ' << *s->simplex << '-' << s->opposite_vertex;
	 cout << " >>" << endl;
      } else {
	 cout << ' ' << f.neighbors() << endl;
      }
   }
}

template <typename E> inline
void beneath_beyond_algo<E>::check_fp(int f_index, const facet_info& f, int p, std::ostringstream& errors) const
{
   E prod= points[p] * f.normal;
   if (f.vertices.contains(p)) {
      if (prod)
	 errors << "facet(" << f_index << ") * incident vertex(" << p << ")=" << prod << endl;
   } else {
      if (prod<=0)
	 errors << "facet(" << f_index << ") * non-incident vertex(" << p << ")=" << prod << endl;
   }
}

// various consistency checks
template <typename E>
void beneath_beyond_algo<E>::check_p(int p) const
{
   if (!AH.rows() || facet_nullspace.rows()) {
      std::ostringstream errors;

      for (typename Entire<ActiveFacets>::const_iterator f=entire(select(nodes(dual_graph), active_facets)); !f.at_end(); ++f)
	 check_fp(f.index(), *f, p, errors);

      if (!errors.str().empty())
	 throw std::runtime_error("beneath_beyond_algo - consistency checks failed:\n" + errors.str());
   }
}

template <typename E>
void beneath_beyond_algo<E>::check_f(int f, int last_p) const
{
   std::ostringstream errors;
   const facet_info& fi=dual_graph.node(f);

   if (points_so_far.empty()) {
      for (typename Entire<sequence>::const_iterator p=entire(range(0,last_p)); !p.at_end(); ++p)
	 check_fp(f, fi, *p, errors);
   } else {
      for (typename Entire<Bitset>::const_iterator p=entire(points_so_far); !p.at_end(); ++p)
	 check_fp(f, fi, *p, errors);
   }

   if (!errors.str().empty())
      throw std::runtime_error("beneath_beyond_algo - consistency checks failed:\n" + errors.str());
}

template <typename E>
void beneath_beyond_algo<E>::dump_p(int p) const
{
   if (!AH.rows() || facet_nullspace.rows()) {
      for (typename Entire<ActiveFacets>::const_iterator f=entire(select(nodes(dual_graph), active_facets)); !f.at_end(); ++f)
	 if (f.degree()) {
	    E prod= points[p] * f->normal;
	    cout << "facet(" << f.index() << "): prod=" << prod << ", sqr_dist=" << double(prod*prod/f->sqr_normal) << '\n';
	 }
   }
}
#endif // POLY_DEBUG

} }

#endif // _POLYMAKE_BENEATH_BEYOND_TCC

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