#include "tropicalbasis.h"

#include <iostream>

#include "buchberger.h"
#include "groebnerengine.h"
#include "tropical.h"
#include "tropical2.h"
#include "division.h"
#include "wallideal.h"
#include "halfopencone.h"
#include "log.h"

#include "timer.h"

static Timer iterativeTropicalBasisTimer("Iterative tropical basis",1);

typedef set<IntegerVector> IntegerVectorSet;

static Polynomial cleverSaturation(Polynomial const &p, IntegerVector const &w)
{
  PolynomialRing theRing=p.getRing();
  if(p.isZero())return p;
  if(w.size()==0)return p;
  Polynomial f=initialForm(p,w);
  debug<<"P:"<<p<<" w: "<<w<<" f: "<<f<<"\n";
  IntegerVector gcd=f.greatestCommonMonomialDivisor();
  if(!gcd.isZero())
    {
      debug<<"OLD"<<p<<"\n";
      debug<<"initialForm"<<initialForm(p,w)<<"\n";
      PolynomialSet reducer(p.getRing());
      reducer.push_back(theRing.monomialFromExponentVector(IntegerVector::allOnes(theRing.getNumberOfVariables()))-theRing.one());
      reducer.markAndScale(LexicographicTermOrder());
      debug<<gcd;
      int s=gcd.max();
      Polynomial all=theRing.monomialFromExponentVector(s*IntegerVector::allOnes(theRing.getNumberOfVariables())-gcd);
      Polynomial g= division(all*p,reducer,LexicographicTermOrder());
      g.saturate();
      debug<<"NEW"<<g<<"\n";
      return g;
    }
  return p;
}

static void initialSaturatingBuchberger(PolynomialSet *g, TermOrder const &termOrder_, IntegerVector const &w_)
{
  int n=w_.size();
  IntegerVectorList temp;
  temp.push_back(IntegerVector::standardVector(n+1,n));
  IntegerVector w=concatenation(w_,IntegerVector(1));w[w.size()-1]=-w_.sum();
  temp.push_back(w);
  MatrixTermOrder termOrder(temp);
  PolynomialRing theRing=g->getRing().withVariablesAppended("Z");
  PolynomialSet sPolynomials(theRing);

  for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
    if(!i->isZero())sPolynomials.push_back(cleverSaturation(i->embeddedInto(theRing),w)); // It is safe and useful to ignore the 0 polynomial

  sPolynomials.push_back(theRing.monomialFromExponentVector(IntegerVector::allOnes(theRing.getNumberOfVariables()))-theRing.one());

  sPolynomials.saturate();
  sPolynomials.markAndScale(termOrder);

  *g=PolynomialSet(theRing);

  while(!sPolynomials.empty())
    {
      Polynomial p=*sPolynomials.begin();
      sPolynomials.pop_front();

      p=division(p,*g,termOrder);
      if(!p.isZero())
        {
          p=cleverSaturation(p,w);
          p.mark(termOrder);
          p.scaleMarkedCoefficientToOne();
          bool isMonomial=p.isMonomial();
          for(PolynomialSet::const_iterator i=g->begin();i!=g->end();i++)
            if((!isMonomial) || (!i->isMonomial())) // 2 % speed up!
            {
              if(!relativelyPrime(i->getMarked().m.exponent,p.getMarked().m.exponent))
                {
                  Polynomial s=sPolynomial(*i,p);
                  s.mark(termOrder); // with respect to some termorder
                  s.scaleMarkedCoefficientToOne();
                  sPolynomials.push_back(s);
                }
            }
          g->push_back(p);
          {
            static int t;
            t++;
            //      if((t&31)==0)fprintf(Stderr," gsize %i  spolys:%i\n",g->size(),sPolynomials.size());
          }
        }
    }
  minimize(g);
  autoReduce(g,termOrder);
}



