#include "polyhedralfan.h"
#include "reversesearch.h"
#include "wallideal.h"
#include "buchberger.h"
#include "printer.h"
#include "timer.h"
#include "symmetry.h"

static Timer polyhedralFanRefinementTimer("Polyhedral fan refinement",1);

PolyhedralFan::PolyhedralFan(int ambientDimension):
  n(ambientDimension)
{
}


PolyhedralFan PolyhedralFan::fullSpace(int n)
{
  PolyhedralFan ret(n);

  ret.cones.insert(PolyhedralCone(n));

  return ret;
}


PolyhedralFan PolyhedralFan::facetsOfCone(PolyhedralCone const &c)
{
  PolyhedralCone C(c);
  C.canonicalize();
  PolyhedralFan ret(C.ambientDimension());

  IntegerVectorList halfSpaces=C.getHalfSpaces();

  for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
    {
      IntegerVectorList a;
      IntegerVectorList b;
      b.push_back(*i);
      PolyhedralCone n=intersection(PolyhedralCone(a,b),c);
      n.canonicalize();
      ret.cones.insert(n);
    }
  return ret;
}


PolyhedralFan PolyhedralFan::bergmanOfPrincipalIdeal(Polynomial const &p1)
{
  Polynomial p=p1.homogenization();
  PolyhedralFan ret(p1.getNumberOfVariables());
  PolynomialSet g;
  g.push_back(p);
  buchberger(&g,LexicographicTermOrder());

  EnumerationTargetCollector gfan;
  LexicographicTermOrder myTermOrder;
  ReverseSearch rs(myTermOrder);

  rs.setEnumerationTarget(&gfan);

  rs.enumerate(g);

  PolynomialSetList theList=gfan.getList();
  for(PolynomialSetList::const_iterator i=theList.begin();i!=theList.end();i++)
    {
      //      AsciiPrinter(Stderr).printPolynomialSet(*i);
      IntegerVectorList inequalities=wallInequalities(*i);
      IntegerVectorList f=wallFlipableNormals(*i,true);
      for(IntegerVectorList::const_iterator j=f.begin();j!=f.end();j++)
	{
	  //	  AsciiPrinter(Stderr).printVector(*j);
	  if(myTermOrder(*j,*j-*j))
	    {
	      IntegerVectorList equalities;
	      equalities.push_back(*j);
	      PolyhedralCone c=PolyhedralCone(inequalities,equalities).withLastCoordinateRemoved();
	      c.canonicalize();
	      ret.cones.insert(c);
	    }
	}
    }

  return ret;
}


PolyhedralFan PolyhedralFan::normalFanOfNewtonPolytope(Polynomial const &p1)
{
  Polynomial p=p1.homogenization();
  PolyhedralFan ret(p1.getNumberOfVariables());
  PolynomialSet g;
  g.push_back(p);
  buchberger(&g,LexicographicTermOrder());

  EnumerationTargetCollector gfan;
  LexicographicTermOrder myTermOrder;
  ReverseSearch rs(myTermOrder);

  rs.setEnumerationTarget(&gfan);

  rs.enumerate(g);

  PolynomialSetList theList=gfan.getList();
  for(PolynomialSetList::const_iterator i=theList.begin();i!=theList.end();i++)
    {
      IntegerVectorList inequalities=wallInequalities(*i);
      IntegerVectorList equalities;

      PolyhedralCone c=PolyhedralCone(inequalities,equalities).withLastCoordinateRemoved();
      c.canonicalize();
      ret.cones.insert(c);
    }

  return ret;
}


void PolyhedralFan::print(class Printer *p)const
{
  p->printString("Printing PolyhedralFan");
  p->printNewLine();
  p->printString("Ambient dimension: ");
  p->printInteger(n);
  p->printNewLine();
  p->printString("Number of cones: ");
  p->printInteger(cones.size());
  p->printNewLine();
  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
    {
      p->printNewLine();
      p->printPolyhedralCone(*i);
    }
  p->printString("Done printing PolyhedralFan.");
  p->printNewLine(); 
}

int PolyhedralFan::getAmbientDimension()const
{
  return n;
}

bool PolyhedralFan::isEmpty()const
{
  return cones.empty();
}

int PolyhedralFan::getMaxDimension()const
{
  assert(!cones.empty());

  return cones.begin()->dimension();
}

int PolyhedralFan::getMinDimension()const
{
  assert(!cones.empty());

  return cones.rbegin()->dimension();
}
 
