#include "wallideal.h"

#include "division.h"
#include "buchberger.h"
#include "timer.h"
#include "printer.h"
#include "lp.h"
#include "polyhedralcone.h"
#include "log.h"

#include <set>

static Timer flipTimer("Flip",10);
static Timer flipTimer1("Flip1",10);
static Timer flipTimer2("Flip2",10);
static Timer flipTimer3("Flip3",10);
static Timer flipTimer4("Flip4",10);
static Timer flipTimer5("Flip5",10);
static Timer coneTimer("fajskfda",10);

Polynomial wallPolynomial(Polynomial const &p, IntegerVector const &wallNormal)
{
  Polynomial r(p.getRing());
  IntegerVector markedExponent=p.getMarked().m.exponent;

  for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
    {
      IntegerVector dif=markedExponent-i->first.exponent;
      if(dependent(dif,wallNormal))
        r+=Polynomial(Term(i->second,i->first));
    }

  r.mark(Monomial(r.getRing(),markedExponent));

  return r;
}

static Polynomial wallPolynomial(Polynomial const &p, IntegerVectorList const &wallEqualities)
{
  Polynomial r(p.getRing());
  IntegerVector markedExponent=p.getMarked().m.exponent;

  for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
    {
      IntegerVector dif=markedExponent-i->first.exponent;
      bool dep=false;
      for(IntegerVectorList::const_iterator j=wallEqualities.begin();j!=wallEqualities.end();j++)
	{
	  if(dependent(dif,*j))
	    {
	      dep=true;
	      break;
	    }
	}

      if(dep || dif.isZero())
	r+=Polynomial(Term(i->second,i->first));
    }

  r.mark(Monomial(p.getRing(),markedExponent));

  return r;
}

PolynomialSet wallIdeal(PolynomialSet const &groebnerBasis, IntegerVector const &wallNormal)
{
  PolynomialRing theRing=groebnerBasis.getRing();
  PolynomialSet r(theRing);

  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
    {
      r.push_back(wallPolynomial(*i,wallNormal));
    }
  return r;
}

PolynomialSet lowerDimensionalWallIdeal(PolynomialSet const &groebnerBasis, IntegerVectorList const &wallEqualities)
{
  PolynomialRing theRing=groebnerBasis.getRing();
  PolynomialSet r(theRing);

  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
    {
      r.push_back(wallPolynomial(*i,wallEqualities));
    }
  return r;
}

IntegerVectorList wallRemoveScaledInequalities(IntegerVectorList const &l)
{
  IntegerVectorList ret;

  for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
    {
      bool found=false;
      
      assert(!i->isZero());

      for(IntegerVectorList::const_iterator k=ret.begin();k!=ret.end();k++)
	if(dependent(*i,*k)&&dotLong(*i,*k)>0)found=true;
      
      if(!found)ret.push_back(*i);
    }
  return ret;
}


IntegerVectorList algebraicTest(IntegerVectorList const &l, PolynomialSet const &groebnerBasis) // TO DO: FIGURE OUT IF THIS TEST WORKS IN THE NON-HOMOGENEOUS CASE
{
  PolynomialRing theRing=groebnerBasis.getRing();
  LexicographicTermOrder T;
  PolynomialSet LT(theRing);
  
  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
    {
      LT.push_back(Polynomial(i->getMarked()));
    }
  

  //fprintf(Stderr,"In Size:%i\n",l.size());
  IntegerVectorList ret;
  for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
    {
      bool accept=true;
      PolynomialSet g2=wallIdeal(groebnerBasis,*i);
      for(PolynomialSet::const_iterator j=g2.begin();j!=g2.end() && accept;j++)
	for(PolynomialSet::const_iterator k=j;k!=g2.end();k++)
	  {
	    //	    fprintf(Stderr,"test\n");
	    if((!j->isMonomial()) || (!k->isMonomial()))
	    if(!relativelyPrime(j->getMarked().m.exponent,k->getMarked().m.exponent))
	      {
		Polynomial s=sPolynomial(*j, *k);
		if(!s.isZero())
		  {
		    bool witness=true;
		    for(TermMap::const_iterator i=s.terms.begin();i!=s.terms.end();i++)
		      {
			bool inideal=false;
			for(PolynomialSet::const_iterator j=LT.begin();j!=LT.end();j++)
			  if(j->getMarked().m.exponent.divides(i->first.exponent))
			    {
			      inideal=true;
			      break;
			    }
			if(inideal)
			  {
			    witness=false;
			    break;
			  }
		      }
		    if(witness)
		      {
			accept=false;
			break;
		      }

/*
		    s.mark(T); // with respect to some termorder
		    s.scaleMarkedCoefficientToOne();

		    if(1)
		      {
			Polynomial t=division(s,LT,T);
			if((t-s).isZero())
			  {
			    accept=false;
			    break;
			  }
		      }
		    	    else
		      {
			s=division(s,g2,T);
			if(!s.isZero())
			  {
			    accept=false;
			    break;
			  }
			  }*/
		  }
	      }
	    
	  }
      if(accept)ret.push_back(*i);
    }
  //  fprintf(Stderr,"Out Size:%i\n",ret.size());
  return ret;
}