PolynomialSet tropicalBasisOfCurve(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension) //Assuming g is homogeneous
{

  /*
   * TODO: Introduce the notion of a tropical prebasis:
   *
   * Definition. A set of polynomials f_1,...,f_m is called a tropical prebasis for the ideal they
   * generate if for every w not in the tropical variety of that ideal there exists a monomial in
   * the ideal generated by the initial forms of f_1,...,f_m w.r.t. w.
   *
   * Computing a tropical prebasis could be faster than computing a tropical basis since fewer
   * groebner bases for the originial ideal might be needed. Still, however, it is relatively easy
   * to determine the tropical variety given a tropical prebasis.
   */
//  bool prebasis=true;
//  debug<<"PREBASIS="<<prebasis<<"\n";
  log2 debug<<"TropicalBasis begin\n";
	log2 debug<<g;
	int homog=linealitySpaceDimension;
  assert(homog>0 || n==0);
  TimerScope ts(&iterativeTropicalBasisTimer);
  PolyhedralFan f(n);
  if(!intersectionFan)intersectionFan=&f;

  //  *intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);
//	log1 fprintf(Stderr,"WARINING USING EXPERIMENTAL TROPICAL HYPERSURFACE INTERSECTION ROUTINE!!\n");
  *intersectionFan=tropicalHyperSurfaceIntersectionClosed(n, g);

  IntegerVectorSet containsNoMonomialCache;

  while(1)
    {
      PolyhedralFan::coneIterator i;

//      {AsciiPrinter P(Stderr);intersectionFan->printWithIndices(&P);}
restart:
//      {AsciiPrinter P(Stderr);intersectionFan->printWithIndices(&P);}
      for(i=intersectionFan->conesBegin();i!=intersectionFan->conesEnd();i++)
	{
//	  log1 cerr<<"!@#$";
	  IntegerVector relativeInteriorPoint=i->getRelativeInteriorPoint();
//	  log1 cerr<<"1234/n";

	  if(containsNoMonomialCache.count(relativeInteriorPoint)>0)
	    {
	      log2 fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
	    }
	  else
	    {
/*	      if(prebasis)
	        {
	          if(containsMonomial(initialForms(g,relativeInteriorPoint)))
	            {
	              intersectionFan->insertFacetsOfCone(*i);
	              intersectionFan->remove(*i);
	              debug<<"LOWERING DIMENSION OF CONE\n";//TODO: checking cones in order of dimension could avoid this.
	              goto restart;
	            }
	        }*/
	      WeightReverseLexicographicTermOrder t(relativeInteriorPoint);
	      log2 fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
	      log2 AsciiPrinter(Stderr).printVector(relativeInteriorPoint);
	      log2 fprintf(Stderr,"\n");
	      PolynomialSet h2=g;
  //            debug<<"g"<<g;

/*	      {Adjust these lines somehow to enable the saturating buchberger.
	        debug<<h2;
	        initialSaturatingBuchberger(&h2, t, relativeInteriorPoint);
	        //buchberger(&h2,t,true);
	        debug<<"SATURATING BUCHBERGER DONE\n";
	      }*/
	     // buchberger(&h2,t);
	      h2=GE_groebnerBasis(h2,t,true);
        //      debug<<"h2"<<h2;
	      log2 fprintf(Stderr,"Done computing Gr\"obner basis.\n");

	  //    debug<<h2;
//	      log3 AsciiPrinter(Stderr).printPolynomialSet(h2);

	      PolynomialSet wall=initialFormsAssumeMarked(h2,relativeInteriorPoint);

	      if(containsMonomial(wall))
		{
		  log2 fprintf(Stderr,"Initial ideal contains a monomial.\n");
		  Polynomial m(computeTermInIdeal(wall));
		  log2 fprintf(Stderr,"Done computing term in ideal\n");

		  Polynomial temp=m-division(m,h2,LexicographicTermOrder());
		  g.push_back(temp);

		  log2 fprintf(Stderr,"Adding element to basis:\n");
		  log2 AsciiPrinter(Stderr).printPolynomial(temp);
		  log2 fprintf(Stderr,"\n");

		  *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
		  break;
		}
	      else
		{
		  if(i->dimension()<=1+homog)
		    //if(!containsMonomial(wall) && i->dimension()<=1+homog)//line for testing perturbation code
		    {
		      log2 fprintf(Stderr,"Initial ideal contains no monomial... caching weight vector.\n");
		      containsNoMonomialCache.insert(relativeInteriorPoint);
		    }
		  else
		    {
		      /* We need to compute the initial ideal with
			 respect to "relativeInteriorPoint" perturbed
			 with a basis of the span of the cone. Instead
			 of perturbing we may as well compute initial
			 ideal successively. We have already computed
			 the initial ideal with respect to
			 "relativeInteriorPoint". To get the perturbed
			 initial ideal we take initial ideal with
			 respect to each vector in the basis of the
			 span.*/
		      IntegerVectorList empty;
		      PolyhedralCone dual=PolyhedralCone(empty,i->getEquations(),i->ambientDimension()).dualCone();
		      dual.canonicalize();
		      IntegerVectorList basis=dual.getEquations();
		      PolynomialSet witnessLiftBasis=h2;//basis with respect to relativeInteriorPoint
		      log2 debug<<"basis"<<basis<<"\n";
		      for(IntegerVectorList::const_iterator j=basis.begin();j!=basis.end();j++)
			{
			  log2 debug<<"wall"<<wall<<"\n";
			  WeightReverseLexicographicTermOrder t(*j);
			  PolynomialSet h3=wall;
			  buchberger(&h3,t);
			  wall=initialFormsAssumeMarked(h3,*j);
			  witnessLiftBasis=liftBasis(h3,witnessLiftBasis);
			}
                      log2 debug<<"wall"<<wall<<"\n";
		      if(containsMonomial(wall))
			{
			  Polynomial m(computeTermInIdeal(wall));
			  Polynomial temp=m-division(m,witnessLiftBasis,LexicographicTermOrder());
			  g.push_back(temp);
			  *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
			  break;
			}
		      else
			{
			  assert(0);
			}
		    }
		}
	    }
	}
      if(i==intersectionFan->conesEnd())break;
    }

  log2 debug<<"TropicalBasis end\n";
  log2 cerr <<"RETURNING";
  return g;
}

