#include "parser.h"
#include "printer.h"
#include "polynomial.h"
#include "division.h"
#include "lp.h"
#include "gfanapplication.h"
#include "polyhedralcone.h"

#include "polymakefile.h"
#include "determinant.h"

class TriangulateApplication : public GFanApplication
{
  FieldOption theFieldOption;
  StringOption inputOption;
  StringOption outputOption;
public:
  bool includeInDefaultInstallation()
  {
    return false;
  }
  const char *helpText()
  {
    return "This program ........\n";
  }
  TriangulateApplication():
    inputOption("-i","Specify the name of the input file.","polymake.out"),
    outputOption("-o","Specify the name of the output file.","polymake.out2")
  {
    registerOptions();
  }    

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

  void printIntList(list<int> const &v)
  {
    FILE *f=Stderr;
    fprintf(f,"{");
    for(list<int>::const_iterator i=v.begin();i!=v.end();i++)
      {
	if(i!=v.begin())fprintf(f," ");
	fprintf(f,"%i",*i);
      }
    fprintf(f,"}\n");
  }
  void printIntListList(list<list<int> > const &l)
  {
    for(list<list<int> >::const_iterator i=l.begin();i!=l.end();i++)
      printIntList(*i);
  }

  typedef list<int> Cone;
  void coneSort(Cone &c)
  {
    c.sort();
  }

  IntegerVectorList coneToVectorList(Cone const &c, IntegerMatrix const &rays)
  {
    IntegerVectorList ret;

    for(Cone::const_iterator i=c.begin();i!=c.end();i++)
      ret.push_back(rays[*i]);

    return ret;
  }

  int coneDim(Cone const &c, IntegerMatrix const &rays)
  {
    return rankOfMatrix(coneToVectorList(c,rays));
  }

  Cone firstSimplex(Cone const &c, IntegerMatrix const &rays)
  {
    Cone ret;
    int d=0;

    for(Cone::const_iterator i=c.begin();i!=c.end();i++)
      {
	ret.push_back(*i);
	if(coneDim(ret,rays)!=ret.size())ret.pop_back();
      }
    return ret;
  }

  IntegerVectorList coneComplement(Cone c, IntegerMatrix const &rays)//returns generators of orth. complement.
  {
    IntegerVectorList equations=coneToVectorList(c,rays);
    IntegerVectorList empty;
    return PolyhedralCone(empty,equations,rays.getWidth()).dualCone().getEquations();
  }

  bool isVisible(int v, Cone const &c, IntegerVectorList const &complement, IntegerMatrix const &rays)
  {
    IntegerVectorList l1=coneToVectorList(c,rays);
    l1.push_back(rays[v]);
    for(IntegerVectorList::const_iterator i=complement.begin();i!=complement.end();i++)
      {
	l1.push_back(*i);
      }
    //    AsciiPrinter(Stderr).printVectorList(l1);
    // fprintf(Stderr,"sgn=%i\n",determinantSign(l1));
    return determinantSign(l1)>0;
  }

