#include "polyhedralcone.h"

#include "lp.h"
#include "subspace.h"
#include "symmetry.h"

//--------------------------------
// PolyhedralCone
//--------------------------------

static bool compareIntegerLists(IntegerVectorList const &a, IntegerVectorList const &b)
{
  assert(a.size()==b.size());

  IntegerVectorList::const_iterator B=b.begin();
  for(IntegerVectorList::const_iterator A=a.begin();A!=a.end();A++)
    {
      if(LexicographicTermOrder()(*A,*B))return true;
      if(LexicographicTermOrder()(*B,*A))return false;
      B++;
    }
  return false;
}

bool operator<(PolyhedralCone const &a, PolyhedralCone const &b)
{
  assert(a.state>=3);
  assert(b.state>=3);

  if(a.linearitySpace.size()<b.linearitySpace.size())return true;
  if(a.linearitySpace.size()>b.linearitySpace.size())return false;

  if(a.n<b.n)return true;
  if(a.n>b.n)return false;

  if(a.halfSpaces.size()<b.halfSpaces.size())return true;
  if(a.halfSpaces.size()>b.halfSpaces.size())return false;

  if(compareIntegerLists(a.linearitySpace,b.linearitySpace))return true;
  if(compareIntegerLists(b.linearitySpace,a.linearitySpace))return false;

  if(compareIntegerLists(a.halfSpaces,b.halfSpaces))return true;
  if(compareIntegerLists(b.halfSpaces,a.halfSpaces))return false;

  return false;
}

bool operator!=(PolyhedralCone const &a, PolyhedralCone const &b)
{
  return (a<b)||(b<a);
}


