// To do: The facet detection in this file may not work for non-homogeneous ideals since facets pass the test even if they do not contain a positive vector. Check if this is the case and fix the bug.

#include "breadthfirstsearch.h"

#include "buchberger.h"
#include "wallideal.h"
#include "printer.h"
#include "lp.h"


BreadthFirstSearch::BreadthFirstSearch(const SymmetryGroup &symmetryGroup_, bool minkowski_):
  numberOfVertices(0),
  numberOfEdges(0),
  symmetryGroup(symmetryGroup_),
  minkowski(minkowski_)
{
}

class Orbit{
  const SymmetryGroup &s;
  PolynomialSet g;
  PolynomialSet initialIdeal;
  IntegerVectorList markedFacets;
  list<int> polynomialLengths(const PolynomialSet &s)
  {
    list<int> l;
    for(PolynomialSet::const_iterator i=s.begin();i!=s.end();i++)
      l.push_back(i->numberOfTerms());
    l.sort();
    return l;
  }
public:
  Orbit(const PolynomialSet &g_, const SymmetryGroup &s_):
    s(s_),
    g(g_),
    initialIdeal(g_.markedTermIdeal())
    //g(g_.markedTermIdeal())    //monomialideal
  {
    g.sort_();
  }
  const PolynomialSet &representative()const
  {
    return g;
  }
  const IntegerVectorList &getMarkedFacets()const
  {
    return markedFacets;
  }
  void markFacet(const IntegerVector &facet)
  {
    markedFacets.push_back(facet);
  }
  bool containsAndMark(const PolynomialSet &p_, const IntegerVector &v)
  {
    PolynomialSet p=p_.markedTermIdeal();

{
  PolynomialSet &g=initialIdeal;
    // simple tests
    if(p.size()!=g.size())return false;

    {
      list<int> P=polynomialLengths(p);
      list<int> G=polynomialLengths(g);
      list<int>::const_iterator i=P.begin();
      for(list<int>::const_iterator j=G.begin();j!=G.end();j++)
	{
	  if(*j!=*i)return false;
	  i++;
	}
    }
    {
      vector<int> a;

      //      sort(a.begin(),a.end());

      IntegerVector P=p.exponentsSum();
      IntegerVector G=g.exponentsSum();
      P.sort();
      G.sort();
      if(P!=G)return false;
    }
    // checking all elements
    for(SymmetryGroup::ElementContainer::const_iterator j=s.elements.begin();j!=s.elements.end();j++)
      {
	PolynomialSet q;
	for(PolynomialSet::const_iterator i=p.begin();i!=p.end();i++)
	  {
	    q.push_back(s.permutePolynomial(*i,*j));
	  }
	q.sort_();

	/*	fprintf(Stderr,"COMPARE\n");
	AsciiPrinter(Stderr).printPolynomialSet(q);
	AsciiPrinter(Stderr).printPolynomialSet(g);
	*/if(q==g)
	  {
	    //    markFacet(SymmetryGroup::composeInverse(*j,v));
	    markFacet(SymmetryGroup::compose(*j,v));
	    //markFacet(v);
	    return true;	
	  }
      }
}
    return false;
  }
  bool permutationFixesInitialIdeal(IntegerVector const &v)
  {
    PolynomialSet G=g.markedTermIdeal();
    G.sort_();
    
    PolynomialSet q;
    for(PolynomialSet::const_iterator i=G.begin();i!=G.end();i++)
      {
	q.push_back(SymmetryGroup::permutePolynomial(*i,v));
      }
    q.sort_();
    return (q==G);
  }
  int orbitSize()
  {
    int groupSize=s.elements.size();

    int numFixed=0;
    for(SymmetryGroup::ElementContainer::const_iterator j=s.elements.begin();j!=s.elements.end();j++)
      {
	//	fprintf(Stderr,"1a\n");
	if(permutationFixesInitialIdeal(*j))numFixed++;
	//	fprintf(Stderr,"1b\n");
      }

    fprintf(Stderr,"numFixed = %i\n",numFixed);
    return groupSize/numFixed;
  }
};

typedef list<Orbit> OrbitList;

class OrbitContainer{
public:
  OrbitList l;
  void push_back(const Orbit &orbit)
  {
        l.push_back(orbit); // breadth first
	//    l.push_front(orbit); // depth first
  }
  bool empty()
  {
    return l.empty();
  }
  Orbit &front()
  {
    return *l.begin();
  }
  void pop_front()
  {
    l.pop_front();
  }
  bool containsAndMark(const PolynomialSet &g, const IntegerVector &v)
  {
    for(OrbitList::iterator i=l.begin();i!=l.end();i++)
      {
	if(i->containsAndMark(g,v))return true;
      }
    return false;
  }
  void print(FILE *f)
  {
    AsciiPrinter p(f);
    for(OrbitList::const_iterator i=l.begin();i!=l.end();i++)
      {
	p.printPolynomialSet(i->representative());
	p.printVectorList(i->getMarkedFacets());
      }
  }
  int size()
  {
    return l.size();
  }
};


