#include "tropical2.h"

#include <stdlib.h>
#include "buchberger.h"
#include "division.h"
#include "tropical.h"
#include "wallideal.h"
#include "dimension.h"

#include "timer.h"

static Timer iterativeTropicalBasisTimer("Iterative tropical basis",1);
static Timer tropicalPrincipalIntersectionTimer("Tropical principal intersection",1);

Polynomial initialForm(Polynomial const &p, IntegerVector const &weight)
{
  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(dot(dif,weight)==0)
	r+=Polynomial(Term(i->second,i->first));
    }
  r.mark(Monomial(markedExponent));

  return r;
}


PolynomialSet initialForms(PolynomialSet const &groebnerBasis, IntegerVector const &weight)
{
  PolynomialSet r;

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



PolyhedralFan tropicalPrincipalIntersection(int n, PolynomialSet const &g, int linealitySpaceDimension)
{
  TimerScope ts(&tropicalPrincipalIntersectionTimer);
  PolyhedralFan ret=PolyhedralFan::fullSpace(n);

  for(PolynomialSet::const_iterator i=g.begin();i!=g.end();i++)
    {
      ret=refinement(ret,PolyhedralFan::bergmanOfPrincipalIdeal(*i),linealitySpaceDimension,true);
    }
  return ret;
}


typedef set<IntegerVector> IntegerVectorSet;

PolynomialSet iterativeTropicalBasis(int n, PolynomialSet g, PolyhedralFan *intersectionFan, int linealitySpaceDimension, bool doPrint) //Assuming g is homogeneous
{
  int homog=linealitySpaceDimension;
  assert(homog>0);
  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;
      PolyhedralFan::coneIterator i;
      for(i=intersectionFan->conesBegin();i!=intersectionFan->conesEnd();i++)
	//     for(i=l.begin();i!=l.end();i++)
	{
	  IntegerVector relativeInteriorPoint=i->getRelativeInteriorPoint();
	  if(containsNoMonomialCache.count(relativeInteriorPoint)>0)
	    {
	      if(doPrint)fprintf(Stderr,"Weight vector found in cache.... contains no monomial.\n");
	    }
	  else
	    {
	      //  fprintf(Stderr,"\n");
	      //  AsciiPrinter(Stderr).printVector(relativeInteriorPoint);
	      WeightReverseLexicographicTermOrder t(relativeInteriorPoint);
	      if(doPrint)fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
	      if(doPrint)AsciiPrinter(Stderr).printVector(relativeInteriorPoint);
	      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=initialForms(h2,relativeInteriorPoint);
	      //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(i->dimension()<=1+homog)
		    //if(!containsMonomial(wall) && i->dimension()<=1+homog)//line for testing perturbation code
		    {
		      if(doPrint)fprintf(Stderr,"Initial ideal contains no monomial... caching weight vector.\n");
		      containsNoMonomialCache.insert(relativeInteriorPoint);
		    }
		  else
		    {
		      /*		      fprintf(Stderr,"Perturbation needed in computation of tropical basis of curve!!!!!!!!!!!!!!!!\n");
		      fprintf(Stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
		      fprintf(Stderr,"THE TROPICAL BASIS PERTURBATION CODE HAS NOT BEEN TESTED.\n");
		      fprintf(Stderr,"PLEASE SEND YOUR EXAMPLE TO THE AUTHOR.\n");
		      fprintf(Stderr,"RUNNING TROPICAL BASIS PERTURBATION CODE.\n");
		      fprintf(Stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
		      fprintf(Stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
		      fprintf(Stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
		      fprintf(Stderr,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
		      */
		      /* 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.*/
		      PolyhedralCone dual=PolyhedralCone(IntegerVectorList(),i->getLinealitySpace(),i->ambientDimension()).dualCone();
		      dual.canonicalize();
		      IntegerVectorList basis=dual.getLinealitySpace();
		      PolynomialSet witnessLiftBasis=h2;//basis with respect to relativeInteriorPoint
		      //  AsciiPrinter(Stderr).printPolynomialSet(wall);
		      //  fprintf(Stderr,"Lifted");
		      //  AsciiPrinter(Stderr).printPolynomialSet(witnessLiftBasis);
		      for(IntegerVectorList::const_iterator j=basis.begin();j!=basis.end();j++)
			{
			  //  AsciiPrinter(Stderr).printVector(*j);
			  WeightReverseLexicographicTermOrder t(*j);
			  PolynomialSet h3=wall;
			  buchberger(&h3,t);
			  wall=initialForms(h3,*j);
			  //  AsciiPrinter(Stderr).printPolynomialSet(wall);
			  witnessLiftBasis=liftBasis(h3,witnessLiftBasis);
			  //  fprintf(Stderr,"Lifted");
			  //  AsciiPrinter(Stderr).printPolynomialSet(witnessLiftBasis);
			}
		      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;
    }
  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=initialForms(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;
}


static PolynomialSet checkList(IntegerVectorList const &l, PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, int h, bool &result, bool onlyCheckRays)
{
  for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
    {
      WeightReverseLexicographicTermOrder t(*i);
      fprintf(Stderr,"Computing Gr\"obner basis with respect to:");
      AsciiPrinter(Stderr).printVector(*i);
      fprintf(Stderr,"\n");
      PolynomialSet h2=groebnerBasis;
      buchberger(&h2,t);
      fprintf(Stderr,"Done computing Gr\"obner basis.\n");
      
      //AsciiPrinter(Stderr).printPolynomialSet(h2);
      PolynomialSet wall=initialForms(h2,*i);
      
      int hdim2=dimensionOfHomogeneitySpace(wall);
      if(hdim2>h)
	{
	  if(!containsMonomial(wall))
	    {
	      fprintf(Stderr,"Iterating recursively.\n");
	      //PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,0);
	      PolynomialSet initialIdeal=guessInitialIdealWithoutMonomial(wall,fullNeighbourBasis,onlyCheckRays);
	      
	      if(fullNeighbourBasis)
		{
		  //*fullNeighbourBasis=liftBasis(initialIdeal,h2);
		  *fullNeighbourBasis=liftBasis(*fullNeighbourBasis,h2);
		}
		  
	      
	      result=true;
	      return initialIdeal;
	    }
	}
    }
  result=false;
  return groebnerBasis;
}

PolynomialSet guessInitialIdealWithoutMonomial(PolynomialSet const &groebnerBasis, PolynomialSet *fullNeighbourBasis, bool onlyCheckRays) //ideal must be homogeneous
  // fullNeighbourBasis is set to a Groebner basis of the full ideal. The returned basis and fullNeighbourBasis have at least one termorder in common
{
  assert(groebnerBasis.isValid());
  if(fullNeighbourBasis)
    {
      assert(fullNeighbourBasis->isValid());
    }

  int n=groebnerBasis.numberOfVariablesInRing();
  int h=dimensionOfHomogeneitySpace(groebnerBasis);
  int d=krullDimension(groebnerBasis);

  if(d==h)
    {
      if(fullNeighbourBasis)*fullNeighbourBasis=groebnerBasis;
      return groebnerBasis;
    }

  {
    IntegerVectorList a;
    PolyhedralCone p=PolyhedralCone(wallInequalities(groebnerBasis),a);
    IntegerVectorList extreme=p.extremeRays();
    fprintf(Stderr,"Extreme rays of Groebner cone:\n");
    AsciiPrinter(Stderr).printVectorList(extreme);

    bool result;
    PolynomialSet r=checkList(extreme,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
    if(result)return r;
  }
  if(onlyCheckRays)
    {
      fprintf(Stderr,"UNABLE TO COMPUTE STARTING CONE.\n");
      assert(0); // We could not guess an intial ideal
    }

  PolyhedralFan f=PolyhedralFan::fullSpace(n);
  /*  for(int i=0;i<d-1;i++)
    {
      IntegerVector v(n);
      for(int j=0;j<n;j++)v[j]=rand()&1;
      IntegerVectorList a,b;
      b.push_back(v);
      PolyhedralFan F(n);
      F.insert(PolyhedralCone(a,b));
      f=refinement(f,F);
    }
  AsciiPrinter P(Stderr);
  f.print(&P);
  */
  int hypersurfacesToGo=groebnerBasis.size();
  for(PolynomialSet::const_iterator i=groebnerBasis.begin();i!=groebnerBasis.end();i++)
    {
      fprintf(Stderr,"Hypersurfaces to go:%i\n",hypersurfacesToGo--);
      fprintf(Stderr,"Max dimension: %i\n",f.getMaxDimension());
      f=refinement(f,PolyhedralFan::bergmanOfPrincipalIdeal(*i));
      f.removeAllExcept(3);
      
      IntegerVectorList l=f.getRelativeInteriorPoints();
      
      bool result;
      PolynomialSet r=checkList(l,groebnerBasis,fullNeighbourBasis,h,result, onlyCheckRays);
      if(result)return r;
    }
  fprintf(Stderr,"UNABLE TO COMPUTE STARTING CONE.\n");
  assert(0); // We could not guess an intial ideal
  return groebnerBasis;
}