void PolyhedralCone::ensureStateAsMinimum(int s)
{
  if((state<1) && (s==1))
    {
      removeRedundantRows(&halfSpaces,&linearitySpace,false);
    }
  if((state<2) && (s>=2))
    {
      //      halfSpaces.sort();
      //      AsciiPrinter(Stderr).printVectorList(halfSpaces);

       if(halfSpaces.size()>25)
	 //	if(0)
	 {
	  IntegerVectorList h1;
	  IntegerVectorList h2;
	  bool a=false;
	  for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
	    {
	      if(a)
		h1.push_back(*i);
	      else
		h2.push_back(*i);
	      a=!a;
	    }
	  PolyhedralCone c1(h1,linearitySpace);
	  PolyhedralCone c2(h2,linearitySpace);
	  c1.ensureStateAsMinimum(2);
	  c2.ensureStateAsMinimum(2);
	  halfSpaces=c1.halfSpaces;
	  for(IntegerVectorList::const_iterator i=c2.halfSpaces.begin();i!=c2.halfSpaces.end();i++)
	    halfSpaces.push_back(*i);
	}
	

       //      fprintf(Stderr,"Number half spaces: %i, number of equations: %i\n",halfSpaces.size(),linearitySpace.size());
      removeRedundantRows(&halfSpaces,&linearitySpace,true);
      //   fprintf(Stderr,"done\n");
    }
  if((state<3) && (s>=3))
    {
      //      fprintf(Stderr,"Number half spaces: %i, number of equations: %i\n",halfSpaces.size(),linearitySpace.size());
      Subspace v(linearitySpace,n);
      
      linearitySpace=v.getRepresentation();
      
      for(IntegerVectorList::iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
	{
	  *i=v.canonicalizeVector(*i);
	}
      halfSpaces.sort();
      //fprintf(Stderr,"done\n");
    }
  state=s;
}

void PolyhedralCone::canonicalize()
{
  ensureStateAsMinimum(3);
}


void PolyhedralCone::findFacets()
{
  ensureStateAsMinimum(2);
}


PolyhedralCone::PolyhedralCone(int ambientDimension):
  n(ambientDimension),
  state(1)
{
}


PolyhedralCone::PolyhedralCone(IntegerVectorList const &halfSpaces_, IntegerVectorList const &linearitySpace_, int ambientDimension):
  halfSpaces(halfSpaces_),
  linearitySpace(linearitySpace_),
  state(0)
{
  n=ambientDimension;
  if(n==-1)
    {
      if(!halfSpaces_.empty())
	n=halfSpaces_.begin()->size();
      else if(!linearitySpace_.empty())
	n=linearitySpace_.begin()->size();
      else
	{
	  assert(0);
	}
    }
  for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
    {
      assert(i->size()==n);
    }
  for(IntegerVectorList::const_iterator i=linearitySpace.begin();i!=linearitySpace.end();i++)
    {
      assert(i->size()==n);
    }
  ensureStateAsMinimum(1);
  //  computeAndReduceLinearitySpace();
}


IntegerVector PolyhedralCone::getRelativeInteriorPoint()const
{
  //  ensureStateAsMinimum(1);
  assert(state>=1);
  IntegerVector ret;
  IntegerVectorList g=linearitySpace;
  int numberOfEqualities=g.size();
  g.insert(g.end(),halfSpaces.begin(),halfSpaces.end());
  IntegerVector equalitySet(g.size());
  for(int i=0;i<numberOfEqualities;i++)equalitySet[i]=1;
  ret=relativeInteriorPoint(g,&equalitySet);

  for(IntegerVectorList::const_iterator i=linearitySpace.begin();i!=linearitySpace.end();i++)
    {
      assert(dotLong(*i,ret)==0);
    }

  for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
    {
      if(!(dotLong(*i,ret)>0))
	{
	  fprintf(Stderr,"PolyhedralCone::relativeInteriorPoint() : halfSpaces not reduced or mistake in cdd interface!!!\n");
	}
    }

  return ret;
}


IntegerVectorList PolyhedralCone::getHalfSpaces()const
{
  return halfSpaces;
}


const IntegerVectorList &PolyhedralCone::getLinealitySpace()const
{
  assert(state>=1);
  return linearitySpace;
}


int PolyhedralCone::ambientDimension()const
{
  return n;
}


int PolyhedralCone::dimensionOfLinealitySpace()const
{
  return getLinealitySpace().size();
}


int PolyhedralCone::dimension()const
{
  assert(state>=1);
  //  ensureStateAsMinimum(1);
  return ambientDimension()-linearitySpace.size();
}


bool PolyhedralCone::isZero()const
{
  return dimension()==0;
}


PolyhedralCone intersection(const PolyhedralCone &a, const PolyhedralCone &b)
{
  assert(a.ambientDimension()==b.ambientDimension());
  IntegerVectorList halfSpaces=a.halfSpaces;
  halfSpaces.insert(halfSpaces.end(),b.halfSpaces.begin(),b.halfSpaces.end());
  IntegerVectorList linearitySpace=a.linearitySpace;
  linearitySpace.insert(linearitySpace.end(),b.linearitySpace.begin(),b.linearitySpace.end());
  return PolyhedralCone(halfSpaces,linearitySpace,a.ambientDimension());
}


PolyhedralCone PolyhedralCone::positiveOrthant(int dimension)
{
  IntegerVectorList halfSpaces;

  for(int i=0;i<dimension;i++)halfSpaces.push_back(IntegerVector::standardVector(dimension,i));

  return PolyhedralCone(halfSpaces,IntegerVectorList(),dimension);
}


bool PolyhedralCone::isInStateMinimum(int s)const
{
  return state>=s;
}

void PolyhedralCone::print(class Printer *p)const
{
  p->printString("Printing PolyhedralCone");
  p->printNewLine();
  p->printString("Ambient dimension: ");
  p->printInteger(n);
  p->printNewLine();
  if(isInStateMinimum(1))
    {
      p->printString("Dimension: ");
      p->printInteger(dimension());
      p->printNewLine();
    }
  p->printString("Linearity space:");
  //  p->printNewLine();
  p->printVectorList(linearitySpace);
  p->printString("Inequalities:");
  p->printVectorList(halfSpaces);
  p->printString("Relative interior point:\n");
  p->printVector(getRelativeInteriorPoint());
  p->printNewLine();
  p->printString("Done printing PolyhedralCone.");
  p->printNewLine();
}


static IntegerVector dehomogenize(IntegerVector const &v)
{
  assert(v.size()>0);

  IntegerVector ret(v.size()-1);

  for(int i=0;i<v.size()-1;i++)ret[i]=v[i];
  return ret;
}

static IntegerVectorList dehomogenize(IntegerVectorList const &l)
{
  IntegerVectorList ret;

  for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
    {
      ret.push_back(dehomogenize(*i));
    }
  return ret;
}

PolyhedralCone PolyhedralCone::withLastCoordinateRemoved()const
{
  assert(n>0);

  return PolyhedralCone(dehomogenize(halfSpaces),dehomogenize(linearitySpace));
}

bool PolyhedralCone::containsPositiveVector()const
{
  PolyhedralCone temp=intersection(*this,PolyhedralCone::positiveOrthant(n));
  IntegerVector v=temp.getRelativeInteriorPoint();
  return v.isPositive();
}


int PolyhedralCone::dimensionOfLargestContainedSubspace()const
{
  if(halfSpaces.empty())return dimension();
  IntegerVectorList a;
  PolyhedralCone temp(a,halfSpaces);
  temp=intersection(temp,*this);

  return temp.dimension();
}


bool PolyhedralCone::contains(IntegerVector const &v)const
{
  for(IntegerVectorList::const_iterator i=linearitySpace.begin();i!=linearitySpace.end();i++)
    {
      if(dotLong(*i,v)!=0)return false;
    }

  for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
    {
      if(dotLong(*i,v)<0)return false;
    }

  return true;
}


bool PolyhedralCone::contains(PolyhedralCone const &c)const
{
  PolyhedralCone c2=intersection(*this,c);
  PolyhedralCone c3=c;
  c2.canonicalize();
  c3.canonicalize();
  return !(c2!=c3);
}


bool PolyhedralCone::containsRelatively(IntegerVector const &v)const
{
  assert(state>=1);
  for(IntegerVectorList::const_iterator i=linearitySpace.begin();i!=linearitySpace.end();i++)
    {
      if(dotLong(*i,v)!=0)return false;
    }

  for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
    {
      if(dotLong(*i,v)<=0)return false;
    }

  return true;
}


PolyhedralCone PolyhedralCone::permuted(IntegerVector const &v)const
{
  PolyhedralCone ret(SymmetryGroup::permuteIntegerVectorList(halfSpaces,v),SymmetryGroup::permuteIntegerVectorList(linearitySpace,v),n);
  if(state>=1)ret.state=1;
  if(state>=2)ret.state=2;

  ret.ensureStateAsMinimum(state);
  return ret;
}


IntegerVectorList PolyhedralCone::extremeRays()const
{
  assert(dimension()==ambientDimension());

  IntegerVectorList ret;
  IntegerVectorList indices=extremeRaysInequalityIndices(halfSpaces);
  for(IntegerVectorList::const_iterator i=indices.begin();i!=indices.end();i++)
    {
      IntegerVectorList equations;
      for(int j=0;j<i->size();j++)
	{
	  IntegerVectorList::const_iterator a=halfSpaces.begin();
	  for(int k=0;k<(*i)[j];k++)
	    {
	      assert(a!=halfSpaces.end());
	      a++;
	    }
	  assert(a!=halfSpaces.end());
	  equations.push_back(*a);
	}
      ret.push_back(PolyhedralCone(halfSpaces,equations).getRelativeInteriorPoint());
    }
  return ret;
}


bool PolyhedralCone::isSimplicial()const
{
  assert(state>=2);
  
  //  ensureStateAsMinimum(2);
  //  AsciiPrinter P(Stderr);
  //  print(&P);
  return dimensionOfLinealitySpace()+getHalfSpaces().size()+dimensionOfLargestContainedSubspace()==n;
}


bool PolyhedralCone::checkDual(PolyhedralCone const &c)const
{
  assert(dimensionOfLargestContainedSubspace()+c.dimension()==ambientDimension());

  for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
    {
      assert(c.contains(*i));
    }
  for(IntegerVectorList::const_iterator i=linearitySpace.begin();i!=linearitySpace.end();i++)
    {
      assert(c.contains(*i));
    }
  return true;
}


PolyhedralCone PolyhedralCone::dualCone()const
{
  assert(state>=1);

  IntegerVectorList dualInequalities,dualEquations;

  dual(ambientDimension(),halfSpaces,linearitySpace,&dualInequalities,&dualEquations);

  PolyhedralCone ret(dualInequalities,dualEquations);

  ret.ensureStateAsMinimum(state);
  //  ret.canonicalize();


  assert(checkDual(ret));
  assert(ret.checkDual(*this));

  return ret;
}