/*
PolynomialSet iterativeTropicalBasisNoPerturbation(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension, bool doPrint) //Assuming g is homogeneous
{
  TimerScope ts(&iterativeTropicalBasisTimer);
  PolyhedralFan f(n);
  if(!intersectionFan)intersectionFan=&f;

  *intersectionFan=tropicalPrincipalIntersection(n,g,linealitySpaceDimension);

  IntegerVectorSet containsNoMonomialCache;

  while(1)
    {
      //      AsciiPrinter(Stderr).printPolyhedralFan(*intersectionFan);
      //      assert(f.getMaxDimension()==1);

      IntegerVectorList l=intersectionFan->getRelativeInteriorPoints();

      IntegerVectorList::const_iterator i;
      for(i=l.begin();i!=l.end();i++)
	{
	  if(containsNoMonomialCache.count(*i)>0)
	    {
	      if(doPrint)fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
	    }
	  else
	    {
	      WeightReverseLexicographicTermOrder t(*i);
	      if(doPrint)fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
	      if(doPrint)AsciiPrinter(Stderr).printVector(*i);
	      if(doPrint)fprintf(Stderr,"\n");
	      PolynomialSet h2=g;
	      buchberger(&h2,t);
	      if(doPrint)fprintf(Stderr,"Done computing Gr\"obner basis.\n");

	      //	      AsciiPrinter(Stderr).printPolynomialSet(h2);
	      PolynomialSet wall=initialFormsAssumeMarked(h2,*i);
	      //fprintf(Stderr,"Wall ideal:\n");
	      //AsciiPrinter(Stderr).printPolynomialSet(wall);

	      if(containsMonomial(wall))
		{
		  if(doPrint)fprintf(Stderr,"Initial ideal contains a monomial.\n");
		  Polynomial m(computeTermInIdeal(wall));
		  if(doPrint)fprintf(Stderr,"Done computing term in ideal\n");

		  Polynomial temp=m-division(m,h2,LexicographicTermOrder());
		  g.push_back(temp);

		  if(doPrint)fprintf(Stderr,"Adding element to basis:\n");
		  if(doPrint)AsciiPrinter(Stderr).printPolynomial(temp);
		  if(doPrint)fprintf(Stderr,"\n");

		  *intersectionFan=refinement(*intersectionFan,PolyhedralFan::bergmanOfPrincipalIdeal(temp),linealitySpaceDimension,true);
		  break;
		}
	      else
		{
		  if(doPrint)fprintf(Stderr,"Initial ideal contains no monomial... caching weight vector.\n");
		  containsNoMonomialCache.insert(*i);
		}
	    }
	}
      if(i==l.end())break;
    }
  return g;
}

*/