1 #include "polyhedralcone.h"
2
3 #include "lp.h"
4 #include "subspace.h"
5 #include "symmetry.h"
6 #include "polymakefile.h"
7 #include "wallideal.h"
8 #include <sstream>
9 #include "triangulation.h"
10 #include "triangulation2.h"
11 #include "linalg.h"
12 #include "linalgfloat.h"
13 #include "continuedfractions.h"
14 #include "polyhedralfan.h"
15 #include "log.h"
16
hasFullRowRank(IntegerVectorList const & l)17 static bool hasFullRowRank(IntegerVectorList const &l)//The tests where this is used can actually be improved by also using the known equations
18 {
19 int s=l.size();
20 if(s==0)return true;
21 if(s>l.front().size())return false;
22 return ::rank(rowsToIntegerMatrix(l))==s;
23 }
24
25 //--------------------------------
26 // PolyhedralCone
27 //--------------------------------
28
compareIntegerLists(IntegerVectorList const & a,IntegerVectorList const & b)29 static bool compareIntegerLists(IntegerVectorList const &a, IntegerVectorList const &b)
30 {
31 assert(a.size()==b.size());
32
33 IntegerVectorList::const_iterator B=b.begin();
34 for(IntegerVectorList::const_iterator A=a.begin();A!=a.end();A++)
35 {
36 if(LexicographicTermOrder()(*A,*B))return true;
37 if(LexicographicTermOrder()(*B,*A))return false;
38 B++;
39 }
40 return false;
41 }
42
operator <(PolyhedralCone const & a,PolyhedralCone const & b)43 bool operator<(PolyhedralCone const &a, PolyhedralCone const &b)
44 {
45 assert(a.state>=3);
46 assert(b.state>=3);
47
48 if(a.dimension()>b.dimension())return true;//INVERTED
49 if(a.dimension()<b.dimension())return false;//INVERTED
50
51 if(a.equations.size()<b.equations.size())return true;
52 if(a.equations.size()>b.equations.size())return false;
53
54 if(a.n<b.n)return true;
55 if(a.n>b.n)return false;
56
57 if(a.inequalities.size()<b.inequalities.size())return true;
58 if(a.inequalities.size()>b.inequalities.size())return false;
59
60 if(compareIntegerLists(a.equations,b.equations))return true;
61 if(compareIntegerLists(b.equations,a.equations))return false;
62
63 if(compareIntegerLists(a.inequalities,b.inequalities))return true;
64 if(compareIntegerLists(b.inequalities,a.inequalities))return false;
65
66 return false;
67 }
68
operator !=(PolyhedralCone const & a,PolyhedralCone const & b)69 bool operator!=(PolyhedralCone const &a, PolyhedralCone const &b)
70 {
71 return (a<b)||(b<a);
72 }
73
74
ensureStateAsMinimum(int s)75 void PolyhedralCone::ensureStateAsMinimum(int s)
76 {
77 if((state<1) && (s==1))
78 {
79 //log0 cerr<<"Number of inequalities: "<<halfSpaces.size()<<" Number of equations: "<<equations.size()<<endl;
80
81 /* if(inequalities.size()==0)
82 {
83 equations=subsetBasis(equations);
84 }
85 else*/
86 {
87 IntegerMatrix temp=rowsToIntegerMatrix(equations,n);
88 FieldMatrix m=integerMatrixToFieldMatrix(temp,Q);
89 m.reduce();
90 m.removeZeroRows();
91
92 IntegerVectorList newInequalities;
93 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
94 {
95 FieldVector w=integerVectorToFieldVector(*i,Q);
96 w=m.canonicalize(w);
97 if(!w.isZero())
98 newInequalities.push_back(w.primitive());
99 }
100
101 inequalities=newInequalities;
102 removeDuplicates(inequalities);
103 equations=fieldMatrixToIntegerMatrixPrimitive(m).getRows();
104 }
105
106 // log0 cerr<<"Number of inequalities: "<<halfSpaces.size()<<" Number of equations: "<<equations.size()<<endl;
107
108 // log0 fprintf(Stderr,"+\n");
109 if(!(preassumptions&PCP_impliedEquationsKnown))
110 if(inequalities.size()>1)//there can be no implied equation unless we have at least two inequalities
111 if(!hasFullRowRank(inequalities))//the inequalities have been moved to an appropriate subspace above, so now we can easily check if the dual cone is simplicial
112 removeRedundantRows(&inequalities,&equations,false);
113 //log0 cerr<<"Number of inequalities: "<<halfSpaces.size()<<" Number of equations: "<<equations.size()<<endl;
114 //log0 fprintf(Stderr,"+\n");
115 }
116 if((state<2) && (s>=2) && !(preassumptions&PCP_facetsKnown))
117 {
118 // halfSpaces.sort();
119 // AsciiPrinter(Stderr).printVectorList(halfSpaces);
120
121 if(inequalities.size()>25)
122 // if(0)
123 {
124 IntegerVectorList h1;
125 IntegerVectorList h2;
126 bool a=false;
127 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
128 {
129 if(a)
130 h1.push_back(*i);
131 else
132 h2.push_back(*i);
133 a=!a;
134 }
135 PolyhedralCone c1(h1,equations);
136 PolyhedralCone c2(h2,equations);
137 c1.ensureStateAsMinimum(2);
138 c2.ensureStateAsMinimum(2);
139 inequalities=c1.inequalities;
140 for(IntegerVectorList::const_iterator i=c2.inequalities.begin();i!=c2.inequalities.end();i++)
141 inequalities.push_back(*i);
142 }
143
144
145 // fprintf(Stderr,"Number half spaces: %i, number of equations: %i\n",halfSpaces.size(),equations.size());
146 if(equations.size())
147 {
148 FieldMatrix equationSpace=integerMatrixToFieldMatrix(rowsToIntegerMatrix(equations,n),Q);
149 equationSpace.reduce();
150 IntegerVectorList halfSpaces2;
151 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
152 {
153 halfSpaces2.push_back(equationSpace.canonicalize(integerVectorToFieldVector(*i,Q)).primitive());
154 }
155 if(hasFullRowRank(halfSpaces2))goto noFallBack;//the inequalities have been moved to an appropriate subspace above, so now we can easily check if the dual cone is simplicial
156 for(IntegerVectorList::const_iterator i=halfSpaces2.begin();i!=halfSpaces2.end();i++)
157 if((i->max()>1000) || (i->min()<-1000))goto fallBack;//more overflows caught in lp_cdd.
158
159 inequalities=fastNormals(halfSpaces2);
160 goto noFallBack;
161 fallBack:
162 removeRedundantRows(&inequalities,&equations,true);
163 noFallBack:;
164 }
165 else
166 inequalities=fastNormals(inequalities);
167 // fprintf(Stderr,"done\n");
168 }
169 if((state<3) && (s>=3))
170 {
171 // fprintf(Stderr,"Number half spaces: %i, number of equations: %i\n",halfSpaces.size(),equations.size());
172 // log1 fprintf(Stderr,"Computing subspace...\n");
173 Subspace v(equations,n);
174
175 equations=v.getRepresentation();
176 // log1 fprintf(Stderr,"...done computing subspace.\n");
177
178
179 for(IntegerVectorList::iterator i=inequalities.begin();i!=inequalities.end();i++)
180 {
181 *i=v.canonicalizeVector(*i);
182 }
183 inequalities.sort();
184 //fprintf(Stderr,"done\n");
185 }
186 state=s;
187 }
188
canonicalize()189 void PolyhedralCone::canonicalize()
190 {
191 ensureStateAsMinimum(3);
192 }
193
194
findFacets()195 void PolyhedralCone::findFacets()
196 {
197 ensureStateAsMinimum(2);
198 }
199
findImpliedEquations()200 void PolyhedralCone::findImpliedEquations()
201 {
202 ensureStateAsMinimum(1);
203 }
204
205 /*PolyhedralCone::PolyhedralCone(IntegerVectorList const &halfSpaces_, IntegerVectorList const &equations_, int ambientDimension, int state):
206 inequalities(halfSpaces_),
207 equations(equations_),
208 state(0),
209 multiplicity(1),
210 n(ambientDimension)
211 {
212 this->state=state;
213 }
214 */
215
PolyhedralCone(int ambientDimension)216 PolyhedralCone::PolyhedralCone(int ambientDimension):
217 n(ambientDimension),
218 state(1),
219 preassumptions(PCP_impliedEquationsKnown|PCP_facetsKnown),
220 multiplicity(1),
221 haveExtremeRaysBeenCached(false),
222 haveGeneratorsOfLinealitySpaceBeenCached(false)
223 {
224 }
225
226
PolyhedralCone(IntegerVectorList const & halfSpaces_,IntegerVectorList const & equations_,int ambientDimension,int preassumptions_)227 PolyhedralCone::PolyhedralCone(IntegerVectorList const &halfSpaces_, IntegerVectorList const &equations_, int ambientDimension, int preassumptions_):
228 //PolyhedralCone::PolyhedralCone(IntegerVectorList const &halfSpaces_, IntegerVectorList const &equations_, int ambientDimension):
229 inequalities(halfSpaces_),
230 equations(equations_),
231 // equations(subsetBasis(equations_)),
232 state(0),
233 preassumptions(preassumptions_),
234 multiplicity(1),
235 haveExtremeRaysBeenCached(false),
236 haveGeneratorsOfLinealitySpaceBeenCached(false)
237 {
238 n=ambientDimension;
239 if(n==-1)
240 {
241 if(!halfSpaces_.empty())
242 n=halfSpaces_.begin()->size();
243 else if(!equations_.empty())
244 n=equations_.begin()->size();
245 else
246 {
247 assert(0);
248 }
249 }
250 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
251 {
252 assert(i->size()==n);
253 }
254 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
255 {
256 /*AsciiPrinter P(Stderr);
257 P.printString("FJDSKFJA\n");
258 P.printVector(*i);
259 */
260 assert(i->size()==n);
261 }
262 ensureStateAsMinimum(1);
263 // computeAndReduceLinearitySpace();
264 }
265
266
getRelativeInteriorPoint() const267 IntegerVector PolyhedralCone::getRelativeInteriorPoint()const
268 {
269 // ensureStateAsMinimum(1);
270 assert(state>=1);
271 IntegerVector ret;
272 IntegerVectorList g=equations;
273 int numberOfEqualities=g.size();
274 g.insert(g.end(),inequalities.begin(),inequalities.end());
275 IntegerVector equalitySet(g.size());
276 for(int i=0;i<numberOfEqualities;i++)equalitySet[i]=1;
277
278 if(!g.empty())
279 ret=relativeInteriorPoint(ambientDimension(),g,&equalitySet);
280 else
281 ret=IntegerVector(n);//cone is the full space, lp code would fail, since the dimension is unknown
282
283 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
284 {
285 assert(dotLong(*i,ret)==0);
286 }
287
288 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
289 {
290 if(!(dotLong(*i,ret)>0))
291 {
292 fprintf(Stderr,"PolyhedralCone::relativeInteriorPoint() : halfSpaces not reduced or mistake in cdd interface!!!\n");
293 }
294 }
295
296 return ret;
297 }
298
299
getHalfSpaces() const300 IntegerVectorList const &PolyhedralCone::getHalfSpaces()const
301 {
302 return inequalities;
303 }
304
305
getEquations() const306 const IntegerVectorList &PolyhedralCone::getEquations()const
307 {
308 assert(state>=1);
309 return equations;
310 }
311
312
getImpliedEquations() const313 const IntegerVectorList &PolyhedralCone::getImpliedEquations()const
314 {
315 const_cast<PolyhedralCone*>(this)->findImpliedEquations();
316 return equations;
317 }
318
319
generatorsOfSpan() const320 IntegerVectorList PolyhedralCone::generatorsOfSpan()const
321 {
322 assert(isInStateMinimum(1));
323 IntegerVectorList empty;
324 PolyhedralCone temp(empty,getEquations(),n);
325 return temp.dualCone().getEquations();
326 }
327
328
generatorsOfLinealitySpace() const329 IntegerVectorList PolyhedralCone::generatorsOfLinealitySpace()const
330 {
331 if(haveGeneratorsOfLinealitySpaceBeenCached)return cachedGeneratorsOfLinealitySpace;
332 IntegerVectorList l=equations;
333 l.insert(l.end(),inequalities.begin(),inequalities.end());
334 FieldMatrix L=integerMatrixToFieldMatrix(rowsToIntegerMatrix(l,ambientDimension()),Q);
335 cachedGeneratorsOfLinealitySpace=fieldMatrixToIntegerMatrixPrimitive(L.reduceAndComputeKernel()).getRows();
336 // return linealitySpace().generatorsOfSpan();
337 haveGeneratorsOfLinealitySpaceBeenCached=true;
338 return cachedGeneratorsOfLinealitySpace;
339 }
340
341
ambientDimension() const342 int PolyhedralCone::ambientDimension()const
343 {
344 return n;
345 }
346
347
codimension() const348 int PolyhedralCone::codimension()const
349 {
350 return ambientDimension()-dimension();
351 // return getEquations().size();
352 }
353
354
dimension() const355 int PolyhedralCone::dimension()const
356 {
357 assert(state>=1);
358 // ensureStateAsMinimum(1);
359 return ambientDimension()-equations.size();
360 }
361
362
isZero() const363 bool PolyhedralCone::isZero()const
364 {
365 return dimension()==0;
366 }
367
368
isFullSpace() const369 bool PolyhedralCone::isFullSpace()const
370 {
371 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
372 if(!i->isZero())return false;
373 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
374 if(!i->isZero())return false;
375 return true;
376 }
377
378
intersection(const PolyhedralCone & a,const PolyhedralCone & b)379 PolyhedralCone intersection(const PolyhedralCone &a, const PolyhedralCone &b)
380 {
381 if(a.ambientDimension()!=b.ambientDimension())
382 {
383 debug<<a.ambientDimension()<<"=="<<b.ambientDimension();
384 assert(a.ambientDimension()==b.ambientDimension());
385 }
386 IntegerVectorList inequalities=a.inequalities;
387 inequalities.insert(inequalities.end(),b.inequalities.begin(),b.inequalities.end());
388 IntegerVectorList equations=a.equations;
389
390
391 equations.insert(equations.end(),b.equations.begin(),b.equations.end());
392 {
393 removeDuplicates(equations);
394 removeDuplicates(inequalities);
395 IntegerVectorList Aequations=a.equations;
396 IntegerVectorList Ainequalities=a.inequalities;
397 removeDuplicates(Aequations);
398 removeDuplicates(Ainequalities);
399 if((Ainequalities.size()==inequalities.size()) && (Aequations.size()==equations.size()))return a;
400 IntegerVectorList Bequations=b.equations;
401 IntegerVectorList Binequalities=b.inequalities;
402 removeDuplicates(Bequations);
403 removeDuplicates(Binequalities);
404 if((Binequalities.size()==inequalities.size()) && (Bequations.size()==equations.size()))return b;
405 }
406
407 return PolyhedralCone(inequalities,equations,a.ambientDimension());
408 }
409
product(const PolyhedralCone & a,const PolyhedralCone & b)410 PolyhedralCone product(const PolyhedralCone &a, const PolyhedralCone &b)
411 {
412 IntegerVectorList equations2;
413 IntegerVectorList inequalities2;
414
415 int n1=a.n;
416 int n2=b.n;
417
418 for(IntegerVectorList::const_iterator i=a.equations.begin();i!=a.equations.end();i++)
419 equations2.push_back(concatenation(*i,IntegerVector(n2)));
420 for(IntegerVectorList::const_iterator i=b.equations.begin();i!=b.equations.end();i++)
421 equations2.push_back(concatenation(IntegerVector(n1),*i));
422 for(IntegerVectorList::const_iterator i=a.inequalities.begin();i!=a.inequalities.end();i++)
423 inequalities2.push_back(concatenation(*i,IntegerVector(n2)));
424 for(IntegerVectorList::const_iterator i=b.inequalities.begin();i!=b.inequalities.end();i++)
425 inequalities2.push_back(concatenation(IntegerVector(n1),*i));
426
427 PolyhedralCone ret(inequalities2,equations2,n1+n2);
428 ret.setMultiplicity(a.getMultiplicity()*b.getMultiplicity());
429 ret.setLinearForm(concatenation(a.getLinearForm(),b.getLinearForm()));
430
431 ret.ensureStateAsMinimum(a.state);
432 ret.ensureStateAsMinimum(b.state);
433
434 return ret;
435 }
436
positiveOrthant(int dimension)437 PolyhedralCone PolyhedralCone::positiveOrthant(int dimension)
438 {
439 IntegerVectorList halfSpaces;
440
441 for(int i=0;i<dimension;i++)halfSpaces.push_back(IntegerVector::standardVector(dimension,i));
442
443 IntegerVectorList empty;
444 return PolyhedralCone(halfSpaces,empty,dimension);
445 }
446
447
isInStateMinimum(int s) const448 bool PolyhedralCone::isInStateMinimum(int s)const
449 {
450 return state>=s;
451 }
452
getState() const453 int PolyhedralCone::getState()const
454 {
455 return state;
456 }
457
458
print(class Printer * p,bool xml) const459 void PolyhedralCone::print(class Printer *p, bool xml)const
460 {
461 if(0)
462 {
463 p->printString("Printing PolyhedralCone");
464 p->printNewLine();
465 p->printString("Ambient dimension: ");
466 p->printInteger(n);
467 p->printNewLine();
468 if(isInStateMinimum(1))
469 {
470 p->printString("Dimension: ");
471 p->printInteger(dimension());
472 p->printNewLine();
473 }
474 p->printString("Linearity space:");
475 // p->printNewLine();
476 p->printVectorList(equations);
477 p->printString("Inequalities:");
478 p->printVectorList(inequalities);
479 p->printString("Relative interior point:\n");
480 p->printVector(getRelativeInteriorPoint());
481 p->printNewLine();
482 p->printString("Done printing PolyhedralCone.");
483 p->printNewLine();
484 }
485 else
486 {
487 PolymakeFile polymakeFile;
488 polymakeFile.create("NONAME","PolyhedralCone","PolyhedralCone",xml);
489 polymakeFile.writeCardinalProperty("AMBIENT_DIM",n);
490 if(isInStateMinimum(1))
491 {
492 polymakeFile.writeCardinalProperty("DIM",dimension());
493 //need to check that the following is done correctly
494 // polymakeFile.writeCardinalProperty("LINEALITY_DIM",dimensionOfLinealitySpace());
495
496 polymakeFile.writeMatrixProperty("IMPLIED_EQUATIONS",rowsToIntegerMatrix(equations,n));
497 }
498 polymakeFile.writeCardinalProperty("LINEALITY_DIM",dimensionOfLinealitySpace());
499 polymakeFile.writeMatrixProperty("LINEALITY_SPACE",rowsToIntegerMatrix(linealitySpace().dualCone().getEquations(),n));
500
501
502 if(isInStateMinimum(2))
503 polymakeFile.writeMatrixProperty("FACETS",rowsToIntegerMatrix(inequalities,n));
504 else
505 polymakeFile.writeMatrixProperty("INEQUALITIES",rowsToIntegerMatrix(inequalities,n));
506
507 polymakeFile.writeCardinalVectorProperty("RELATIVE_INTERIOR_POINT",getRelativeInteriorPoint());
508
509
510 stringstream s;
511 polymakeFile.writeStream(s);
512 string S=s.str();
513 p->printString(S.c_str());
514 }
515 }
516
517
printAsFan(class Printer * p) const518 void PolyhedralCone::printAsFan(class Printer *p)const
519 {
520 IntegerVectorList empty;
521 PolyhedralCone C2=*this;
522 C2.canonicalize();
523 PolyhedralFan F(n);F.insert(C2);
524 F.printWithIndices(p,FPF_default);
525 }
526
527
dehomogenize(IntegerVector const & v)528 static IntegerVector dehomogenize(IntegerVector const &v)
529 {
530 assert(v.size()>0);
531
532 IntegerVector ret(v.size()-1);
533
534 for(int i=0;i<v.size()-1;i++)ret[i]=v[i];
535 return ret;
536 }
537
dehomogenize(IntegerVectorList const & l)538 static IntegerVectorList dehomogenize(IntegerVectorList const &l)
539 {
540 IntegerVectorList ret;
541
542 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
543 {
544 ret.push_back(dehomogenize(*i));
545 }
546 return ret;
547 }
548
withLastCoordinateRemoved() const549 PolyhedralCone PolyhedralCone::withLastCoordinateRemoved()const
550 {
551 assert(n>0);
552
553 return PolyhedralCone(dehomogenize(inequalities),dehomogenize(equations));
554 }
555
containsPositiveVector() const556 bool PolyhedralCone::containsPositiveVector()const
557 {
558 PolyhedralCone temp=intersection(*this,PolyhedralCone::positiveOrthant(n));
559 IntegerVector v=temp.getRelativeInteriorPoint();
560 return v.isPositive();
561 }
562
563
dimensionOfLinealitySpace() const564 int PolyhedralCone::dimensionOfLinealitySpace()const
565 {
566 if(inequalities.empty())return dimension();
567 IntegerVectorList a;
568 PolyhedralCone temp(a,inequalities);
569 temp=intersection(temp,*this);
570
571 return temp.dimension();
572 }
573
574
contains(IntegerVector const & v) const575 bool PolyhedralCone::contains(IntegerVector const &v)const
576 {
577 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
578 {
579 if(dotLong(*i,v)!=0)return false;
580 }
581
582 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
583 {
584 if(dotLong(*i,v)<0)return false;
585 }
586
587 return true;
588 }
589
590
contains(IntegerVectorList const & l) const591 bool PolyhedralCone::contains(IntegerVectorList const &l)const
592 {
593 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
594 if(!contains(*i))return false;
595 return true;
596 }
597
contains(PolyhedralCone const & c) const598 bool PolyhedralCone::contains(PolyhedralCone const &c)const
599 {
600 PolyhedralCone c2=intersection(*this,c);
601 PolyhedralCone c3=c;
602 c2.canonicalize();
603 c3.canonicalize();
604 return !(c2!=c3);
605 }
606
607
containsPerturbed(IntegerVectorList const & l) const608 bool PolyhedralCone::containsPerturbed(IntegerVectorList const &l)const
609 {
610 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
611 {
612 for(IntegerVectorList::const_iterator j=l.begin();j!=l.end();j++)
613 if(dotLong(*i,*j)!=0)return false;
614 }
615
616 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
617 {
618 for(IntegerVectorList::const_iterator j=l.begin();j!=l.end();j++)
619 if(dotLong(*i,*j)<0)return false;
620 else if(dotLong(*i,*j)>0) break;
621 }
622 return true;
623 }
624
containsRelatively(IntegerVector const & v) const625 bool PolyhedralCone::containsRelatively(IntegerVector const &v)const
626 {
627 assert(state>=1);
628 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
629 {
630 if(dotLong(*i,v)!=0)return false;
631 }
632
633 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
634 {
635 if(dotLong(*i,v)<=0)return false;
636 }
637
638 return true;
639 }
640
641
permuted(IntegerVector const & v) const642 PolyhedralCone PolyhedralCone::permuted(IntegerVector const &v)const
643 {
644 PolyhedralCone ret(SymmetryGroup::permuteIntegerVectorList(inequalities,v),SymmetryGroup::permuteIntegerVectorList(equations,v),n,(areImpliedEquationsKnown()?PCP_impliedEquationsKnown:0)+(areFacetsKnown()?PCP_facetsKnown:0));
645 if(state>=1)ret.state=1;
646 if(state>=2)ret.state=2;
647
648 ret.ensureStateAsMinimum(state);
649 ret.setMultiplicity(multiplicity);
650 return ret;
651 }
652
getUniquePoint() const653 IntegerVector PolyhedralCone::getUniquePoint()const
654 {
655 IntegerVectorList rays=extremeRays();
656 IntegerVector ret(n);
657 for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
658 ret+=*i;
659
660 assert(containsRelatively(ret));//remove this check later
661 return ret;
662 }
663
getUniquePointFromExtremeRays(IntegerVectorList const & extremeRays) const664 IntegerVector PolyhedralCone::getUniquePointFromExtremeRays(IntegerVectorList const &extremeRays)const
665 {
666 IntegerVector ret(n);
667 for(IntegerVectorList::const_iterator i=extremeRays.begin();i!=extremeRays.end();i++)
668 if(contains(*i))ret+=*i;
669 return ret;
670 }
671
672 /**
673 * Returns a primitive vector parallel to the projection of v onto the orthogonal complement of E.
674 */
primitiveProjection(IntegerVectorList const & E,IntegerVector & v,bool useFloat)675 static IntegerVector primitiveProjection(IntegerVectorList const &E, IntegerVector &v, bool useFloat)
676 {
677 int n=v.size();
678 if(useFloat)
679 {
680 IntegerMatrix E2=rowsToIntegerMatrix(E,n);
681 linalgfloat::Matrix E3(E.size(),n);
682 for(int i=0;i<E3.getHeight();i++)
683 for(int j=0;j<E3.getWidth();j++)
684 E3[i][j]=E2[i][j];
685 cerr<<E3;
686 E3.orthogonalize();
687 linalgfloat::Vector v2(n);
688 for(int i=0;i<n;i++)v2[i]=v[i];
689 cerr<<E3;
690 linalgfloat::Vector coefficients=E3.projectionCoefficients(v2);
691
692 cerr<<coefficients<<"\n";
693
694 vector<double> temp(coefficients.size());
695 for(int i=0;i<coefficients.size();i++)temp[i]=coefficients[i];
696 vector<int> tempInt(coefficients.size());
697 int denominator;
698 doubleVectorToFractions(temp, tempInt, denominator);
699 IntegerVector tempInt2(tempInt.size());for(int i=0;i<tempInt.size();i++)tempInt2[i]=tempInt[i];
700 if(tempInt2.max()>1000)goto fallback;
701 if(tempInt2.min()<-1000)goto fallback;
702 if(denominator>1000)goto fallback;
703 if(denominator<-1000)goto fallback;
704 IntegerVector ret=denominator*v;
705 int I=0;
706 for(IntegerVectorList::const_iterator i=E.begin();i!=E.end();i++,I++)ret-=coefficients[I]*(*i);
707 if(!E2.inKernel(ret))goto fallback;
708 return normalized(ret);
709 }
710 fallback:
711 debug<<"FLOATCONELINALG FALLBACK\n";
712
713 IntegerVectorList inequalities;
714 inequalities.push_back(v);
715 FieldMatrix linealitySpaceOrth=combineOnTop(integerMatrixToFieldMatrix(rowsToIntegerMatrix(E,n),Q),integerMatrixToFieldMatrix(rowsToIntegerMatrix(inequalities,n),Q));
716 FieldMatrix temp=combineOnTop(linealitySpaceOrth.reduceAndComputeKernel(),integerMatrixToFieldMatrix(rowsToIntegerMatrix(E,n),Q));
717 FieldMatrix temp2=temp.reduceAndComputeKernel();
718 assert(temp2.getHeight()==1);
719 return temp2[0].primitive();
720 }
721
extremeRays(IntegerVectorList const * generatorsOfLinealitySpace) const722 IntegerVectorList PolyhedralCone::extremeRays(IntegerVectorList const *generatorsOfLinealitySpace)const
723 {
724 assert((dimension()==ambientDimension()) || (state>=3));
725
726 if(haveExtremeRaysBeenCached)return cachedExtremeRays;
727
728 //The generators of the lineality space are not known, we may as well compute them. Since we typycally will use them at least once.
729 IntegerVectorList generatorsOfLinealitySpace2;
730 if(!generatorsOfLinealitySpace)
731 {
732 FieldMatrix linealitySpaceOrth=combineOnTop(integerMatrixToFieldMatrix(rowsToIntegerMatrix(this->equations,n),Q),integerMatrixToFieldMatrix(rowsToIntegerMatrix(this->inequalities,n),Q));
733 generatorsOfLinealitySpace2=fieldMatrixToIntegerMatrixPrimitive(linealitySpaceOrth.reduceAndComputeKernel()).getRows();
734 generatorsOfLinealitySpace=&generatorsOfLinealitySpace2;
735 }
736 // this->print(&debug);
737 /* This code actually works even for lower dimensional cones if they have been canonicalized. */
738
739 // AsciiPrinter(Stderr).printVectorList(halfSpaces);
740 // AsciiPrinter(Stderr).printVectorList(equations);
741
742 IntegerVectorList ret;
743
744 // log0 fprintf(Stderr,"calling double description (cddlib)\n");
745 IntegerVectorList indices=extremeRaysInequalityIndices(inequalities);
746 // log0 fprintf(Stderr,"returning\n");
747 // log0 fprintf(Stderr,"computing relative interior points\n");
748
749 // debug<<"INDICES"<<indices;
750 for(IntegerVectorList::const_iterator i=indices.begin();i!=indices.end();i++)
751 {
752 // log0 AsciiPrinter(Stderr)<<*i;
753 if(1)
754 {
755 /* At this point we know lineality space, implied equations and
756 also inequalities for the ray. To construct a vector on the
757 ray which is stable under (or indendent of) angle and
758 linarity preserving transformation we find the dimension 1
759 subspace orthorgonal to the implied equations and the
760 lineality space and pick a suitable primitive generator */
761
762 /* To be more precise,
763 * let E be the set of equations, and v the inequality defining a ray R.
764 * We wish to find a vector satisfying these, but it must also be orthogonal
765 * to the lineality space of the cone, that is, in the span of {E,v}.
766 * One way to get such a vector is to project v to E an get a vector p.
767 * Then v-p is in the span of {E,v} by construction.
768 * The vector v-p is also in the orthogonal complement to E by construction,
769 * that is, the span of R.
770 * We wish to argue that it is not zero.
771 * That would imply that v=p, meaning that v is in the span of the equations.
772 * However, that would contradict that R is a ray.
773 * In case v-p does not satisfy the inequality v (is this possible?), we change the sign.
774 *
775 * As a consequence we need the following procedure
776 * primitiveProjection():
777 * Input: E,v
778 * Output: A primitive representation of the vector v-p, where p is the projection of v onto E
779 *
780 * Notice that the output is a Q linear combination of the input and that p is
781 * a linear combination of E. The check that p has been computed correctly,
782 * it suffices to check that v-p satisfies the equations E.
783 * The routine will actually first compute a multiple of v-p.
784 * It will do this using floating point arithmetics. It will then transform
785 * the coefficients to get the multiple of v-p into integers. Then it
786 * verifies in exact arithmetics, that with these coefficients we get a point
787 * satisfying E. It then returns the primitive vector on the ray v-p.
788 * In case of a failure it falls back to an implementation using rational arithmetics.
789 */
790
791
792 IntegerVector asVector(inequalities.size());
793 // log0 AsciiPrinter(Stderr)<<asVector;
794 for(int j=0;j<i->size();j++)asVector[(*i)[j]]=1;
795 // log0 AsciiPrinter(Stderr)<<asVector;
796
797 IntegerVectorList equations=this->equations;
798 IntegerVectorList inequalities;
799
800 IntegerVector theInequality;
801
802 // log0 AsciiPrinter(Stderr)<<inequalities;
803
804 IntegerVectorList::const_iterator a=this->inequalities.begin();
805 for(int j=0;j<asVector.size();j++,a++)
806 if(asVector[j])
807 equations.push_back(*a);
808 else
809 theInequality=*a;
810
811 assert(!theInequality.isZero());
812
813 IntegerVector thePrimitiveVector;
814 if(generatorsOfLinealitySpace)
815 {
816 IntegerMatrix temp=rowsToIntegerMatrix(equations,n);
817 temp.append(rowsToIntegerMatrix(*generatorsOfLinealitySpace,n));
818
819 // debug<<*generatorsOfLinealitySpace;
820 // debug.printVectorList(temp.getRows());
821 log3 debug<<temp.getRows();
822 thePrimitiveVector=vectorInKernel(temp);
823 }
824 else
825 {
826 //log0 AsciiPrinter(Stderr)<<equations;
827 /* {
828 IntegerVectorList equations2=this->equations;
829 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)equations2.push_back(*i);
830 debug<<primitiveProjection(equations2,theInequality,true)<<"\n";
831 debug<<primitiveProjection(equations2,theInequality,false)<<"\n";
832 }
833 */
834 /** TODO: These calls are slow, but used often in the symmetric traversal. Maybe they should be done in floating point somehow.*/
835 FieldMatrix linealitySpaceOrth=combineOnTop(integerMatrixToFieldMatrix(rowsToIntegerMatrix(this->equations,n),Q),integerMatrixToFieldMatrix(rowsToIntegerMatrix(this->inequalities,n),Q));
836 FieldMatrix temp=combineOnTop(linealitySpaceOrth.reduceAndComputeKernel(),integerMatrixToFieldMatrix(rowsToIntegerMatrix(equations,n),Q));
837 FieldMatrix temp2=temp.reduceAndComputeKernel();
838
839 assert(temp2.getHeight()==1);
840 thePrimitiveVector=temp2[0].primitive();
841
842 // debug<<thePrimitiveVector<<"\n";
843 }
844 if(!contains(thePrimitiveVector))thePrimitiveVector=-thePrimitiveVector;
845 ret.push_back(thePrimitiveVector);
846 }
847 else
848 {
849 /* IntegerVectorList equations;
850 for(int j=0;j<i->size();j++)
851 {
852 IntegerVectorList::const_iterator a=halfSpaces.begin();
853 for(int k=0;k<(*i)[j];k++)
854 {
855 assert(a!=halfSpaces.end());
856 a++;
857 }
858 assert(a!=halfSpaces.end());
859 equations.push_back(*a);
860 }*/
861
862 IntegerVector asVector(inequalities.size());
863 for(int j=0;j<i->size();j++)asVector[(*i)[j]]=1;
864
865 IntegerVectorList equations=this->equations;
866 IntegerVectorList inequalities;
867
868 IntegerVectorList::const_iterator a=inequalities.begin();
869 for(int j=0;j<asVector.size();j++,a++)
870 if(asVector[j])
871 equations.push_back(*a);
872 else
873 inequalities.push_back(*a);
874
875
876 //log0 fprintf(Stderr,"Equations %i, HalfSpaces: %i\n",equations.size(),inequalities.size());
877 IntegerVector u=PolyhedralCone(inequalities,equations).getRelativeInteriorPoint();
878 if(!u.isZero())
879 ret.push_back(u);
880 else
881 {
882 log2 fprintf(Stderr,"Remember to fix cdd double description interface\n");
883 }
884 }
885 }
886 // log0 fprintf(Stderr,"done computing relative interior points\n");
887
888
889 /* //triangulation method. Keep this code.
890 {
891 IntegerMatrix temp=rowsToIntegerMatrix(halfSpaces);
892 log0 fprintf(Stderr,"calling double description (triangulation)\n");
893 IntegerVectorList ret2=Triangulation::normals(temp);
894 log0 fprintf(Stderr,"returning\n");
895 return ret2;
896
897 AsciiPrinter(Stderr).printVectorList(halfSpaces);
898 AsciiPrinter(Stderr).printVectorList(equations);
899 fprintf(Stderr,"dim:%i\n",dimension());
900 //ret.sort();
901 AsciiPrinter(Stderr).printVectorList(ret);
902 // AsciiPrinter(Stderr).printVectorList(fastNormals(ret2));
903 AsciiPrinter(Stderr).printVectorList(ret2);
904
905 fprintf(stderr,"-----------------------\n");
906 }
907 */
908
909 cachedExtremeRays=ret;
910 haveExtremeRaysBeenCached=true;
911
912 return ret;
913 }
914
915
isSimplicial() const916 bool PolyhedralCone::isSimplicial()const
917 {
918 assert(state>=2);
919
920 // ensureStateAsMinimum(2);
921 // AsciiPrinter P(Stderr);
922 // print(&P);
923 return codimension()+getHalfSpaces().size()+dimensionOfLinealitySpace()==n;
924 }
925
926
checkDual(PolyhedralCone const & c) const927 bool PolyhedralCone::checkDual(PolyhedralCone const &c)const
928 {
929 assert(dimensionOfLinealitySpace()+c.dimension()==ambientDimension());
930
931 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
932 {
933 assert(c.contains(*i));
934 }
935 for(IntegerVectorList::const_iterator i=equations.begin();i!=equations.end();i++)
936 {
937 assert(c.contains(*i));
938 }
939 return true;
940 }
941
942
dualCone() const943 PolyhedralCone PolyhedralCone::dualCone()const
944 {
945 assert(state>=1);
946
947 IntegerVectorList dualInequalities,dualEquations;
948
949 dual(ambientDimension(),inequalities,equations,&dualInequalities,&dualEquations);
950
951 PolyhedralCone ret(dualInequalities,dualEquations);
952
953 ret.ensureStateAsMinimum(state);
954 // ret.canonicalize();
955
956
957 assert(checkDual(ret));
958 assert(ret.checkDual(*this));
959
960 return ret;
961 }
962
963
negated() const964 PolyhedralCone PolyhedralCone::negated()const
965 {
966 IntegerVectorList inequalities2;
967 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)inequalities2.push_back(-*i);
968 // PolyhedralCone ret(inequalities2,equations,n);
969 PolyhedralCone ret(inequalities2,equations,n,(areFacetsKnown()?PCP_facetsKnown:0)|(areImpliedEquationsKnown()?PCP_impliedEquationsKnown:0));
970 ret.ensureStateAsMinimum(state);
971 return ret;
972 }
973
linealitySpace() const974 PolyhedralCone PolyhedralCone::linealitySpace()const
975 {
976 IntegerVectorList l1=getEquations();
977 IntegerVectorList l2=getHalfSpaces();
978
979 l1.splice(l1.begin(),l2);
980
981 IntegerVectorList temp;
982 PolyhedralCone ret(temp,l1,n);
983 ret.ensureStateAsMinimum(state);
984 return ret;
985 }
986
987
span() const988 PolyhedralCone PolyhedralCone::span()const
989 {
990 return PolyhedralCone(IntegerVectorList(),getImpliedEquations(),n);
991 }
992
993
getMultiplicity() const994 int PolyhedralCone::getMultiplicity()const
995 {
996 return multiplicity;
997 }
998
999
setMultiplicity(int m)1000 void PolyhedralCone::setMultiplicity(int m)
1001 {
1002 multiplicity=m;
1003 }
1004
1005
quotientLatticeBasis() const1006 IntegerVectorList PolyhedralCone::quotientLatticeBasis()const
1007 {
1008 assert(isInStateMinimum(1));// Implied equations must have been computed in order to know the span of the cone
1009
1010 return basisOfQuotientLattice(equations,inequalities,n);
1011 }
1012
1013
semiGroupGeneratorOfRay() const1014 IntegerVector PolyhedralCone::semiGroupGeneratorOfRay()const
1015 {
1016 IntegerVectorList temp=quotientLatticeBasis();
1017 assert(temp.size()==1);
1018 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
1019 if(dotLong(temp.front(),*i)<0)
1020 {
1021 temp.front()=-temp.front();
1022 break;
1023 }
1024 return temp.front();
1025 }
1026
1027
getLinearForm() const1028 IntegerVector const &PolyhedralCone::getLinearForm()const
1029 {
1030 return linearForm;
1031 }
1032
1033
setLinearForm(IntegerVector const & linearForm_)1034 void PolyhedralCone::setLinearForm(IntegerVector const &linearForm_)
1035 {
1036 linearForm=linearForm_;
1037 }
1038
1039
link(IntegerVector const & w) const1040 PolyhedralCone PolyhedralCone::link(IntegerVector const &w)const
1041 {
1042 /*
1043 * Observe that the inequalities giving rise to facets
1044 * also give facets in the link, if they are kept as
1045 * inequalities. This means that the state cannot decrease when taking links.
1046 *
1047 */
1048 // assert(state>=3);
1049 IntegerVectorList inequalities2;
1050 for(IntegerVectorList::const_iterator j=inequalities.begin();j!=inequalities.end();j++)
1051 if(dotLong(w,*j)==0)inequalities2.push_back(*j);
1052 // PolyhedralCone C(inequalities2,getEquations(),n);
1053 // C.canonicalize();
1054 // PolyhedralCone C(inequalities2,getEquations(),n,state);//STATE-----------------------------------------------------
1055 PolyhedralCone C(inequalities2,getEquations(),n,(areImpliedEquationsKnown()?PCP_impliedEquationsKnown:0)|(areFacetsKnown()?PCP_facetsKnown:0));
1056 C.ensureStateAsMinimum(state);
1057
1058 C.setLinearForm(getLinearForm());
1059 C.setMultiplicity(getMultiplicity());
1060
1061 return C;
1062 }
1063
1064
givenByRays(IntegerVectorList const & generators,IntegerVectorList const & linealitySpace,int n)1065 PolyhedralCone PolyhedralCone::givenByRays(IntegerVectorList const &generators, IntegerVectorList const &linealitySpace, int n)
1066 {
1067 //rewrite modulo lineality space
1068 IntegerVectorList newGenerators=generators;
1069 {
1070 Subspace l(linealitySpace,n);
1071 for(IntegerVectorList::const_iterator i=generators.begin();i!=generators.end();i++)
1072 newGenerators.push_back(l.canonicalizeVector(*i));
1073 }
1074
1075 PolyhedralCone dual(newGenerators,linealitySpace,n);
1076 dual.findFacets();
1077 dual.canonicalize();
1078 IntegerVectorList inequalities=dual.extremeRays();
1079
1080 IntegerVectorList span=generators;
1081 for(IntegerVectorList::const_iterator i=linealitySpace.begin();i!=linealitySpace.end();i++)span.push_back(*i);
1082 FieldMatrix m2Q=integerMatrixToFieldMatrix(rowsToIntegerMatrix(span,n),Q);
1083 IntegerVectorList equations=fieldMatrixToIntegerMatrixPrimitive(m2Q.reduceAndComputeKernel()).getRows();
1084
1085 return PolyhedralCone(inequalities,equations,n);
1086 }
1087
1088
volume() const1089 FieldElement PolyhedralCone::volume()const
1090 {
1091 AsciiPrinter P(Stderr);
1092
1093 PolyhedralCone Ctemp=intersection(*this,this->linealitySpace().dualCone());
1094
1095 Ctemp.canonicalize();
1096 cerr<<"testestests";
1097 IntegerVectorList eq2=rowsToIntegerMatrix(Ctemp.equations,n).transposed().getRows();
1098 eq2.push_front(IntegerVector(Ctemp.equations.size()));
1099 eq2=rowsToIntegerMatrix(eq2).transposed().getRows();
1100 cerr<<"testestests2";
1101
1102 IntegerVectorList in2=rowsToIntegerMatrix(Ctemp.inequalities,n).transposed().getRows();
1103 in2.push_front(IntegerVector(Ctemp.inequalities.size()));
1104 in2=rowsToIntegerMatrix(in2).transposed().getRows();
1105
1106 for(int i=0;i<n;i++)
1107 {
1108 in2.push_back(IntegerVector::standardVector(n+1,i+1)+IntegerVector::standardVector(n+1,0));
1109 in2.push_back(-IntegerVector::standardVector(n+1,i+1)+IntegerVector::standardVector(n+1,0));
1110 }
1111
1112 debug<<in2;
1113 debug<<eq2;
1114
1115 PolyhedralCone lifted(in2,eq2,n+1);
1116
1117
1118 lifted.print(&debug);
1119
1120 lifted.canonicalize();
1121 IntegerMatrix A=rowsToIntegerMatrix(lifted.extremeRays()).transposed();//points are columns
1122
1123 cerr << "height " << A.getHeight() << " width" <<A.getWidth()<<endl;
1124
1125 P<<A.transposed().getRows();
1126
1127 FieldMatrix A2(Q,A.getHeight()-1,A.getWidth());
1128 for(int i=0;i<A.getHeight()-1;i++)
1129 {
1130 A2[i]=integerVectorToFieldVector(A[i+1],Q)/integerVectorToFieldVector(A[0],Q);
1131 }
1132 A2=A2.transposed();//Now points are rows
1133
1134 P<<"Triangulating\n";
1135 list<Triangulation::Cone> T=Triangulation::triangulate(A.transposed(),true);//revlex
1136 P<<"Done triangulating\n";
1137
1138
1139 P<<(int)T.size();
1140
1141 FieldElement ret(Q);
1142
1143 for(list<Triangulation::Cone>::const_iterator i=T.begin();i!=T.end();i++)
1144 {
1145 FieldMatrix S=A2.submatrixRows(*i);
1146
1147 debug<<"M"<<S<<"M";
1148
1149
1150 FieldMatrix S1(Q,S.getHeight()-1,S.getWidth());
1151 for(int j=0;j<S1.getHeight();j++)
1152 S1[j]=S[j+1]-S[0];
1153
1154 FieldMatrix S2=S1*(S1.transposed());
1155 FieldElement square=S2.reduceAndComputeDeterminant();
1156
1157 // ret+=square.squareroot();
1158 }
1159
1160 return ret;
1161 }
1162
1163
1164 #include "halfopencone.h"
hasFace(PolyhedralCone const & f) const1165 bool PolyhedralCone::hasFace(PolyhedralCone const &f)const
1166 {
1167 if(!contains(f))return false;
1168 IntegerVectorList linealitySpace=this->linealitySpace().dualCone().getEquations();
1169 IntegerVectorList rays=extremeRays();
1170
1171 for(IntegerVectorList::const_iterator i=linealitySpace.begin();i!=linealitySpace.end();i++)
1172 if(!f.contains(*i))return false;
1173
1174 IntegerVectorList strict;
1175 for(IntegerVectorList::const_iterator i=rays.begin();i!=rays.end();i++)
1176 if(!f.contains(*i))
1177 strict.push_back(*i);
1178
1179 IntegerVectorList linealitySpace2=f.linealitySpace().dualCone().getEquations();
1180 IntegerVectorList rays2=f.extremeRays();
1181
1182 for(IntegerVectorList::const_iterator i=rays2.begin();i!=rays2.end();i++)
1183 {
1184 linealitySpace2.push_back(*i);
1185 }
1186 IntegerVectorList empty;
1187 HalfOpenCone C(n,linealitySpace2, empty, strict);
1188 return !C.isEmpty();
1189 }
1190
faceContaining(IntegerVector const & v) const1191 PolyhedralCone PolyhedralCone::faceContaining(IntegerVector const &v)const
1192 {
1193 assert(n==v.size());
1194 if(!contains(v))
1195 {
1196 debug<<"Not contained in cone:\n"<<v<<"\n";
1197 assert(contains(v));
1198 }
1199 IntegerVectorList newEquations=equations;
1200 IntegerVectorList newInequalities;
1201 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
1202 if(dotLong(*i,v))
1203 newInequalities.push_back(*i);
1204 else
1205 newEquations.push_back(*i);
1206
1207 PolyhedralCone ret(newInequalities,newEquations,n,(state>=1)?PCP_impliedEquationsKnown:0);
1208 ret.ensureStateAsMinimum(state);
1209 return ret;
1210 }
1211
1212
faceContainingPerturbed(IntegerVectorList const & l) const1213 PolyhedralCone PolyhedralCone::faceContainingPerturbed(IntegerVectorList const &l)const
1214 {
1215 IntegerVectorList newEquations=equations;
1216 IntegerVectorList newInequalities;
1217 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
1218 {
1219 bool onFacet=true;
1220 for(IntegerVectorList::const_iterator j=l.begin();j!=l.end();j++)
1221 if(dotLong(*i,*j)){onFacet=false;break;}
1222 if(!onFacet)
1223 newInequalities.push_back(*i);
1224 else
1225 newEquations.push_back(*i);
1226 }
1227 PolyhedralCone ret(newInequalities,newEquations,n,(state>=1)?PCP_impliedEquationsKnown:0);
1228 ret.ensureStateAsMinimum(state);
1229 return ret;
1230 }
1231
1232
projection(int newn) const1233 PolyhedralCone PolyhedralCone::projection(int newn)const
1234 {
1235 assert(newn<=n);
1236 assert(newn>=0);
1237 IntegerVectorList rays=extremeRays();
1238 IntegerVectorList lines=linealitySpace().generatorsOfSpan();
1239 rays=rowsToIntegerMatrix(rays,n).submatrix(0,0,rays.size(),newn).getRows();
1240 lines=rowsToIntegerMatrix(lines,n).submatrix(0,0,lines.size(),newn).getRows();
1241 return givenByRays(rays,lines,newn);
1242 }
1243
sum(PolyhedralCone const & A,PolyhedralCone const & B)1244 PolyhedralCone sum(PolyhedralCone const &A, PolyhedralCone const &B)
1245 {
1246 IntegerVectorList raysA=A.extremeRays();
1247 IntegerVectorList raysB=B.extremeRays();
1248 for(IntegerVectorList::const_iterator i=raysB.begin();i!=raysB.end();i++)raysA.push_back(*i);
1249 IntegerVectorList gensA=A.generatorsOfLinealitySpace();
1250 IntegerVectorList gensB=B.generatorsOfLinealitySpace();
1251 for(IntegerVectorList::const_iterator i=gensB.begin();i!=gensB.end();i++)gensA.push_back(*i);
1252
1253
1254 return PolyhedralCone::givenByRays(raysA,gensA,A.ambientDimension());
1255 }
1256
inducedLink(PolyhedralCone const & c,PolyhedralCone const & l)1257 IntegerVectorList inducedLink(PolyhedralCone const &c, PolyhedralCone const &l)
1258 {
1259 PolyhedralCone A=sum(c,l);
1260 A.canonicalize();
1261 IntegerVectorList ret;
1262 if(A.dimension()==l.dimension())return ret;
1263 assert(l.dimension()+1==A.dimension());
1264 int nineq=A.getHalfSpaces().size();
1265 assert(nineq<=1);
1266 if(nineq==1)
1267 {
1268 ret.push_back(A.getUniquePoint());
1269 }
1270 else
1271 {
1272 PolyhedralCone ADual=A.dualCone();
1273 IntegerVectorList L=l.getEquations();
1274 for(IntegerVectorList::const_iterator i=L.begin();i!=L.end();i++)
1275 if(!ADual.contains(*i))
1276 {
1277 IntegerVectorList temp;temp.push_back(*i);
1278 PolyhedralCone ray(temp,A.getImpliedEquations());
1279 ray.canonicalize();
1280 IntegerVector v=ray.getUniquePoint();
1281 ret.push_back(v);
1282 ret.push_back(-v);
1283 break;
1284 }
1285 assert(!ret.empty());
1286 }
1287 return ret;
1288 }
1289
doesSatisfyInequalityExpensive(IntegerVector const & ineq) const1290 bool PolyhedralCone::doesSatisfyInequalityExpensive(IntegerVector const &ineq)const
1291 {
1292 if(!haveExtremeRaysBeenCached)this->extremeRays();
1293 if(!haveGeneratorsOfLinealitySpaceBeenCached)this->generatorsOfLinealitySpace();
1294 assert(haveGeneratorsOfLinealitySpaceBeenCached);
1295 assert(haveExtremeRaysBeenCached);
1296
1297 for(IntegerVectorList::const_iterator i=cachedGeneratorsOfLinealitySpace.begin();i!=cachedGeneratorsOfLinealitySpace.end();i++)
1298 if(dotLong(*i,ineq)!=0)return false;
1299 for(IntegerVectorList::const_iterator i=cachedExtremeRays.begin();i!=cachedExtremeRays.end();i++)
1300 if(dotLong(*i,ineq)<0)return false;
1301
1302 return true;
1303 }
1304
1305
triangulation() const1306 list<PolyhedralCone> PolyhedralCone::triangulation()const
1307 {
1308 IntegerVectorList rays=this->extremeRays();
1309 rays.sort();
1310 IntegerVectorList lines=this->generatorsOfLinealitySpace();
1311 IntegerMatrix M=rowsToIntegerMatrix(rays,n);
1312 IntegerMatrix M2=M.transposed();
1313 Triangulation2 T(M2);
1314 T.triangulate();
1315 T.changeToTriangulationInducedBy(LexicographicTermOrder());
1316 list<PolyhedralCone> ret;
1317 for(set<IntegerVector>::const_iterator i=T.bases.begin();i!=T.bases.end();i++)
1318 {
1319 IntegerVectorList gens;
1320 for(int j=0;j<i->size();j++)gens.push_back(M[(*i)[j]]);
1321 ret.push_back(PolyhedralCone::givenByRays(gens,lines,n));
1322 }
1323 return ret;
1324 }
1325