PolyhedralFan refinement(const PolyhedralFan &a, const PolyhedralFan &b, int cutOffDimension, bool allowASingleConeOfCutOffDimension)
{
  TimerScope ts(&polyhedralFanRefinementTimer);
  //  fprintf(Stderr,"PolyhedralFan refinement: #A=%i #B=%i\n",a.cones.size(),b.cones.size());
  int conesSkipped=0;
  int numberOfComputedCones=0;
  bool linealityConeFound=!allowASingleConeOfCutOffDimension;
  assert(a.getAmbientDimension()==b.getAmbientDimension());

  PolyhedralFan ret(a.getAmbientDimension());

  for(PolyhedralConeList::const_iterator A=a.cones.begin();A!=a.cones.end();A++)
    for(PolyhedralConeList::const_iterator B=b.cones.begin();B!=b.cones.end();B++)
      {
	PolyhedralCone c=intersection(*A,*B);
	int cdim=c.dimension();
	//	assert(cdim>=linealitySpaceDimension);
	bool thisIsLinealityCone=(cutOffDimension>=cdim);
	if((!thisIsLinealityCone)||(!linealityConeFound))
	  {
	    numberOfComputedCones++;
	    c.canonicalize();
	    ret.cones.insert(c);
	    linealityConeFound=linealityConeFound || thisIsLinealityCone;
	  }
	else
	  {
	    conesSkipped++;
	  }
      }
  //  fprintf(Stderr,"Number of skipped cones: %i, lineality cone found: %i\n",conesSkipped,linealityConeFound);
  //  fprintf(Stderr,"Number of computed cones: %i, number of unique cones: %i\n",numberOfComputedCones,ret.cones.size());

  return ret;
}

IntegerVectorList PolyhedralFan::getRays(int dim)
{
  IntegerVectorList ret;
  for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
    {
      if(i->dimension()==dim)
	ret.push_back(i->getRelativeInteriorPoint());
    }
  return ret;
}


IntegerVectorList PolyhedralFan::getRelativeInteriorPoints()
{
  IntegerVectorList ret;
  for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
    {
      ret.push_back(i->getRelativeInteriorPoint());
    }
  return ret;
}


PolyhedralCone const& PolyhedralFan::highestDimensionalCone()
{
  return *cones.rbegin();
}

void PolyhedralFan::insert(PolyhedralCone const &c)
{
  cones.insert(c);
}

void PolyhedralFan::remove(PolyhedralCone const &c)
{
  cones.erase(c);
}

void PolyhedralFan::removeAllExcept(int a)
{
  PolyhedralConeList::iterator i=cones.begin();
  while(a>0)
    {
      if(i==cones.end())return;
      i++;
      a--;
    }
  cones.erase(i,cones.end());
}

void PolyhedralFan::removeAllLowerDimensional()
{
  if(!cones.empty())
    {
      int d=getMaxDimension();
      PolyhedralConeList::iterator i=cones.begin();
      while(i!=cones.end() && i->dimension()==d)i++;
      cones.erase(i,cones.end());
    }
}


PolyhedralFan PolyhedralFan::facetComplex()const
{
  //  fprintf(Stderr,"Computing facet complex...\n");
  PolyhedralFan ret(n);

  for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
    {
      PolyhedralFan a=facetsOfCone(*i);
      for(PolyhedralConeList::const_iterator j=a.cones.begin();j!=a.cones.end();j++)
	ret.insert(*j);
    }
  //  fprintf(Stderr,"Done computing facet complex.\n");
  return ret;
}


PolyhedralFan PolyhedralFan::facetComplexSymmetry(SymmetryGroup const &sym)const
{
  fprintf(Stderr,"Computing facet complex...\n");
  PolyhedralFan ret(n);

  for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
    {
      PolyhedralFan a=facetsOfCone(*i);
      for(PolyhedralConeList::const_iterator j=a.cones.begin();j!=a.cones.end();j++)
	{
	  IntegerVector v=j->getRelativeInteriorPoint();
	  bool alreadyInRet=false;
	  for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
	    {
	      IntegerVector u=SymmetryGroup::compose(*k,v);
	      if(!j->containsRelatively(u))
		{
		  for(PolyhedralConeList::const_iterator l=ret.cones.begin();l!=ret.cones.end();l++)
		    if(l->containsRelatively(u))alreadyInRet=true;
		}
	    }
	  if(!alreadyInRet)ret.insert(*j);
	}
    }
  fprintf(Stderr,"Done computing facet complex.\n");
  return ret;
}

/*
IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup *sym)const
{
  assert(!cones.empty());
  int h=cones.begin()->dimensionOfLargestContainedSubspace();
  PolyhedralFan f=*this;
  while(f.getMaxDimension()!=h+1)
    {
      f=f.facetComplex();
    }

  IntegerVectorList rays;

  PolyhedralFan done(n);
  for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
    if(!done.contains(*i))
      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
	{
	  PolyhedralCone cone=i->permuted(*k);
	  if(!done.contains(cone))
	    {
	      rays.push_back(cone.getRelativeInteriorPoint());
	      done.insert(cone);
	    }
	}
  return rays;
}
*/

IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup *sym)const
{
  assert(!cones.empty());
  int h=cones.begin()->dimensionOfLargestContainedSubspace();
  PolyhedralFan f=*this;
  while(f.getMaxDimension()!=h+1)
    {
      f=f.facetComplexSymmetry(*sym);
    }

  IntegerVectorList rays;

  for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
    {
      PolyhedralFan done(n);
      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
	{
	  PolyhedralCone cone=i->permuted(*k);
	  if(!done.contains(cone))
	    {
	      rays.push_back(cone.getRelativeInteriorPoint());
	      done.insert(cone);
	    }
	}
    }
  return rays;
}


/*void PolyhedralFan::printWithIndices(class Printer *p, SymmetryGroup *sym)const //fan must be pure
{
  //  print(p);
  SymmetryGroup symm(n);
  if(!sym)sym=&symm;
  assert(!cones.empty());
  int h=cones.begin()->dimensionOfLargestContainedSubspace();
  fprintf(Stdout,"Rays:\n");
  //  IntegerVectorList rays;//=f.getRelativeInteriorPoints();
  fprintf(Stderr,"Computing rays...\n");
  IntegerVectorList rays=getRaysInPrintingOrder(sym);
  fprintf(Stderr,"Done computing rays.\n");
  p->printVectorList(rays,true);
  PolyhedralFan f=*this;

  //  while(f.getMaxDimension()>=h)
  IntegerVector fvector(f.getMaxDimension()-h);
  while(f.getMaxDimension()!=h)
    {
      int currentDimension=f.getMaxDimension()-h;
      fprintf(Stderr,"Processing dimension %i cones...\n",currentDimension);
      PolyhedralFan done(n);
      IntegerVector rayIncidenceCounter(rays.size());
      p->printString("Printing index list for dimension ");
      p->printInteger(currentDimension);
      p->printString(" cones:\n");
      p->printString("{\n");
      int faceIndex =0;
      bool split=false;
      bool addComma=false;
      for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
	{
	  if(!done.contains(*i))
	    {
	      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
		{
		  PolyhedralCone cone=i->permuted(*k);
		  if(!done.contains(cone))
		    {
		      // p->printString("Face ");
		      // p->printInteger(faceIndex);
		      // p->printString(": ");
		      int rayIndex=0;
		      IntegerVector indices(0);
		      for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
			{
			  if(cone.contains(*j))
			    {
			      indices.grow(indices.size()+1);
			      indices[indices.size()-1]=rayIndex;
			      rayIncidenceCounter[rayIndex]++;
			    }
			  rayIndex++;
			}
		      if(addComma)
			{
			  p->printString(",");
			  p->printNewLine();
			  if(split)
			    p->printNewLine();
			  split=false;
			}
		      p->printVector(indices,true);
		      addComma=true;
		      faceIndex++;
		      done.insert(cone);
		    }
		}
	      split=true;//p->printNewLine();
	    }
	}
      p->printString("}\n");
      p->printString("Number of dimension ");
      p->printInteger(f.getMaxDimension()-h);
      p->printString(" cones incident to each ray:\n");
      p->printVector(rayIncidenceCounter);
      p->printNewLine();
      fvector[f.getMaxDimension()-h-1]=faceIndex;
      f=f.facetComplex();
      fprintf(Stderr,"Done processing dimension %i cones.\n",currentDimension);
      //      fvector.grow(fvector.size()+1);fvector[fvector.size()-1]=faceIndex;
    }
  p->printString("F-vector:\n");
  p->printVector(fvector);
  p->printNewLine();
}
*/

