#include "parser.h"
#include "printer.h"
#include "polynomial.h"
#include "division.h"
#include "buchberger.h"
#include "wallideal.h"
#include "lp.h"
#include "reversesearch.h"
#include "breadthfirstsearch.h"
#include "termorder.h"
#include "ep_standard.h"
#include "ep_xfig.h"
#include "gfanapplication.h"
#include "timer.h"
#include "matrix.h"

static Timer globalTimer("Global-timer",1);


class GraphEdge{
public:
  int source;
  int destination;
  IntegerVector direction;
  int edgeIndex;
  GraphEdge(int source_,int destination_, IntegerVector direction_):
    source(source_),
    destination(destination_),
    direction(direction_),
    edgeIndex(-1)
  {
  }
};

class GraphVertex{
public:
  bool marked;
  list<GraphEdge> edges;
  int parent;
  IntegerVector coordinates;
  GraphVertex():
    marked(false)
  {
  }
  void print(FILE *f)const
  {
    for(list<GraphEdge>::const_iterator i=edges.begin();i!=edges.end();i++)
      {
	fprintf(f,"\t[%i] %i\t",i->edgeIndex,i->destination);
	AsciiPrinter(f).printVector(i->direction);
	//	AsciiPrinter(f).printIntegerVector(i->direction);
	fprintf(f,"\n");
      }
  }
};

class MyTarget: public EnumerationTarget
{
  PolynomialSetList &l;
 public:
  MyTarget(PolynomialSetList &l_):
    l(l_)
  {
  }
  void beginEnumeration(const PolynomialSet &groebnerBasis)
  {
  }
  void endEnumeration()
  {
  }
  bool basis(const PolynomialSet &groebnerBasis)
  {
    PolynomialSet temp(groebnerBasis);;
    temp.sort_();
    l.push_back(temp);
    return true;
  }
};

typedef vector<GraphVertex> VertexVector;

class Graph
{
  VertexVector vertices;
  int dimension;
  int numberOfEdges;
  vector<GraphEdge*> edgePointers;
public:
  Graph(PolynomialSet const &g)
  {
    dimension=g.numberOfVariablesInRing();
    fprintf(Stderr,"dimension%i\n",dimension);
    PolynomialSetList bases;
    MyTarget t(bases);
    LexicographicTermOrder termOrder;
    //    StandardGradedLexicographicTermOrder termOrder;
	ReverseSearch ea(termOrder);
    ea.setEnumerationTarget(&t);
    ea.enumerate(g);

    //   AsciiPrinter(Stderr).printPolynomialSetList(bases);


    vertices=VertexVector(bases.size());
    int j=0;
    for(PolynomialSetList::const_iterator i=bases.begin();i!=bases.end();i++)
      {
	fprintf(Stderr,"%i\n",j);
	IntegerVectorList flipable=wallFlipableNormals(*i,false);
	for(IntegerVectorList::const_iterator k=flipable.begin();k!=flipable.end();k++)
	  {
	    {
	      PolynomialSet g2;
	      /*fprintf(Stderr,"flipping\n");
	      AsciiPrinter(Stderr).printPolynomialSet(*i);
	      AsciiPrinter(Stderr).printVector(*k);*/
	      g2=flip(*i,*k);
	      g2.sort_();
	      /*fprintf(Stderr,"done flipping\n");
		AsciiPrinter(Stderr).printPolynomialSet(g2);*/
	      int index=0;
	      PolynomialSetList::const_iterator destination;
	      for(destination=bases.begin();destination!=bases.end();destination++)
		{
		  if(destination->markedTermIdeal()==g2.markedTermIdeal())break;
		  index++;
		}
	      assert(destination!=bases.end());
	      /*	      fprintf(Stderr,"Matching: %i :\n",index);
	      AsciiPrinter(Stderr).printPolynomialSet(g2);
	      AsciiPrinter(Stderr).printPolynomialSet(*destination);*/
	      vertices[j].edges.push_back(GraphEdge(j,index,*k));
	    }
	  }
	j++;
      }

    print();

    int edgeIndex=0;
    for(int i=0;i<vertices.size();i++)
      {
	for(list<GraphEdge>::iterator j=vertices[i].edges.begin();j!=vertices[i].edges.end();j++)
	  {
	    if(j->edgeIndex==-1)
	      {
		list<GraphEdge>::iterator k;
		for(k=vertices[j->destination].edges.begin();k!=vertices[j->destination].edges.end();k++)
		  {
		    if(k->destination==i)break;
		  }
		if(k==vertices[j->destination].edges.end())
		  print();
		assert(k!=vertices[j->destination].edges.end());
		assert(k->edgeIndex==-1);
		k->edgeIndex=j->edgeIndex=edgeIndex++;
		k->direction=-j->direction;
		edgePointers.push_back(&(*j));
	      }
	  }
      }
    numberOfEdges=edgeIndex;
  }
  void markVertices(bool b)
  {
    for(VertexVector::iterator i=vertices.begin();i!=vertices.end();i++)i->marked=b;
  }
  void print(FILE *f=Stderr)
  {
    int index=0;
    for(VertexVector::const_iterator i=vertices.begin();i!=vertices.end();i++)
      {
	fprintf(f,"%i:\t",index);
	AsciiPrinter(f).printVector(i->coordinates);
	fprintf(f,"\n");
	i->print(f);
	index++;
      }
  }
  void applyCoefficients(IntegerVector const &e)
  {
    markVertices(false);
    list<int> active;
    vertices[0].coordinates=IntegerVector(dimension);
    active.push_back(0);
    vertices[0].parent=0;
    while(!active.empty())
      {
	int c=active.front();
	active.pop_front();
	for(list<GraphEdge>::const_iterator i=vertices[c].edges.begin();i!=vertices[c].edges.end();i++)
	  {
	    if(!vertices[i->destination].marked)
	      {
		vertices[i->destination].coordinates=vertices[c].coordinates+e[i->edgeIndex]*i->direction;
		active.push_back(i->destination);
		vertices[i->destination].marked=true;
		vertices[i->destination].parent=c;
	      }
	  }
      }
  }
  IntegerVector findPath(int v)//, vector<list<GraphEdge>::const_iterator> *iterators)
  {
    IntegerVector ret(numberOfEdges);

    while(v!=0)
      {
	int edge=-1;
	for(list<GraphEdge>::const_iterator i=vertices[v].edges.begin();i!=vertices[v].edges.end();i++)
	  if(i->destination==vertices[v].parent)
	    {
	      edge=i->edgeIndex;
	      //	      if(iterators)(*iterators)[edge]=i;
	      break;
	    }
	assert(edge!=-1);
	ret[edge]+=1;
	v=vertices[v].parent;
      }

    return ret;
  }

