/*
 * app_tropicalprevarietycomponents.cpp
 *
 *  Created on: 7 Feb 2021
 *      Author: anders
 */




#include <iostream>
#include <fstream>
#include <stdlib.h>
#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 "termorder.h"
#include "ep_standard.h"
#include "ep_xfig.h"
#include "polyhedralcone.h"
#include "gfanapplication.h"
#include "saturation.h"
#include "field_rationals.h"
#include "field_zmodpz.h"
#include "field_rationalfunctions.h"
#include "symmetry.h"
#include "linalg.h"
#include "fieldlp.h"
#include "integer.h"
#include "polynomialgcd.h"
#include "packedmonomial.h"
#include "gfanlib_zcone.h"
#include "gfanlib_tableau.h"
#include "gfanlib_hypersurfaceintersection.h"
#include "gfanlib_circuittableint.h"
#include "gfanlib_mixedvolume.h"
#include "lll.h"
#include "gfanlib_circuittableinteger.h"
#include "tropical2.h"

using namespace gfan;

class TropicalPrevarietyComponentsApplication : public GFanApplication
{
	SimpleOption fibersOption;
	StringOption saveAs;
	StringOption loadFrom;
	IntegerOption optionNThreads;
	IntegerOption optionHalfspace;
  IntegerOption optionNumberOfDimensionsToClear;
public:
  bool includeInDefaultInstallation()
  {
    return false;
  }
  const char *helpText()
  {
    return "This program computes the tropical prevariety of a set of polynomials given in a polynomialring with first variable being regarded as a field element with valuation 1.\n";
  }
  TropicalPrevarietyComponentsApplication():
	  fibersOption("--fibers","Iterate over fibers of a coordinate projection. In this case the projection happens on the last variable of the ring."),
	  optionNThreads("-j","Number of threads",40),
	  optionNumberOfDimensionsToClear("-c","Number of dimensions to remove in global pointedness test.",0),
	  saveAs("--saveas","Specify filename for saving intermediate result."),
	  loadFrom("--loadfrom","Specify filename for loading intermediate result."),
	  optionHalfspace("--half","specicfy right hand side of a strict inequality")
  {
	  optionHalfspace.hide();
	  registerOptions();
  }

  const char *name()
  {
    return "_tropicalprevarietycomponents";
  }

  template<typename typ>
  gfan::Matrix<typ> convertMatrix(IntMatrix const &m)
  {
	  gfan::Matrix<typ> ret(m.getHeight(),m.getWidth());
	  for(int i=0;i<m.getHeight();i++)
		  for(int j=0;j<m.getWidth();j++)
			  ret[i][j]=typ(m[i][j]);
	  return ret;
  }

  IntMatrix rowsToIntegerMatrix2(IntegerVectorList const &l)//taken from app_mixedvolume.cpp
  {
	  assert(l.size());
	  IntMatrix ret(l.size(),l.front().size());
	  IntegerVectorList::const_iterator I=l.begin();
	  for(int i=0;i<ret.getHeight();i++,I++)
		  for(int j=0;j<ret.getWidth();j++)
			  ret[i][j]=(*I)[j];
	  return ret;
  }

/*template<typename typ>
int numberOfEdges(gfan::Matrix<typ> const &A)
{
	int ret=0;
	for(int i=0;i<A.getHeight();i++)
		for(int j=0;j<i;j++)
		{
			gfan::Matrix<typ> B(0,A.getWidth());
			for(int k=0;k<A.getHeight();k++)B.appendRow(A[i].toVector()-A[k].toVector());
			for(int k=0;k<A.getHeight();k++)B.appendRow(A[j].toVector()-A[k].toVector());
			auto C=Cone<typ>(B.transposed());
			ret+=C.getDimension()==C.getAmbientDimension()-1;
		}
	std::cerr<<ret<<"\n";
	return ret;
}*/

	class EquivalenceRelation{
		vector<int> parentList;
		vector<int> classNamesOfRoots;//only valid to look up with root indices
	public:
		void grow(int M)
		{
			for(int a=parentList.size();a<M;a++)
			{
				parentList.push_back(a);
				classNamesOfRoots.push_back(-1);
			}
		}
		void rek(int i, int j)//let all in iths path point to j
		{
			int old=parentList[i];
			parentList[i]=j;
			if(old!=i)
				rek(old,j);
		}
		void makeEquivalent(int i, int j)
		{
			grow(i+1);
			grow(j+1);
			rek(i,j);
		}
		bool isRoot(int i)
		{
			return i==parentList[i];
		}
		int getRepresentative(int i)
		{
			if(isRoot(i))return i;
			return getRepresentative(parentList[i]);
		}
		int getClassName(int i)
		{
			assert(isRoot(i));
			return classNamesOfRoots[i];
		}
		int buildClassNames()
		{
			int numberOfClasses=0;

			for(int i=0;i<parentList.size();i++)
				if(isRoot(i))
					classNamesOfRoots[i]=numberOfClasses++;

			return numberOfClasses;
		}
	};