void PolyhedralFan::printWithIndices(class Printer *p, SymmetryGroup *sym)const //fan must be pure
{
  //  print(p);
  SymmetryGroup symm(n);
  if(!sym)sym=&symm;
  assert(!cones.empty());
  int h=cones.begin()->dimensionOfLargestContainedSubspace();
  fprintf(Stdout,"Rays:\n");
  //  IntegerVectorList rays;//=f.getRelativeInteriorPoints();
  fprintf(Stderr,"Computing rays...\n");
  IntegerVectorList rays=getRaysInPrintingOrder(sym);
  fprintf(Stderr,"Done computing rays.\n");
  p->printVectorList(rays,true);
  PolyhedralFan f=*this;

  //  while(f.getMaxDimension()>=h)
  IntegerVector fvector(f.getMaxDimension()-h);
  while(f.getMaxDimension()!=h)
    {
      int currentDimension=f.getMaxDimension()-h;
      fprintf(Stderr,"Processing dimension %i cones...\n",currentDimension);
      IntegerVector rayIncidenceCounter(rays.size());
      p->printString("Printing index list for dimension ");
      p->printInteger(currentDimension);
      p->printString(" cones:\n");
      p->printString("{\n");
      int faceIndex =0;
      bool split=false;
      bool addComma=false;
      for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
	{
	  PolyhedralFan done(n);
	  for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
	    {
	      PolyhedralCone cone=i->permuted(*k);
	      if(!done.contains(cone))
		{
		  // p->printString("Face ");
		  // p->printInteger(faceIndex);
		  // p->printString(": ");
		  int rayIndex=0;
		  IntegerVector indices(0);
		  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
		    {
		      if(cone.contains(*j))
			{
			  indices.grow(indices.size()+1);
			  indices[indices.size()-1]=rayIndex;
			  rayIncidenceCounter[rayIndex]++;
			}
		      rayIndex++;
		    }
		  if(addComma)
		    {
		      p->printString(",");
		      p->printNewLine();
		      if(split)
			p->printNewLine();
		      split=false;
		    }
		  p->printVector(indices,true);
		  addComma=true;
		  faceIndex++;
		  done.insert(cone);
		}
	    }
	  split=true;//p->printNewLine();
	}
      p->printString("}\n");
      p->printString("Number of dimension ");
      p->printInteger(f.getMaxDimension()-h);
      p->printString(" cones incident to each ray:\n");
      p->printVector(rayIncidenceCounter);
      p->printNewLine();
      fvector[f.getMaxDimension()-h-1]=faceIndex;
      f=f.facetComplexSymmetry(*sym);
      fprintf(Stderr,"Done processing dimension %i cones.\n",currentDimension);
      //      fvector.grow(fvector.size()+1);fvector[fvector.size()-1]=faceIndex;
    }
  p->printString("F-vector:\n");
  p->printVector(fvector);
  p->printNewLine();
}

IncidenceList PolyhedralFan::getIncidenceList(SymmetryGroup *sym)const //fan must be pure
{
  IncidenceList ret;
  SymmetryGroup symm(n);
  if(!sym)sym=&symm;
  assert(!cones.empty());
  int h=cones.begin()->dimensionOfLargestContainedSubspace();
  IntegerVectorList rays=getRaysInPrintingOrder(sym);
  PolyhedralFan f=*this;

  while(f.getMaxDimension()!=h)
    {
      IntegerVectorList l;
      PolyhedralFan done(n);
      IntegerVector rayIncidenceCounter(rays.size());
      int faceIndex =0;
      for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
	{
	  if(!done.contains(*i))
	    {
	      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
		{
		  PolyhedralCone cone=i->permuted(*k);
		  if(!done.contains(cone))
		    {
		      int rayIndex=0;
		      IntegerVector indices(0);
		      for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
			{
			  if(cone.contains(*j))
			    {
			      indices.grow(indices.size()+1);
			      indices[indices.size()-1]=rayIndex;
			      rayIncidenceCounter[rayIndex]++;
			    }
			  rayIndex++;
			}
		      l.push_back(indices);
		      faceIndex++;
		      done.insert(cone);
		    }
		}
	    }
	}
      ret[f.getMaxDimension()]=l;
      f=f.facetComplex();
    }
  return ret;
}


void PolyhedralFan::makePure()
{
  if(getMaxDimension()!=getMinDimension())removeAllLowerDimensional();
}

bool PolyhedralFan::contains(PolyhedralCone const &c)const
{
  return cones.count(c);
}


PolyhedralFan::coneIterator PolyhedralFan::conesBegin()const
{
  return cones.begin();
}


PolyhedralFan::coneIterator PolyhedralFan::conesEnd()const
{
  return cones.end();
}


bool PolyhedralFan::isRefinementOf(PolyhedralFan const &f)const
{
  /*  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
    {
      static int t;
      fprintf(Stderr,"%i\n",t++);
      for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
	{
	  PolyhedralCone c=intersection(*i,*j);
	  if(c.dimension()==n)
	    {
	     
	      if(!j->contains(*i))
		return false;
	    }
	}
	}*/

  int t=0;
  for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
    {
      int t1=0;
      for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
	{
	  PolyhedralCone c=intersection(*i,*j);
	  if(c.dimension()==n)
	    {
	      t1++;
	      if(!j->contains(*i))
		return false;
	    }
	}
      fprintf(Stderr,"%i\n",t1);
      t+=t1;
    }
  fprintf(Stderr,"%i\n",t);
  
      /*      bool found=false;
      for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
	  if(j->contains(*i)){
	    found=true;
	    break;
	  }
      if(!found)
	{
	  for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
	    {
	      PolyhedralCone c=intersection(*i,*j);
	      if(c.dimension()==n){
		AsciiPrinter(Stderr).printVector(c.getRelativeInteriorPoint());
	      }
	    }
	  
	  return false;
	  }*/

  return true;
}