  void computeCoordinazationOld()
  {
    IntegerVectorList cycles;

    IntegerVectorList lp;
    IntegerVector e(numberOfEdges);
    for(int i=0;i<numberOfEdges;i++)e[i]=1;
    do
      {
	int troubleEdgeIndex=0;
	list<GraphEdge>::const_iterator troubleEdgeIterator;

	for(int i=0;i<e.size();i++)
	  {
	    assert(e[i]>0);
	  }
	for(IntegerVectorList::const_iterator i=lp.begin();i!=lp.end();i++)
	  {
	    assert(dot(*i,e)==0);
	  }

	applyCoefficients(e);
	print();

	bool ok=true;
	int v1,v2;
	for(int v=0;v<vertices.size();v++)
	  {
	    for(list<GraphEdge>::const_iterator edge=vertices[v].edges.begin();edge!=vertices[v].edges.end();edge++)
	      {
		v1=v;
		v2=edge->destination;
		troubleEdgeIndex=edge->edgeIndex;
		troubleEdgeIterator=edge;
		if(!((vertices[v2].coordinates-vertices[v1].coordinates)==e[edge->edgeIndex]*edge->direction))
		  {
		    ok=false;
		    break;
		  }
	      }
	    if(!ok)break;
	  }
	if(ok)
	  {
	    fprintf(Stderr,"Coordinazation found!!\n");
	    return;
	  }
	fprintf(Stderr,"Coordinazation problem v1: %i  v2: %i\n",v1,v2);
	
	//	vector<list<GraphEdge>::const_iterator> i1(numberOfEdges),i2(numberOfEdges);
	IntegerVector p1=findPath(v1);//,&i1);
	IntegerVector p2=findPath(v2);//,&i2);
	fprintf(Stderr,"path1:");
	AsciiPrinter(Stderr).printVector(p1);
	fprintf(Stderr,"\n");
	fprintf(Stderr,"path2:");
	AsciiPrinter(Stderr).printVector(p2);
	fprintf(Stderr,"\n");

	IntegerVector ps=p1-p2;
	ps[troubleEdgeIndex]-=1;
	
	cycles.push_back(ps);

	for(int i=0;i<dimension;i++)
	  {
	    //	    vector<list<GraphEdge>::const_iterator>::const_iterator a1=i1.begin();
	    //	    vector<list<GraphEdge>::const_iterator>::const_iterator a2=i2.begin();

	    IntegerVector v(numberOfEdges);
	    for(int j=0;j<numberOfEdges;j++)
	      {
		//	int oldv=v[j];
		//	int s1=0,s2=0;
		//		if(p1[j]){s1=(*a1)->direction[i];}
		//		if(p2[j]){s2=(*a2)->direction[i];}
		//		v[j]=s1-s2;
		//if(j==troubleEdgeIndex)v[j]-=troubleEdgeIterator->direction[i];
		
 
		//oldv-=v[j];
		//	fprintf(Stderr,"directionsize:%i\n",((*a1)->direction.size()));
		v[j]+=ps[j]*(edgePointers[j]->direction[i]);
	    //		a1++;
	    //	a2++;
	      }
	    lp.push_back(v);
	  }
	fprintf(Stderr,"LP program (size=%i):\n",lp.size());
	AsciiPrinter(Stderr).printVectorList(lp);
	fprintf(Stderr,"cycles:\n");
	AsciiPrinter(Stderr).printVectorList(cycles);
	{
	  IntegerVector used(cycles.begin()->size());
	  for(IntegerVectorList::const_iterator i=cycles.begin();i!=cycles.end();i++)
	    for(int j=0;j<used.size();j++)
	      if((*i)[j])used[j]=1;
	  fprintf(Stderr,"UsedEges:\n");
	  AsciiPrinter(Stderr).printVector(used);	
	  fprintf(Stderr,"\n");
	  for(int j=0;j<used.size();j++)
	    if(used[j])fprintf(Stderr,"%i\n",j);
	}
      }
    while(positiveVectorInKernel(lp,&e));
  }
  void computeCoordinazation()
  {
    //Compute the a basis for the flow space
    IntegerVectorList incidenceMatrix;
    for(vector<GraphEdge*>::const_iterator i=edgePointers.begin();i!=edgePointers.end();i++)
      {
	IntegerVector e(vertices.size());
	e[(*i)->source]=-1;
	e[(*i)->destination]=1;
	incidenceMatrix.push_back(e);
      }

    AsciiPrinter(Stderr).printVectorList(incidenceMatrix);

    PolyhedralCone c(IntegerVectorList(),rowsToIntegerMatrix(incidenceMatrix).transposed().getRows(),edgePointers.size());
    PolyhedralCone d=c.dualCone();
    IntegerVectorList cycles2=d.getLinealitySpace();
    AsciiPrinter(Stderr).printVectorList(cycles2);

    IntegerVectorList A;
    for(IntegerVectorList::const_iterator i=cycles2.begin();i!=cycles2.end();i++)
      {
	for(int j=0;j<dimension;j++)
	  {
	    IntegerVector row(i->size());
	    int K=0;
	    for(vector<GraphEdge*>::const_iterator k=edgePointers.begin();k!=edgePointers.end();k++)
	      {
		row[K]=(*i)[K]*(*k)->direction[j];
		K++;
	      }
	    A.push_back(row);
	  }
      }
    AsciiPrinter(Stderr).printVectorList(A);

    PolyhedralCone ker(IntegerVectorList(),A,numberOfEdges);
    if(ker.containsPositiveVector())
      {
	fprintf(Stderr,"Fan is regular\n");
      }
    else
      {
	fprintf(Stderr,"Fan is non-regular\n");
      }



    IntegerVectorList cycles;
    IntegerVectorList lp;
    IntegerVector e(numberOfEdges);
    for(int i=0;i<numberOfEdges;i++)e[i]=1;
    do
      {
	int troubleEdgeIndex=0;
	list<GraphEdge>::const_iterator troubleEdgeIterator;

	for(int i=0;i<e.size();i++)
	  {
	    assert(e[i]>0);
	  }
	for(IntegerVectorList::const_iterator i=lp.begin();i!=lp.end();i++)
	  {
	    assert(dot(*i,e)==0);
	  }

	applyCoefficients(e);
	print();

	bool ok=true;
	int v1,v2;
	for(int v=0;v<vertices.size();v++)
	  {
	    for(list<GraphEdge>::const_iterator edge=vertices[v].edges.begin();edge!=vertices[v].edges.end();edge++)
	      {
		v1=v;
		v2=edge->destination;
		troubleEdgeIndex=edge->edgeIndex;
		troubleEdgeIterator=edge;
		if(!((vertices[v2].coordinates-vertices[v1].coordinates)==e[edge->edgeIndex]*edge->direction))
		  {
		    ok=false;
		    break;
		  }
	      }
	    if(!ok)break;
	  }
	if(ok)
	  {
	    fprintf(Stderr,"Coordinazation found!!\n");
	    return;
	  }
	fprintf(Stderr,"Coordinazation problem v1: %i  v2: %i\n",v1,v2);
	
	//	vector<list<GraphEdge>::const_iterator> i1(numberOfEdges),i2(numberOfEdges);
	IntegerVector p1=findPath(v1);//,&i1);
	IntegerVector p2=findPath(v2);//,&i2);
	fprintf(Stderr,"path1:");
	AsciiPrinter(Stderr).printVector(p1);
	fprintf(Stderr,"\n");
	fprintf(Stderr,"path2:");
	AsciiPrinter(Stderr).printVector(p2);
	fprintf(Stderr,"\n");

	IntegerVector ps=p1-p2;
	ps[troubleEdgeIndex]-=1;
	
	cycles.push_back(ps);

	for(int i=0;i<dimension;i++)
	  {
	    //	    vector<list<GraphEdge>::const_iterator>::const_iterator a1=i1.begin();
	    //	    vector<list<GraphEdge>::const_iterator>::const_iterator a2=i2.begin();

	    IntegerVector v(numberOfEdges);
	    for(int j=0;j<numberOfEdges;j++)
	      {
		//	int oldv=v[j];
		//	int s1=0,s2=0;
		//		if(p1[j]){s1=(*a1)->direction[i];}
		//		if(p2[j]){s2=(*a2)->direction[i];}
		//		v[j]=s1-s2;
		//if(j==troubleEdgeIndex)v[j]-=troubleEdgeIterator->direction[i];
		
 
		//oldv-=v[j];
		//	fprintf(Stderr,"directionsize:%i\n",((*a1)->direction.size()));
		v[j]+=ps[j]*(edgePointers[j]->direction[i]);
	    //		a1++;
	    //	a2++;
	      }
	    lp.push_back(v);
	  }
	fprintf(Stderr,"LP program (size=%i):\n",lp.size());
	AsciiPrinter(Stderr).printVectorList(lp);
	fprintf(Stderr,"cycles:\n");
	AsciiPrinter(Stderr).printVectorList(cycles);
	{
	  IntegerVector used(cycles.begin()->size());
	  for(IntegerVectorList::const_iterator i=cycles.begin();i!=cycles.end();i++)
	    for(int j=0;j<used.size();j++)
	      if((*i)[j])used[j]=1;
	  fprintf(Stderr,"UsedEges:\n");
	  AsciiPrinter(Stderr).printVector(used);	
	  fprintf(Stderr,"\n");
	  for(int j=0;j<used.size();j++)
	    if(used[j])fprintf(Stderr,"%i\n",j);
	}
      }
    while(positiveVectorInKernel(lp,&e));
  }
};

