#include "wallideal.h"

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

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;
  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(markedExponent));

  return r;
}

static Polynomial wallPolynomial(Polynomial const &p, IntegerVectorList const &wallEqualities)
{
  Polynomial r;
  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(markedExponent));

  return r;
}

PolynomialSet wallIdeal(PolynomialSet const &groebnerBasis, IntegerVector const &wallNormal)
{
  PolynomialSet r;

  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)
{
  PolynomialSet r;

  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))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
{
  LexicographicTermOrder T;
  PolynomialSet LT;
  
  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);
      //      AsciiPrinter(Stderr).printPolynomialSet(g2);
      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())
    {
      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)
{
  //  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;
  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)),PolyhedralCone::positiveOrthant(n));
  
  return p.containsPositiveVector();
}


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);
    PolyhedralCone p=PolyhedralCone(normals,a);
    p.findFacets();
    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(best);

  return p;
}


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

  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);
  c.canonicalize();
  return c;
}

PolyhedralCone groebnerCone(PolynomialSet const &reducedGroebnerBasis)
{
  //  IntegerVectorList l=wallInequalities(reducedGroebnerBasis);
  IntegerVectorList l=algebraicTest(wallInequalities(reducedGroebnerBasis),reducedGroebnerBasis);
  IntegerVectorList a;
  PolyhedralCone c(l,a);
  c.canonicalize();
  return c;
}

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


PolynomialSet liftBasis(PolynomialSet const &toBeLifted, PolynomialSet const &originalBasisForFullIdeal)
{
  assert(toBeLifted.isValid());
  assert(originalBasisForFullIdeal.isValid());

  PolynomialSet newBasis;

  //  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;
}