1 #include <sstream>
2 #include "polyhedralfan.h"
3 #include "reversesearch.h"
4 #include "wallideal.h"
5 #include "buchberger.h"
6 #include "printer.h"
7 #include "timer.h"
8 #include "symmetry.h"
9 #include "polymakefile.h"
10 #include "symmetriccomplex.h"
11 #include "linalg.h"
12 #include "lp.h"
13 #include "codimoneconnectedness.h"
14 #include "symmetrictraversal.h"
15 #include "traverser_groebnerfan.h"
16 #include "tropical2.h"
17 #include "log.h"
18 
19 static Timer polyhedralFanRefinementTimer("Polyhedral fan refinement",1);
20 
PolyhedralFan(int ambientDimension)21 PolyhedralFan::PolyhedralFan(int ambientDimension):
22   n(ambientDimension)
23 {
24 }
25 
26 
fullSpace(int n)27 PolyhedralFan PolyhedralFan::fullSpace(int n)
28 {
29   PolyhedralFan ret(n);
30 
31   PolyhedralCone temp(n);
32   temp.canonicalize();
33   ret.cones.insert(temp);
34 
35   return ret;
36 }
37 
38 
halfSpace(int n,int i)39 PolyhedralFan PolyhedralFan::halfSpace(int n, int i)
40 {
41   assert(0<=i);
42   assert(i<n);
43   PolyhedralFan ret(n);
44 
45   IntegerVector v(n);
46   v[i]=-1;
47   IntegerVectorList l;
48   IntegerVectorList empty;
49   l.push_back(v);
50   ret.cones.insert(PolyhedralCone(l,empty,n));
51 
52   return ret;
53 }
54 
55 
facetsOfCone(PolyhedralCone const & c)56 PolyhedralFan PolyhedralFan::facetsOfCone(PolyhedralCone const &c)
57 {
58   PolyhedralCone C(c);
59   C.canonicalize();
60   PolyhedralFan ret(C.ambientDimension());
61 
62   IntegerVectorList halfSpaces=C.getHalfSpaces();
63 
64   for(IntegerVectorList::const_iterator i=halfSpaces.begin();i!=halfSpaces.end();i++)
65     {
66       IntegerVectorList a;
67       IntegerVectorList b;
68       b.push_back(*i);
69       PolyhedralCone n=intersection(PolyhedralCone(a,b),c);
70       n.canonicalize();
71       ret.cones.insert(n);
72     }
73   return ret;
74 }
75 
complementOfCone(PolyhedralCone const & c,bool includec)76 PolyhedralFan PolyhedralFan::complementOfCone(PolyhedralCone const &c, bool includec)
77 {
78   PolyhedralCone C=c;
79   C.canonicalize();
80   IntegerVectorList inequalities=C.getHalfSpaces();
81   IntegerVectorList equations=C.getEquations();
82 
83   for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
84     inequalities.push_back(*i);
85 
86   int n=C.ambientDimension();
87 
88   PolyhedralFan ret=PolyhedralFan::fullSpace(n);
89   for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
90     {
91       PolyhedralFan temp(C.ambientDimension());
92       for(int j=-1;j<2;j+=2)
93 	{
94 	  IntegerVectorList inequality;
95 	  inequality.push_back(j* *i);
96 	  IntegerVectorList empty;
97 	  PolyhedralCone tempC(inequality,empty,n);
98 	  tempC.canonicalize();
99 	  temp.insert(tempC);
100 	}
101       ret=refinement(ret,temp);
102     }
103   if(!includec)ret.remove(C);
104   return ret;
105 }
106 
bergmanOfPrincipalIdeal(Polynomial const & p1)107 PolyhedralFan PolyhedralFan::bergmanOfPrincipalIdeal(Polynomial const &p1)
108 {
109   PolynomialRing theRing=p1.getRing();
110   PolynomialRing theSecondRing=theRing.withVariablesAppended("H");
111   Polynomial p=p1.homogenization(theSecondRing);
112   PolyhedralFan ret(p1.getNumberOfVariables());
113   PolynomialSet g(theSecondRing);
114   g.push_back(p);
115   buchberger(&g,LexicographicTermOrder());
116 
117   EnumerationTargetCollector gfan;
118   LexicographicTermOrder myTermOrder;
119   ReverseSearch rs(myTermOrder);
120 
121   rs.setEnumerationTarget(&gfan);
122 
123   rs.enumerate(g);
124 
125   PolynomialSetList theList=gfan.getList();
126   for(PolynomialSetList::const_iterator i=theList.begin();i!=theList.end();i++)
127     {
128       //      AsciiPrinter(Stderr).printPolynomialSet(*i);
129       IntegerVectorList inequalities=wallInequalities(*i);
130       IntegerVectorList f=wallFlipableNormals(*i,true);
131       for(IntegerVectorList::const_iterator j=f.begin();j!=f.end();j++)
132 	{
133 	  //	  AsciiPrinter(Stderr).printVector(*j);
134 	  if(myTermOrder(*j,*j-*j))
135 	    {
136 	      IntegerVectorList equalities;
137 	      equalities.push_back(*j);
138 	      PolyhedralCone c=PolyhedralCone(inequalities,equalities).withLastCoordinateRemoved();
139 	      c.canonicalize();
140 	      c.setLinearForm(i->begin()->getMarked().m.exponent.subvector(0,p1.getNumberOfVariables()));
141 //	      c.setMultiplicity(gcdOfVector(j->subvector(0,j->size()-1)));
142 
143 	      Polynomial temp=initialForm(p1,c.getRelativeInteriorPoint());
144 	      Polynomial temp1=initialForm(temp,j->subvector(0,j->size()-1));
145 	      Polynomial temp2=initialForm(temp,-j->subvector(0,j->size()-1));
146 	      temp1.mark(LexicographicTermOrder());
147 	      temp2.mark(LexicographicTermOrder());
148 	      c.setMultiplicity(gcdOfVector(temp1.getMarked().m.exponent-temp2.getMarked().m.exponent));
149 
150 	      ret.cones.insert(c);
151 	    }
152 	}
153     }
154 
155   return ret;
156 }
157 
158 
normalFanOfNewtonPolytope(Polynomial const & p1)159 PolyhedralFan PolyhedralFan::normalFanOfNewtonPolytope(Polynomial const &p1)
160 {
161   PolynomialRing theRing=p1.getRing();
162   PolynomialRing theSecondRing=theRing.withVariablesAppended("H");
163   Polynomial p=p1.homogenization(theSecondRing);
164   PolyhedralFan ret(p1.getNumberOfVariables());
165   PolynomialSet g(theSecondRing);
166   //  PolynomialSet g(theRing);//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167   g.push_back(p);
168   buchberger(&g,LexicographicTermOrder());
169 
170   EnumerationTargetCollector gfan;
171 /*  {//old enumeration strategy
172   LexicographicTermOrder myTermOrder;
173   ReverseSearch rs(myTermOrder);
174 
175   rs.setEnumerationTarget(&gfan);
176 
177   rs.enumerate(g);
178   }
179 */
180   {//new enumeration strategy
181     GroebnerFanTraverser traverser(g);
182     traverser.setIsKnownToBeComplete(true);
183     TargetGlue target(gfan);
184     symmetricTraverse(traverser,target);
185   }
186 
187   PolynomialSetList theList=gfan.getList();
188   for(PolynomialSetList::const_iterator i=theList.begin();i!=theList.end();i++)
189     {
190       IntegerVectorList inequalities=wallInequalities(*i);
191       IntegerVectorList equalities;
192 
193       PolyhedralCone c=PolyhedralCone(inequalities,equalities).withLastCoordinateRemoved();
194       c.canonicalize();
195       c.setLinearForm(i->begin()->getMarked().m.exponent.subvector(0,c.ambientDimension()));
196       ret.cones.insert(c);
197     }
198 
199   return ret;
200 }
201 
202 
print(class Printer * p) const203 void PolyhedralFan::print(class Printer *p)const
204 {
205   p->printString("Printing PolyhedralFan");
206   p->printNewLine();
207   p->printString("Ambient dimension: ");
208   p->printInteger(n);
209   p->printNewLine();
210   p->printString("Number of cones: ");
211   p->printInteger(cones.size());
212   p->printNewLine();
213   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
214     {
215       p->printNewLine();
216       p->printPolyhedralCone(*i);
217     }
218   p->printString("Done printing PolyhedralFan.");
219   p->printNewLine();
220 }
221 
getAmbientDimension() const222 int PolyhedralFan::getAmbientDimension()const
223 {
224   return n;
225 }
226 
isEmpty() const227 bool PolyhedralFan::isEmpty()const
228 {
229   return cones.empty();
230 }
231 
getMaxDimension() const232 int PolyhedralFan::getMaxDimension()const
233 {
234   assert(!cones.empty());
235 
236   return cones.begin()->dimension();
237 }
238 
getMinDimension() const239 int PolyhedralFan::getMinDimension()const
240 {
241   assert(!cones.empty());
242 
243   return cones.rbegin()->dimension();
244 }
245 
refinement(const PolyhedralFan & a,const PolyhedralFan & b,int cutOffDimension,bool allowASingleConeOfCutOffDimension,bool stable)246 PolyhedralFan refinement(const PolyhedralFan &a, const PolyhedralFan &b, int cutOffDimension, bool allowASingleConeOfCutOffDimension, bool stable)
247 {
248   //Josephine Yu contributed to this function.
249 
250   TimerScope ts(&polyhedralFanRefinementTimer);
251   //  fprintf(Stderr,"PolyhedralFan refinement: #A=%i #B=%i\n",a.cones.size(),b.cones.size());
252   int conesSkipped=0;
253   int numberOfComputedCones=0;
254   bool linealityConeFound=!allowASingleConeOfCutOffDimension;
255   assert(a.getAmbientDimension()==b.getAmbientDimension());
256 
257   PolyhedralFan ret(a.getAmbientDimension());
258 
259   for(PolyhedralConeList::const_iterator A=a.cones.begin();A!=a.cones.end();A++)
260     for(PolyhedralConeList::const_iterator B=b.cones.begin();B!=b.cones.end();B++)
261       {
262 	PolyhedralCone c=intersection(*A,*B);
263 	int cdim=c.dimension();
264 	//	assert(cdim>=linealitySpaceDimension);
265 	bool thisIsLinealityCone=(cutOffDimension>=cdim);
266 //        if(((!thisIsLinealityCone)||(!linealityConeFound))&&((!stable)||(sum(*A,*B).dimension()==c.ambientDimension())))
267 /*          if(((!thisIsLinealityCone)||(!linealityConeFound))&&((!stable)||
268               (rank(combineOnTop(rowsToIntegerMatrix(A->generatorsOfSpan(),c.ambientDimension()),
269                             rowsToIntegerMatrix(B->generatorsOfSpan(),c.ambientDimension())))==c.ambientDimension())))
270 */            if(((!thisIsLinealityCone)||(!linealityConeFound))&&((!stable)||
271                     (A->dimension()+B->dimension()+::rank(combineOnTop(rowsToIntegerMatrix(A->getImpliedEquations(),c.ambientDimension()),
272                                             rowsToIntegerMatrix(B->getImpliedEquations(),c.ambientDimension())))==2*c.ambientDimension())))
273 	  {
274 	    numberOfComputedCones++;
275 	    c.canonicalize();
276 	    ret.cones.insert(c);
277 	    linealityConeFound=linealityConeFound || thisIsLinealityCone;
278 	  }
279 	else
280 	  {
281 	    conesSkipped++;
282 	  }
283       }
284   //  fprintf(Stderr,"Number of skipped cones: %i, lineality cone found: %i\n",conesSkipped,linealityConeFound);
285   //  fprintf(Stderr,"Number of computed cones: %i, number of unique cones: %i\n",numberOfComputedCones,ret.cones.size());
286 
287   return ret;
288 }
289 
290 
product(const PolyhedralFan & a,const PolyhedralFan & b)291 PolyhedralFan product(const PolyhedralFan &a, const PolyhedralFan &b)
292 {
293   PolyhedralFan ret(a.getAmbientDimension()+b.getAmbientDimension());
294 
295   for(PolyhedralConeList::const_iterator A=a.cones.begin();A!=a.cones.end();A++)
296     for(PolyhedralConeList::const_iterator B=b.cones.begin();B!=b.cones.end();B++)
297       ret.insert(product(*A,*B));
298 
299   return ret;
300 }
301 
302 
getRays(int dim)303 IntegerVectorList PolyhedralFan::getRays(int dim)
304 {
305   IntegerVectorList ret;
306   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
307     {
308       if(i->dimension()==dim)
309 	ret.push_back(i->getRelativeInteriorPoint());
310     }
311   return ret;
312 }
313 
314 
getRelativeInteriorPoints()315 IntegerVectorList PolyhedralFan::getRelativeInteriorPoints()
316 {
317   IntegerVectorList ret;
318   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
319     {
320       ret.push_back(i->getRelativeInteriorPoint());
321     }
322   return ret;
323 }
324 
325 
highestDimensionalCone() const326 PolyhedralCone const& PolyhedralFan::highestDimensionalCone()const
327 {
328   return *cones.rbegin();
329 }
330 
331 
convexHull() const332 PolyhedralCone PolyhedralFan::convexHull()const
333 {
334   if(cones.empty())return PolyhedralCone::givenByRays(IntegerVectorList(),IntegerVectorList(),n);
335   return PolyhedralCone::givenByRays(getRaysInPrintingOrder(0,0),cones.begin()->generatorsOfLinealitySpace(),n);
336 }
337 
338 
insert(PolyhedralCone const & c)339 void PolyhedralFan::insert(PolyhedralCone const &c)
340 {
341   cones.insert(c);
342 }
343 
344 
insertFacetsOfCone(PolyhedralCone const & c)345 void PolyhedralFan::insertFacetsOfCone(PolyhedralCone const &c)
346 {
347   PolyhedralFan facets=facetsOfCone(c);
348   for(PolyhedralConeList::const_iterator i=facets.cones.begin();i!=facets.cones.end();i++)insert(*i);
349 }
350 
351 
remove(PolyhedralCone const & c)352 void PolyhedralFan::remove(PolyhedralCone const &c)
353 {
354   cones.erase(c);
355 }
356 
removeAllExcept(int a)357 void PolyhedralFan::removeAllExcept(int a)
358 {
359   PolyhedralConeList::iterator i=cones.begin();
360   while(a>0)
361     {
362       if(i==cones.end())return;
363       i++;
364       a--;
365     }
366   cones.erase(i,cones.end());
367 }
368 
removeAllLowerDimensional()369 void PolyhedralFan::removeAllLowerDimensional()
370 {
371   if(!cones.empty())
372     {
373       int d=getMaxDimension();
374       PolyhedralConeList::iterator i=cones.begin();
375       while(i!=cones.end() && i->dimension()==d)i++;
376       cones.erase(i,cones.end());
377     }
378 }
379 
380 
facetComplex() const381 PolyhedralFan PolyhedralFan::facetComplex()const
382 {
383   //  fprintf(Stderr,"Computing facet complex...\n");
384   PolyhedralFan ret(n);
385 
386   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
387     {
388       PolyhedralFan a=facetsOfCone(*i);
389       for(PolyhedralConeList::const_iterator j=a.cones.begin();j!=a.cones.end();j++)
390 	ret.insert(*j);
391     }
392   //  fprintf(Stderr,"Done computing facet complex.\n");
393   return ret;
394 }
395 
396 
fullComplex() const397 PolyhedralFan PolyhedralFan::fullComplex()const
398 {
399   PolyhedralFan ret=*this;
400 
401   while(1)
402     {
403       log2 debug<<"looping";
404       bool doLoop=false;
405       PolyhedralFan facets=ret.facetComplex();
406       log2 debug<<"number of facets"<<facets.size()<<"\n";
407       for(PolyhedralConeList::const_iterator i=facets.cones.begin();i!=facets.cones.end();i++)
408 	if(!ret.contains(*i))
409 	  {
410 	    ret.insert(*i);
411 	    doLoop=true;
412 	  }
413       if(!doLoop)break;
414     }
415   return ret;
416 }
417 
418 
419 /*
420 PolyhedralFan PolyhedralFan::facetComplexSymmetry(SymmetryGroup const &sym, bool keepRays, bool dropLinealitySpace)const
421 {
422   log1 fprintf(Stderr,"Computing facet complex...\n");
423   PolyhedralFan ret(n);
424 
425   if(keepRays)
426     for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
427       if(i->dimension()==i->dimensionOfLinealitySpace()+1)ret.insert(*i);
428 
429   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
430     {
431       PolyhedralFan a=facetsOfCone(*i);
432       for(PolyhedralConeList::const_iterator j=a.cones.begin();j!=a.cones.end();j++)
433 	{
434 	  if((!dropLinealitySpace) || j->dimension()!=j->dimensionOfLinealitySpace())
435 	    {
436 	      IntegerVector v=j->getRelativeInteriorPoint();
437 	      bool alreadyInRet=false;
438 	      for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
439 		{
440 		  IntegerVector u=SymmetryGroup::compose(*k,v);
441 		  if(!j->containsRelatively(u))
442 		    {
443 		      for(PolyhedralConeList::const_iterator l=ret.cones.begin();l!=ret.cones.end();l++)
444 			if(l->containsRelatively(u))alreadyInRet=true;
445 		    }
446 		}
447 	      if(!alreadyInRet)ret.insert(*j);
448 	    }
449 	}
450     }
451   log1 fprintf(Stderr,"Done computing facet complex.\n");
452   return ret;
453 }
454 */
455 
facetComplexSymmetry(SymmetryGroup const & sym,bool keepRays,bool dropLinealitySpace) const456 PolyhedralFan PolyhedralFan::facetComplexSymmetry(SymmetryGroup const &sym, bool keepRays, bool dropLinealitySpace)const
457 {
458   log1 fprintf(Stderr,"Computing facet complex...\n");
459   PolyhedralFan ret(n);
460 
461   vector<IntegerVector> relIntTable;
462   vector<int> dimensionTable;
463 
464   if(keepRays)
465     for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
466       if(i->dimension()==i->dimensionOfLinealitySpace()+1)
467 	{
468 	  relIntTable.push_back(i->getRelativeInteriorPoint());
469 	  dimensionTable.push_back(i->dimension());
470 	  ret.insert(*i);
471 	}
472 
473   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();i++)
474     {
475       int iDim=i->dimension();
476       if(dropLinealitySpace && (i->dimension()==i->dimensionOfLinealitySpace()+1))continue;
477 
478       //      i->findFacets();
479       IntegerVectorList normals=i->getHalfSpaces();
480       for(IntegerVectorList::const_iterator j=normals.begin();j!=normals.end();j++)
481 	{
482 	  bool alreadyInRet=false;
483 	  for(int l=0;l<relIntTable.size();l++)
484 	    {
485 	      if(dimensionTable[l]==iDim-1)
486 		for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
487 		  {
488 		    IntegerVector u=SymmetryGroup::compose(*k,relIntTable[l]);
489 		    if((dotLong(*j,u)==0) && (i->contains(u)))
490 		      {
491 			alreadyInRet=true;
492 			goto exitLoop;
493 		      }
494 		  }
495 	    }
496 	exitLoop:
497 	  if(!alreadyInRet)
498 	    {
499 	      IntegerVectorList equations=i->getEquations();
500 	      IntegerVectorList inequalities=i->getHalfSpaces();
501 	      equations.push_back(*j);
502 	      PolyhedralCone c(inequalities,equations,n);
503 	      c.canonicalize();
504 	      ret.insert(c);
505 	      relIntTable.push_back(c.getRelativeInteriorPoint());
506 	      dimensionTable.push_back(c.dimension());
507 	    }
508 	}
509     }
510   log1 fprintf(Stderr,"Done computing facet complex.\n");
511   return ret;
512 }
513 
514 /*
515 IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup *sym)const
516 {
517   assert(!cones.empty());
518   int h=cones.begin()->dimensionOfLargestContainedSubspace();
519   PolyhedralFan f=*this;
520   while(f.getMaxDimension()!=h+1)
521     {
522       f=f.facetComplex();
523     }
524 
525   IntegerVectorList rays;
526 
527   PolyhedralFan done(n);
528   for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
529     if(!done.contains(*i))
530       for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
531 	{
532 	  PolyhedralCone cone=i->permuted(*k);
533 	  if(!done.contains(cone))
534 	    {
535 	      rays.push_back(cone.getRelativeInteriorPoint());
536 	      done.insert(cone);
537 	    }
538 	}
539   return rays;
540 }
541 */
542 
543 
stableRay(PolyhedralCone const & c,SymmetryGroup const * sym)544 IntegerVector PolyhedralFan::stableRay(PolyhedralCone const &c, SymmetryGroup const *sym)
545 {
546   PolyhedralCone C=c;//cast away const instead?
547 
548   IntegerVector v=C.getRelativeInteriorPoint();
549 
550   IntegerVector ret(v.size());
551   for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
552     {
553       IntegerVector v2=SymmetryGroup::compose(*k,v);
554       if(c.contains(v2))ret+=v2;
555     }
556   return normalized(ret);
557 }
558 
559 
stableRay(IntegerVector const & v,IntegerVectorList const & equations,IntegerVectorList const & inequalities,SymmetryGroup const * sym)560 IntegerVector PolyhedralFan::stableRay(IntegerVector const &v, IntegerVectorList const &equations, IntegerVectorList const &inequalities, SymmetryGroup const *sym)
561 {
562   IntegerVector ret(v.size());
563   for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
564     {
565       IntegerVector v2=SymmetryGroup::compose(*k,v);
566       bool containsV2=true;
567 
568       for(IntegerVectorList::const_iterator l=equations.begin();l!=equations.end();l++)
569 	if(dotLong(*l,v2)!=0)
570 	  {
571 	    containsV2=false;
572 	    goto leave;
573 	  }
574       for(IntegerVectorList::const_iterator l=inequalities.begin();l!=inequalities.end();l++)
575 	if(dotLong(*l,v2)<0)
576 	  {
577 	    containsV2=false;
578 	    goto leave;
579 	  }
580     leave:
581       if(containsV2)ret+=v2;
582     }
583   return normalized(ret);
584 }
585 
586 /* Slow version using facetComplexSymmetry()
587 IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup *sym)const
588 {
589   assert(!cones.empty());
590   int h=cones.begin()->dimensionOfLinealitySpace();
591   PolyhedralFan f=*this;
592   if(f.getMaxDimension()==h)return IntegerVectorList();
593   while(f.getMaxDimension()>h+1)
594     {
595       f=f.facetComplexSymmetry(*sym,true,true);
596     }
597   IntegerVectorList rays;
598 
599   for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
600     {
601       if(i->dimension()!=i->dimensionOfLinealitySpace())//This check is needed since the above while loop may not be run and therefore the lineality space may not have been removed.
602 	{
603 	  bool found=false;
604 	  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
605 	    if(i->contains(*j))found=true;
606 
607 	  if(!found)
608 	    {
609 	      //IntegerVector interiorPointForOrbit=i->getRelativeInteriorPoint();
610 	      IntegerVector interiorPointForOrbit=stableRay(*i,sym);
611 	      PolyhedralFan done(n);
612 	      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
613 		{
614 		  PolyhedralCone cone=i->permuted(*k);
615 
616 		  if(!done.contains(cone))
617 		    {
618 		      rays.push_back(SymmetryGroup::compose(*k,interiorPointForOrbit));
619 		      done.insert(cone);
620 		    }
621 		}
622 	    }
623 	}
624     }
625   return rays;
626   }*/
627 
rayComplexSymmetry(SymmetryGroup const & sym) const628 PolyhedralFan PolyhedralFan::rayComplexSymmetry(SymmetryGroup const &sym)const
629 {
630   //  log0 fprintf(Stderr,"rayComplexSymmetry - begin\n");
631   PolyhedralFan ret(n);
632   log1 fprintf(Stderr,"Computing rays of %i cones\n",(int)cones.size());
633   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
634     {
635       {
636 	static int t;
637 	if(!((t++)%10))log1 fprintf(Stderr,"%i\n",t);
638       }
639       //  log0 fprintf(Stderr,"calling\n");
640       IntegerVectorList rays=i->extremeRays();
641       //log0 fprintf(Stderr,"returning\n");
642       for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
643 	{
644 	  bool alreadyInRet=false;
645 	  for(PolyhedralConeList::const_iterator I=ret.cones.begin();I!=ret.cones.end();I++)
646 	    for(SymmetryGroup::ElementContainer::const_iterator k=sym.elements.begin();k!=sym.elements.end();k++)
647 	      {
648 		IntegerVector u=SymmetryGroup::compose(*k,*j);
649 		if((I->contains(u)))
650 		  {
651 		    alreadyInRet=true;
652 		    goto exitLoop;
653 		  }
654 	      }
655 	exitLoop:
656 	  IntegerVectorList equations=i->getEquations();
657 	  IntegerVectorList inequalities1=i->getHalfSpaces();
658 	  IntegerVectorList inequalities2;
659 	  for(IntegerVectorList::const_iterator k=inequalities1.begin();k!=inequalities1.end();k++)
660 	    {
661 	      if(dotLong(*j,*k))
662 		inequalities2.push_back(*k);
663 	      else
664 		equations.push_back(*k);
665 	    }
666 	  PolyhedralCone C(inequalities2,equations,n);
667 	  C.canonicalize();
668 	  ret.insert(C);
669 	}
670     }
671   //  log0 fprintf(Stderr,"rayComplexSymmetry - end\n");
672   return ret;
673 }
674 
675 
676 #if 0
677 IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup const *sym)const
678 {
679   assert(!cones.empty());
680   int h=cones.begin()->dimensionOfLinealitySpace();
681 
682   /*
683   PolyhedralFan f=*this;
684   if(f.getMaxDimension()==h)return IntegerVectorList();
685   while(f.getMaxDimension()>h+1)
686     {
687       f=f.facetComplexSymmetry(*sym,true,true);
688     }
689   */
690   PolyhedralFan f=rayComplexSymmetry(*sym);
691   IntegerVectorList rays;
692 
693   log1 fprintf(Stderr,"Number of cones in RayComplex: %i\n",f.cones.size());
694 
695   for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
696     {
697       static int t;
698       log1 fprintf(Stderr,"%i\n",t++);
699       if(i->dimension()!=i->dimensionOfLinealitySpace())//This check is needed since the above while loop may not be run and therefore the lineality space may not have been removed.
700 	{
701 	  bool found=false;
702 	  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
703 	    if(i->contains(*j))
704 	      {
705 		found=true;
706 		break;
707 	      }
708 	  if(!found)
709 	    {
710 	      //IntegerVector interiorPointForOrbit=i->getRelativeInteriorPoint();
711 	      IntegerVector interiorPointForOrbit=stableRay(*i,sym);
712 	      //    PolyhedralFan done(n);
713 
714 	      //Check that this works:
715 	      set<IntegerVector> thisOrbitsRays;
716 
717 	      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
718 		{
719 		  IntegerVector temp=SymmetryGroup::compose(*k,interiorPointForOrbit);
720 		  thisOrbitsRays.insert(temp);
721 		}
722 	      for(set<IntegerVector>::const_iterator i=thisOrbitsRays.begin();i!=thisOrbitsRays.end();i++)rays.push_back(*i);
723 	      //Instead of this:
724 	      /*	      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
725 		{
726 		  bool found=false;
727 		  IntegerVector temp=SymmetryGroup::compose(*k,interiorPointForOrbit);
728 		  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)//REWRITE WITH LOGARITHMIC SEARCH
729 		    if(*j==temp)
730 		      {
731 			found=true;
732 			break;
733 		      }
734 		  if(!found)
735 		    {
736 		      PolyhedralCone cone=i->permuted(*k);
737 		      rays.push_back(SymmetryGroup::compose(*k,interiorPointForOrbit));
738 		      //      done.insert(cone);
739 		    }
740 		}
741 	      */
742 	    }
743 	}
744     }
745   return rays;
746 }
747 
748 #elseif 0
749 //version used until Sep 2010
750 IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup const *sym, bool upToSymmetry)const
751 {
752 	SymmetryGroup localsym(n);
753 	if(!sym)sym=&localsym;
754   IntegerVectorList rays;
755   log1 fprintf(Stderr,"Computing rays of %i cones\n",cones.size());
756   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
757     {
758       {
759 	static int t;
760 	if(!((t++)%10))log1 fprintf(Stderr,"%i\n",t);
761       }
762       IntegerVectorList temp=i->extremeRays();
763       //      AsciiPrinter(Stderr).printVectorList(temp);
764 
765       for(IntegerVectorList::const_iterator j=temp.begin();j!=temp.end();j++)
766 	{
767 	  IntegerVectorList equations=i->getEquations();
768 	  IntegerVectorList inequalities1=i->getHalfSpaces();
769 	  IntegerVectorList inequalities2;
770 	  for(IntegerVectorList::const_iterator k=inequalities1.begin();k!=inequalities1.end();k++)
771 	    {
772 	      if(dotLong(*j,*k))
773 		inequalities2.push_back(*k);
774 	      else
775 		equations.push_back(*k);
776 	    }
777 	  bool isFound=false;
778 	  for(IntegerVectorList::const_iterator j2=rays.begin();j2!=rays.end();j2++)
779 	    for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
780 	      {
781 		bool isInCone=true;
782 		IntegerVector v=SymmetryGroup::compose(*k,*j2);
783 		for(IntegerVectorList::const_iterator l=equations.begin();l!=equations.end();l++)
784 		  if(dotLong(*l,v)!=0)
785 		    {
786 		      isInCone=false;
787 		      goto leave;
788 		    }
789 		for(IntegerVectorList::const_iterator l=inequalities2.begin();l!=inequalities2.end();l++)
790 		  if(dotLong(*l,v)<0)
791 		    {
792 		      isInCone=false;
793 		      goto leave;
794 		    }
795 	      leave:
796 		if(isInCone)
797 		  {
798 		    isFound=true;
799 		    goto leave2;
800 		  }
801 	      }
802 	leave2:
803 	  if(!isFound)
804 	    {
805 	      IntegerVector ray=stableRay(*j,equations,inequalities2,sym);
806 	      rays.push_back(ray);
807 	    }
808 	}
809     }
810   rays.sort();
811   if(upToSymmetry)return rays;
812   IntegerVectorList ret;
813   for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
814     {
815       set<IntegerVector> thisOrbitsRays;
816       for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
817 	{
818 	  IntegerVector temp=SymmetryGroup::compose(*k,*i);
819 	  thisOrbitsRays.insert(temp);
820 	}
821       for(set<IntegerVector>::const_iterator i=thisOrbitsRays.begin();i!=thisOrbitsRays.end();i++)ret.push_back(*i);
822     }
823   return ret;
824 }
825 #else
getRaysInPrintingOrder(SymmetryGroup const * sym,bool upToSymmetry) const826 IntegerVectorList PolyhedralFan::getRaysInPrintingOrder(SymmetryGroup const *sym, bool upToSymmetry)const
827 {
828   /*
829    * The ordering changed in this version. Previously the orbit representatives stored in "rays" were
830    * just the first extreme ray from the orbit that appeared. Now it is gotten using "orbitRepresentative"
831    * which causes the ordering in which the different orbits appear to change.
832    */
833 
834   if(cones.empty())return IntegerVectorList();
835   IntegerVectorList generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();//all cones have the same lineality space
836 
837         SymmetryGroup localsym(n);
838         if(!sym)sym=&localsym;
839   set<IntegerVector> rays;
840   log1 fprintf(Stderr,"Computing rays of %i cones\n",(int)cones.size());
841   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
842     {
843       {
844         static int t;
845         if(!((t++)%10))log1 fprintf(Stderr,"%i\n",t);
846       }
847       IntegerVectorList temp=i->extremeRays(&generatorsOfLinealitySpace);
848       for(IntegerVectorList::const_iterator j=temp.begin();j!=temp.end();j++)
849         rays.insert(sym->orbitRepresentative(*j));
850     }
851   IntegerVectorList ret;
852   if(upToSymmetry)
853     for(set<IntegerVector>::const_iterator i=rays.begin();i!=rays.end();i++)ret.push_back(*i);
854   else
855     for(set<IntegerVector>::const_iterator i=rays.begin();i!=rays.end();i++)
856       {
857         set<IntegerVector> thisOrbitsRays;
858         for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
859           {
860             IntegerVector temp=SymmetryGroup::compose(*k,*i);
861             thisOrbitsRays.insert(temp);
862           }
863         for(set<IntegerVector>::const_iterator i=thisOrbitsRays.begin();i!=thisOrbitsRays.end();i++)ret.push_back(*i);
864       }
865   return ret;
866 }
867 
868 
869 #endif
870 
871 
872 /*void PolyhedralFan::printWithIndices(class Printer *p, SymmetryGroup *sym)const //fan must be pure
873 {
874   //  print(p);
875   SymmetryGroup symm(n);
876   if(!sym)sym=&symm;
877   assert(!cones.empty());
878   int h=cones.begin()->dimensionOfLargestContainedSubspace();
879   fprintf(Stdout,"Rays:\n");
880   //  IntegerVectorList rays;//=f.getRelativeInteriorPoints();
881   fprintf(Stderr,"Computing rays...\n");
882   IntegerVectorList rays=getRaysInPrintingOrder(sym);
883   fprintf(Stderr,"Done computing rays.\n");
884   p->printVectorList(rays,true);
885   PolyhedralFan f=*this;
886 
887   //  while(f.getMaxDimension()>=h)
888   IntegerVector fvector(f.getMaxDimension()-h);
889   while(f.getMaxDimension()!=h)
890     {
891       int currentDimension=f.getMaxDimension()-h;
892       fprintf(Stderr,"Processing dimension %i cones...\n",currentDimension);
893       PolyhedralFan done(n);
894       IntegerVector rayIncidenceCounter(rays.size());
895       p->printString("Printing index list for dimension ");
896       p->printInteger(currentDimension);
897       p->printString(" cones:\n");
898       p->printString("{\n");
899       int faceIndex =0;
900       bool split=false;
901       bool addComma=false;
902       for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
903 	{
904 	  if(!done.contains(*i))
905 	    {
906 	      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
907 		{
908 		  PolyhedralCone cone=i->permuted(*k);
909 		  if(!done.contains(cone))
910 		    {
911 		      // p->printString("Face ");
912 		      // p->printInteger(faceIndex);
913 		      // p->printString(": ");
914 		      int rayIndex=0;
915 		      IntegerVector indices(0);
916 		      for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
917 			{
918 			  if(cone.contains(*j))
919 			    {
920 			      indices.grow(indices.size()+1);
921 			      indices[indices.size()-1]=rayIndex;
922 			      rayIncidenceCounter[rayIndex]++;
923 			    }
924 			  rayIndex++;
925 			}
926 		      if(addComma)
927 			{
928 			  p->printString(",");
929 			  p->printNewLine();
930 			  if(split)
931 			    p->printNewLine();
932 			  split=false;
933 			}
934 		      p->printVector(indices,true);
935 		      addComma=true;
936 		      faceIndex++;
937 		      done.insert(cone);
938 		    }
939 		}
940 	      split=true;//p->printNewLine();
941 	    }
942 	}
943       p->printString("}\n");
944       p->printString("Number of dimension ");
945       p->printInteger(f.getMaxDimension()-h);
946       p->printString(" cones incident to each ray:\n");
947       p->printVector(rayIncidenceCounter);
948       p->printNewLine();
949       fvector[f.getMaxDimension()-h-1]=faceIndex;
950       f=f.facetComplex();
951       fprintf(Stderr,"Done processing dimension %i cones.\n",currentDimension);
952       //      fvector.grow(fvector.size()+1);fvector[fvector.size()-1]=faceIndex;
953     }
954   p->printString("F-vector:\n");
955   p->printVector(fvector);
956   p->printNewLine();
957 }
958 */
959 
960  /*
961 void PolyhedralFan::printWithIndices(class Printer *p, SymmetryGroup *sym, const char *polymakeFileName)const //fan must be pure
962 {
963   IntegerVector multiplicities;
964   SymmetryGroup symm(n);
965   PolymakeFile polymakeFile;
966 
967 
968 
969 
970 
971 
972   if(polymakeFileName)
973     {
974       polymakeFile.create(polymakeFileName,"PolyhedralFan","PolyhedralFan");
975     }
976 
977   if(!sym)sym=&symm;
978   assert(!cones.empty());
979   int h=cones.begin()->dimensionOfLinealitySpace();
980   fprintf(Stdout,"Rays:\n");
981   //  IntegerVectorList rays;//=f.getRelativeInteriorPoints();
982   fprintf(Stderr,"Computing rays...\n");
983   IntegerVectorList rays=getRaysInPrintingOrder(sym);
984   fprintf(Stderr,"Done computing rays.\n");
985 
986 
987   SymmetricComplex symCom(n,rays,*sym);
988 
989 
990   if(p)
991     p->printVectorList(rays,true);
992   if(polymakeFileName)
993     {
994       polymakeFile.writeCardinalProperty("AMBIENT_DIM",n);
995       polymakeFile.writeCardinalProperty("DIM",getMaxDimension());
996       polymakeFile.writeCardinalProperty("LINEALITY_DIM",h);
997       polymakeFile.writeMatrixProperty("RAYS",rowsToIntegerMatrix(rays,n));
998       polymakeFile.writeCardinalProperty("N_RAYS",rays.size());
999       polymakeFile.writeMatrixProperty("LINEALITY_SPACE",rowsToIntegerMatrix(highestDimensionalCone().linealitySpace().dualCone().getEquations()));
1000       polymakeFile.writeMatrixProperty("ORTH_LINEALITY_SPACE",rowsToIntegerMatrix(highestDimensionalCone().linealitySpace().getEquations()));
1001       polymakeFile.writeCardinalProperty("PURE",1);
1002     }
1003 
1004 
1005   PolyhedralFan f=*this;
1006   stringstream s;
1007 
1008   //  while(f.getMaxDimension()>=h)
1009   IntegerVector fvector(f.getMaxDimension()-h);
1010   bool isHighestDimension=true;
1011   while(f.getMaxDimension()!=h)
1012     {
1013       int currentDimension=f.getMaxDimension()-h;
1014       fprintf(Stderr,"Processing dimension %i cones...\n",currentDimension);
1015       IntegerVector rayIncidenceCounter(rays.size());
1016       if(p)
1017 	{
1018 	  p->printString("Printing index list for dimension ");
1019 	  p->printInteger(currentDimension);
1020 	  p->printString(" cones:\n");
1021 	  p->printString("{\n");
1022 	}
1023       int faceIndex =0;
1024       bool split=false;
1025       bool addComma=false;
1026       for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
1027 	{
1028 	  {
1029 	    SymmetricComplex::Cone c;
1030 	    c.dimension=i->dimension();
1031 
1032 	    int rayIndex=0;
1033 	    for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
1034 	      {
1035 		if(i->contains(*j))c.insert(rayIndex);
1036 		rayIndex++;
1037 	      }
1038 	    symCom.insert(c);
1039 	  }
1040 
1041 
1042 	  PolyhedralFan done(n);
1043 	  for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
1044 	    {
1045 	      PolyhedralCone cone=i->permuted(*k);
1046 	      if(!done.contains(cone))
1047 		{
1048 		  if(isHighestDimension)
1049 		    {
1050 		      multiplicities.grow(multiplicities.size()+1);
1051 		      multiplicities[multiplicities.size()-1]=i->getMultiplicity();
1052 		    }
1053 		  // p->printString("Face ");
1054 		  // p->printInteger(faceIndex);
1055 		  // p->printString(": ");
1056 		  int rayIndex=0;
1057 		  IntegerVector indices(0);
1058 		  for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
1059 		    {
1060 		      if(cone.contains(*j))
1061 			{
1062 			  indices.grow(indices.size()+1);
1063 			  indices[indices.size()-1]=rayIndex;
1064 			  rayIncidenceCounter[rayIndex]++;
1065 			}
1066 		      rayIndex++;
1067 		    }
1068 		  if(addComma)
1069 		    {
1070 		      if(p)
1071 			{
1072 			  p->printString(",");
1073 			  p->printNewLine();
1074 			}
1075 		      if(isHighestDimension)s<<endl;
1076 		      if(split)
1077 			{
1078 			  if(p)p->printNewLine();
1079 			  //	  s<<endl;
1080 			}
1081 		      split=false;
1082 		    }
1083 		  if(p)
1084 		    p->printVector(indices,true);
1085 		  if(isHighestDimension)
1086 		    {
1087 		      s << '{';
1088 		      for(int i=0;i<indices.size();i++)
1089 			{
1090 			  if(i)s<<" ";
1091 			  s<<indices[i];
1092 			}
1093 		      s << '}';
1094 		    }
1095 		  addComma=true;
1096 		  faceIndex++;
1097 		  done.insert(cone);
1098 		}
1099 	    }
1100 	  split=true;//p->printNewLine();
1101 	}
1102       if(p)
1103 	{
1104 	  p->printString("}\n");
1105 	  p->printString("Number of dimension ");
1106 	  p->printInteger(f.getMaxDimension()-h);
1107 	  p->printString(" cones incident to each ray:\n");
1108 	  p->printVector(rayIncidenceCounter);
1109 	  p->printNewLine();
1110 	}
1111       fvector[f.getMaxDimension()-h-1]=faceIndex;
1112       fprintf(Stderr,"TESTESTSETST\n");
1113       f=f.facetComplexSymmetry(*sym);
1114       fprintf(Stderr,"TESTESTSETST\n");
1115       fprintf(Stderr,"Done processing dimension %i cones.\n",currentDimension);
1116       //      fvector.grow(fvector.size()+1);fvector[fvector.size()-1]=faceIndex;
1117       isHighestDimension=false;
1118     }
1119   if(p)
1120     {
1121       p->printString("Multiplicities:\n");
1122       p->printVector(multiplicities);
1123       p->printNewLine();
1124       p->printString("F-vector:\n");
1125       p->printVector(fvector);
1126       p->printNewLine();
1127     }
1128 
1129   if(polymakeFileName)
1130     {
1131       polymakeFile.writeCardinalVectorProperty("F_VECTOR",fvector);
1132       s<<endl;
1133       polymakeFile.writeStringProperty("MAXIMAL_CONES",s.str());
1134 
1135       IntegerVectorList m;
1136       m.push_back(multiplicities);
1137       polymakeFile.writeMatrixProperty("MULTIPLICITIES",rowsToIntegerMatrix(m).transposed());
1138     }
1139   if(polymakeFileName)
1140     polymakeFile.close();
1141 
1142 
1143   AsciiPrinter(Stdout).printString(symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,true));
1144 }
1145 */
1146 
1147 
1148 
1149   /*MARKS CONES AS NONMAXIMAL IN THE SYMMETRIC COMPLEX IN WHICH THEY WILL BE INSERTED -not to be confused with the facet testing in the code
1150    */
computeFacets(SymmetricComplex::Cone const & theCone,IntegerMatrix const & rays,IntegerVectorList const & facetCandidates,SymmetricComplex const & theComplex)1151  static list<SymmetricComplex::Cone> computeFacets(SymmetricComplex::Cone const &theCone, IntegerMatrix const &rays, IntegerVectorList const &facetCandidates, SymmetricComplex const &theComplex/*, int linealityDim*/)
1152 {
1153   set<SymmetricComplex::Cone> ret;
1154 
1155   for(IntegerVectorList::const_iterator i=facetCandidates.begin();i!=facetCandidates.end();i++)
1156     {
1157       set<int> indices;
1158 
1159       bool notAll=false;
1160       for(vector<int>::const_iterator j=theCone.indices.begin();j!=theCone.indices.end();j++)
1161 	if(dotLong(rays[*j],*i)==0)
1162 	  indices.insert(*j);
1163 	else
1164 	  notAll=true;
1165 
1166       SymmetricComplex::Cone temp(indices,theCone.dimension-1,0,false,theComplex);
1167       /*      temp.multiplicity=0;
1168       temp.dimension=theCone.dimension-1;
1169       temp.setIgnoreSymmetry(true);
1170       */
1171       if(notAll)ret.insert(temp);
1172 
1173     }
1174   //  fprintf(Stderr,"HEJ!!!!\n");
1175 
1176   list<SymmetricComplex::Cone> ret2;
1177   for(set<SymmetricComplex::Cone>::const_iterator i=ret.begin();i!=ret.end();i++)
1178     {
1179       bool isMaximal=true;
1180 
1181       /*      if(i->indices.size()+linealityDim<i->dimension)//#3
1182 	isMaximal=false;
1183 	else*/
1184 	for(set<SymmetricComplex::Cone>::const_iterator j=ret.begin();j!=ret.end();j++)
1185 	  if(i!=j && i->isSubsetOf(*j))
1186 	    {
1187 	      isMaximal=false;
1188 	      break;
1189 	    }
1190       if(isMaximal)
1191 	{
1192 	  SymmetricComplex::Cone temp(i->indexSet(),i->dimension,i->multiplicity,true,theComplex);
1193 	  temp.setKnownToBeNonMaximal(); // THIS IS WHERE WE SET THE CONES NON-MAXIMAL FLAG
1194 	  //	  temp.setIgnoreSymmetry(false);
1195 	  ret2.push_back(temp);
1196 	}
1197     }
1198   return ret2;
1199 }
1200 
1201 
addFacesToSymmetricComplex(SymmetricComplex & c,PolyhedralCone const & cone,IntegerVectorList const & facetCandidates,IntegerVectorList const & generatorsOfLinealitySpace)1202 void addFacesToSymmetricComplex(SymmetricComplex &c, PolyhedralCone const &cone, IntegerVectorList const &facetCandidates, IntegerVectorList const &generatorsOfLinealitySpace)
1203 {
1204   IntegerMatrix const &rays=c.getVertices();
1205   set<int> indices;
1206 
1207 //  for(int j=0;j<rays.getHeight();j++)if(cone.contains(rays[j]))indices.insert(j);
1208 
1209   IntegerVectorList l=cone.extremeRays(&generatorsOfLinealitySpace);
1210   for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)indices.insert(c.indexOfVertex(*i));
1211 
1212   addFacesToSymmetricComplex(c,indices,facetCandidates,cone.dimension(),cone.getMultiplicity());
1213 }
1214 
addFacesToSymmetricComplex(SymmetricComplex & c,set<int> const & indices,IntegerVectorList const & facetCandidates,int dimension,int multiplicity)1215 void addFacesToSymmetricComplex(SymmetricComplex &c, set<int> const &indices, IntegerVectorList const &facetCandidates, int dimension, int multiplicity)
1216 {
1217   IntegerMatrix const &rays=c.getVertices();
1218   list<SymmetricComplex::Cone> clist;
1219   {
1220 
1221     SymmetricComplex::Cone temp(indices,dimension,multiplicity,true,c);
1222     //    temp.dimension=cone.dimension();
1223     //   temp.multiplicity=cone.getMultiplicity();
1224     clist.push_back(temp);
1225   }
1226 
1227   //  int linealityDim=cone.dimensionOfLinealitySpace();
1228 
1229   while(!clist.empty())
1230     {
1231       SymmetricComplex::Cone &theCone=clist.front();
1232 
1233       if(!c.contains(theCone))
1234 	{
1235 	  log2
1236 	  {
1237 	    static int t;
1238 	    if((t&1023)==0)
1239 	      {
1240 	    	fprintf(Stderr,"clist size:%i\n",(int)clist.size());
1241 	      }
1242 	    t++;
1243 	  }
1244 
1245 	  c.insert(theCone);
1246 	  //	  log0 fprintf(Stderr,"ADDING\n");
1247 	  list<SymmetricComplex::Cone> facets=computeFacets(theCone,rays,facetCandidates,c/*,linealityDim*/);
1248 	  clist.splice(clist.end(),facets);
1249 	}
1250       clist.pop_front();
1251     }
1252 
1253 }
1254 
1255 
1256 /**
1257    Produce strings that express the vectors in terms of rays of the fan modulo the lineality space. Symmetry is ignored??
1258  */
renamingStrings(IntegerVectorList const & theVectors,IntegerVectorList const & originalRays,IntegerVectorList const & linealitySpace,SymmetryGroup * sym) const1259 vector<string> PolyhedralFan::renamingStrings(IntegerVectorList const &theVectors, IntegerVectorList const &originalRays, IntegerVectorList const &linealitySpace, SymmetryGroup *sym)const
1260 {
1261   vector<string> ret;
1262   for(IntegerVectorList::const_iterator i=theVectors.begin();i!=theVectors.end();i++)
1263     {
1264       for(PolyhedralConeList::const_iterator j=cones.begin();j!=cones.end();j++)
1265 	{
1266 	  if(j->contains(*i))
1267 	    {
1268 	      vector<int> relevantIndices;
1269 	      IntegerVectorList relevantRays=linealitySpace;
1270 	      int K=0;
1271 	      for(IntegerVectorList::const_iterator k=originalRays.begin();k!=originalRays.end();k++,K++)
1272 		if(j->contains(*k))
1273 		  {
1274 		    relevantIndices.push_back(K);
1275 		    relevantRays.push_back(*k);
1276 		  }
1277 
1278 	      FieldMatrix LFA(Q,relevantRays.size(),n);
1279 	      int J=0;
1280 	      for(IntegerVectorList::const_iterator j=relevantRays.begin();j!=relevantRays.end();j++,J++)
1281 		LFA[J]=integerVectorToFieldVector(*j,Q);
1282 	      FieldVector LFB=concatenation(integerVectorToFieldVector(*i,Q),FieldVector(Q,relevantRays.size()));
1283 	      LFA=LFA.transposed();
1284 	      FieldVector LFX=LFA.solver().canonicalize(LFB);
1285 	      stringstream s;
1286 	      if(LFX.subvector(0,n).isZero())
1287 	        {
1288 		  s<<"Was:";
1289 	          FieldVector S=LFX.subvector(n+linealitySpace.size(),LFX.size());
1290 		  for(int k=0;k<S.size();k++)
1291 		    if(!S[k].isZero())
1292 		      s<<"+"<<S[k].toString()<<"*["<<relevantIndices[k]<<"] ";
1293 	        }
1294 	      ret.push_back(s.str());
1295 	      break;
1296 	    }
1297 	}
1298     }
1299   return ret;
1300 }
1301 
toSymmetricComplex(SymmetryGroup * sym)1302 SymmetricComplex PolyhedralFan::toSymmetricComplex(SymmetryGroup *sym)
1303 {
1304   SymmetryGroup symm(n);
1305 	  if(!sym)sym=&symm;
1306 
1307 	  IntegerVectorList rays=getRaysInPrintingOrder(sym);
1308 
1309 	  SymmetricComplex symCom(n,rays,*sym);
1310 
1311 	  if(cones.empty())return symCom;
1312 	  IntegerVectorList generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();
1313 
1314 	  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1315 	    {
1316 	      {
1317 		static int t;
1318 		log1 fprintf(Stderr,"Adding faces of cone %i\n",t++);
1319 	      }
1320 	      log2 fprintf(Stderr,"Dim: %i\n",i->dimension());
1321 
1322 	      addFacesToSymmetricComplex(symCom,*i,i->getHalfSpaces(),generatorsOfLinealitySpace);
1323 	    }
1324 
1325 	  log1 cerr<<"Remapping";
1326 	  symCom.remap();
1327 	  log1 cerr<<"Done remapping";
1328 	  return symCom;
1329 }
1330 
printWithIndices(class Printer * p,int flags,SymmetryGroup * sym,vector<string> const * comments) const1331 void PolyhedralFan::printWithIndices(class Printer *p, int flags, SymmetryGroup *sym, vector<string> const *comments)const
1332 //void PolyhedralFan::printWithIndices(class Printer *p, bool printMultiplicities, SymmetryGroup *sym, bool group, bool ignoreCones, bool xml, bool tPlaneSort, vector<string> const *comments)const
1333 {
1334   assert(p);
1335   //  IntegerVector multiplicities;
1336   SymmetryGroup symm(n);
1337 
1338   PolymakeFile polymakeFile;
1339 //  polymakeFile.create("NONAME","fan","PolyhedralFan",flags&FPF_xml);
1340   polymakeFile.create("NONAME","fan","SymmetricFan",flags&FPF_xml);
1341 
1342   bool produceXml=polymakeFile.isXmlFormat();
1343 
1344 
1345   if(!sym)sym=&symm;
1346 
1347   if(cones.empty())
1348     {
1349       p->printString("Polyhedral fan is empty. Printing not supported.\n");
1350       return;
1351     }
1352 
1353   int h=cones.begin()->dimensionOfLinealitySpace();
1354 
1355   log1 fprintf(Stderr,"Computing rays.\n");
1356   IntegerVectorList rays=getRaysInPrintingOrder(sym);
1357 
1358   SymmetricComplex symCom(n,rays,*sym);
1359 
1360   polymakeFile.writeCardinalProperty("AMBIENT_DIM",n);
1361   polymakeFile.writeCardinalProperty("DIM",getMaxDimension());
1362   polymakeFile.writeCardinalProperty("LINEALITY_DIM",h);
1363   polymakeFile.writeMatrixProperty("RAYS",rowsToIntegerMatrix(rays,n),true,comments);
1364   polymakeFile.writeCardinalProperty("N_RAYS",rays.size());
1365   IntegerVectorList linealitySpaceGenerators=highestDimensionalCone().linealitySpace().dualCone().getEquations();
1366   polymakeFile.writeMatrixProperty("LINEALITY_SPACE",rowsToIntegerMatrix(linealitySpaceGenerators,n));
1367   polymakeFile.writeMatrixProperty("ORTH_LINEALITY_SPACE",rowsToIntegerMatrix(highestDimensionalCone().linealitySpace().getEquations(),n));
1368 
1369   if(flags & FPF_primitiveRays)
1370   {
1371 	 IntegerVectorList primitiveRays;
1372 	 for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
1373 		 for(PolyhedralConeList::const_iterator j=cones.begin();j!=cones.end();j++)
1374 			 if(j->contains(*i)&&(j->dimensionOfLinealitySpace()+1==j->dimension()))
1375 					 primitiveRays.push_back(j->semiGroupGeneratorOfRay());
1376 
1377 	  polymakeFile.writeMatrixProperty("PRIMITIVE_RAYS",rowsToIntegerMatrix(primitiveRays,n));
1378   }
1379 
1380 
1381   IntegerVectorList generatorsOfLinealitySpace=cones.begin()->generatorsOfLinealitySpace();
1382 
1383   log1 fprintf(Stderr,"Building symmetric complex.\n");
1384   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1385     {
1386       {
1387 	static int t;
1388 	log1 fprintf(Stderr,"Adding faces of cone %i\n",t++);
1389       }
1390       log2 fprintf(Stderr,"Dim: %i\n",i->dimension());
1391 
1392       addFacesToSymmetricComplex(symCom,*i,i->getHalfSpaces(),generatorsOfLinealitySpace);
1393     }
1394 
1395   log1 cerr<<"Remapping";
1396   symCom.remap();
1397   log1 cerr<<"Done remapping";
1398 
1399 
1400   PolyhedralFan f=*this;
1401 
1402   //  IntegerVector fvector(f.getMaxDimension()-h);
1403 
1404 
1405 
1406   //fprintf(Stderr,"maxdim %i h %i\n",f.getMaxDimension(),h);
1407   /*  while(!f.cones.empty())
1408     {
1409       int currentDimension=f.getMaxDimension()-h;
1410       IntegerVector rayIncidenceCounter(rays.size());
1411       int faceIndex =0;
1412       bool split=false;
1413       bool addComma=false;
1414       for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
1415 	{
1416 	  {
1417 	    SymmetricComplex::Cone c;
1418 	    c.dimension=i->dimension();
1419 	    c.multiplicity=i->getMultiplicity();
1420 
1421 	    int rayIndex=0;
1422 	    for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
1423 	      {
1424 		if(i->contains(*j))c.insert(rayIndex);
1425 		rayIndex++;
1426 	      }
1427 	    symCom.insert(c);
1428 	  }
1429 	}
1430       //      fvector[f.getMaxDimension()-h-1]=faceIndex;
1431       f=f.facetComplexSymmetry(*sym);
1432     }
1433   */
1434   log1 fprintf(Stderr,"Computing f-vector.\n");
1435   IntegerVector fvector=symCom.fvector();
1436   polymakeFile.writeCardinalVectorProperty("F_VECTOR",fvector);
1437   log1 fprintf(Stderr,"Done computing f-vector.\n");
1438 
1439   if(flags&FPF_boundedInfo)
1440     {
1441       log1 fprintf(Stderr,"Computing bounded f-vector.\n");
1442       IntegerVector fvectorBounded=symCom.fvector(true);
1443       polymakeFile.writeCardinalVectorProperty("F_VECTOR_BOUNDED",fvectorBounded);
1444       log1 fprintf(Stderr,"Done computing bounded f-vector.\n");
1445     }
1446 /*
1447  * Removed to make the Polymake people happy.
1448  *    {
1449     int euler=0;
1450     int mul=-1;
1451     for(int i=0;i<fvector.size();i++,mul*=-1)euler+=mul*fvector[i];
1452     polymakeFile.writeCardinalProperty("MY_EULER",euler);
1453   }
1454 */
1455   log1 fprintf(Stderr,"Checking if complex is simplicial and pure.\n");
1456   polymakeFile.writeBooleanProperty("SIMPLICIAL",symCom.isSimplicial());
1457   polymakeFile.writeBooleanProperty("PURE",symCom.isPure());
1458   log1 fprintf(Stderr,"Done checking.\n");
1459 
1460 
1461   if(flags&FPF_conesCompressed)
1462   {
1463     polymakeFile.writeArrayArrayIntProperty("SYMMETRY_GENERATORS",rowsToIntegerMatrix(sym->getUniqueGenerators(),n));
1464     log1 fprintf(Stderr,"Producing list of cones up to symmetry.\n");
1465     polymakeFile.writeStringProperty("CONES_ORBITS",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),false,flags&FPF_group,0,true,flags&FPF_tPlaneSort,produceXml));
1466     log1 fprintf(Stderr,"Done producing list of cones up to symmetry.\n");
1467     log1 fprintf(Stderr,"Producing list of maximal cones up to symmetry.\n");
1468     stringstream multiplicities;
1469     polymakeFile.writeStringProperty("MAXIMAL_CONES_ORBITS",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,flags&FPF_group, &multiplicities,true,flags&FPF_tPlaneSort,produceXml));
1470     if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES_ORBITS",multiplicities.str());
1471     log1 fprintf(Stderr,"Done producing list of maximal cones up to symmetry.\n");
1472   }
1473 
1474   if(flags&FPF_conesExpanded)
1475     {
1476       if(flags&FPF_cones)
1477 	{
1478 	  log1 fprintf(Stderr,"Producing list of cones.\n");
1479 	  polymakeFile.writeStringProperty("CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),false,flags&FPF_group,0,false,flags&FPF_tPlaneSort,produceXml));
1480 	  log1 fprintf(Stderr,"Done producing list of cones.\n");
1481 	}
1482       if(flags&FPF_maximalCones)
1483 	{
1484 	  log1 fprintf(Stderr,"Producing list of maximal cones.\n");
1485 	  stringstream multiplicities;
1486 	  polymakeFile.writeStringProperty("MAXIMAL_CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,flags&FPF_group, &multiplicities,false,flags&FPF_tPlaneSort,produceXml));
1487 	  if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES",multiplicities.str());
1488 	  log1 fprintf(Stderr,"Done producing list of maximal cones.\n");
1489 	}
1490     }
1491 
1492   if(flags&FPF_values)
1493     {
1494       {
1495 	IntegerVectorList values;
1496 	for(IntegerVectorList::const_iterator i=linealitySpaceGenerators.begin();i!=linealitySpaceGenerators.end();i++)
1497 	  {
1498 	    IntegerVector v(1);
1499 	    v[0]=evaluatePiecewiseLinearFunction(*i);
1500 	    values.push_back(v);
1501 	  }
1502 	polymakeFile.writeMatrixProperty("LINEALITY_VALUES",rowsToIntegerMatrix(values,1));
1503       }
1504       {
1505 	IntegerVectorList values;
1506 	for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
1507 	  {
1508 	    IntegerVector v(1);
1509 	    v[0]=evaluatePiecewiseLinearFunction(*i);
1510 	    values.push_back(v);
1511 	  }
1512 	polymakeFile.writeMatrixProperty("RAY_VALUES",rowsToIntegerMatrix(values,1));
1513       }
1514     }
1515 
1516 
1517   log1 fprintf(Stderr,"Producing final string for output.\n");
1518   stringstream s;
1519   polymakeFile.writeStream(s);
1520   string S=s.str();
1521   log1 fprintf(Stderr,"Printing string.\n");
1522   p->printString(S.c_str());
1523   log1 fprintf(Stderr,"Done printing string.\n");
1524 }
1525 
1526 
readFan(string const & filename,bool onlyMaximal,IntegerVector * w,set<int> const * coneIndices,SymmetryGroup const * sym,bool readCompressedIfNotSym)1527 PolyhedralFan PolyhedralFan::readFan(string const &filename, bool onlyMaximal, IntegerVector *w, set<int> const *coneIndices, SymmetryGroup const *sym, bool readCompressedIfNotSym)
1528 {
1529     PolymakeFile inFile;
1530     inFile.open(filename.c_str());
1531 
1532     int n=inFile.readCardinalProperty("AMBIENT_DIM");
1533     int nRays=inFile.readCardinalProperty("N_RAYS");
1534     IntegerMatrix rays=inFile.readMatrixProperty("RAYS",nRays,n);
1535     int linealityDim=inFile.readCardinalProperty("LINEALITY_DIM");
1536     IntegerMatrix linealitySpace=inFile.readMatrixProperty("LINEALITY_SPACE",linealityDim,n);
1537 
1538 
1539     const char *sectionName=0;
1540     const char *sectionNameMultiplicities=0;
1541     if(sym || readCompressedIfNotSym)
1542     {
1543       sectionName=(onlyMaximal)?"MAXIMAL_CONES_ORBITS":"CONES_ORBITS";
1544       sectionNameMultiplicities=(onlyMaximal)?"MULTIPLICITIES_ORBITS":"DUMMY123";
1545     }
1546       else
1547       {  sectionName=(onlyMaximal)?"MAXIMAL_CONES":"CONES";
1548       sectionNameMultiplicities=(onlyMaximal)?"MULTIPLICITIES":"DUMMY123";
1549       }
1550 
1551 
1552     IntegerVector w2(n);
1553     if(w==0)w=&w2;
1554 
1555     SymmetryGroup sym2(n);
1556     if(sym==0)sym=&sym2;
1557 
1558     vector<list<int> > cones=inFile.readMatrixIncidenceProperty(sectionName);
1559     IntegerVectorList r;
1560 
1561     bool hasMultiplicities=inFile.hasProperty(sectionNameMultiplicities);
1562     IntegerMatrix multiplicities(0,0);
1563     if(hasMultiplicities)multiplicities=inFile.readMatrixProperty(sectionNameMultiplicities,cones.size(),1);
1564 
1565 
1566     PolyhedralFan ret(n);
1567 
1568     log2 cerr<< "Number of orbits to expand "<<cones.size()<<endl;
1569     for(int i=0;i<cones.size();i++)
1570       if(coneIndices==0 || coneIndices->count(i))
1571 	{
1572 	  log2 cerr<<"Expanding symmetries of cone"<<endl;
1573 	  /*	  for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
1574 	    {
1575 	      IntegerVectorList coneRays;
1576 	      for(list<int>::const_iterator j=cones[i].begin();j!=cones[i].end();j++)
1577 		coneRays.push_back(SymmetryGroup::compose(*perm,rays[*j]));
1578 	      if(isInNonNegativeSpan(*w,coneRays,linealitySpace.getRows()))
1579 		{
1580 		  PolyhedralCone C=PolyhedralCone::givenByRays(coneRays,linealitySpace.getRows(),n);
1581 		  C.canonicalize();
1582 		  ret.insert(C);
1583 		}
1584 	    }
1585 	  */
1586 	  {
1587 	    IntegerVectorList coneRays;
1588 	    for(list<int>::const_iterator j=cones[i].begin();j!=cones[i].end();j++)
1589 	      coneRays.push_back((rays[*j]));
1590 	    PolyhedralCone C=PolyhedralCone::givenByRays(coneRays,linealitySpace.getRows(),n);
1591 	    if(hasMultiplicities)C.setMultiplicity(multiplicities[i][0]);
1592 	    for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
1593 	      {
1594 		if(C.contains(SymmetryGroup::composeInverse(*perm,*w)))
1595 		  {
1596 		    PolyhedralCone C2=C.permuted(*perm);
1597 		    C2.canonicalize();
1598 		    ret.insert(C2);
1599 		  }
1600 	      }
1601 	  }
1602 	}
1603     return ret;
1604 }
1605 
1606 
getIncidenceList(SymmetryGroup * sym) const1607 IncidenceList PolyhedralFan::getIncidenceList(SymmetryGroup *sym)const //fan must be pure
1608 {
1609   IncidenceList ret;
1610   SymmetryGroup symm(n);
1611   if(!sym)sym=&symm;
1612   assert(!cones.empty());
1613   int h=cones.begin()->dimensionOfLinealitySpace();
1614   IntegerVectorList rays=getRaysInPrintingOrder(sym);
1615   PolyhedralFan f=*this;
1616 
1617   while(f.getMaxDimension()!=h)
1618     {
1619       IntegerVectorList l;
1620       PolyhedralFan done(n);
1621       IntegerVector rayIncidenceCounter(rays.size());
1622       int faceIndex =0;
1623       for(PolyhedralConeList::const_iterator i=f.cones.begin();i!=f.cones.end();i++)
1624 	{
1625 	  if(!done.contains(*i))
1626 	    {
1627 	      for(SymmetryGroup::ElementContainer::const_iterator k=sym->elements.begin();k!=sym->elements.end();k++)
1628 		{
1629 		  PolyhedralCone cone=i->permuted(*k);
1630 		  if(!done.contains(cone))
1631 		    {
1632 		      int rayIndex=0;
1633 		      IntegerVector indices(0);
1634 		      for(IntegerVectorList::const_iterator j=rays.begin();j!=rays.end();j++)
1635 			{
1636 			  if(cone.contains(*j))
1637 			    {
1638 			      indices.grow(indices.size()+1);
1639 			      indices[indices.size()-1]=rayIndex;
1640 			      rayIncidenceCounter[rayIndex]++;
1641 			    }
1642 			  rayIndex++;
1643 			}
1644 		      l.push_back(indices);
1645 		      faceIndex++;
1646 		      done.insert(cone);
1647 		    }
1648 		}
1649 	    }
1650 	}
1651       ret[f.getMaxDimension()]=l;
1652       f=f.facetComplex();
1653     }
1654   return ret;
1655 }
1656 
1657 
makePure()1658 void PolyhedralFan::makePure()
1659 {
1660   if(getMaxDimension()!=getMinDimension())removeAllLowerDimensional();
1661 }
1662 
contains(PolyhedralCone const & c) const1663 bool PolyhedralFan::contains(PolyhedralCone const &c)const
1664 {
1665   return cones.count(c);
1666 }
1667 
1668 
coneContaining(IntegerVector const & v) const1669 PolyhedralCone PolyhedralFan::coneContaining(IntegerVector const &v)const
1670 {
1671   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1672     if(i->contains(v))return i->faceContaining(v);
1673   debug<<"Vector "<<v<<" not contained in support of fan\n";
1674   assert(0);
1675 }
1676 
1677 
conesBegin() const1678 PolyhedralFan::coneIterator PolyhedralFan::conesBegin()const
1679 {
1680   return cones.begin();
1681 }
1682 
1683 
conesEnd() const1684 PolyhedralFan::coneIterator PolyhedralFan::conesEnd()const
1685 {
1686   return cones.end();
1687 }
1688 
1689 
isRefinementOf(PolyhedralFan const & f) const1690 bool PolyhedralFan::isRefinementOf(PolyhedralFan const &f)const
1691 {
1692   /*  for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1693     {
1694       static int t;
1695       fprintf(Stderr,"%i\n",t++);
1696       for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
1697 	{
1698 	  PolyhedralCone c=intersection(*i,*j);
1699 	  if(c.dimension()==n)
1700 	    {
1701 
1702 	      if(!j->contains(*i))
1703 		return false;
1704 	    }
1705 	}
1706 	}*/
1707 
1708   int t=0;
1709   for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
1710     {
1711       int t1=0;
1712       for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1713 	{
1714 	  PolyhedralCone c=intersection(*i,*j);
1715 	  if(c.dimension()==n)
1716 	    {
1717 	      t1++;
1718 	      if(!j->contains(*i))
1719 		return false;
1720 	    }
1721 	}
1722       fprintf(Stderr,"%i\n",t1);
1723       t+=t1;
1724     }
1725   fprintf(Stderr,"%i\n",t);
1726 
1727       /*      bool found=false;
1728       for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
1729 	  if(j->contains(*i)){
1730 	    found=true;
1731 	    break;
1732 	  }
1733       if(!found)
1734 	{
1735 	  for(PolyhedralConeList::const_iterator j=f.cones.begin();j!=f.cones.end();j++)
1736 	    {
1737 	      PolyhedralCone c=intersection(*i,*j);
1738 	      if(c.dimension()==n){
1739 		AsciiPrinter(Stderr).printVector(c.getRelativeInteriorPoint());
1740 	      }
1741 	    }
1742 
1743 	  return false;
1744 	  }*/
1745 
1746   return true;
1747 }
1748 
1749 /*	    for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
1750 	      {
1751 		if(C.contains(SymmetryGroup::composeInverse(*perm,*w)))
1752 		  {
1753 		    PolyhedralCone C2=C.permuted(*perm);
1754 		    C2.canonicalize();
1755 		    ret.insert(C2);
1756 		  }
1757 	      }
1758 */
link(IntegerVector const & w,SymmetryGroup * sym) const1759 PolyhedralFan PolyhedralFan::link(IntegerVector const &w, SymmetryGroup *sym)const
1760 {
1761   SymmetryGroup symL(n);
1762   if(!sym)sym=&symL;
1763 
1764   PolyhedralFan ret(n);
1765 
1766   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1767     {
1768       for(SymmetryGroup::ElementContainer::const_iterator perm=sym->elements.begin();perm!=sym->elements.end();perm++)
1769 	{
1770 	  IntegerVector w2=SymmetryGroup::composeInverse(*perm,w);
1771 	  if(i->contains(w2))
1772 	    {
1773 	      IntegerVectorList equations=i->getEquations();
1774 	      IntegerVectorList inequalities1=i->getHalfSpaces();
1775 	      IntegerVectorList inequalities2;
1776 	      for(IntegerVectorList::const_iterator j=inequalities1.begin();j!=inequalities1.end();j++)
1777 		if(dotLong(w2,*j)==0)inequalities2.push_back(*j);
1778 	      PolyhedralCone C(inequalities2,equations,n);
1779 	      C.canonicalize();
1780 	      C.setLinearForm(i->getLinearForm());
1781 	      PolyhedralCone C2=C.permuted(*perm);
1782 	      C2.canonicalize();
1783 	      C2.setMultiplicity(i->getMultiplicity());
1784 	      ret.insert(C2);
1785 	    }
1786 	}
1787     }
1788   return ret;
1789 }
1790 
link(IntegerVector const & w) const1791 PolyhedralFan PolyhedralFan::link(IntegerVector const &w)const
1792 {
1793   PolyhedralFan ret(n);
1794 
1795   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1796     {
1797       if(i->contains(w))
1798 	{
1799 	  IntegerVectorList equations=i->getEquations();
1800 	  IntegerVectorList inequalities1=i->getHalfSpaces();
1801 	  IntegerVectorList inequalities2;
1802 	  for(IntegerVectorList::const_iterator j=inequalities1.begin();j!=inequalities1.end();j++)
1803 	    if(dotLong(w,*j)==0)inequalities2.push_back(*j);
1804 	  PolyhedralCone C(inequalities2,equations,n);
1805 	  C.canonicalize();
1806 	  C.setLinearForm(i->getLinearForm());
1807       C.setMultiplicity(i->getMultiplicity());
1808 	  ret.insert(C);
1809 	}
1810     }
1811   return ret;
1812 }
1813 
1814 
evaluatePiecewiseLinearFunction(IntegerVector const & x) const1815 int64 PolyhedralFan::evaluatePiecewiseLinearFunction(IntegerVector const &x)const
1816 {
1817   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1818     {
1819       if(i->contains(x))return dotLong(i->getLinearForm(),x);
1820     }
1821   assert(0);
1822   return 0;
1823 }
1824 
1825 
volume(int d,SymmetryGroup * sym) const1826 FieldElement PolyhedralFan::volume(int d, SymmetryGroup *sym)const
1827 {
1828   SymmetryGroup symL(n);
1829   if(!sym)sym=&symL;
1830 
1831   FieldElement ret(Q);
1832 
1833   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1834     {
1835       if(i->dimension()==d)
1836 	{
1837 	  IntegerVector w=stableRay(*i,sym);
1838 	  ret=ret+Q.zHomomorphism(sym->orbitSize(w))*i->volume();
1839 	}
1840     }
1841   return ret;
1842 }
1843 
1844 
isConnected(SymmetryGroup * sym) const1845 bool PolyhedralFan::isConnected(SymmetryGroup *sym)const
1846 {
1847   SymmetryGroup symL(n);
1848   if(!sym)sym=&symL;
1849 
1850   CodimOneConnectednessTester ct;
1851 
1852   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1853     {
1854       log2 cerr<<"Computing ridges of facet." << endl;
1855       PolyhedralFan ridges=facetsOfCone(*i);
1856       IntegerVectorList interiorPoints;
1857       for(PolyhedralConeList::const_iterator j=ridges.cones.begin();j!=ridges.cones.end();j++)
1858 	interiorPoints.push_back(sym->orbitRepresentative(j->getUniquePoint()));
1859       ct.insertFacetOrbit(interiorPoints);
1860     }
1861   return ct.isConnected();
1862 }
1863 
1864 
size() const1865 int PolyhedralFan::size()const
1866 {
1867   return cones.size();
1868 }
1869 
dimensionOfLinealitySpace() const1870 int PolyhedralFan::dimensionOfLinealitySpace()const
1871 {
1872   assert(!cones.empty());
1873   return cones.begin()->dimensionOfLinealitySpace();
1874 }
1875 
negated() const1876 PolyhedralFan PolyhedralFan::negated()const
1877 {
1878   PolyhedralFan ret(n);
1879 
1880   for(coneIterator i=conesBegin();i!=conesEnd();i++)
1881     ret.insert(i->negated());
1882   return ret;
1883 }
1884 
1885 
isCompatible(PolyhedralCone const & c) const1886 bool PolyhedralFan::isCompatible(PolyhedralCone const &c)const
1887 {
1888   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1889     {
1890       PolyhedralCone C=intersection(c,*i);
1891       C.canonicalize();
1892       if(!c.hasFace(C))return false;
1893       if(!i->hasFace(C))return false;
1894     }
1895   return true;
1896 }
1897 
merge(PolyhedralCone const & c)1898 void PolyhedralFan::merge(PolyhedralCone const &c)
1899 {
1900   AsciiPrinter P(Stderr);
1901 
1902   if(isCompatible(c))insert(c);
1903   else
1904     {
1905       assert(0);//Does not work in general.
1906     }
1907   /*
1908   //  P<<"BEFORE MERGE-------------------------" <<*this;
1909 
1910   PolyhedralFan ret=complementOfCone(c,false);
1911   //  P<<"COMPLEMENT OF CONE-------------------------" <<ret;
1912   ret=refinement(*this,ret);
1913 
1914   PolyhedralFan C(c.ambientDimension());
1915   C.insert(c);
1916 
1917   for(PolyhedralConeList::const_iterator i=cones.begin();i!=cones.end();i++)
1918     C=refinement(complementOfCone(*i,true),C);
1919 
1920   for(PolyhedralConeList::const_iterator i=C.cones.begin();i!=C.cones.end();i++)
1921     ret.insert(*i);
1922   *this=ret;
1923 
1924   // P<<"AFTER MERGE" <<*this;
1925   */
1926 }
1927 
1928 
1929 
removeNonMaximal()1930 void PolyhedralFan::removeNonMaximal()
1931 {
1932   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();)
1933     {
1934       IntegerVector w=i->getRelativeInteriorPoint();
1935       bool containedInOther=false;
1936       for(PolyhedralConeList::iterator j=cones.begin();j!=cones.end();j++)
1937 	if(j!=i)
1938 	  {
1939 	    if(j->contains(w)){containedInOther=true;break;}
1940 	  }
1941       if(containedInOther)
1942 	{
1943 	  PolyhedralConeList::iterator k=i;
1944 	  i++;
1945 	  cones.erase(k);
1946 	}
1947       else i++;
1948     }
1949 }
1950 
1951 
moveConesToVector(vector<PolyhedralCone> & v)1952 void PolyhedralFan::moveConesToVector(vector<PolyhedralCone> &v)
1953 {
1954   for(PolyhedralConeList::iterator i=cones.begin();i!=cones.end();)
1955     {
1956       v.push_back(*i);
1957       PolyhedralConeList::iterator k=i;
1958       i++;
1959       cones.erase(k);
1960     }
1961 }
1962 
1963 
1964 
balancingEquations()1965 IntegerMatrix PolyhedralFan::balancingEquations()
1966 {
1967   PolyhedralFan &f=*this;
1968   int n=f.getAmbientDimension();
1969   PolyhedralCone conv=f.convexHull();
1970 //debug<<"convexhull">>conv;
1971 f.removeAllLowerDimensional();
1972 PolyhedralFan skeleton=f.facetComplex();
1973 int d=f.getMaxDimension();
1974 int numberOfFacets=0;
1975 for(PolyhedralFan::coneIterator i=f.conesBegin();i!=f.conesEnd();i++)numberOfFacets+=(i->dimension()==d);
1976 int numberOfRidges=0;
1977 vector<list<pair<IntegerVector,int> > > ridges;
1978 map<IntegerVector,int> ridgeIndices;
1979 vector<IntegerVectorList> ridgesGeneratorOfSpan;
1980 vector<bool> ridgesIgnore;
1981 for(PolyhedralFan::coneIterator i=skeleton.conesBegin();i!=skeleton.conesEnd();i++)
1982   {
1983     ridgesIgnore.push_back(!conv.containsRelatively(i->getRelativeInteriorPoint()));
1984     assert(i->dimension()==d-1);
1985     ridgeIndices[i->getUniquePoint()]=ridges.size();
1986     ridges.push_back(list<pair<IntegerVector,int> >());
1987     ridgesGeneratorOfSpan.push_back(i->generatorsOfSpan());
1988   }
1989 int I=0;
1990 for(PolyhedralFan::coneIterator i=f.conesBegin();i!=f.conesEnd();i++,I++)
1991   {
1992     PolyhedralFan temp=PolyhedralFan::facetsOfCone(*i);
1993     for(PolyhedralFan::coneIterator j=temp.conesBegin();j!=temp.conesEnd();j++)
1994       {
1995         IntegerVector ridgeVector=j->getUniquePoint();
1996         int rIndex=ridgeIndices[ridgeVector];
1997         ridges[rIndex].push_back(pair<IntegerVector,int>(i->link(ridgeVector).semiGroupGeneratorOfRay(),I));
1998       }
1999   }
2000 IntegerMatrix equations(0,numberOfFacets/*+ridges.size()*(d-1)*/);
2001 //      int K=numberOfFacets;
2002 for(int I=0;I<ridges.size();I++)
2003   {
2004     if(ridgesIgnore[I])continue;
2005     IntegerMatrix equationsNewTransposed(numberOfFacets/*+ridges.size()*(d-1)*/,n);
2006     for(list<pair<IntegerVector,int> >::const_iterator j=ridges[I].begin();j!=ridges[I].end();j++)
2007       {
2008         equationsNewTransposed[j->second]=j->first;
2009       }
2010     FieldMatrix tem=integerMatrixToFieldMatrix(rowsToIntegerMatrix(ridgesGeneratorOfSpan[I],n),Q);
2011     IntegerMatrix equationsOfSpan=fieldMatrixToIntegerMatrixPrimitive(tem.reduceAndComputeKernel());
2012 //            for(IntegerVectorList::const_iterator i=ridgesGeneratorOfLinSpace[I].begin();i!=ridgesGeneratorOfLinSpace[I].end();i++)
2013 //              equationsNewTransposed[K++]=*i;
2014 
2015     equations.append(equationsOfSpan*equationsNewTransposed.transposed());
2016   }
2017 //        debug<<equations.getRows();
2018 return equations;
2019 }
2020 
2021 
containsInSupport(IntegerVector const & w)2022 bool PolyhedralFan::containsInSupport(IntegerVector const &w)
2023 {
2024   for(coneIterator i=conesBegin();i!=cones.end();i++)
2025     if(i->contains(w))return true;
2026   return false;
2027 }
2028 
triangulation() const2029 PolyhedralFan PolyhedralFan::triangulation()const
2030 {
2031   PolyhedralFan ret(n);
2032   for(coneIterator i=conesBegin();i!=conesEnd();i++)
2033     {
2034       list<PolyhedralCone> temp=i->triangulation();
2035       for(list<PolyhedralCone>::const_iterator j=temp.begin();j!=temp.end();j++)
2036         {
2037           PolyhedralCone c=*j;
2038           c.canonicalize();
2039           ret.insert(c);
2040         }
2041     }
2042   return ret;
2043 }
2044 
2045