  list<Cone> triangulateRek(int d, Cone const &c, IntegerMatrix const &rays)
  {
    list<Cone> ret;

    assert(d>1);
    if(d==2)
      {
	Cone c2;
	
	c2.push_back(*c.begin());
	c2.push_back(*(++(c.begin())));
	ret.push_back(c2);


	//fprintf(Stderr,"Base\n");
	//printIntList(c2);
	return ret;
      }
    list<Cone> boundary=triangulateRek(d-1,c,rays);

    Cone l;
    Cone::const_iterator i;
    for(i=c.begin();i!=c.end();i++)
      {
	l.push_back(*i);
	if(coneDim(l,rays)==d)break;
      }
    assert(i!=c.end());
    //    fprintf(Stderr,"A\n");

    IntegerVectorList complement;
    {
      Cone temp=*boundary.begin();
      temp.push_back(*i);
      complement=coneComplement(temp,rays);

      //      AsciiPrinter(Stderr).printVectorList(complement);

    }
    {
      int N=boundary.size();
      for(list<Cone>::const_iterator i=boundary.begin();--N>=0;i++)
	{
	  Cone temp=*i;
	  int a=temp.back();
	  temp.pop_back();
	  int b=temp.back();
	  temp.pop_back();
	  temp.push_back(a);
	  temp.push_back(b);
	  boundary.push_back(temp);
	  //  fprintf(Stderr,"Boundary lifted from lower dimension:\n");
	  //  printIntListList(boundary);
	}
    }


    for(;i!=c.end();i++)
      {
	{//We are done if we leave the subspace we are working in:
	  bool done=false;
	  for(IntegerVectorList::const_iterator j=complement.begin();j!=complement.end();j++)
	    if(dotLong(*j,rays[*i])!=0)
	      {
		done=true;
		break;
	      }
	  if(done)break;
	}

	for(list<Cone>::iterator j=boundary.begin();j!=boundary.end();j++)
	  if(isVisible(*i,*j,complement,rays))
	    {
	      //     fprintf(Stderr,"Found visible\n");
	      Cone b=*j;
	      list<Cone>::iterator tempj=j;
	      j++;
	      boundary.erase(tempj);
	      j--;
	      
	      
	      for(Cone::const_iterator k=b.begin();k!=b.end();k++)
		{
		  Cone temp;
		  for(Cone::const_iterator l=b.begin();l!=b.end();l++)
		    if(l!=k)
		      temp.push_back(*l);
		    else
		      temp.push_back(*i);
		  {
		    Cone temp2=temp;
		    coneSort(temp2);
		    bool found=false;
		    for(list<Cone>::iterator l=boundary.begin();l!=boundary.end();l++)
		      {
			Cone temp3=*l;
			coneSort(temp3);
			//fprintf(Stderr,"comparing:");
			//printIntList(temp3);
			//printIntList(temp2);
			if(temp3==temp2)
			  {
			    //  fprintf(Stderr,"erasing.");

			    boundary.erase(l);
			    found=true;
			    break;
			  }
		      }
		    if(!found)boundary.push_back(temp);


		    //   fprintf(Stderr,"New boundary:\n");
		    // printIntListList(boundary);
		    // fprintf(Stderr,".\n");

		  }
		}
	      b.push_back(*i);
	      ret.push_back(b);
	    }
      }
    return ret;
  }
  list<Cone> triangulate(Cone c, IntegerMatrix const &rays) //computes a lexicographic triangulation
  {
    coneSort(c);
    return triangulateRek(coneDim(c,rays),c,rays);
  } 

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

    PolymakeFile inFile;
    fprintf(Stderr,"Test\n");
    inFile.open(inputOption.getValue());
    fprintf(Stderr,"Test\n");


    int n=inFile.readCardinalProperty("AMBIENT_DIM");
    //int d=inFile.readCardinalProperty("DIM");
    int nRays=inFile.readCardinalProperty("N_RAYS");

    fprintf(Stderr,"%i %i\n",n,nRays);

    IntegerMatrix rays=inFile.readMatrixProperty("RAYS",nRays,n);

    vector<list<int> > cones=inFile.readMatrixIncidenceProperty("MAXIMAL_CONES");

    for(vector<list<int> >::const_iterator i=cones.begin();i!=cones.end();i++)
      printIntList(*i);

    AsciiPrinter(Stderr).printVectorList(rays.getRows());

    vector<Cone> simplicialComplex;

    for(vector<list<int> >::const_iterator i=cones.begin();i!=cones.end();i++)
      {
	//	fprintf(Stderr,"Triangulating:\n");
	//	printIntList(*i);
	list<Cone> coneList=triangulate(*i,rays);
	// fprintf(Stderr,"Result:\n");
	for(list<Cone>::const_iterator j=coneList.begin();j!=coneList.end();j++)
	  {
	    // printIntList(*j);
	    simplicialComplex.push_back(*j);
	  }
      }

    PolymakeFile outFile;
    outFile.create(outputOption.getValue(),"topaz","SimplicialComplex");

    outFile.writeCardinalProperty("N_VERTICES",nRays);
    //    outFile.writeCardinalProperty("DIM",d);
    outFile.writeIncidenceMatrixProperty("INPUT_FACES",simplicialComplex);

    outFile.close();

    return 0;
  }
};

static TriangulateApplication theApplication;