	template<typename typ>
	bool nonEmptyIntersectionOfDualsInteriors(vector<Cone<typ>> &cones)
	{
		if(cones.size()==0)return true;
		int ambientDimension=cones[0].getAmbientDimension();
		// The duals are generated by the equations as line generators and inequalities as half line generators.

		// We write up a system of linear inequalities
		// -x=A_i\lambda_i with an i for each cone in cones
		// with some lambdas being strictly positive
		// There is one lambda entry for each generator of a dual cone.
		// The variables -x_1,...,-x_n represent the supposed common point in the interior of the dual points

		// We first
		//  * collect all generators of the duals,
		//  * count the number of lambdas
		//  * and the number of lambdas that have positivity constraints
		vector<pair<gfan::Matrix<typ>, gfan::Matrix<typ> > > temp;

		int numberOfVariables=ambientDimension;
		int numberOfPositiveVariables=0;
		for(auto &c:cones)
		{
			auto eq=c.getEquations().transposed();
			auto ineq=c.getInequalities().transposed();
			temp.push_back(make_pair(ineq,eq));
			numberOfVariables+=eq.getHeight();
			numberOfVariables+=ineq.getHeight();
			numberOfPositiveVariables+=ineq.getHeight();
		}

		// We then set up the equations and strict inequalities, while there are no non-strict

	  std::cerr<<"ambientDimension:"<<ambientDimension<<"\n";
	  std::cerr<<"cones.size():"<<(int)cones.size()<<"\n";
	  std::cerr<<"numberOfVar"<<numberOfVariables<<"\n";

	  int numberOfEquations=ambientDimension*cones.size();
	  std::cerr<<"numerOfEq"<<numberOfEquations<<"\n";
	  std::cerr<<"numberOfPos"<<numberOfPositiveVariables<<"\n";

		gfan::Matrix<typ> equations(numberOfVariables,numberOfEquations);
		gfan::Matrix<typ> nonStrict(numberOfVariables,0);
		gfan::Matrix<typ> strict(numberOfVariables,numberOfPositiveVariables);
		int i=0;
		int equationsIndex=ambientDimension;
		int l=0;
		for(auto &c:cones)
		{
			equations.setSubMatrix(0,i*ambientDimension,ambientDimension,(i+1)*ambientDimension,gfan::Matrix<typ>::identity(ambientDimension));
			equations.setSubMatrix(equationsIndex,i*ambientDimension,equationsIndex+temp[i].second.getHeight(),(i+1)*ambientDimension,temp[i].second);
			equationsIndex+=temp[i].second.getHeight();
			equations.setSubMatrix(equationsIndex,i*ambientDimension,equationsIndex+temp[i].first.getHeight(),(i+1)*ambientDimension,temp[i].first);
			for(int k=0;k<temp[i].first.getHeight();k++)
				strict[equationsIndex+k][l++]=typ(1);
			equationsIndex+=temp[i].first.getHeight();
			i++;
		}

/*		std::cerr<<nonStrict<<"\n";
		std::cerr<<M<<"\n";
		std::cerr<<Strict<<"\n";
*/

		// We use a HalfOpenCone to answer the answer the question

		HalfOpenCone<typ> solutionSet(nonStrict,equations,strict);

		return !solutionSet.isEmpty();
	}
	template<typename typ>
	bool sumOfConesContains(vector<Cone<typ>> v, Vector<typ> const &w)
	{
		int n=w.size();
		vector<gfan::Matrix<typ>> raysVector;
		vector<gfan::Matrix<typ>> linesVector;
		int nLines=0;
		int nRays=0;
		for(auto &i:v)
		{
			assert(n==i.getAmbientDimension());
			auto R=i.getRays();
			R.normalizeRows();
			raysVector.push_back(R);
			auto L=i.getLinealitySpace();
			L.normalizeRows();
			linesVector.push_back(L);
			nLines+=linesVector.back().getHeight();
			nRays+=raysVector.back().getHeight();
		}

//		std::cerr<<"RAYS\n";
//		for(auto &a:raysVector)std::cerr<<a;
//		std::cerr<<"LINES\n";
//		for(auto &a:linesVector)std::cerr<<a;
		gfan::Matrix<typ> nonStrict(1+nRays+nLines,nRays);
		for(int i=0;i<nRays;i++)nonStrict[i+1][i]=1;
		gfan::Matrix<typ> equations(1+nRays+nLines,n);
		equations.setSubMatrix(0,0,1,n,gfan::Matrix<typ>::rowVectorMatrix(w));
		int pos=1;
		for(int i=0;i<v.size();i++)
		{
			equations.setSubMatrix(pos,0,pos+raysVector[i].getHeight(),n,raysVector[i]);
			pos+=raysVector[i].getHeight();
		}
		for(int i=0;i<v.size();i++)
		{
			equations.setSubMatrix(pos,0,pos+linesVector[i].getHeight(),n,linesVector[i]);
			pos+=linesVector[i].getHeight();
		}
		gfan::Matrix<typ> strict(1+nRays+nLines,1);
		strict[0][0]=-1;
	  //	  std::cerr<<nonStrict.toString()<<"\n";
	  //	  std::cerr<<equations.toString()<<"\n";
	  //	  std::cerr<<strict.toString()<<"\n";
		HalfOpenCone<typ> C(nonStrict,equations,strict);
	  //std::cerr<<C.toString();
		return !C.isEmpty();
	}
	template<typename typ>
	bool sumOfHalfOpenConesContainsOrigin(vector<HalfOpenCone<typ>> const &v)
	{
		if(v.size()==0)return true;// Note that while this is correct it may not be the desired behavior.
		std::vector<Cone<typ>> l;
		for(auto &c:v)
			l.push_back(c.lifted);
		return sumOfConesContains(l,-gfan::Vector<typ>::standardVector(l[0].getAmbientDimension(),l[0].getAmbientDimension()-1));
	}
	template<typename typ>
	class Fraction{
	public:
		typ numerator;
		typ denominator;
		Fraction():
			numerator(0),
			denominator(1)
		{
		}
		Fraction(typ const &numerator_,typ const &denominator_):
			numerator(numerator_),
			denominator(denominator_)
		{
			assert(denominator.isPositive());
			Vector<typ> v(2);
			v[0]=denominator;
			v[1]=numerator;
			auto d=gcd(v);
			assert(d.isPositive());
			denominator=denominator/d;
			numerator=numerator/d;
		}
		bool operator<(Fraction<typ> const &b)const
		{
			return numerator*b.denominator<b.numerator*denominator;
		}
		string toString()const
		{
			return numerator.toString()+string("/")+denominator.toString();
		}
	};

  
	template <typename typ>
	class IntervalBound
	{
		bool isFinite;
		bool isPlusInfinity;
		bool isMinusInfinity;
		Fraction<typ> value;
		IntervalBound():
			isFinite(false),
			isPlusInfinity(false),
			isMinusInfinity(false)
		{
		}
	public:
	  bool getIsFinite()const
	  {
	    return isFinite;
	  }
	  Fraction<typ> getValue()const
	  {
	    assert(isFinite);
	    return value;
	  }
		static IntervalBound plusInfinity()
		{
			IntervalBound ret;
			ret.isPlusInfinity=true;
			return ret;
		}
		static IntervalBound minusInfinity()
		{
			IntervalBound ret;
			ret.isMinusInfinity=true;
			return ret;
		}
		IntervalBound(Fraction<typ> const &value_):
			value(value_),
			isFinite(true),
			isPlusInfinity(false),
			isMinusInfinity(false)
		{
		}
		bool operator<(IntervalBound const &b)const//comparing infinities give false
		{
			if(isPlusInfinity)
				return false;
			if(b.isPlusInfinity)
				return true;
			if(b.isMinusInfinity)
				return false;
			if(isMinusInfinity)
				return true;
			return value<b.value;
		}
		friend IntervalBound min(IntervalBound const &a, IntervalBound const &b)
		{
			if(a<b)return a;
			return b;
		}
		friend IntervalBound max(IntervalBound const &a, IntervalBound const &b)
		{
			if(b<a)return a;
			return b;
		}
		string toString()const
		{
			if(isPlusInfinity)
				return string("infinity");
			if(isMinusInfinity)
				return string("-infinity");
			return value.toString();
		}
	};

	// For simplicity intervals are open unless start==end
	template<typename typ>
	class Interval{
	public:
		IntervalBound<typ> start;
		IntervalBound<typ> end;
		bool empty;
	public:
		Interval():
			start(IntervalBound<typ>::minusInfinity()),
			end(IntervalBound<typ>::plusInfinity()),
			empty(false)
	{
	}
		Interval(IntervalBound<typ> const &start_,IntervalBound<typ> const &end_):
			start(start_),
			end(end_)
		{
			empty=!(start<end);
		}
		static Interval emptyInterval()
		{
			Interval ret(IntervalBound<typ>::plusInfinity(),IntervalBound<typ>::minusInfinity());
			ret.empty=true;
			return ret;
		}
		Interval intersectedWith(Interval const &b)const
		{
			return Interval(max(start,b.start),min(end,b.end));
		}
		string	toString()const
		{
			return string("[")+start.toString()+string(",")+end.toString()+string("]");
		}
	  Fraction<typ> getInteriorPoint()const
	  {
	    if(start.getIsFinite())
	      {
		if(end.getIsFinite())
		  {
		    auto A=start.getValue();
		    auto B=end.getValue();
		    return Fraction<typ>(A.numerator*B.denominator+A.denominator*B.numerator,A.denominator*B.denominator+A.denominator*B.denominator);
		  }
		else
		  {
		    auto temp=start.getValue();
		    return Fraction<typ>(temp.numerator+temp.denominator,temp.denominator);
		  }
	      }
	    else
	      {
		if(end.getIsFinite())
		  {
		    auto temp=end.getValue();
		    return Fraction<typ>(temp.numerator-temp.denominator,temp.denominator);
		  }
		else
		  {
		    return Fraction<typ>(0,1);
		  }
	      }
	  }
	};