IntegerVectorList wallInequalities(PolynomialSet const &groebnerBasis)
{
  IntegerVectorList ret;

  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
    {
      IntegerVector markedExponent=i->getMarked().m.exponent;

      for(TermMap::const_iterator j=i->terms.begin();j!=i->terms.end();j++)
        {
          IntegerVector dif=markedExponent-j->first.exponent;

          if(!dif.isZero())
            {
	      /* bool found=false; //These four lines are now done in wallRemoveScaledInequalities()

              for(IntegerVectorList::const_iterator k=ret.begin();k!=ret.end();k++)
                if(dependent(dif,*k))found=true;
             
              if(!found)
	      */
	      ret.push_back(dif);
            }
        }
    }

  if(ret.empty())
    {
      log3 fprintf(Stderr,"WARNING: wallideal.cpp - No defining equalities of the Groebner cone! Is your ideal equal to the ring?\n");
    }

  return wallRemoveScaledInequalities(ret);
  //  return algebraicTest(wallRemoveScaledInequalities(ret),groebnerBasis);
}


PolynomialSet flip(PolynomialSet const &groebnerBasis, IntegerVector const &wallNormal)
{
  PolynomialRing theRing=groebnerBasis.getRing();
  //  fprintf(Stderr,"ENTERING flip\n");
  //  fprintf(Stderr,"flip - start\n");
  //  AsciiPrinter(Stderr).printPolynomialSet(groebnerBasis);
  //  AsciiPrinter(Stderr).printVector(wallNormal);

  TimerScope ts(&flipTimer);
  //Subroutine 3.7 in [Sturmfels]

  // Step 1
  //  fprintf(Stderr,"flip - step1\n");  
  flipTimer1.on();
  PolynomialSet wall=wallIdeal(groebnerBasis,wallNormal);
  wall.markAndScale(WeightTermOrder(wallNormal));// This marking will be used later on when we lift
  //  fprintf(Stderr,"Changed order:\n");
  // AsciiPrinter(Stderr).printPolynomialSet(wall);
  flipTimer1.off();
  
  // Step 2
  //  fprintf(Stderr,"flip - step2\n"); 
  // AsciiPrinter(Stderr).printPolynomialSet(wall);
  flipTimer2.on();
  PolynomialSet oldWall=wall;
  WeightTermOrder flipOrder(-wallNormal);
  buchberger(&wall,flipOrder);
  flipTimer2.off();

  // Step 3
  //  fprintf(Stderr,"flip - step3\n");  
  flipTimer3.on();
  PolynomialSet newBasis(theRing);
  flipTimer3.off();

  // Step 4 Lift
  //  fprintf(Stderr,"flip - lifting\n");
  //  fprintf(Stderr,"flip - step4\n");  

  {
    // liftBasis() could be used for this!!!!

    TimerScope ts(&flipTimer4);
    for(PolynomialSet::const_iterator j=wall.begin();j!=wall.end();j++)
      {
	newBasis.push_back(divisionLift(*j, oldWall, groebnerBasis, LexicographicTermOrder()));
	/*{
	    // The following should also work:
	    Polynomial q=*j-division(*j,groebnerBasis,LexicographicTermOrder());
	    assert(!q.isZero());
	    
	    newBasis.push_back(q);
	    }*/
      }
  }
  
  // Step 5 Autoreduce
  //  fprintf(Stderr,"flip - autoreduce\n");
  //  AsciiPrinter(Stderr).printPolynomialSet(wall);
  //  AsciiPrinter(Stderr).printPolynomialSet(newBasis);
  //  fprintf(Stderr,"flip - step5\n");  
  {
    TimerScope ts(&flipTimer5);
    PolynomialSet::const_iterator k=wall.begin();
    for(PolynomialSet::iterator j=newBasis.begin();j!=newBasis.end();j++)
      {
	j->mark(k->getMarked().m);
	k++;
      }
    
    //  fprintf(Stderr,"Marked order:\n");
    //  AsciiPrinter(Stderr).printPolynomialSet(wall);
    //  AsciiPrinter(Stderr).printPolynomialSet(newBasis);
    //  fprintf(Stderr,"Not reduced lifted basis:\n");
    //  AsciiPrinter(Stderr).printPolynomialSet(newBasis);

    static int t;
    t++;
    t&=0;
    if(t==0)
      {
	// fprintf(Stderr,"autoreducing ..");
	autoReduce(&newBasis,LexicographicTermOrder());
	// fprintf(Stderr,".. done\n");
      }
    
    //autoReduce(&newBasis,StandardGradedLexicographicTermOrder());
  }
  
  //  fprintf(Stderr,"flip - done\n");
  //  fprintf(Stderr,"LEAVING flip\n");

  //fprintf(Stderr,"%i",newBasis.size());

  return newBasis;
}