class ConstructionApplication : public GFanApplication
{
  FieldOption theFieldOption;
public:
  bool includeInDefaultInstallation()
  {
    return false;
  }
  ConstructionApplication()  {
    registerOptions();
  }    

  char *name()
  {
    return "_construction";
  }

  int main()
  {
    Field::printList(Stderr);
    LpSolver::printList(Stderr);
    lpSetSolver("cddgmp");
     

    PolynomialSet g;
    PolynomialSet gOld;

    g=FileParser(Stdin).parsePolynomialSet();
    
    AsciiPrinter(Stderr).printPolynomialSet(g);

    globalTimer.on();

    fprintf(Stderr,"Computing Groebner Basis...\n");
    //        buchberger(&g,myOrder);
    buchberger(&g,StandardGradedLexicographicTermOrder());
    AsciiPrinter(Stderr).printPolynomialSet(g);

    fprintf(Stderr,"A reduced Groebner basis has been computed\n");
      
    Graph edgeGraph(g);
    edgeGraph.computeCoordinazation();
    //    edgeGraph.print();

    globalTimer.off();
    //    if(optionPerformanceTimer.getValue())Timer::printList();
 
    return 0;
  }
  const char *helpText()
  {
    return "helppp";
  }
};

static ConstructionApplication theApplication;