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