bool wallContainsPositiveVector(IntegerVector const &wallNormal)
{
  //This is not right I think
  int n=wallNormal.size();
  for(int i=0;i<n;i++)if(wallNormal[i]<0)return true;

  return false;
}

void wallAddCoordinateWalls(IntegerVectorList &normals)
{
  assert(!normals.empty());
  int n=normals.begin()->size();
  for(int i=0;i<n;i++)normals.push_back(IntegerVector::standardVector(n,i));
}


bool isIdealHomogeneous(PolynomialSet const &groebnerBasis)
{
  int n=groebnerBasis.numberOfVariablesInRing();
  IntegerVectorList a;
  PolyhedralCone p=intersection(PolyhedralCone(a,wallInequalities(groebnerBasis),n),PolyhedralCone::positiveOrthant(n));
  
  return p.containsPositiveVector();
}


/* This routine is a preprocessing step for redudancy removal.
   The routine normalizes the list of vectors in gcd sense.
   It removes duplicates.
   It removes direction that are sums of other directions.
   The input must satisfy:
     Input must be pointed, meaning that there must exist a codimension one subspace with all the input vectors strictly on one side.
     It particular, there is no zero-vector in the list (It would be easy to change the routine to handle this case)
   These requirements guarantee that *i and *j are not removed in the line with the comment.
 */
IntegerVectorList normalizedWithSumsAndDuplicatesRemoved(IntegerVectorList const &a)
{
  IntegerVectorList ret;
  set<IntegerVector> b;

  for(IntegerVectorList::const_iterator i=a.begin();i!=a.end();i++)
    {
      assert(!(i->isZero()));
      b.insert(normalized(*i));
    }

  for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
    for(set<IntegerVector>::const_iterator j=i;j!=b.end();j++)
	if(i!=j)b.erase(normalized(*i+*j));//this can never remove *i or *j

  for(set<IntegerVector>::const_iterator i=b.begin();i!=b.end();i++)
    ret.push_back(*i);

  return ret;
}

IntegerVectorList wallFlipableNormals(PolynomialSet const &groebnerBasis, bool isKnownToBeHomogeneous)
{
  // New optimised version using PolyhedralCone
  int n=groebnerBasis.numberOfVariablesInRing();
  IntegerVectorList a;
  //AsciiPrinter(Stderr).printVectorList(wallInequalities(groebnerBasis));
  //PolyhedralCone p=intersection(PolyhedralCone(wallInequalities(groebnerBasis),a),PolyhedralCone::positiveOrthant(n));

  IntegerVectorList normals=algebraicTest(wallInequalities(groebnerBasis),groebnerBasis);
    coneTimer.on();
    //    PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);

    /*    AsciiPrinter(Stderr).printVectorList(normals);
    AsciiPrinter(Stderr).printInteger(normals.size());
    */
    normals=normalizedWithSumsAndDuplicatesRemoved(normals);

    //    AsciiPrinter(Stderr).printVectorList(normals);
    //   AsciiPrinter(Stderr).printInteger(normals.size());
    
    PolyhedralCone p=PolyhedralCone(normals,a,n);
    //  fprintf(Stderr,"Finding facets\n");
    p.findFacets();
    //  fprintf(Stderr,"Done Finding facets\n");
    coneTimer.off();

  IntegerVectorList ilist=p.getHalfSpaces();
  IntegerVectorList ret;
  for(IntegerVectorList::iterator i=ilist.begin();i!=ilist.end();i++)
    {
      bool facetOK=true;
      if(!isKnownToBeHomogeneous)
	{
	  *i=-1 * (*i);
	  PolyhedralCone temp=intersection(PolyhedralCone(ilist,a),PolyhedralCone::positiveOrthant(n));
	  facetOK=(temp.dimension()==n);
	  *i=-1 * (*i);
	}

      if(facetOK)
	ret.push_back(*i);
    }

  // New version using PolyhedralCone
  /*  int n=groebnerBasis.numberOfVariablesInRing();
  IntegerVectorList a;
  //  AsciiPrinter(Stderr).printVectorList(wallInequalities(groebnerBasis));
  PolyhedralCone p=intersection(PolyhedralCone(wallInequalities(groebnerBasis),a),PolyhedralCone::positiveOrthant(n));
  //  PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
  p.findFacets();
  IntegerVectorList ilist=p.getHalfSpaces();
  IntegerVectorList ret;
  for(IntegerVectorList::const_iterator i=ilist.begin();i!=ilist.end();i++)
    if(wallContainsPositiveVector(*i))ret.push_back(*i);
  */

  /* Old code not using PolyhedralCone
  IntegerVectorList ilist=wallInequalities(groebnerBasis);
  wallAddCoordinateWalls(ilist);
  IntegerVectorList ret;
  for(IntegerVectorList::const_iterator i=ilist.begin();i!=ilist.end();i++)
    if(isFacet(ilist,i)&&wallContainsPositiveVector(*i))
      ret.push_back(*i);
  */
  return ret;
}


