1 /*
2 * lib_cone.cpp
3 *
4 * Created on: Sep 29, 2010
5 * Author: anders
6 */
7
8 #include "gfanlib_zcone.h"
9
10 #include <vector>
11 #include <set>
12 #include <sstream>
13
14 //extern "C"{
15 #ifdef NOCDDPREFIX
16 #include "setoper.h"
17 #include "cdd.h"
18 #else
19 #include "cdd/setoper.h"
20 #include "cdd/cdd.h"
21 #endif
22 //}
23
24 namespace gfan{
isCddlibRequired()25 bool isCddlibRequired()
26 {
27 return true;
28 }
initializeCddlibIfRequired()29 void initializeCddlibIfRequired() // calling this frequently will cause memory leaks because deinitialisation is not possible with old versions of cddlib.
30 {
31 dd_set_global_constants();
32 }
deinitializeCddlibIfRequired()33 void deinitializeCddlibIfRequired()
34 {
35 // dd_free_global_constants();
36 }
ensureCddInitialisation()37 static void ensureCddInitialisation()
38 {
39 // A more complicated initialisation than the following (meaning attempts to count the number of times
40 // cddlib was requested to be initialised) would require cddlib to be thread aware.
41 // The error below is implemented with an assert(0) because throwing an exception may leave the impression that
42 // it is possible to recover from this error. While that may be true, it would not work in full generality,
43 // as the following if statement cannot test whether dd_free_global_constants() has also been called.
44 // Moverover, in multithreaded environments it would be quite difficult to decide if cddlib was initialised.
45 if(!dd_one[0]._mp_num._mp_d)
46 {
47 std::cerr<<"CDDLIB HAS NOT BEEN INITIALISED!\n"
48 "\n"
49 "Fix this problem by calling the following function in your initialisation code:\n"
50 "dd_set_global_constants();\n"
51 "(after possibly setting the gmp allocators) and\n"
52 "dd_free_global_constants()\n"
53 "in your deinitialisation code (only available for cddlib version>=094d).\n"
54 "This requires the header includes:\n"
55 "#include \"cdd/setoper.h\"\n"
56 "#include \"cdd/cdd.h\"\n"
57 "\n"
58 "Alternatively, you may call gfan:initializeCddlibIfRequired() and deinitializeCddlibIfRequired()\n"
59 "if gfanlib is the only code using cddlib. If at some point cddlib is no longer required by gfanlib\n"
60 "these functions may do nothing.\n"
61 "Because deinitialisation is not possible in cddlib <094d, the functions may leak memory and should not be called often.\n"
62 "\n"
63 "This error message will never appear if the initialisation was done properly, and therefore never appear in a shipping version of your software.\n";
64 assert(0);
65 }
66 /*
67 static bool initialized;
68 if(!initialized)
69 {
70 dd_set_global_constants();
71 initialized=true;
72 }*/
73 }
74
75
76 class LpSolver
77 {
ZMatrix2MatrixGmp(ZMatrix const & g,dd_ErrorType * Error)78 static dd_MatrixPtr ZMatrix2MatrixGmp(ZMatrix const &g, dd_ErrorType *Error)
79 {
80 int n=g.getWidth();
81 dd_MatrixPtr M=NULL;
82 dd_rowrange m_input,i;
83 dd_colrange d_input,j;
84 dd_RepresentationType rep=dd_Inequality;
85 // dd_boolean found=dd_FALSE, newformat=dd_FALSE, successful=dd_FALSE;
86 char command[dd_linelenmax], comsave[dd_linelenmax];
87 dd_NumberType NT;
88
89 (*Error)=dd_NoError;
90
91 rep=dd_Inequality; // newformat=dd_TRUE;
92
93 m_input=g.getHeight();
94 d_input=n+1;
95
96 NT=dd_Rational;
97
98 M=dd_CreateMatrix(m_input, d_input);
99 M->representation=rep;
100 M->numbtype=NT;
101
102 for (i = 0; i < m_input; i++) {
103 dd_set_si(M->matrix[i][0],0);
104 for (j = 1; j < d_input; j++) {
105 g[i][j-1].setGmp(mpq_numref(M->matrix[i][j]));
106 mpz_set_ui(mpq_denref(M->matrix[i][j]), 1);
107 mpq_canonicalize(M->matrix[i][j]);
108 }
109 }
110
111 // successful=dd_TRUE;
112
113 return M;
114 }
ZMatrix2MatrixGmp(ZMatrix const & inequalities,ZMatrix const & equations,dd_ErrorType * err)115 static dd_MatrixPtr ZMatrix2MatrixGmp(ZMatrix const &inequalities, ZMatrix const &equations, dd_ErrorType *err)
116 {
117 ZMatrix g=inequalities;
118 g.append(equations);
119 int numberOfInequalities=inequalities.getHeight();
120 int numberOfRows=g.getHeight();
121 dd_MatrixPtr A=NULL;
122 ensureCddInitialisation();
123 A=ZMatrix2MatrixGmp(g, err);
124 for(int i=numberOfInequalities;i<numberOfRows;i++)
125 set_addelem(A->linset,i+1);
126 return A;
127 }
getConstraints(dd_MatrixPtr A,bool returnEquations)128 static ZMatrix getConstraints(dd_MatrixPtr A, bool returnEquations)
129 {
130 int rowsize=A->rowsize;
131 int n=A->colsize-1;
132
133 ZMatrix ret(0,n);
134 for(int i=0;i<rowsize;i++)
135 {
136 bool isEquation=set_member(i+1,A->linset);
137 if(isEquation==returnEquations)
138 {
139 QVector v(n);
140 for(int j=0;j<n;j++)v[j]=Rational(A->matrix[i][j+1]);
141 ret.appendRow(QToZVectorPrimitive(v));
142 }
143 }
144 return ret;
145 }
isFacet(ZMatrix const & g,int index)146 static bool isFacet(ZMatrix const &g, int index)
147 {
148 bool ret;
149 // dd_MatrixPtr M=NULL,M2=NULL,M3=NULL;
150 dd_MatrixPtr M=NULL;
151 // dd_colrange d;
152 dd_ErrorType err=dd_NoError;
153 // dd_rowset redrows,linrows,ignoredrows, basisrows;
154 // dd_colset ignoredcols, basiscols;
155 // dd_DataFileType inputfile;
156 FILE *reading=NULL;
157
158 ensureCddInitialisation();
159
160 M=ZMatrix2MatrixGmp(g, &err);
161 if (err!=dd_NoError) goto _L99;
162
163 // d=M->colsize;
164
165 //static dd_Arow temp;
166 dd_Arow temp;
167 dd_InitializeArow(g.getWidth()+1,&temp);
168
169 ret= !dd_Redundant(M,index+1,temp,&err);
170
171 dd_FreeMatrix(M);
172 dd_FreeArow(g.getWidth()+1,temp);
173
174 if (err!=dd_NoError) goto _L99;
175 return ret;
176 _L99:
177 assert(0);
178 return false;
179 }
180
181 /*
182 Heuristic for checking if inequality of full dimensional cone is a
183 facet. If the routine returns true then the inequality is a
184 facet. If it returns false it is unknown.
185 */
fastIsFacetCriterion(ZMatrix const & normals,int i)186 static bool fastIsFacetCriterion(ZMatrix const &normals, int i)
187 {
188 int n=normals.getWidth();
189 for(int j=0;j<n;j++)
190 if(normals[i][j].sign()!=0)
191 {
192 int sign=normals[i][j].sign();
193 bool isTheOnly=true;
194 for(int k=0;k<normals.getHeight();k++)
195 if(k!=i)
196 {
197 if(normals[i][j].sign()==sign)
198 {
199 isTheOnly=false;
200 break;
201 }
202 }
203 if(isTheOnly)return true;
204 }
205 return false;
206 }
207
fastIsFacet(ZMatrix const & normals,int i)208 static bool fastIsFacet(ZMatrix const &normals, int i)
209 {
210 if(fastIsFacetCriterion(normals,i))return true;
211 return isFacet(normals,i);
212 }
213
214 class MyHashMap
215 {
216 typedef std::vector<std::set<ZVector> > Container;
217 Container container;
218 int tableSize;
219 public:
220 class iterator
221 {
222 class MyHashMap &hashMap;
223 int index; // having index==-1 means that we are before/after the elements.
224 std::set<ZVector>::iterator i;
225 public:
operator ++()226 bool operator++()
227 {
228 if(index==-1)goto search;
229 i++;
230 while(i==hashMap.container[index].end())
231 {
232 search:
233 index++;
234 if(index>=hashMap.tableSize){
235 index=-1;
236 return false;
237 }
238 i=hashMap.container[index].begin();
239 }
240 return true;
241 }
operator *() const242 ZVector const & operator*()const
243 {
244 return *i;
245 }
operator *()246 ZVector operator*()
247 {
248 return *i;
249 }
iterator(MyHashMap & hashMap_)250 iterator(MyHashMap &hashMap_):
251 hashMap(hashMap_)
252 {
253 index=-1;
254 }
255 };
function(const ZVector & v)256 unsigned int function(const ZVector &v)
257 {
258 unsigned int ret=0;
259 int n=v.size();
260 for(int i=0;i<n;i++)
261 ret=(ret<<3)+(ret>>29)+v.UNCHECKEDACCESS(i).hashValue();
262 return ret%tableSize;
263 }
MyHashMap(int tableSize_)264 MyHashMap(int tableSize_):
265 container(tableSize_),
266 tableSize(tableSize_)
267 {
268 assert(tableSize_>0);
269 }
insert(const ZVector & v)270 void insert(const ZVector &v)
271 {
272 container[function(v)].insert(v);
273 }
erase(ZVector const & v)274 void erase(ZVector const &v)
275 {
276 container[function(v)].erase(v);
277 }
begin()278 iterator begin()
279 {
280 iterator ret(*this);
281 ++ ret;
282 return ret;
283 }
size()284 int size()
285 {
286 iterator i=begin();
287 int ret=0;
288 do{ret++;}while(++i);
289 return ret;
290 }
291 };
292
293
normalizedWithSumsAndDuplicatesRemoved(ZMatrix const & a)294 static ZMatrix normalizedWithSumsAndDuplicatesRemoved(ZMatrix const &a)
295 {
296 // TODO: write a version of this function which will work faster if the entries fit in 32bit
297 if(a.getHeight()==0)return a;
298 int n=a.getWidth();
299 ZVector temp1(n);
300 // ZVector temp2(n);
301 ZMatrix ret(0,n);
302 MyHashMap b(a.getHeight());
303
304 for(int i=0;i<a.getHeight();i++)
305 {
306 assert(!(a[i].toVector().isZero()));
307 b.insert(a[i].toVector().normalized());
308 }
309
310 {
311 MyHashMap::iterator i=b.begin();
312
313 do
314 {
315 MyHashMap::iterator j=i;
316 while(++j)
317 {
318 ZVector const &I=*i;
319 ZVector const &J=*j;
320 for(int k=0;k<n;k++)temp1[k]=I.UNCHECKEDACCESS(k)+J.UNCHECKEDACCESS(k);
321 // normalizedLowLevel(temp1,temp2);
322 // b.erase(temp2);//this can never remove *i or *j
323 b.erase(temp1.normalized());//this can never remove *i or *j
324 }
325 }
326 while(++i);
327 }
328 ZMatrix original(0,n);
329 {
330 MyHashMap::iterator i=b.begin();
331 do
332 {
333 original.appendRow(*i);
334 }
335 while(++i);
336 }
337
338 for(int i=0;i!=original.getHeight();i++)
339 for(int j=0;j!=a.getHeight();j++)
340 if(!dependent(original[i].toVector(),a[j].toVector()))
341 {
342 ZVector const &I=original[i];
343 ZVector const &J=a[j];
344 for(int k=0;k<n;k++)temp1[k]=I.UNCHECKEDACCESS(k)+J.UNCHECKEDACCESS(k);
345 // normalizedLowLevel(temp1,temp2);
346 // b.erase(temp2);//this can never remove *i or *j
347 b.erase(temp1.normalized());//this can never remove *i or *j
348 }
349 {
350 MyHashMap::iterator i=b.begin();
351 do
352 {
353 ZVector temp=*i;
354 ret.appendRow(temp);
355 }
356 while(++i);
357 }
358 return ret;
359 }
360 public:
fastNormals(ZMatrix const & inequalities)361 static ZMatrix fastNormals(ZMatrix const &inequalities)
362 {
363 ZMatrix normals=normalizedWithSumsAndDuplicatesRemoved(inequalities);
364 for(int i=0;i!=normals.getHeight();i++)
365 if(!fastIsFacet(normals,i))
366 {
367 normals[i]=normals[normals.getHeight()-1];
368 normals.eraseLastRow();
369 i--;
370 }
371 return normals;
372 }
removeRedundantRows(ZMatrix & inequalities,ZMatrix & equations,bool removeInequalityRedundancies)373 void removeRedundantRows(ZMatrix &inequalities, ZMatrix &equations, bool removeInequalityRedundancies)
374 {
375 ensureCddInitialisation();
376
377 int numberOfEqualities=equations.getHeight();
378 int numberOfInequalities=inequalities.getHeight();
379 int numberOfRows=numberOfEqualities+numberOfInequalities;
380
381 if(numberOfRows==0)return;//the full space, so description is already irredundant
382
383 // dd_rowset r=NULL;
384 ZMatrix g=inequalities;
385 g.append(equations);
386
387 // dd_LPSolverType solver=dd_DualSimplex;
388 dd_MatrixPtr A=NULL;
389 dd_ErrorType err=dd_NoError;
390
391 A=ZMatrix2MatrixGmp(g,&err);
392 if (err!=dd_NoError) goto _L99;
393
394 for(int i=numberOfInequalities;i<numberOfRows;i++)
395 set_addelem(A->linset,i+1);
396
397 A->objective=dd_LPmax;
398
399 dd_rowset impl_linset;
400 dd_rowset redset;
401 dd_rowindex newpos;
402
403 if(removeInequalityRedundancies)
404 dd_MatrixCanonicalize(&A, &impl_linset, &redset, &newpos, &err);
405 else
406 dd_MatrixCanonicalizeLinearity(&A, &impl_linset, &newpos, &err);
407
408 if (err!=dd_NoError) goto _L99;
409
410 {
411 int n=A->colsize-1;
412 equations=ZMatrix(0,n); //TODO: the number of rows needed is actually known
413 inequalities=ZMatrix(0,n); //is known by set_card(). That might save some copying.
414
415 {
416 int rowsize=A->rowsize;
417 QVector point(n);
418 for(int i=0;i<rowsize;i++)
419 {
420 for(int j=0;j<n;j++)point[j]=Rational(A->matrix[i][j+1]);
421 ((set_member(i+1,A->linset))?equations:inequalities).appendRow(QToZVectorPrimitive(point));
422 }
423 }
424 assert(set_card(A->linset)==equations.getHeight());
425 assert(A->rowsize==equations.getHeight()+inequalities.getHeight());
426
427 set_free(impl_linset);
428 if(removeInequalityRedundancies)
429 set_free(redset);
430 free(newpos);
431
432 dd_FreeMatrix(A);
433 return;
434 }
435 _L99:
436 assert(!"Cddlib reported error when called by Gfanlib.");
437 }
relativeInteriorPoint(const ZMatrix & inequalities,const ZMatrix & equations)438 ZVector relativeInteriorPoint(const ZMatrix &inequalities, const ZMatrix &equations)
439 {
440 QVector retUnscaled(inequalities.getWidth());
441 ensureCddInitialisation();
442 int numberOfEqualities=equations.getHeight();
443 int numberOfInequalities=inequalities.getHeight();
444 int numberOfRows=numberOfEqualities+numberOfInequalities;
445
446 // dd_rowset r=NULL;
447 ZMatrix g=inequalities;
448 g.append(equations);
449
450 dd_LPSolverType solver=dd_DualSimplex;
451 dd_MatrixPtr A=NULL;
452 dd_ErrorType err=dd_NoError;
453
454 A=ZMatrix2MatrixGmp(g,&err);
455 if (err!=dd_NoError) goto _L99;
456 {
457 dd_LPSolutionPtr lps1;
458 dd_LPPtr lp,lp1;
459
460 for(int i=0;i<numberOfInequalities;i++)
461 dd_set_si(A->matrix[i][0],-1);
462 for(int i=numberOfInequalities;i<numberOfRows;i++)
463 set_addelem(A->linset,i+1);
464
465 A->objective=dd_LPmax;
466 lp=dd_Matrix2LP(A, &err);
467 if (err!=dd_NoError) goto _L99;
468
469 lp1=dd_MakeLPforInteriorFinding(lp);
470 dd_LPSolve(lp1,solver,&err);
471 if (err!=dd_NoError) goto _L99;
472
473 lps1=dd_CopyLPSolution(lp1);
474
475 assert(!dd_Negative(lps1->optvalue));
476
477 for (int j=1; j <(lps1->d)-1; j++)
478 retUnscaled[j-1]=Rational(lps1->sol[j]);
479
480 dd_FreeLPData(lp);
481 dd_FreeLPSolution(lps1);
482 dd_FreeLPData(lp1);
483 dd_FreeMatrix(A);
484 return QToZVectorPrimitive(retUnscaled);
485 }
486 _L99:
487 assert(0);
488 return QToZVectorPrimitive(retUnscaled);
489 }
dual(ZMatrix const & inequalities,ZMatrix const & equations,ZMatrix & dualInequalities,ZMatrix & dualEquations)490 void dual(ZMatrix const &inequalities, ZMatrix const &equations, ZMatrix &dualInequalities, ZMatrix &dualEquations)
491 {
492 int result;
493
494 dd_MatrixPtr A=NULL;
495 dd_ErrorType err=dd_NoError;
496
497 ensureCddInitialisation();
498
499 A=ZMatrix2MatrixGmp(inequalities, equations, &err);
500
501 dd_PolyhedraPtr poly;
502 poly=dd_DDMatrix2Poly2(A, dd_LexMin, &err);
503
504 if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) assert(0);
505
506 dd_MatrixPtr A2=dd_CopyGenerators(poly);
507
508 dualInequalities=getConstraints(A2,false);
509 dualEquations=getConstraints(A2,true);
510
511 dd_FreeMatrix(A2);
512 dd_FreeMatrix(A);
513 dd_FreePolyhedra(poly);
514
515 return;
516 _L99:
517 assert(0);
518 }
519 // this procedure is take from cddio.c.
dd_ComputeAinc(dd_PolyhedraPtr poly)520 static void dd_ComputeAinc(dd_PolyhedraPtr poly)
521 {
522 /* This generates the input incidence array poly->Ainc, and
523 two sets: poly->Ared, poly->Adom.
524 */
525 dd_bigrange k;
526 dd_rowrange i,m1;
527 dd_colrange j;
528 dd_boolean redundant;
529 dd_MatrixPtr M=NULL;
530 mytype sum,temp;
531
532 dd_init(sum); dd_init(temp);
533 if (poly->AincGenerated==dd_TRUE) goto _L99;
534
535 M=dd_CopyOutput(poly);
536 poly->n=M->rowsize;
537 m1=poly->m1;
538
539 /* this number is same as poly->m, except when
540 poly is given by nonhomogeneous inequalty:
541 !(poly->homogeneous) && poly->representation==Inequality,
542 it is poly->m+1. See dd_ConeDataLoad.
543 */
544 poly->Ainc=(set_type*)calloc(m1, sizeof(set_type));
545 for(i=1; i<=m1; i++) set_initialize(&(poly->Ainc[i-1]),poly->n);
546 set_initialize(&(poly->Ared), m1);
547 set_initialize(&(poly->Adom), m1);
548
549 for (k=1; k<=poly->n; k++){
550 for (i=1; i<=poly->m; i++){
551 dd_set(sum,dd_purezero);
552 for (j=1; j<=poly->d; j++){
553 dd_mul(temp,poly->A[i-1][j-1],M->matrix[k-1][j-1]);
554 dd_add(sum,sum,temp);
555 }
556 if (dd_EqualToZero(sum)) {
557 set_addelem(poly->Ainc[i-1], k);
558 }
559 }
560 if (!(poly->homogeneous) && poly->representation==dd_Inequality){
561 if (dd_EqualToZero(M->matrix[k-1][0])) {
562 set_addelem(poly->Ainc[m1-1], k); /* added infinity inequality (1,0,0,...,0) */
563 }
564 }
565 }
566
567 for (i=1; i<=m1; i++){
568 if (set_card(poly->Ainc[i-1])==M->rowsize){
569 set_addelem(poly->Adom, i);
570 }
571 }
572 for (i=m1; i>=1; i--){
573 if (set_card(poly->Ainc[i-1])==0){
574 redundant=dd_TRUE;
575 set_addelem(poly->Ared, i);
576 }else {
577 redundant=dd_FALSE;
578 for (k=1; k<=m1; k++) {
579 if (k!=i && !set_member(k, poly->Ared) && !set_member(k, poly->Adom) &&
580 set_subset(poly->Ainc[i-1], poly->Ainc[k-1])){
581 if (!redundant){
582 redundant=dd_TRUE;
583 }
584 set_addelem(poly->Ared, i);
585 }
586 }
587 }
588 }
589 dd_FreeMatrix(M);
590 poly->AincGenerated=dd_TRUE;
591 _L99:;
592 dd_clear(sum); dd_clear(temp);
593 }
594
595
extremeRaysInequalityIndices(const ZMatrix & inequalities)596 std::vector<std::vector<int> > extremeRaysInequalityIndices(const ZMatrix &inequalities)
597 {
598 int dim2=inequalities.getHeight();
599 if(dim2==0)return std::vector<std::vector<int> >();
600 // int dimension=inequalities.getWidth();
601
602 dd_MatrixPtr A=NULL;
603 dd_ErrorType err=dd_NoError;
604
605 ensureCddInitialisation();
606 A=ZMatrix2MatrixGmp(inequalities, &err);
607
608 dd_PolyhedraPtr poly;
609 poly=dd_DDMatrix2Poly2(A, dd_LexMin, &err);
610
611 if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) assert(0);
612 if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
613
614 std::vector<std::vector<int> > ret;
615
616 /*
617 How do we interpret the cddlib output? For a long ting gfan has
618 been using poly->n as the number of rays of the cone and thus
619 returned sets of indices that actually gave the lineality
620 space. The mistake was then caught later in PolyhedralCone. On Feb
621 17 2009 gfan was changed to check the length of each set to make
622 sure that it does not specify the lineality space and only return
623 those sets giving rise to rays. This does not seem to be the best
624 strategy and might even be wrong.
625 */
626
627
628 for (int k=1; k<=poly->n; k++)
629 {
630 int length=0;
631 for (int i=1; i<=poly->m1; i++)
632 if(set_member(k,poly->Ainc[i-1]))length++;
633 if(length!=dim2)
634 {
635 std::vector<int> v(length);
636 int j=0;
637 for (int i=1; i<=poly->m1; i++)
638 if(set_member(k,poly->Ainc[i-1]))v[j++]=i-1;
639 ret.push_back(v);
640 }
641 }
642
643 dd_FreeMatrix(A);
644 dd_FreePolyhedra(poly);
645
646 return ret;
647 _L99:
648 assert(0);
649 return std::vector<std::vector<int> >();
650 }
651
652 };
653
654 LpSolver lpSolver;
655
isInStateMinimum(int s) const656 bool ZCone::isInStateMinimum(int s)const
657 {
658 return state>=s;
659 }
660
661
operator <(ZCone const & a,ZCone const & b)662 bool operator<(ZCone const &a, ZCone const &b)
663 {
664 assert(a.state>=3);
665 assert(b.state>=3);
666
667 if(a.n<b.n)return true;
668 if(a.n>b.n)return false;
669
670 if(a.equations<b.equations)return true;
671 if(b.equations<a.equations)return false;
672
673 if(a.inequalities<b.inequalities)return true;
674 if(b.inequalities<a.inequalities)return false;
675
676 return false;
677 }
678
679
operator !=(ZCone const & a,ZCone const & b)680 bool operator!=(ZCone const &a, ZCone const &b)
681 {
682 return (a<b)||(b<a);
683 }
684
685
ensureStateAsMinimum(int s) const686 void ZCone::ensureStateAsMinimum(int s)const
687 {
688 if((state<1) && (s==1))
689 {
690 {
691 QMatrix m=ZToQMatrix(equations);
692 m.reduce();
693 m.removeZeroRows();
694
695 ZMatrix newInequalities(0,inequalities.getWidth());
696 for(int i=0;i<inequalities.getHeight();i++)
697 {
698 QVector w=ZToQVector(inequalities[i]);
699 w=m.canonicalize(w);
700 if(!w.isZero())
701 newInequalities.appendRow(QToZVectorPrimitive(w));
702 }
703
704 inequalities=newInequalities;
705 inequalities.sortAndRemoveDuplicateRows();
706 equations=QToZMatrixPrimitive(m);
707 }
708
709 if(!(preassumptions&PCP_impliedEquationsKnown))
710 if(inequalities.getHeight()>1)//there can be no implied equation unless we have at least two inequalities
711 lpSolver.removeRedundantRows(inequalities,equations,false);
712
713 assert(inequalities.getWidth()==equations.getWidth());
714 }
715 if((state<2) && (s>=2) && !(preassumptions&PCP_facetsKnown))
716 {
717 /* if(inequalities.size()>25)
718 {
719 IntegerVectorList h1;
720 IntegerVectorList h2;
721 bool a=false;
722 for(IntegerVectorList::const_iterator i=inequalities.begin();i!=inequalities.end();i++)
723 {
724 if(a)
725 h1.push_back(*i);
726 else
727 h2.push_back(*i);
728 a=!a;
729 }
730 PolyhedralCone c1(h1,equations);
731 PolyhedralCone c2(h2,equations);
732 c1.ensureStateAsMinimum(2);
733 c2.ensureStateAsMinimum(2);
734 inequalities=c1.inequalities;
735 for(IntegerVectorList::const_iterator i=c2.inequalities.begin();i!=c2.inequalities.end();i++)
736 inequalities.push_back(*i);
737 }
738 */
739 if(equations.getHeight())
740 {
741 QMatrix m=ZToQMatrix(equations);
742 m.reduce();
743 m.REformToRREform();
744 ZMatrix inequalities2(0,equations.getWidth());
745 for(int i=0;i<inequalities.getHeight();i++)
746 {
747 inequalities2.appendRow(QToZVectorPrimitive(m.canonicalize(ZToQVector(inequalities[i]))));
748 }
749 inequalities=LpSolver::fastNormals(inequalities2);
750 goto noFallBack;
751 fallBack://alternativ (disabled)
752 lpSolver.removeRedundantRows(inequalities,equations,true);
753 noFallBack:;
754 }
755 else
756 inequalities=LpSolver::fastNormals(inequalities);
757 }
758 if((state<3) && (s>=3))
759 {
760 QMatrix equations2=ZToQMatrix(equations);
761 equations2.reduce(false,false,true);
762 equations2.REformToRREform(true);
763 for(int i=0;i<inequalities.getHeight();i++)
764 {
765 inequalities[i]=QToZVectorPrimitive(equations2.canonicalize(ZToQVector(inequalities[i])));
766 }
767 inequalities.sortRows();
768 equations=QToZMatrixPrimitive(equations2);
769 }
770 if(state<s)
771 state=s;
772 }
773
operator <<(std::ostream & f,ZCone const & c)774 std::ostream &operator<<(std::ostream &f, ZCone const &c)
775 {
776 f<<"Ambient dimension:"<<c.n<<std::endl;
777 f<<"Inequalities:"<<std::endl;
778 f<<c.inequalities<<std::endl;
779 f<<"Equations:"<<std::endl;
780 f<<c.equations<<std::endl;
781 return f;
782 }
783
toString() const784 std::string ZCone::toString()const
785 {
786 std::stringstream f;
787 f<<*this;
788 return f.str();
789 }
790
ZCone(int ambientDimension)791 ZCone::ZCone(int ambientDimension):
792 preassumptions(PCP_impliedEquationsKnown|PCP_facetsKnown),
793 state(1),
794 n(ambientDimension),
795 multiplicity(1),
796 linearForms(ZMatrix(0,ambientDimension)),
797 inequalities(ZMatrix(0,ambientDimension)),
798 equations(ZMatrix(0,ambientDimension)),
799 haveExtremeRaysBeenCached(false)
800 {
801 }
802
803
ZCone(ZMatrix const & inequalities_,ZMatrix const & equations_,int preassumptions_)804 ZCone::ZCone(ZMatrix const &inequalities_, ZMatrix const &equations_, int preassumptions_):
805 preassumptions(preassumptions_),
806 state(0),
807 n(inequalities_.getWidth()),
808 multiplicity(1),
809 linearForms(ZMatrix(0,inequalities_.getWidth())),
810 inequalities(inequalities_),
811 equations(equations_),
812 haveExtremeRaysBeenCached(false)
813 {
814 assert(preassumptions_<4);//OTHERWISE WE ARE DOING SOMETHING STUPID LIKE SPECIFYING AMBIENT DIMENSION
815 assert(equations_.getWidth()==n);
816 ensureStateAsMinimum(1);
817 }
818
canonicalize()819 void ZCone::canonicalize()
820 {
821 ensureStateAsMinimum(3);
822 }
823
findFacets()824 void ZCone::findFacets()
825 {
826 ensureStateAsMinimum(2);
827 }
828
getFacets() const829 ZMatrix ZCone::getFacets()const
830 {
831 ensureStateAsMinimum(2);
832 return inequalities;
833 }
834
findImpliedEquations()835 void ZCone::findImpliedEquations()
836 {
837 ensureStateAsMinimum(1);
838 }
839
getImpliedEquations() const840 ZMatrix ZCone::getImpliedEquations()const
841 {
842 ensureStateAsMinimum(1);
843 return equations;
844 }
845
getRelativeInteriorPoint() const846 ZVector ZCone::getRelativeInteriorPoint()const
847 {
848 ensureStateAsMinimum(1);
849 // assert(state>=1);
850
851 return lpSolver.relativeInteriorPoint(inequalities,equations);
852 }
853
getUniquePoint() const854 ZVector ZCone::getUniquePoint()const
855 {
856 ZMatrix rays=extremeRays();
857 ZVector ret(n);
858 for(int i=0;i<rays.getHeight();i++)
859 ret+=rays[i];
860
861 return ret;
862 }
863
getUniquePointFromExtremeRays(ZMatrix const & extremeRays) const864 ZVector ZCone::getUniquePointFromExtremeRays(ZMatrix const &extremeRays)const
865 {
866 ZVector ret(n);
867 for(int i=0;i<extremeRays.getHeight();i++)
868 if(contains(extremeRays[i]))ret+=extremeRays[i];
869 return ret;
870 }
871
872
ambientDimension() const873 int ZCone::ambientDimension()const
874 {
875 return n;
876 }
877
878
codimension() const879 int ZCone::codimension()const
880 {
881 return ambientDimension()-dimension();
882 }
883
884
dimension() const885 int ZCone::dimension()const
886 {
887 // assert(state>=1);
888 ensureStateAsMinimum(1);
889 return ambientDimension()-equations.getHeight();
890 }
891
892
dimensionOfLinealitySpace() const893 int ZCone::dimensionOfLinealitySpace()const
894 {
895 ZMatrix temp=inequalities;
896 temp.append(equations);
897 ZCone temp2(ZMatrix(0,n),temp);
898 return temp2.dimension();
899 }
900
901
isOrigin() const902 bool ZCone::isOrigin()const
903 {
904 return dimension()==0;
905 }
906
907
isFullSpace() const908 bool ZCone::isFullSpace()const
909 {
910 for(int i=0;i<inequalities.getHeight();i++)
911 if(!inequalities[i].isZero())return false;
912 for(int i=0;i<equations.getHeight();i++)
913 if(!equations[i].isZero())return false;
914 return true;
915 }
916
917
intersection(const ZCone & a,const ZCone & b)918 ZCone intersection(const ZCone &a, const ZCone &b)
919 {
920 assert(a.ambientDimension()==b.ambientDimension());
921 ZMatrix inequalities=a.inequalities;
922 inequalities.append(b.inequalities);
923 ZMatrix equations=a.equations;
924 equations.append(b.equations);
925
926 equations.sortAndRemoveDuplicateRows();
927 inequalities.sortAndRemoveDuplicateRows();
928
929 {
930 ZMatrix Aequations=a.equations;
931 ZMatrix Ainequalities=a.inequalities;
932 Aequations.sortAndRemoveDuplicateRows();
933 Ainequalities.sortAndRemoveDuplicateRows();
934 if((Ainequalities.getHeight()==inequalities.getHeight()) && (Aequations.getHeight()==equations.getHeight()))return a;
935 ZMatrix Bequations=b.equations;
936 ZMatrix Binequalities=b.inequalities;
937 Bequations.sortAndRemoveDuplicateRows();
938 Binequalities.sortAndRemoveDuplicateRows();
939 if((Binequalities.getHeight()==inequalities.getHeight()) && (Bequations.getHeight()==equations.getHeight()))return b;
940 }
941
942 return ZCone(inequalities,equations);
943 }
944
945 /*
946 PolyhedralCone product(const PolyhedralCone &a, const PolyhedralCone &b)
947 {
948 IntegerVectorList equations2;
949 IntegerVectorList inequalities2;
950
951 int n1=a.n;
952 int n2=b.n;
953
954 for(IntegerVectorList::const_iterator i=a.equations.begin();i!=a.equations.end();i++)
955 equations2.push_back(concatenation(*i,IntegerVector(n2)));
956 for(IntegerVectorList::const_iterator i=b.equations.begin();i!=b.equations.end();i++)
957 equations2.push_back(concatenation(IntegerVector(n1),*i));
958 for(IntegerVectorList::const_iterator i=a.inequalities.begin();i!=a.inequalities.end();i++)
959 inequalities2.push_back(concatenation(*i,IntegerVector(n2)));
960 for(IntegerVectorList::const_iterator i=b.inequalities.begin();i!=b.inequalities.end();i++)
961 inequalities2.push_back(concatenation(IntegerVector(n1),*i));
962
963 PolyhedralCone ret(inequalities2,equations2,n1+n2);
964 ret.setMultiplicity(a.getMultiplicity()*b.getMultiplicity());
965 ret.setLinearForm(concatenation(a.getLinearForm(),b.getLinearForm()));
966
967 ret.ensureStateAsMinimum(a.state);
968 ret.ensureStateAsMinimum(b.state);
969
970 return ret;
971 }*/
972
973
positiveOrthant(int dimension)974 ZCone ZCone::positiveOrthant(int dimension)
975 {
976 return ZCone(ZMatrix::identity(dimension),ZMatrix(0,dimension));
977 }
978
979
givenByRays(ZMatrix const & generators,ZMatrix const & linealitySpace)980 ZCone ZCone::givenByRays(ZMatrix const &generators, ZMatrix const &linealitySpace)
981 {
982 //rewrite modulo lineality space
983 /* ZMatrix newGenerators(generators.getHeight(),generators.getWidth());
984 {
985 QMatrix l=ZToQMatrix(linealitySpace);
986 l.reduce();
987 for(int i=0;i<generators.getHeight();i++)
988 newGenerators[i]=QToZVectorPrimitive(l.canonicalize(ZToQVector(generators[i])));
989 }
990 */
991 // ZCone dual(newGenerators,linealitySpace);
992 ZCone dual(generators,linealitySpace);
993 // dual.findFacets();
994 // dual.canonicalize();
995 ZMatrix inequalities=dual.extremeRays();
996
997 /* ZMatrix span=generators;
998 span.append(linealitySpace);
999 QMatrix m2Q=ZToQMatrix(span);
1000 ZMatrix equations=QToZMatrixPrimitive(m2Q.reduceAndComputeKernel());
1001 */
1002 ZMatrix equations=dual.generatorsOfLinealitySpace();
1003 // equations.reduce();equations.removeZeroRows();
1004
1005
1006 return ZCone(inequalities,equations,PCP_impliedEquationsKnown|PCP_facetsKnown);
1007 }
1008
1009
containsPositiveVector() const1010 bool ZCone::containsPositiveVector()const
1011 {
1012 ZCone temp=intersection(*this,ZCone::positiveOrthant(n));
1013 return temp.getRelativeInteriorPoint().isPositive();
1014 }
1015
1016
contains(ZVector const & v) const1017 bool ZCone::contains(ZVector const &v)const
1018 {
1019 for(int i=0;i<equations.getHeight();i++)
1020 {
1021 if(!dot(equations[i],v).isZero())return false;
1022 }
1023 for(int i=0;i<inequalities.getHeight();i++)
1024 {
1025 if(dot(inequalities[i],v).sign()<0)return false;
1026 }
1027 return true;
1028 }
1029
1030
containsRowsOf(ZMatrix const & m) const1031 bool ZCone::containsRowsOf(ZMatrix const &m)const
1032 {
1033 for(int i=0;i<m.getHeight();i++)
1034 if(!contains(m[i]))return false;
1035 return true;
1036 }
1037
1038
contains(ZCone const & c) const1039 bool ZCone::contains(ZCone const &c)const
1040 {
1041 ZCone c2=intersection(*this,c);
1042 ZCone c3=c;
1043 c2.canonicalize();
1044 c3.canonicalize();
1045 return !(c2!=c3);
1046 }
1047
1048
containsRelatively(ZVector const & v) const1049 bool ZCone::containsRelatively(ZVector const &v)const
1050 {
1051 ensureStateAsMinimum(1);
1052 // assert(state>=1);
1053 for(int i=0;i<equations.getHeight();i++)
1054 {
1055 if(!dot(equations[i],v).isZero())return false;
1056 }
1057 for(int i=0;i<inequalities.getHeight();i++)
1058 {
1059 if(dot(inequalities[i],v).sign()<=0)return false;
1060 }
1061 return true;
1062 }
1063
1064
isSimplicial() const1065 bool ZCone::isSimplicial()const
1066 {
1067 // assert(state>=2);
1068 ensureStateAsMinimum(2);
1069 return codimension()+inequalities.getHeight()+dimensionOfLinealitySpace()==n;
1070 }
1071
1072
linealitySpace() const1073 ZCone ZCone::linealitySpace()const
1074 {
1075 ZCone ret(ZMatrix(0,n),combineOnTop(equations,inequalities));
1076 // ret.ensureStateAsMinimum(state);
1077 return ret;
1078 }
1079
1080
dualCone() const1081 ZCone ZCone::dualCone()const
1082 {
1083 ensureStateAsMinimum(1);
1084 // assert(state>=1);
1085
1086 ZMatrix dualInequalities,dualEquations;
1087 lpSolver.dual(inequalities,equations,dualInequalities,dualEquations);
1088 ZCone ret(dualInequalities,dualEquations);
1089 ret.ensureStateAsMinimum(state);
1090
1091 return ret;
1092 }
1093
1094
negated() const1095 ZCone ZCone::negated()const
1096 {
1097 ZCone ret(-inequalities,equations,(areFacetsKnown()?PCP_facetsKnown:0)|(areImpliedEquationsKnown()?PCP_impliedEquationsKnown:0));
1098 // ret.ensureStateAsMinimum(state);
1099 return ret;
1100 }
1101
1102
extremeRays(ZMatrix const * generatorsOfLinealitySpace) const1103 ZMatrix ZCone::extremeRays(ZMatrix const *generatorsOfLinealitySpace)const
1104 {
1105 // assert((dimension()==ambientDimension()) || (state>=3));
1106 if(dimension()!=ambientDimension())
1107 ensureStateAsMinimum(3);
1108
1109 if(haveExtremeRaysBeenCached)return cachedExtremeRays;
1110 ZMatrix ret(0,n);
1111 std::vector<std::vector<int> > indices=lpSolver.extremeRaysInequalityIndices(inequalities);
1112
1113 for(unsigned i=0;i<indices.size();i++)
1114 {
1115 /* At this point we know lineality space, implied equations and
1116 also inequalities for the ray. To construct a vector on the
1117 ray which is stable under (or indendent of) angle and
1118 linarity preserving transformation we find the dimension 1
1119 subspace orthorgonal to the implied equations and the
1120 lineality space and pick a suitable primitive generator */
1121
1122 /* To be more precise,
1123 * let E be the set of equations, and v the inequality defining a ray R.
1124 * We wish to find a vector satisfying these, but it must also be orthogonal
1125 * to the lineality space of the cone, that is, in the span of {E,v}.
1126 * One way to get such a vector is to project v to E an get a vector p.
1127 * Then v-p is in the span of {E,v} by construction.
1128 * The vector v-p is also in the orthogonal complement to E by construction,
1129 * that is, the span of R.
1130 * We wish to argue that it is not zero.
1131 * That would imply that v=p, meaning that v is in the span of the equations.
1132 * However, that would contradict that R is a ray.
1133 * In case v-p does not satisfy the inequality v (is this possible?), we change the sign.
1134 *
1135 * As a consequence we need the following procedure
1136 * primitiveProjection():
1137 * Input: E,v
1138 * Output: A primitive representation of the vector v-p, where p is the projection of v onto E
1139 *
1140 * Notice that the output is a Q linear combination of the input and that p is
1141 * a linear combination of E. The check that p has been computed correctly,
1142 * it suffices to check that v-p satisfies the equations E.
1143 * The routine will actually first compute a multiple of v-p.
1144 * It will do this using floating point arithmetics. It will then transform
1145 * the coefficients to get the multiple of v-p into integers. Then it
1146 * verifies in exact arithmetics, that with these coefficients we get a point
1147 * satisfying E. It then returns the primitive vector on the ray v-p.
1148 * In case of a failure it falls back to an implementation using rational arithmetics.
1149 */
1150
1151
1152 std::vector<int> asVector(inequalities.getHeight());
1153 for(unsigned j=0;j<indices[i].size();j++){asVector[indices[i][j]]=1;}
1154 ZMatrix equations=this->equations;
1155 ZVector theInequality;
1156
1157 for(unsigned j=0;j<asVector.size();j++)
1158 if(asVector[j])
1159 equations.appendRow(inequalities[j]);
1160 else
1161 theInequality=inequalities[j];
1162
1163 assert(!theInequality.isZero());
1164
1165 ZVector thePrimitiveVector;
1166 if(generatorsOfLinealitySpace)
1167 {
1168 QMatrix temp=ZToQMatrix(combineOnTop(equations,*generatorsOfLinealitySpace));
1169 thePrimitiveVector=QToZVectorPrimitive(temp.reduceAndComputeVectorInKernel());
1170 }
1171 else
1172 {
1173 QMatrix linealitySpaceOrth=ZToQMatrix(combineOnTop(this->equations,inequalities));
1174
1175
1176 QMatrix temp=combineOnTop(linealitySpaceOrth.reduceAndComputeKernel(),ZToQMatrix(equations));
1177 thePrimitiveVector=QToZVectorPrimitive(temp.reduceAndComputeVectorInKernel());
1178 }
1179 if(!contains(thePrimitiveVector))thePrimitiveVector=-thePrimitiveVector;
1180 ret.appendRow(thePrimitiveVector);
1181 }
1182
1183 cachedExtremeRays=ret;
1184 haveExtremeRaysBeenCached=true;
1185
1186 return ret;
1187 }
1188
1189
getMultiplicity() const1190 Integer ZCone::getMultiplicity()const
1191 {
1192 return multiplicity;
1193 }
1194
1195
setMultiplicity(Integer const & m)1196 void ZCone::setMultiplicity(Integer const &m)
1197 {
1198 multiplicity=m;
1199 }
1200
1201
getLinearForms() const1202 ZMatrix ZCone::getLinearForms()const
1203 {
1204 return linearForms;
1205 }
1206
1207
setLinearForms(ZMatrix const & linearForms_)1208 void ZCone::setLinearForms(ZMatrix const &linearForms_)
1209 {
1210 linearForms=linearForms_;
1211 }
1212
1213
quotientLatticeBasis() const1214 ZMatrix ZCone::quotientLatticeBasis()const
1215 {
1216 // assert(isInStateMinimum(1));// Implied equations must have been computed in order to know the span of the cone
1217 ensureStateAsMinimum(1);
1218
1219
1220 int a=equations.getHeight();
1221 int b=inequalities.getHeight();
1222
1223 // Implementation below could be moved to nonLP part of code.
1224
1225 // small vector space defined by a+b equations.... big by a equations.
1226
1227 ZMatrix M=combineLeftRight(combineLeftRight(
1228 equations.transposed(),
1229 inequalities.transposed()),
1230 ZMatrix::identity(n));
1231 M.reduce(false,true);
1232 /*
1233 [A|B|I] is reduced to [A'|B'|C'] meaning [A'|B']=C'[A|B] and A'=C'A.
1234
1235 [A'|B'] is in row echelon form, implying that the rows of C' corresponding to zero rows
1236 of [A'|B'] generate the lattice cokernel of [A|B] - that is the linealityspace intersected with Z^n.
1237
1238 [A'] is in row echelon form, implying that the rows of C' corresponding to zero rows of [A'] generate
1239 the lattice cokernel of [A] - that is the span of the cone intersected with Z^n.
1240
1241 It is clear that the second row set is a superset of the first. Their difference is a basis for the quotient.
1242 */
1243 ZMatrix ret(0,n);
1244
1245 for(int i=0;i<M.getHeight();i++)
1246 if(M[i].toVector().subvector(0,a).isZero()&&!M[i].toVector().subvector(a,a+b).isZero())
1247 {
1248 ret.appendRow(M[i].toVector().subvector(a+b,a+b+n));
1249 }
1250
1251 return ret;
1252 }
1253
1254
semiGroupGeneratorOfRay() const1255 ZVector ZCone::semiGroupGeneratorOfRay()const
1256 {
1257 ZMatrix temp=quotientLatticeBasis();
1258 assert(temp.getHeight()==1);
1259 for(int i=0;i<inequalities.getHeight();i++)
1260 if(dot(temp[0].toVector(),inequalities[i].toVector()).sign()<0)
1261 {
1262 temp[0]=-temp[0].toVector();
1263 break;
1264 }
1265 return temp[0];
1266 }
1267
1268
link(ZVector const & w) const1269 ZCone ZCone::link(ZVector const &w)const
1270 {
1271 /* Observe that the inequalities giving rise to facets
1272 * also give facets in the link, if they are kept as
1273 * inequalities. This means that the state cannot decrease
1274 * when taking links - that is why we specify the PCP flags.
1275 */
1276 ZMatrix inequalities2(0,n);
1277 for(int j=0;j<inequalities.getHeight();j++)
1278 if(dot(w,inequalities[j]).sign()==0)inequalities2.appendRow(inequalities[j]);
1279 ZCone C(inequalities2,equations,(areImpliedEquationsKnown()?PCP_impliedEquationsKnown:0)|(areFacetsKnown()?PCP_facetsKnown:0));
1280 C.ensureStateAsMinimum(state);
1281
1282 C.setLinearForms(getLinearForms());
1283 C.setMultiplicity(getMultiplicity());
1284
1285 return C;
1286 }
1287
hasFace(ZCone const & f) const1288 bool ZCone::hasFace(ZCone const &f)const
1289 {
1290 if(!contains(f.getRelativeInteriorPoint()))return false;
1291 ZCone temp=faceContaining(f.getRelativeInteriorPoint());
1292 temp.canonicalize();
1293 // ZCone temp2=*this;
1294 ZCone temp2=f;
1295 temp2.canonicalize();
1296 // std::cout << temp << std::endl;
1297 // std::cout << temp2 << std::endl;
1298
1299 return !(temp2!=temp);
1300 }
1301
faceContaining(ZVector const & v) const1302 ZCone ZCone::faceContaining(ZVector const &v)const
1303 {
1304 assert(n==(int)v.size());
1305 assert(contains(v));
1306 ZMatrix newEquations=equations;
1307 ZMatrix newInequalities(0,n);
1308 for(int i=0;i<inequalities.getHeight();i++)
1309 if(dot(inequalities[i],v).sign()!=0)
1310 newInequalities.appendRow(inequalities[i]);
1311 else
1312 newEquations.appendRow(inequalities[i]);
1313
1314 ZCone ret(newInequalities,newEquations,(state>=1)?PCP_impliedEquationsKnown:0);
1315 ret.ensureStateAsMinimum(state);
1316 return ret;
1317 }
1318
1319
getInequalities() const1320 ZMatrix ZCone::getInequalities()const
1321 {
1322 return inequalities;
1323 }
1324
1325
getEquations() const1326 ZMatrix ZCone::getEquations()const
1327 {
1328 return equations;
1329 }
1330
1331
generatorsOfSpan() const1332 ZMatrix ZCone::generatorsOfSpan()const
1333 {
1334 ensureStateAsMinimum(1);
1335 QMatrix l=ZToQMatrix(equations);
1336 return QToZMatrixPrimitive(l.reduceAndComputeKernel());
1337 }
1338
1339
generatorsOfLinealitySpace() const1340 ZMatrix ZCone::generatorsOfLinealitySpace()const
1341 {
1342 QMatrix l=ZToQMatrix(combineOnTop(inequalities,equations));
1343 return QToZMatrixPrimitive(l.reduceAndComputeKernel());
1344 }
1345
1346 }
1347