	template<typename typ>
	class ProjectionIterator{
		Cone<typ> recessionOfNonTrivialFiber;
	  gfan::Matrix<typ> generatorsOfRecession;
	public:
		HalfOpenCone<typ> coneOverThis;
		bool plusInfinityCovered;
		bool minusInfinityCovered;
		int parameterIndex; //This should maybe get renamed to projectionCoordinate
		int separatorIterator; //runs from zero to separators.size()
		vector<Fraction<typ>> separators;
		ProjectionIterator(HalfOpenCone<typ> const &coneOverThis_, int parameterIndex_):
			coneOverThis(coneOverThis_),
			parameterIndex(parameterIndex_),
			plusInfinityCovered(false),
			minusInfinityCovered(false)
		{
			auto closure=coneOverThis.closure();
			{
				auto temp=closure.intersectedWithIthCoordinateHyperplaneEmbedded(parameterIndex);
				recessionOfNonTrivialFiber=temp.intersectedWithIthCoordinateHyperplaneEmbedded(0);
				generatorsOfRecession=recessionOfNonTrivialFiber.getRays();
				//std::cerr<<"Number of ray generators"<<generatorsOfRecession.getHeight()<<"\n";
			}
//			std::cerr<<"Linality dim of closure:"<<closure.getDimensionOfLinealitySpace()<<"\n";
			bool linealityIsNonZero=closure.getDimensionOfLinealitySpace()!=0;
			bool doesHaveSubdivision=true;
			if(linealityIsNonZero)
			{
				// It is possible that the lineality space projects to a single point - in which case the image of the projection indeed gets subdivided.
				auto hyperplane=Cone<typ>::hyperPlane(Vector<typ>::standardVector(closure.getAmbientDimension(),0));
				auto temp=intersection(hyperplane,closure);
				doesHaveSubdivision=temp.isImplied(Vector<typ>::standardVector(closure.getAmbientDimension(),parameterIndex));
				if(!doesHaveSubdivision)
					plusInfinityCovered=minusInfinityCovered=true;
			}

			if(doesHaveSubdivision)
			{
				auto rays=closure.getRays();
//				std::cerr<<"RAYS:\n"<<rays<<"\n";
				int infIndex=0;//rays.getWidth()-1;
				set<Fraction<typ>> separators1;
				for(int i=0;i<rays.getHeight();i++)
				{
					if(!rays[i][infIndex].isZero())
						separators1.insert(Fraction<typ>(rays[i][parameterIndex],-rays[i][infIndex]));
					else
					{
						if(rays[i][parameterIndex].isPositive())plusInfinityCovered=true;  // Did we get the signs correct here?
						if(rays[i][parameterIndex].isNegative())minusInfinityCovered=true; // ?
					}
				}
				separators.reserve(separators1.size());
				for(auto f:separators1)
					separators.push_back(f);
			}

			separatorIterator=0;
		}
		IntervalBound<typ> getLowerBound()
		{
			if(separatorIterator==0)
				return IntervalBound<typ>::minusInfinity();
			return separators[separatorIterator-1];
		}
		IntervalBound<typ> getUpperBound()
		{
			if(separatorIterator==separators.size())
				return IntervalBound<typ>::plusInfinity();
			return separators[separatorIterator];
		}
	  	Interval<typ> getInterval()
		{
	  		return Interval<typ>(getLowerBound(),getUpperBound());
		}
		bool isFirstInterval()
		{
			return separatorIterator==0;
		}
		bool isLastInterval()
		{
			return separatorIterator==separators.size();
		}
		void next()
		{
			assert(!isLastInterval());
			separatorIterator++;
		}
		bool isFiberEmpty()//for open interval
		{
			if(isLastInterval())
				return !plusInfinityCovered;
			if(isFirstInterval())
				return !minusInfinityCovered;
			return false;
		}
		Cone<typ> recessionOfFiberClosure()//for open interval
		{
			assert(!isFiberEmpty());
			return recessionOfNonTrivialFiber;
		}
	  gfan::Matrix<typ> recessionGenerators()
	  {
			assert(!isFiberEmpty());
			return generatorsOfRecession;
	  }
	};
	template<typename typ>
	class ProjectionsIterator{
	public:
		vector<ProjectionIterator<typ>> v;
		ProjectionsIterator(std::vector<HalfOpenCone<typ>> const &F, int parameterIndex)
		{
			v.reserve(F.size());
			for(auto &p:F)
				v.emplace_back(p,parameterIndex);
		}
	  	Interval<typ> getInterval()
		{
	  		Interval<typ> ret;
	  		for(auto &a:v)
	  			ret=ret.intersectedWith(a.getInterval());
	  		return ret;
		}
		bool isLastInterval()
		{
			for(auto &a:v)
				if(!a.isLastInterval())return false;
			return true;
		}
		void next()
		{
			IntervalBound<typ> u=getInterval().end;
			for(auto &a:v)
				if(!(u<a.getInterval().end))
					a.next();
		}
	};
	template<typename typ>
	void testFiber(std::vector<HalfOpenCone<typ> > const &collection)
	{
		//compare to  ~/gfan-svn/trunk/anton/RESULTS/n-4---2020-03-04/
	}

#if 0
	//This FaceLattice should be remove as a complete implementation would have too much code overlap with Complex
	template<typename typ>
	class FaceLattice{
		gfan::Matrix<typ> rays;
		using Face=vector<int>;
		vector<Face> faces;
		vector<bool> isFaceDeleted;
	public:
		FaceLattice(vector<HalfOpenCone<typ>> const &v)
		{
			gfan::Matrix<typ> linealitySpace;
			bool coneFound=false;
			for(int i=0;i<v.size();i++)
				if(!v[i].isEmpty())
				{
					linealitySpace=v[i].closure().getLinealitySpace();
					coneFound=true;
					break;
				}
			assert(coneFound);

			RayCollector<typ> collector(linealitySpace);

			int numberOfNonEmptyHalfOpenCones=0;
			for(int i=0;i<v.size();i++)
				if(!v[i].isEmpty())
				{
					numberOfNonEmptyHalfOpenCones++;
					auto rays=v.closure().getRays();
					Face face;
					for(int i=0;i<rays.getHeight();i++)
					{
						face.push_back(collector.lookup(normalize(rays[i].toVector())));
//					if(rays[i][0].isNegative())
//						toBeMadeEquivalent.insert(index);
					}
					faces.push_back(face);
					//NOW FIND ALL FACES
					isFaceDeleted.push_back(false);
				}
			rays=collector.getRays();
		}
/*		class NeighbourIterator
		{
			NeighbourIterator
		}*/
	};
#endif
	//
	HalfOpenCone<gfan::CircuitTableInt64> leakCone(gfan::Complex<gfan::CircuitTableInt64> &C, vector<int> const &Cindices, vector<int> const &Findices)
	{
		vector<int> commonBounded;
		for(auto &i:Findices)
			if(!C.vertices[i][0].isZero())
				commonBounded.push_back(i);

		vector<int> indicesCBounded;
		for(auto &i:Cindices)
			if(!C.vertices[i][0].isZero())
				indicesCBounded.push_back(i);

	    typedef gfan::CircuitTableInt64 typ; // type for integers
		gfan::Matrix<typ> lines=C.linealitySpace;
		for(int a=1;a<commonBounded.size();a++)
		  lines.appendRow(C.vertices[commonBounded[0]][0]*C.vertices[commonBounded[a]].toVector()-
				  C.vertices[commonBounded[a]][0]*C.vertices[commonBounded[0]].toVector());
if(1)
		for(auto &i:Findices)
			if(C.vertices[i][0].isZero())
				lines.appendRow(C.vertices[i]);


	//	std::cerr<<"111\n";
		gfan::Matrix<typ> rays(0,lines.getWidth());
		for(int i:Cindices)
		  if(C.vertices[i][0].isZero())rays.appendRow(C.vertices[i]);
	//	std::cerr<<"222\n";

		auto temp=setDifference(indicesCBounded,commonBounded);
		for(int a:temp)
		  rays.appendRow(-(C.vertices[commonBounded[0]][0]*C.vertices[a].toVector()-
				 C.vertices[a][0]*C.vertices[commonBounded[0]].toVector())); // The sign change is needed because the homogenising coordinate is negative
	//	std::cerr<<"333\n";


		lines.normalizeRows();
		lines.sortAndRemoveDuplicateRows();
		rays.normalizeRows();
		rays.sortAndRemoveDuplicateRows();

		GeneratedCone<typ> K(rays.transposed(),lines.transposed());
	//	std::cerr<<"AAAAAAAAA\n";
		auto Korth=K.getOrthogonalComplement();
		auto Knorm=K.getFacetNormals();
		Korth.normalizeRows();
		Korth.sortAndRemoveDuplicateRows();
		Knorm.normalizeRows();
		Knorm.sortAndRemoveDuplicateRows();

		HalfOpenCone<typ> KRelInterior(gfan::Matrix<typ>(K.getAmbientDimension(),0),Korth.transposed(),Knorm.transposed());//Put in empty matrix instead of the nonstrict inequalities
	//	std::cerr<<"BBBBBBBBBB\n";
		return KRelInterior;
	}
	vector<vector<int>> maximalSubsets(vector<vector<int>> const &C)
			{
		vector<vector<int>> ret;
				for(auto &v:C)
				{
					int nSetsBigger=0;
					for(auto &V:C)
						if(isSubsetOf(v,V))
						{
							nSetsBigger++;
							if(nSetsBigger>1)
								break;
						}
					if(nSetsBigger==1)
						ret.push_back(v);
				}
				return ret;
			}
  bool shouldBeDeleted(gfan::Complex<gfan::CircuitTableInt64> &C, gfan::Complex<gfan::CircuitTableInt64>::Cone const &p)
  {
    std::cerr << "shouldBeDeleted: p.dimension = " << p.dimension << std::endl;  
    typedef gfan::CircuitTableInt64 typ; // type for integers
    vector<HalfOpenCone<typ>> leakCones;

    { // INSERT RECESSION OF p
      gfan::Matrix<typ> unboundedDirections(0,C.ambientDimension);
      for(auto &i:p.indices)
	if(C.vertices[i][0].isZero()) // unbounded
	  unboundedDirections.appendRow(C.vertices[i]);
      GeneratedCone<typ> rC(unboundedDirections.transposed());
      auto rCorth=rC.getOrthogonalComplement();
      auto rCnorm=rC.getFacetNormals();
      rCorth.normalizeRows();
      rCorth.sortAndRemoveDuplicateRows();
      rCnorm.normalizeRows();
      rCnorm.sortAndRemoveDuplicateRows();
      auto rowOfOnes = gfan::Matrix<typ>::rowVectorMatrix(gfan::Vector<typ>::allOnes(rCnorm.getHeight()));
      auto sumOfRows = rowOfOnes * rCnorm;
      HalfOpenCone<typ> rCRelInterior(
				      rCnorm.transposed(), // non-strict inequalities
				      rCorth.transposed(), // non-strict equations
				      sumOfRows.transposed()  // strict inequalities
				      );//Put in empty matrix instead of the nonstrict inequalities
      leakCones.push_back(rCRelInterior);
      std::cerr << "shouldBeDeleted: rCRelInterior.getDimension() = " << rCRelInterior.getDimension() << std::endl;  
    } // end INSERT RECESSION OF p
	
    set<gfan::Complex<typ>::Cone const*> S;
    for(auto &i:p.indices)
      if(!C.vertices[i][0].isZero()) // bounded
	for(auto conePointer:C.rayIncidences[i])
	  S.insert(conePointer);

    std::cerr<<"shouldBeDeleted: Number of adjacent polyhedra (including self):"<<S.size()<<"\n";
    for(auto c:S)
    	if(!c->markedAsDeleted)
      {
	if(setDifference(c->indices,p.indices).size()==0)continue;
							
	// We find generators for the link
	auto indicesP=p.indices;
	auto indicesC=c->indices;
	vector<int> indicesPBounded;
	for(auto a:indicesP)
	  if(!C.vertices[a][0].isZero())
	    {
	      assert(C.vertices[a][0].isNegative());// is this correct?
	      indicesPBounded.push_back(a);
	    }
	vector<int> indicesCBounded;
	for(auto a:indicesC)
	  if(!C.vertices[a][0].isZero())indicesCBounded.push_back(a);
	auto commonBounded=intersection(indicesCBounded,indicesPBounded);
//	std::cerr<<"P:"<<toString2(indicesP)<<"\nPbounded:"<<toString2(indicesPBounded)<<"\n";
//	std::cerr<<"C:"<<toString2(indicesC)<<"\nCbounded:"<<toString2(indicesCBounded)<<"\n";
//	std::cerr<<"Intersection:"<<toString2(commonBounded)<<"\n";

	assert(commonBounded.size());

if(0)
	leakCones.push_back(leakCone(C,c->indices,intersection(c->indices,indicesP)));
else
{
	auto temp=maximalSubsets(C.facesContainedInSubset(intersection(c->indices,indicesP)));
	for(auto &a:temp)
		leakCones.push_back(leakCone(C,c->indices,a));
}



	// To improve the above: if intersection(c->indices,indicesP) was marked as deleted,
	// then we should consider facets of the intersection - these may still be present in the
	// precomplex, but they have smaller leak cones.
	// We need to find all maximal sets among those in the complex contained in the intersection.
      }
 //   std::cerr<<"Number of leak cones:\n"<<leakCones.size()<<"\n";
//    for(auto &a:leakCones)std::cerr<<a.toString();
    bool containsOrigin=sumOfHalfOpenConesContainsOrigin(leakCones);
//    std::cerr<<"CONTAINS ORIGIN:"<<containsOrigin<<"\n";

    //if(!containsOrigin)C.
    return !containsOrigin;
  }
  //  template<typename typ>
  void simplifyComplex(gfan::Complex<gfan::CircuitTableInt64> &C)
  {
    typedef gfan::CircuitTableInt64 typ; // type for integers
    bool somethingChanged=true;
    while(somethingChanged)
      {
	somethingChanged=false;
	for(auto &p:C.cones)
	  if(!p.markedAsDeleted)
	    {
	      //if(p.dimension==C.getMaxDim())
		{
		  if(shouldBeDeleted(C,p))
		    {
		      std::cerr<<C.fvector().toString();
		      p.markAsDeleted();
		      C.getMaxDim(true);
		      std::cerr<<"Deleting\n";//<<C.toString(FPF_rays|gfan::FanPrintingFlags::FPF_cones);
		      std::cerr<<C.fvector().toString();
		      somethingChanged=true;
		    }
		}
	    }
      }
  }
//  typedef gfan::CircuitTableInt64 typ; // type for integers
  typedef gfan::CircuitTableInteger typ; // type for integers
  void simplifyComplexSculpting(gfan::Complex<typ> &C)
  {
	  std::cerr<<"SCULPTING\n";
	  std::cerr<<"FVECTOR:"<<C.fvector().toString()<<"\n";
	  restart:
	  for(int i=0;i<C.vertices.getHeight();i++)
		  if(!C.vertices[i][0].isZero())
			  if(C.rayIncidences[i].size())
	  {
			  std::cerr<<"Checking vertex index:"<<i<<"\n";
			  gfan::Matrix<typ> rays(0,C.ambientDimension);
		  for(auto cone:C.rayIncidences[i])
			  for(auto j:cone->indices)
				  if(!C.vertices[j][0].isZero())
					rays.appendRow(C.vertices[j][0]*C.vertices[i].toVector()-C.vertices[i][0]*C.vertices[j].toVector());
				  else
					  rays.appendRow(C.vertices[j].toVector());
//		  std::cerr<<rays;
		  gfan::GeneratedCone<typ> K(rays.transposed());
//		  std::cerr<<"LINDIM"<<K.getDimensionOfLinealitySpace()<<"\n";
		  if(K.getDimensionOfLinealitySpace()==0)
		  {
			  // find separating hyperplane
			  // ask if vertices[i] is in convex hull of neighbours. Unbounded directions should also be included
			  gfan::Matrix<typ> neighbours(0,C.ambientDimension);
			  for(auto cone:C.rayIncidences[i])
				  for(auto j:cone->indices)
					  if(j!=i)
						  neighbours.appendRow(C.vertices[j].toVector());
			  neighbours.normalizeRows();
			  neighbours.sortAndRemoveDuplicateRows();

			  std::cerr<<C.vertices[i].toVector()<<"Has Neighbours:\n"<<neighbours.toString();
			  gfan::GeneratedCone<typ> convexHull(neighbours.transposed());

			  gfan::Vector<typ> separating(C.ambientDimension);
			  separating[0]=1;
			  separating[1]=2;
			  separating[2]=3;
			  assert(convexHull.getAmbientDimension()==separating.size());
			  std::cerr<<"IF STATEMENT\n";
			  std::cerr<<C.vertices[i].toVector()<<convexHull.getAmbientDimension()<<"\n";
			  if(convexHull.contains(C.vertices[i].toVector(),&separating))
			  {

				  assert(0);
			  }
			  else
			  {
				  separating=separating.normalized();
				  std::cerr<<"Separating normal:"<<separating.toString()<<"\n";
				  auto A=neighbours*gfan::Matrix<typ>::rowVectorMatrix(separating).transposed();
//				  std::cerr<<A;
				  assert(A.transposed()[0].toVector().isNonNegative());
				  auto a=dot(C.vertices[i].toVector(),separating);
//				  std::cerr<<a<<"\n";
				  assert(a.sign()<0);
//				  std::cerr<<"size"<<separating.size()<<"\n";
				  auto hyperplane=gfan::Cone<typ>::hyperPlane(separating);
				  auto halfspace=gfan::Cone<typ>::halfSpace(separating);
				  std::cerr<<"A1\n";
				  assert(!halfspace.contains(C.vertices[i].toVector()));
				  std::cerr<<"A2\n";
				  //halfspace=gfan::Cone<CircuitTableInt64>::halfSpace(separating);
				  //now we iterate over every relatively open cone involving i
				  //--- compute the intersection with the separating half space
				  //--- compute extreme vertices and rays
				  //--- remember this for all intersections
				  //--- remove the cone
				  //update ray/vertex list with new rays/vertices
				  //insert each earlier computed intersection polyhedron
				  //ALSO insert intersection with separating HYPERplane
				  vector<Complex<typ>::Cone*> oldCones;
				  vector<pair<gfan::Matrix<typ>,int>> newCones;
				  for(auto cone:C.rayIncidences[i])
				  {
					  std::cerr<<"lloop\n";
					  oldCones.push_back(const_cast<Complex<typ>::Cone*>(cone));

					  gfan::Matrix<typ> lines(0,C.ambientDimension);
					  // should add lineality space if any
					  gfan::Matrix<typ> rays(0,C.ambientDimension);


					  std::cerr<<"processing:"<<cone->toString()<<"\n";
					  for(auto j:cone->indices)
						  rays.appendRow(C.vertices[j].toVector());
//					  std::cerr<<"RAYS:\n"<<matrixToString(rays)<<"\n";
//					  std::cerr<<"LINES:\n"<<matrixToString(lines)<<"\n";
					  gfan::GeneratedCone<typ> p(rays.transposed(),lines.transposed());
					  std::cerr<<"Original dimension:"<<cone->dimension<<"\n";
					  std::cerr<<"p dimension:"<<p.getDimension()<<"\n";
					  assert(cone->dimension==p.getDimension());
					  //					  std::cerr<<"OKa\n";
//					  std::cerr<<matrixToString(p.getFacetNormals())<<matrixToString(p.getOrthogonalComplement());
					  auto pInHDescription=gfan::Cone<typ>(p.getFacetNormals().normalizedRows().transposed(),p.getOrthogonalComplement().normalizedRows().transposed());
//					  std::cerr<<"OKa2\n";
					  auto newPolyhedron=intersection(pInHDescription,halfspace);
					  auto newFace=intersection(pInHDescription,hyperplane);
					  auto newRays=newPolyhedron.getRays().normalizedRows();
					  auto newRays2=newFace.getRays().normalizedRows();
//					  std::cerr<<"NewRays1"<<matrixToString(newRays)<<"NewRays2"<<matrixToString(newRays2);

					  int dim1=newPolyhedron.getDimension();
					  int dim2=newFace.getDimension();
					  std::cerr<<"new dimension"<<dim1<<"\n";
					  std::cerr<<"new dimension"<<dim2<<"\n";
					  if(dim1>0)
					  {
						  newCones.push_back(std::make_pair(newRays,dim1));
						  std::cerr<<"ADDING intersection.\n";
					  }
					  if(dim2>0)
						  if(dim2<dim1)
						  {
							  newCones.push_back(std::make_pair(newRays2,dim2));
							  std::cerr<<"ADDING face of intersection.\n";
						  }
					  {
						  // check that we can rebuild cones with the right dimension
						  gfan::GeneratedCone<typ> p(newRays.transposed(),lines.transposed());
						  assert(p.getDimension()==dim1);
						  gfan::GeneratedCone<typ> p2(newRays2.transposed(),lines.transposed());
						  assert(p2.getDimension()==dim2);
					  }
				  }
				  C.substitute(oldCones, newCones);
				  std::cerr<<"FVECTOR:"<<C.fvector().toString()<<"\n";
				  std::cerr<<"Restarting\n";
				  goto restart;
			  }
		  }
	  }
  }
  void simplify(vector<gfan::HalfOpenCone<typ>> &halfOpenCones)
  {
//    typedef gfan::CircuitTableInt64 typ; // type for integers
    auto linealitySpace=halfOpenCones[0].closure().getLinealitySpace();
    RayCollector<typ> collector(linealitySpace);
    for(auto &c:halfOpenCones)
      {
	auto rays=c.closure().getRays();
	for(int i=0;i<rays.getHeight();i++)
	  collector.lookup(normalize(rays[i].toVector()));//duplicate code?
      }

    gfan::Complex<typ> C=extractFaceComplex(halfOpenCones,collector,linealitySpace);

    C.buildRayIncidences();

				      
    // Plan:
    // Make ConeContainer a vector instead of a set
    // Write a C.buildInverses method
    // For each polyhedron P do inverse look ups for each bounded ray
    // For each face of P take the correspondig union
    // Perform LP test for P that involves a cone for each face-union pair
    std::cerr<<"Simplifying\n"<<C.toString(FPF_rays|gfan::FanPrintingFlags::FPF_cones);
/*    if(0)
    	simplifyComplex(C);
    else*/
    	simplifyComplexSculpting(C);
    std::cerr<<"Result:\n"<<C.toString(FPF_rays|gfan::FanPrintingFlags::FPF_cones);


  }
  int main()
  {
    //	typedef gfan::CircuitTableInt128 typ; // type for integers
    	    typedef gfan::CircuitTableInt64 typ; // type for integers
	    //	    typedef gfan::CircuitTableInteger typ; // type for integers

	if(0)
	{
		vector<Cone<typ>> L;
		stringstream A("{(1,0),(0,1)}");
		stringstream B("{(0,0),(0,0)}");
		stringstream C("{(1,-1),(1,0)}");
//		stringstream C("{(0,-1),(1,0)}");
		stringstream DD("{(0,0),(0,0)}");
//		std::cerr<<"A"<<gfan::Matrix<typ>::readMatrix(A,2)<<"\n";
//		std::cerr<<"B"<<gfan::Matrix<typ>::readMatrix(B,2)<<"\n";
		L.emplace_back(gfan::Matrix<typ>::readMatrix(A,2),gfan::Matrix<typ>::readMatrix(B,2));
		L.emplace_back(gfan::Matrix<typ>::readMatrix(C,2),gfan::Matrix<typ>::readMatrix(DD,2));
		for(auto &c:L)
		{
			std::cerr<<c.toString();
		}
		std::cerr<<nonEmptyIntersectionOfDualsInteriors(L);
		exit(0);
	}



    FileParser S(Stdin); // allows to read and parse polynomials, etc.
    const PolynomialSet g = S.parsePolynomialSetWithRing(); // system of polynomials
    std::cerr<<"DONE reading input\n";

    /*		auto g=StringParser(
//				  "Q[t,x,y]{x+y+t,xxxtt+yyytt+1tt+xy}"
				"Q[x01,x02,x03,x04,x05,x06,x07,x08,x09,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30]"
"{1+x02*x11^2-x03*x11^2,"
"1+x02*x12^2-x04*x12^2,"
"1+x02*x13^2-x05*x13^2,"
"1+x02*x14^2-x06*x14^2,"
"1+x03*x15^2-x04*x15^2,"
"1+x03*x16^2-x05*x16^2,"
"1+x03*x17^2-x06*x17^2,"
"1+x04*x18^2-x05*x18^2,"
"1+x04*x19^2-x06*x19^2,"
"1+x05*x20^2-x06*x20^2,"
"1+x21^2-x07*x21^2,"
"1+x22^2-x08*x22^2,"
"1+x23^2-x09*x23^2,"
"1+x24^2-x10*x24^2,"
"1+x07*x25^2-x08*x25^2,"
"1+x07*x26^2-x09*x26^2,"
"1+x07*x27^2-x10*x27^2,"
"1+x08*x28^2-x09*x28^2,"
"1+x08*x29^2-x10*x29^2,"
"1+x09*x30^2-x10*x30^2,"
"x02-x01^3*x11*x21^3-x01^5*x12*x22^3-x01^7*x13*x23^3-x01^11*x14*x24^3,"
"x03+x01^2*x11*x21^3-x01^5*x15*x25^3-x01^7*x16*x26^3-x01^11*x17*x27^3,"
"x04+x01^2*x12*x22^3+x01^3*x15*x25^3-x01^7*x18*x28^3-x01^11*x19*x29^3,"
"x05+x01^2*x13*x23^3+x01^3*x16*x26^3+x01^5*x18*x28^3-x01^11*x20*x30^3,"
"x06+x01^2*x14*x24^3+x01^3*x17*x27^3+x01^5*x19*x29^3+x01^7*x20*x30^3,"
"1-x01^3*x11^3*x21-x01^5*x12^3*x22-x01^7*x13^3*x23-x01^11*x14^3*x24,"
"x07+x01^2*x11^3*x21-x01^5*x15^3*x25-x01^7*x16^3*x26-x01^11*x17^3*x27,"
"x08+x01^2*x12^3*x22+x01^3*x15^3*x25-x01^7*x18^3*x28-x01^11*x19^3*x29,"
"x09+x01^2*x13^3*x23+x01^3*x16^3*x26+x01^5*x18^3*x28-x01^11*x20^3*x30,"
"x10+x01^2*x14^3*x24+x01^3*x17^3*x27+x01^5*x19^3*x29+x01^7*x20^3*x30,"
"x01^2*x02+x01^3*x03+x01^5*x04+x01^7*x05+x01^11*x06,"
"x01^2+x01^3*x07+x01^5*x08+x01^7*x09+x01^11*x10,"
"x01^2*x02+x01^3*x03*x07+x01^5*x04*x08-x01^5*x11*x21+x01^7*x05*x09-x01^7*x12*x22-x01^8*x15*x25-x01^9*x13*x23-x01^10*x16*x26+x01^11*x06*x10-x01^12*x18*x28-x01^13*x14*x24-x01^14*x17*x27-x01^16*x19*x29-x01^18*x20*x30}"
				  "Q[t,x,y]{x+y+t,xxxtt+yyytt+1tt+xy}"
			  "Q[t,x,y]{x+y+t}"
	  ).parsePolynomialSetWithRing();
    */	  
          vector<gfan::Matrix<typ> > configurations;
	  for(auto &p:g)
		  configurations.push_back(convertMatrix<typ>(rowsToIntegerMatrix2(p.exponents())));

	  vector<HalfOpenCone<typ> > f;

	  //We pretend that first coordinate is special (that it is the t parameter)
	  int n=configurations[0].getWidth();
	  int nonEmptyIntersections=0;
//	  vector<int> edgeCounts;
//	  for(int i=0;i<configurations.size();i++)edgeCounts.push_back(numberOfEdges(configurations[i]));
	  vector<int> used(configurations.size());
	  PolytopeIntersectionData<typ> data(configurations);

	  gfan::Matrix<typ> strictInequalities(n,1);
	  strictInequalities[0][0]=typ(-1);

	  if(optionHalfspace.getValue())
	  {
		  auto temp=-gfan::Vector<typ>::allOnes(n);
		  temp[0]=typ(-optionHalfspace.getValue());
		  strictInequalities=strictInequalities.transposed();
		  strictInequalities.appendRow(temp);
		  strictInequalities=strictInequalities.transposed();
		  std::cerr<<"STRICTINEQUALITIES:"<<strictInequalities.toString()<<"\n";
	  }

	  ProgressCounter progress;
	  CommonStatistics statistics;
	  
	  if(loadFrom.getValue())
	  {
		  string loadName(loadFrom.getValue());
		  std::cerr<<"Loading vector of HalfOpenCones from file \""<<loadName<<"\"\n";
		  std::ifstream F;
		  F.open(loadName.c_str(),std::ifstream::in);
		  f=loadHalfOpenConeVector<typ>(F);
		  F.close();
	  }
	  else
	  {
	  std::cerr<<"ABOUT TO INTERSECT\n";
	  commonRefinement(
			  gfan::HalfOpenCone<typ>(
					  gfan::Matrix<typ>(n,0),
					  gfan::Matrix<typ>(n,0),
					  strictInequalities),
			  used,configurations,f,nonEmptyIntersections/*,&edgeCounts*/,data,RelationTable(data.layout),progress,&statistics,optionNThreads.getValue()/*32*//*8*/);
	  }
	  if(saveAs.getValue())
	  {
		  string saveName(saveAs.getValue());
		  std::cerr<<"Saving vector of HalfOpenCones in file \""<<saveName<<"\"\n";
		  std::ofstream F;
		  F.open(saveName.c_str(),std::ofstream::out);
		  save(f,F);
		  F.close();
	  }

	  assert(!f.empty());

	  std::cerr<<"INTERSECTION COMPLETED\n";
	  std::cerr<<"NUMBER OF HALFOPEN CONES:"<<f.size()<<"\n";

	  // Convert
	  vector<HalfOpenCone<CircuitTableInteger> > F;
//	  vector<HalfOpenCone<CircuitTableInt64> > F;
	  F.reserve(f.size());
	  for(auto &c:f)F.emplace_back(c);

	  if(0)
	    {
	      simplify(F);
	      return 0;
	    }

	  std::cerr<<"NUMBER OF HALFOPEN CONES:"<<F.size()<<"\n";

	  if(0)
	  {
		  for(auto &c:F)
		  {
			  auto w=c.getRelativeInteriorPoint();
			  IntegerVector W(w.size());
			  for(int i=0;i<W.size();i++)
				  W[i]=w[i].toInt64();

			  auto inF=initialForms(g,W);
			  for(auto &f:inF)
				  assert(f.numberOfTerms()>1);
		  }
	  }



	  if(fibersOption.getValue())
	  {
		  std::cout<<"Printing intervals:\n";
		  ProjectionsIterator<typ> i(f,g.getRing().getNumberOfVariables()-1);
		  int II=0;
		  bool first=true;
		  do
		  {
		    std::cerr<<"Cone number "<<II++<<"\n";
			  if(!first)i.next();
		    i.getInterval();
			  std::cerr<<i.getInterval().toString()<<"\n";

			  int count=std::count_if(i.v.begin(),i.v.end(),[](auto &a){return !a.isFiberEmpty();});


			  int count2=0;
			  {
			    Fraction<typ> fraction=i.getInterval().getInteriorPoint();
			    HalfOpenCone<typ> H=HalfOpenCone<typ>::hyperplane(fraction.numerator*Vector<typ>::standardVector(n,0)+fraction.denominator*Vector<typ>::standardVector(n,g.getRing().getNumberOfVariables()-1));
				std::cerr<<"LOOP. fraction:"<<fraction.toString()<<"\n";
			    for(auto &iter:i.v)
			      {
				auto A=intersection2(iter.coneOverThis,H);
				count2+=!A.isEmpty();
			      }
			  }

			  if(0)
			    {
			      vector<Cone<typ>> recCones;
			      recCones.reserve(count);
			      for(auto &a:i.v)
				if(!a.isFiberEmpty())
				  recCones.push_back(a.recessionOfFiberClosure());

			      std::cerr<<"Count"<<count<<"\n";
			      std::cerr<<"Count2"<<count2<<"\n";
			      if(count<100)
				std::cerr<<"nonEmpty:"<<nonEmptyIntersectionOfDualsInteriors(recCones)<<"\n";
			      else
				std::cerr<<"Skipping.\n";

			    }
			  else if(0)
			    {
			      int N=g.getRing().getNumberOfVariables();
			      gfan::Matrix<gfan::CircuitTableInteger> strict(0,N);
			      for(auto &a:i.v)
				if(!a.isFiberEmpty())
				    strict.append(gfan::Matrix<gfan::CircuitTableInteger>(a.recessionGenerators()));
			      std::cerr<<"rows:"<<strict.getHeight()<<"\n";
			      HalfOpenCone<gfan::CircuitTableInteger> certificatesForPointedness(gfan::Matrix<gfan::CircuitTableInteger>(N,0),gfan::Matrix<gfan::CircuitTableInteger>(N,0),strict.transposed());
			      if(certificatesForPointedness.isEmpty())
				{
				  std::cerr<<"not pointed\n";
				}
			      else
				{
				  std::cerr<<"pointed\n";
				}

			    }

			  else
			    {
			      std::cerr<<"Processing slice\n";

			      {
				      vector<HalfOpenCone<typ>> halfOpenCones;
				      Fraction<typ> fraction=i.getInterval().getInteriorPoint();
				      HalfOpenCone<typ> H=HalfOpenCone<typ>::hyperplane(fraction.numerator*Vector<typ>::standardVector(n,0)+fraction.denominator*Vector<typ>::standardVector(n,g.getRing().getNumberOfVariables()-1));
				      for(auto &c:i.v)
				      {
				    	  auto A=intersection2(c.coneOverThis,H);
				    	  if(!A.isEmpty())
				    		  halfOpenCones.push_back(A);
				      }
				      if(halfOpenCones.size())
				      {
//					simplify(halfOpenCones);
				    	  assert(0/*Since simplify has been instantiated with the wrong type, the call above cannot happen without converting */);
				      }
				      else
				    	  std::cerr<<"Slice is empty\n";
			      }



			      vector<pair<Cone<typ>,int>> V;
			      int LO=1000;
			      int HI=-1;

			      vector<HalfOpenCone<typ>> halfOpenCones;
			    Fraction<typ> fraction=i.getInterval().getInteriorPoint();
			    HalfOpenCone<typ> H=HalfOpenCone<typ>::hyperplane(fraction.numerator*Vector<typ>::standardVector(n,0)+fraction.denominator*Vector<typ>::standardVector(n,g.getRing().getNumberOfVariables()-1));



			    for(auto &c:i.v)
				{
				  {
				    auto A=intersection2(c.coneOverThis,H);
				    if(!A.isEmpty())
				      {
					V.emplace_back(make_pair<Cone<typ>,int>(A.closure(),A.getDimension()));
					halfOpenCones.push_back(A);
					LO=min(LO,A.getDimension());
					HI=max(HI,A.getDimension());
				      }
				  }
				}
			      std::cerr<<"LODIM:"<<LO<<" HIDIM:"<<HI<<" #halfopencones in slice:"<<halfOpenCones.size()<<"\n";
			      std::cerr<<"FVECTOR:"<<fvector(halfOpenCones).toString()<<"\n";
			      for(auto &c:V)
				if(c.second==HI)
				{
				  int adjacentCones=0;
				  auto CClosure=c.first;
				  //				  std::cerr<<"Processing cone";


				  vector<HalfOpenCone<typ>> leakCones;
				  auto H2=Cone<typ>::hyperPlane(Vector<typ>::standardVector(n,0));
				  {
				    // c.first is the closure of a halfopen kegle over snit a polyeder. We need the recession. Intersect with infinity
				    auto C=intersection(c.first,H2);
				    auto T=C.getFacetNormals().transposed();
				    auto O=C.getOrthogonalComplement().transposed();
				    T.normalizeRows();
				    T.sortAndRemoveDuplicateRows();
				    O.normalizeRows();
				    O.sortAndRemoveDuplicateRows();
				    gfan::Vector<typ> v(T.getWidth());
				    for(int i=0;i<T.getHeight();i++)
				      v+=T[i].toVector();
				    leakCones.emplace_back(T.transposed(),O.transposed(),gfan::Matrix<typ>::rowVectorMatrix(v).transposed());
				    //need to check that the above code is correct.
				  }
				  for(auto &k:V)
				    {
				      auto KClosure=k.first;
				      auto temp=intersection(CClosure,KClosure);

				      if(temp.getDimension()!=temp.getDimensionOfLinealitySpace())
				      {
					adjacentCones++;

					auto w=temp.getRelativeInteriorPoint();
					auto t2=KClosure.link(w);
					auto l=intersection(H2,t2);
					
					// For relatively open cone:
					// leakCones.emplace_back(gfan::Matrix<typ>(l.getAmbientDimension(),0),l.getOrthogonalComplement(),l.getFacetNormals());
					// For closed cone missing its lineality space:
					auto T=l.getFacetNormals().transposed();
					T.normalizeRows();
					T.sortAndRemoveDuplicateRows();
					auto O=l.getOrthogonalComplement().transposed();
					O.normalizeRows();
					O.sortAndRemoveDuplicateRows();
					gfan::Vector<typ> v(T.getWidth());
					for(int i=0;i<T.getHeight();i++)
					  v+=T[i].toVector();
					leakCones.emplace_back(T.transposed(),O.transposed(),gfan::Matrix<typ>::rowVectorMatrix(v).transposed());
				      }
				    }
				  std::cerr<<"#adjacent cones "<<adjacentCones<<"\n";
				  std::cerr<<"CONTAINS ORIGIN:"<<sumOfHalfOpenConesContainsOrigin(leakCones)<<"\n";
				}
			    }
			  
			  first=false;
		  }
		  while(!i.isLastInterval());

	  }


      vector<HalfOpenCone<CircuitTableInteger> > f2;
      vector<gfan::Matrix<CircuitTableInteger> > f2Rays;
      f2.reserve(f.size());
      f2Rays.reserve(f.size());
      for(auto &a:f)
    	  if(!a.isEmpty())
    	  {
    		  auto temp=HalfOpenCone<CircuitTableInteger>(a);
    		  f2Rays.push_back(temp.closure().getRays());
    		  f2.push_back(temp);
    	  }
		  else
		  {
			  std::cerr<<"Empty set computed\n"<<a.toString()<<"\nLifted:\n"<<a.lifted.toString()<<"\n";
		  }

      std::cerr<<"DONE CONVERTING FROM INT64 AND COMPUTING RAYS\n";

	  auto linealitySpace=f2.begin()->closure().getLinealitySpace();
	  RayCollector<CircuitTableInteger> collector(linealitySpace);

	  int numberOfNonEmptyHalfOpenCones=0;
	  EquivalenceRelation rel;
	  for(int i=0;i<f2.size();i++)
		  {
			  set<int> toBeMadeEquivalent;
			  if(!(numberOfNonEmptyHalfOpenCones&255))std::cerr<<numberOfNonEmptyHalfOpenCones<<"\n";
			  numberOfNonEmptyHalfOpenCones++;
			  std::cerr<<"Getting rays\n";
			  auto rays=f2Rays[i];
			  std::cerr<<"Preparing"<<rays.getHeight()<<"\n";
			  for(int i=0;i<rays.getHeight();i++)
			  {
				  int index=collector.lookup(normalize(rays[i].toVector()));
				  if(rays[i][0].isNegative())
					  toBeMadeEquivalent.insert(index);
			  }
			  std::cerr<<"Making equivalent\n";
			  if(!toBeMadeEquivalent.empty())
				  for(auto i=toBeMadeEquivalent.begin();i!=toBeMadeEquivalent.end();i++)
					  rel.makeEquivalent(*i,*toBeMadeEquivalent.begin());
		  }
	  std::cerr<<"DONE SETTING UP EQUIVALENCES\n";
	  int N_RAYS=collector.getRays().getHeight();
	  rel.grow(collector.getRays().getHeight());

	  std::cerr<<"N_RAYS"<<N_RAYS<<"\n";

	  // First we build a vector for each class - later we must delete some entries because each half line is its own class

	  int nClasses=rel.buildClassNames();
	  vector<pair<vector<HalfOpenCone<CircuitTableInteger> >,vector<gfan::Matrix<CircuitTableInteger> > > > components(nClasses);

	  for(int I=0;I<f2.size();I++)
		  {
			  auto rays=f2Rays[I];
//			  std::cerr<<rays.toString()<<"\n";
			  int i=0;
			  for(i=0;i<rays.getHeight();i++)
				  if(rays[i][0].isNegative())
					  break;
			  if(i==rays.getHeight())
				  std::cerr<<"THIS SHOULD NOT HAPPEN\n";
			  else
			  {
				  int index=rel.getClassName(rel.getRepresentative(collector.lookup(normalize(rays[i].toVector()))));
				  components[index].first.push_back(f2[I]);
				  components[index].second.push_back(f2Rays[I]);
			  }
		  }

	  std::cerr<<"COMPONENTSVECTOR SIZE:"<<components.size()<<"\n";
	  int numberOfComponents=0;
	  for(auto &F:components)
		  if(F.first.size())
		  {
			  numberOfComponents++;
		  }
	  std::cerr<<"NUMBER OF COMPONENTS:"<<numberOfComponents<<"\n";

	  int componentIndex=0;
	  bool allGood=true;
	  for(auto &F:components)
		  if(F.first.size())
		  {
			  std::cerr<<"Component"<<componentIndex<<"\n";
			  RayCollector<CircuitTableInteger> collector2(linealitySpace);
			  for(int i=0;i<F.first.size();i++)
				  {
					  auto rays=F.second[i];
					  for(int i=0;i<rays.getHeight();i++)
					  {
						  int index=collector2.lookup(normalize(rays[i].toVector()));
					  }
				  }
			  auto componentRays=collector2.getRays();
			  gfan::Matrix<CircuitTableInteger> actualRays(0,componentRays.getWidth());
			  for(int i=0;i<componentRays.getHeight();i++)
				  if(componentRays[i][0].isZero())
					  actualRays.appendRow(componentRays[i]);

			  std::cerr<<"actualRays:\n"<<actualRays<<"\n";
			  std::cerr<<"linealitySpace:\n"<<linealitySpace<<"\n";
	    for(int i=1;i<=optionNumberOfDimensionsToClear.getValue();i++)
	      {
		actualRays.setSubMatrix(0,i,actualRays.getHeight(),i+1,actualRays);
		linealitySpace.setSubMatrix(0,i,linealitySpace.getHeight(),i+1,linealitySpace);
		std::cerr<<"actualRays:\n"<<actualRays<<"\n";
		std::cerr<<"linealitySpace:\n"<<linealitySpace<<"\n";
	      }
			  gfan::GeneratedCone<CircuitTableInteger> temp(actualRays.submatrix(0,1,actualRays.getHeight(),actualRays.getWidth()).transposed(),linealitySpace.submatrix(0,1,linealitySpace.getHeight(),linealitySpace.getWidth()).transposed());
			  std::cerr<<"Dimension of Lineality space of recession:"<<temp.getDimensionOfLinealitySpace()<<"\n";
			  if(temp.getDimensionOfLinealitySpace()==0)
			    std::cerr<<"SUCCESS\n";
			  else
			    {
			      std::cerr<<"FAILURE\n";
			      allGood=false;
			    }

			  componentIndex++;
		  }

	  std::cerr<<"AllGood:"<<allGood<<"\n";

	  {
	    auto allRays=collector.getRays();
	    gfan::Matrix<CircuitTableInteger> actualRays(0,allRays.getWidth());
	    for(int i=0;i<allRays.getHeight();i++)
	      if(allRays[i][0].isZero())
		actualRays.appendRow(allRays[i]);

	    std::cerr<<"actualRays:\n"<<actualRays<<"\n";
	    std::cerr<<"linealitySpace:\n"<<linealitySpace<<"\n";
	    for(int i=1;i<=optionNumberOfDimensionsToClear.getValue();i++)
	      {
		actualRays.setSubMatrix(0,i,actualRays.getHeight(),i+1,actualRays);
		linealitySpace.setSubMatrix(0,i,linealitySpace.getHeight(),i+1,linealitySpace);
		std::cerr<<"actualRays:\n"<<actualRays<<"\n";
		std::cerr<<"linealitySpace:\n"<<linealitySpace<<"\n";
	      }
	    gfan::GeneratedCone<CircuitTableInteger> temp(actualRays.submatrix(0,1,actualRays.getHeight(),actualRays.getWidth()).transposed(),linealitySpace.submatrix(0,1,linealitySpace.getHeight(),linealitySpace.getWidth()).transposed());
	    std::cerr<<"Dimension of Lineality space of recession:"<<temp.getDimensionOfLinealitySpace()<<"\n";
	    if(temp.getDimensionOfLinealitySpace()==0)
	      std::cerr<<"Global recession is pointed:YES\n";
	    else
	      {
		std::cerr<<"Global recession is pointed:NO\n";
		allGood=false;
	      }
	  }

	  vector<Vector<int>> fvectors;

	  
	  componentIndex=0;
	  for(auto &F:components)
		  if(F.first.size())
		  {
			  auto FV=fvector(F.first);
			  std::cerr<<componentIndex<<" "<<FV.toString()<<"\n";
			  fvectors.push_back(FV);
			  componentIndex++;
		  }

	  auto fv=fvector(f2);
	  std::cerr<<fv.toString()<<"\n";

	  for(auto &v:fvectors)
		  fv-=v;

	  std::cerr<<"DIFFERENCE:\n"<<fv.toString()<<"\n";

#if 0

//				std::cerr<<"THE RAYS\n"<<collector.getRays()<<"\n";
	  std::cerr<<"N_RAYS\n"<<collector.getRays().getHeight()<<"\n";
/*				for(auto &c:f)
		{
			std::cerr<<"HALFOPENCONE:\n";
			c.extractFaceComplex();
		}*/
///*				std::cerr<<*/ extractFaceComplex(f,collector,linealitySpace);

	  std::cerr<<fvector(f).toString()<<"\n";
#endif

	  std::cerr<<"Number of computed half open cones:"<<f2.size()<<"\n";
	  std::cerr<<"Number of non-empty half open cones:"<<numberOfNonEmptyHalfOpenCones<<"\n";
	  std::cerr<<"NonEmptyIntersections:"<<statistics.numberOfIntermediateVertices<<"\n";

	  return 0;
  }
};

static TropicalPrevarietyComponentsApplication theApplication;