Polynomial flipMinkowski(Polynomial p, IntegerVector const &wallNormal)
{
  IntegerVector best=p.getMarked().m.exponent;

  for(TermMap::const_iterator i=p.terms.begin();i!=p.terms.end();i++)
    {
      IntegerVector e=i->first.exponent;
      IntegerVector diff=e-best;
      if(dependent(diff,wallNormal)&&dot(diff,wallNormal)<0)best=e;
    }
  p.mark(Monomial(p.getRing(),best));

  return p;
}


PolynomialSet flipMinkowski(PolynomialSet const &groebnerBasis, IntegerVector const &wallNormal)
{
  PolynomialSet r(groebnerBasis.getRing());

  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
    r.push_back(flipMinkowski(*i,wallNormal));

  return r;
}


PolyhedralCone homogeneitySpace(PolynomialSet const &reducedGroebnerBasis)
{
  IntegerVectorList l=wallInequalities(reducedGroebnerBasis);
  IntegerVectorList a;
  PolyhedralCone c(a,l,reducedGroebnerBasis.getRing().getNumberOfVariables());
  c.canonicalize();
  return c;
}

PolyhedralCone groebnerCone(PolynomialSet const &reducedGroebnerBasis, bool useAlgebraicTest)
{
  int n=reducedGroebnerBasis.getRing().getNumberOfVariables();
  IntegerVectorList l=wallInequalities(reducedGroebnerBasis);
  if(useAlgebraicTest)l=algebraicTest(l,reducedGroebnerBasis);
  l=normalizedWithSumsAndDuplicatesRemoved(l);
  IntegerVectorList a;
  PolyhedralCone c(l,a,n);
  c.canonicalize();
  return c;
}

int dimensionOfHomogeneitySpace(PolynomialSet const &reducedGroebnerBasis)
{
  return homogeneitySpace(reducedGroebnerBasis).dimension();
}


PolynomialSet liftBasis(PolynomialSet const &toBeLifted, PolynomialSet const &originalBasisForFullIdeal)
{
  PolynomialRing theRing=toBeLifted.getRing();
  assert(toBeLifted.isValid());
  assert(originalBasisForFullIdeal.isValid());

  PolynomialSet newBasis(theRing);

  //  fprintf(Stderr,"LIFTING:");
  //  AsciiPrinter(Stderr).printPolynomialSet(toBeLifted);

  for(PolynomialSet::const_iterator j=toBeLifted.begin();j!=toBeLifted.end();j++)
    {
      assert(!j->isZero());
      //AsciiPrinter(Stderr).printVector(j->getMarked().m.exponent);
  
      Polynomial q=*j-division(*j,originalBasisForFullIdeal,LexicographicTermOrder());
      assert(!q.isZero());
      q.mark(j->getMarked().m);      	    
      newBasis.push_back(q);
    }
  autoReduce(&newBasis,LexicographicTermOrder());
  //  fprintf(Stderr,"TO:");
  //  AsciiPrinter(Stderr).printPolynomialSet(newBasis);

  return newBasis;
}


bool isMarkingConsistent(PolynomialSet const &g)
{
  IntegerVectorList empty;
  PolyhedralCone c(wallInequalities(g),empty,g.getRing().getNumberOfVariables());
  c=intersection(c,PolyhedralCone::positiveOrthant(c.ambientDimension()));
  log1 AsciiPrinter(Stderr).printPolyhedralCone(c);
  return c.dimension()==c.ambientDimension();
}