void BreadthFirstSearch::setSubspace(IntegerVectorList const &subspacePerp_)
{
  subspacePerp=subspacePerp_;
}


void BreadthFirstSearch::enumerate(const PolynomialSet &groebnerBasis)
{
  int numberOfClosedVertices=0;
  int numberOfVertices=0;

  targetBeginEnumeration(groebnerBasis);

  OrbitContainer active;

  active.push_back(Orbit(groebnerBasis,symmetryGroup));

  while(!active.empty())
    {
      if(!subspacePerp.empty())
	{
	  IntegerVectorList inequalities=wallInequalities(active.front().representative());
	  for(IntegerVectorList::const_iterator i=subspacePerp.begin();i!=subspacePerp.end();i++)
	    inequalities.push_front(*i);
	  IntegerVector equalitySet(inequalities.size());
	  for(int i=0;i<subspacePerp.size();i++)
	    equalitySet[i]=1;
	  if(!hasInteriorPoint(inequalities,true,&equalitySet))
	    {
	      active.pop_front();
	      continue;
	    }
	}

      {
	static int n;
	n++;
	if(!(n%10))
	  fprintf(Stderr,"%i\n",n);
      }
      //  fprintf(Stderr,"Active set:\n");
      // active.print(Stderr);
      // fprintf(Stderr,"end active set\n");
       //       fgetc(Stdin);

      OrbitList::iterator currentIter=active.l.begin();
      Orbit &current=active.front();

      PolynomialSet g=current.representative();
      

      if(!targetBasis(g))break;

      IntegerVectorList normals=wallInequalities(g);

      fprintf(Stderr,"Number of inequalities:%i\n",normals.size());
      for(IntegerVectorList::iterator i=normals.begin();i!=normals.end();i++)
	{
	  if(isFacet(normals,i))
	    {
	      if(wallContainsPositiveVector(*i))
		{
		  /*		  fprintf(Stderr,"Possible flip:\n");
		  AsciiPrinter(Stderr).printVector(*i);
		  AsciiPrinter(Stderr).printPolynomialSet(flip(g,*i));
		  */
		  bool found=false;
		  const IntegerVectorList &marked=current.getMarkedFacets();
		  for(SymmetryGroup::ElementContainer::const_iterator j=symmetryGroup.elements.begin();j!=symmetryGroup.elements.end() &&!found;j++)
		    if(current.permutationFixesInitialIdeal(*j))
		      for(IntegerVectorList::const_iterator k=marked.begin();k!=marked.end();k++)
			if(dependent(*k,SymmetryGroup::compose(*j,*i))){
			  /*  fprintf(Stderr,"Was marked with vector:");
			      AsciiPrinter(Stderr).printVector(*k);*/
			  found=true;break;
			}
		  if(!found)
		    {
		      PolynomialSet neighbour=(minkowski)?flipMinkowski(g,*i):flip(g,*i);

		      /*    fprintf(Stderr,"current:\n");
		      AsciiPrinter(Stderr).printPolynomialSet(g);
		      fprintf(Stderr,"flipped:\n");
		      AsciiPrinter(Stderr).printPolynomialSet(neighbour);
		      */
		      fprintf(Stderr,"--");
		      if(!active.containsAndMark(neighbour,*i))
			{
			  fprintf(Stderr,"-------------------------adding\n");
			  active.push_back(Orbit(neighbour,symmetryGroup));
			  active.containsAndMark(neighbour,*i);
			}
		      else
			fprintf(Stderr,"---------------------------not adding\n");
		    }
		  else fprintf(Stderr,"-- marked\n");
		}
	    }
	  else
	    {
	      IntegerVectorList::iterator temp=i;
	      temp++;
	      normals.erase(i);
	      temp--;
	      i=temp;
	    }
	}
      //      fprintf(Stderr,"New active set:\n");
      //    active.print(Stderr);
      // fprintf(Stderr,"end active set\n");


      numberOfVertices+=current.orbitSize();

      active.pop_front();  //breadth first
      //active.l.erase(currentIter); //depth first

      numberOfClosedVertices++;

      fprintf(Stderr,"Number of active vertices  %i\n",active.size());
      fprintf(Stderr,"Number of closed vertices: %i\n",numberOfClosedVertices);
      fprintf(Stderr,"Number of vertices  %i\n",numberOfVertices);

    }

  targetEndEnumeration();
}