1 #ifndef TRIANGULATION2_H_INCLUDED 2 #define TRIANGULATION2_H_INCLUDED 3 4 #include <set> 5 #include <iostream> 6 #include <algorithm> 7 using namespace std; 8 #include "polyhedralfan.h" 9 #include "matrix.h" 10 #include "field_rationals.h" 11 #include "graph.h" 12 #include "determinant.h" 13 #include "lp.h" 14 #include "wallideal.h" 15 #include "linalg.h" 16 #include "printer.h" 17 #include "log.h" 18 #include "symmetry.h" 19 #include "triangulation.h" 20 21 /** 22 TODO: Merge this class with the Triangulation class. 23 */ 24 25 class Triangulation2 26 { 27 // IntegerVector circuit(IntegerMatrix const &A, IntegerVector const &basis, int j); 28 29 /** 30 A projective vector configuration. The columns are the vectors. 31 */ 32 IntegerMatrix A; 33 IntegerMatrix Atransposed; 34 /** 35 The dual of a subdivision of the cone over A is a polyhedron with facet normals among the columns. 36 The subdivision being a triangulation is equivalent to the polyhedron being simple. 37 Bases is either the sets of facet normals giving rise to a poiont of the polyhedron 38 OR the list of triangles of the triangulation. 39 */ 40 public: Triangulation2(IntegerMatrix const & A_)41 Triangulation2(IntegerMatrix const &A_): 42 A(A_), 43 Atransposed(A_.transposed()) 44 { 45 } 46 set<IntegerVector> bases; getN()47 int getN()const 48 { 49 return A.getWidth(); 50 } getD()51 int getD()const 52 { 53 return A.getHeight(); 54 } complement(IntegerVector const & v,int n)55 IntegerVector complement(IntegerVector const &v, int n)const 56 { 57 IntegerVector ret(0); 58 int j=0; 59 for(int i=0;i<n;i++) 60 { 61 if(j>=v.size()) 62 ret.push_back(i); 63 else 64 if(i==v[j]) 65 j++; 66 else 67 { 68 ret.push_back(i); 69 } 70 } 71 return ret; 72 } 73 /** 74 * Test if the configuration 75 */ 76 #if 0 77 bool isHomogeneous() 78 { 79 /* FieldMatrix temp=integerMatrixToFieldMatrix(A.transposed(),Q); 80 FieldMatrix s=temp.solver(); 81 FieldVector y=s.canonicalize(concatenation(FieldVector::allOnes(Q,temp.getHeight()),FieldVector(Q,temp.getWidth()))); 82 bool ret=y.subvector(0,temp.getHeight()).isZero(); 83 */ 84 IntegerVector temp(A.getWidth()); 85 bool ret=!positiveVectorInKernel(A.getRows(),&temp); 86 //debug<<A.getRows()<<temp; 87 log1 debug<<(ret?"VECTOR CONFIGURATION IS \"HOMOGENEOUS\"\n":"VECTOR CONFIGURATION IS NOT \"HOMOGENEOUS\"\n"); 88 return ret; 89 } 90 #endif 91 /** 92 * The support of the secondary fan is given by {x:{y:y^TA<=x^T}!=emptyset}. 93 * Hence the support is the projection of {(x,y):y^TA<=x^T} onto the X component. 94 * This function computes the support which is a polyhedral cone. 95 */ secondaryFanSupport()96 PolyhedralCone secondaryFanSupport() 97 { 98 // if(!isHomogeneous())<---not working 99 // if(!positiveVectorInKernel(A.getRows(),0)) 100 { 101 FieldMatrix A2=integerMatrixToFieldMatrix(A,Q); 102 IntegerVectorList empty; 103 PolyhedralCone C(fieldMatrixToIntegerMatrix(combineLeftRight(FieldMatrix::identity(Q,A2.getWidth()),A2.transposed())).getRows(),empty,A.getHeight()+A.getWidth()); 104 C.canonicalize(); 105 return C.projection(A2.getWidth()); 106 } 107 /* else 108 return PolyhedralCone(A.getWidth());*/ 109 } secondaryCone()110 PolyhedralCone secondaryCone()const 111 { 112 IntegerVectorList empty; 113 return PolyhedralCone(facets(),empty,A.getWidth(),PCP_impliedEquationsKnown); 114 } interiorPoint()115 IntegerVector interiorPoint()const 116 { 117 // IntegerVectorList empty; 118 //return PolyhedralCone(inequalities(),empty,A.getWidth()).getRelativeInteriorPoint(); 119 //return PolyhedralCone(facets(),empty,A.getWidth()).getRelativeInteriorPoint(); 120 return secondaryCone().getRelativeInteriorPoint(); 121 } 122 123 determinantInequality3(list<int> const & indices)124 IntegerVector determinantInequality3(list<int> const &indices) const 125 { 126 IntegerMatrix AA(A.getHeight(),indices.size()); 127 int i=0; 128 for(list<int>::const_iterator I=indices.begin();I!=indices.end();I++,i++) 129 { 130 for(int j=0;j<A.getHeight();j++) 131 AA[j][i]=A[j][*I]; 132 } 133 FieldMatrix AAA=integerMatrixToFieldMatrix(AA,Q); 134 IntegerVector temp=AAA.reduceAndComputeVectorInKernel().primitive(); 135 IntegerVector u(A.getWidth()); 136 i=0; 137 for(list<int>::const_iterator I=indices.begin();I!=indices.end();I++,i++) 138 u[*I]=-temp[i]; 139 return u; 140 } 141 determinantInequality2(list<int> const & indices)142 IntegerVector determinantInequality2(list<int> const &indices) const 143 { 144 IntegerMatrix AA(A.getHeight(),indices.size()); 145 int i=0; 146 for(list<int>::const_iterator I=indices.begin();I!=indices.end();I++,i++) 147 { 148 for(int j=0;j<A.getHeight();j++) 149 AA[j][i]=A[j][*I]; 150 } 151 IntegerVector v=vectorInKernel(AA); 152 IntegerVectorList V=AA.getRows(); 153 V.push_back(v); 154 IntegerVector u(A.getWidth()); 155 i=0; 156 for(list<int>::const_iterator I=indices.begin();I!=indices.end();I++,i++) 157 u[*I]=v[i]; 158 if(determinantSign(V)==1) 159 return -u; 160 return u; 161 } 162 163 determinantInequality1(list<int> const & indices)164 IntegerVector determinantInequality1(list<int> const &indices) const 165 { 166 // fprintf(stderr,"SET:"); 167 // for(list<int>::const_iterator i=indices.begin();i!=indices.end();i++)fprintf(stderr," %i",*i); 168 // fprintf(stderr,"\n"); 169 int m=indices.size(); 170 FieldVector ret(Q,A.getWidth()); 171 172 int I=m;//HER 173 for(list<int>::const_iterator i=indices.begin();i!=indices.end();i++,I++) 174 { 175 IntegerMatrix B(m-1,A.getHeight()); 176 int J=0; 177 for(list<int>::const_iterator j=indices.begin();j!=indices.end();j++) 178 if(i!=j) 179 { 180 for(int k=0;k<B.getHeight();k++) 181 { 182 B[k][J]=A[k][*j]; 183 } 184 J++; 185 } 186 FieldMatrix B2=integerMatrixToFieldMatrix(B,Q); 187 ret[*i]=(Q.zHomomorphism(1-2*(I&1)))*B2.reduceAndComputeDeterminant(); 188 } 189 // AsciiPrinter T(Stderr); ret.print(T); 190 return ret.primitive(); 191 } determinantInequality(list<int> const & indices)192 IntegerVector determinantInequality(list<int> const &indices) const 193 { 194 return determinantInequality3(indices); 195 196 IntegerVector r1=determinantInequality1(indices); 197 IntegerVector r2=determinantInequality2(indices); 198 IntegerVector r3=determinantInequality3(indices); 199 debug<<r1<<"\n"<<r2<<"\n"<<r3<<"\n"; 200 return r3; 201 } subsetRows(IntegerMatrix const & ATransposed,IntegerVector const & cols)202 IntegerMatrix subsetRows(IntegerMatrix const &ATransposed, IntegerVector const &cols)const 203 { 204 IntegerMatrix ret(cols.size(),ATransposed.getWidth()); 205 206 for(int i=0;i!=cols.size();i++)ret[i]=ATransposed[cols[i]]; 207 return ret; 208 } volume(IntegerVector const & v,IntegerMatrix const & ATransposed)209 FieldElement volume(IntegerVector const &v, IntegerMatrix const &ATransposed)const 210 { 211 IntegerMatrix B(A.getHeight(),A.getHeight()); 212 for(int j=0;j<v.size();j++) 213 B[j]=ATransposed[v[j]]; 214 FieldMatrix B2=integerMatrixToFieldMatrix(B,Q); 215 FieldElement det=B2.reduceAndComputeDeterminant(); 216 if(det.sign()<0)return -det; 217 return det; 218 } coDimensionOneTriangles()219 set<set<int> > coDimensionOneTriangles()const 220 { 221 set<set<int> > codimOne; 222 for(set<IntegerVector>::const_iterator j=bases.begin();j!=bases.end();j++) 223 { 224 set<int> temp; 225 for(int i=0;i<j->size();i++) 226 { 227 temp.insert((*j)[i]); 228 } 229 for(int i=0;i<j->size();i++) 230 { 231 temp.erase((*j)[i]); 232 codimOne.insert(temp); 233 temp.insert((*j)[i]); 234 } 235 } 236 return codimOne; 237 } edgeGraph()238 Graph edgeGraph()const 239 { 240 Graph ret(bases.size()); 241 242 set<set<int> > codimOne=coDimensionOneTriangles(); 243 244 for(set<set<int> >::const_iterator i=codimOne.begin();i!=codimOne.end();i++) 245 { 246 list<int> triangIndices; 247 int J=0; 248 for(set<IntegerVector>::const_iterator j=bases.begin();j!=bases.end();j++,J++) 249 { 250 bool isSuperSet=true; 251 for(set<int>::const_iterator k=i->begin();k!=i->end();k++) 252 { 253 bool hasBeenFound=false; 254 for(int l=0;l<j->size();l++) 255 if((*j)[l]==*k) 256 { 257 hasBeenFound=true; 258 break; 259 } 260 if(!hasBeenFound) 261 { 262 isSuperSet=false; 263 break; 264 } 265 } 266 if(isSuperSet) 267 triangIndices.push_back(J); 268 } 269 if(triangIndices.size()==2) 270 ret.addEdge(triangIndices.front(),triangIndices.back()); 271 } 272 273 // cerr << ret.toString(); 274 275 return ret; 276 } difference(IntegerVector const & v,set<int> const & s)277 set<int> difference(IntegerVector const &v, set<int> const &s)const 278 { 279 set<int> ret; 280 281 for(int i=0;i<v.size();i++)if(s.count(v[i])==0)ret.insert(v[i]); 282 return ret; 283 } 284 /** Computes v setminus s - or actually it returns true if this set has size 1 and false other wise. 285 * In case of true the unique element is stored in theDifferencs. 286 * The routine assumes that the number of elements in v is one larger than that of s. 287 */ differenceOne(IntegerVector const & v,set<int> const & s,int & theDifference)288 bool differenceOne(IntegerVector const &v, set<int> const &s, int &theDifference)const 289 { 290 int diffSize=0; 291 for(int i=0;i<v.size();i++) 292 if(s.count(v[i])==0) 293 { 294 diffSize++; 295 if(diffSize==2)return false; 296 theDifference=v[i]; 297 } 298 return diffSize==1; 299 } 300 /** 301 Seems to compute the outer normals of the secondary cone of the triangulation. 302 */ 303 inequalitiesFast()304 IntegerVectorList inequalitiesFast()const 305 { 306 IntegerVectorList ret; 307 308 // cerr<<"bases.size"<<bases.size()<<"A.height"<<A.getHeight()<<endl; 309 310 // we get an inequality for every A-column not a vertex, i.e. not appearing in the triangulation 311 for(int i=0;i<A.getWidth();i++) 312 { 313 bool appears=false; 314 for(set<IntegerVector>::const_iterator j=bases.begin();j!=bases.end();j++) 315 { 316 for(int k=0;k<j->size();k++) 317 if(i==(*j)[k]){appears=true;goto done;} 318 } 319 done: 320 if(!appears) 321 { 322 // We now find a triangle which covers the ith column 323 set<IntegerVector>::const_iterator j; 324 int lastSign=0; 325 for(j=bases.begin();j!=bases.end();j++) 326 { 327 IntegerMatrix ATransposed=Atransposed; 328 bool containsI=true; 329 IntegerMatrix subMatrix=subsetRows(Atransposed,*j); 330 int sign=determinantSign(subMatrix.getRows()); 331 for(int k=0;k<j->size();k++) 332 { 333 subMatrix[k]=ATransposed[i]; 334 if(sign*determinantSign(subMatrix.getRows())<0) 335 { 336 containsI=false; 337 break; 338 } 339 subMatrix[k]=ATransposed[(*j)[k]]; 340 } 341 lastSign=sign; 342 if(containsI)break; 343 } 344 if(j==bases.end()) 345 { 346 AsciiPrinter(Stderr).printVector(Atransposed[i]); 347 } 348 assert(j!=bases.end()); 349 list<int> temp; 350 for(int k=0;k<j->size();k++)temp.push_back((*j)[k]); 351 temp.push_back(i); 352 ret.push_back(-lastSign* determinantInequality(temp));//Do we need to use the sign here?//HERE//MAY4 353 } 354 } 355 356 357 // we get an inequality for every codim one simplex of the triangulation, i.e. every edge of the polyhedron 358 359 // TODO: the call to determinantSign below can be avoided if we carefully keep track of orientation 360 // - that is, the bases should be stored with an orientation flag just as it is done in the trinaglation class. 361 // debug<<">\n"; 362 #if 0 363 { 364 set<set<int> > codimOne; 365 for(set<IntegerVector>::const_iterator j=bases.begin();j!=bases.end();j++) 366 { 367 set<int> temp; 368 for(int i=0;i<j->size();i++) 369 { 370 temp.insert((*j)[i]); 371 } 372 for(int i=0;i<j->size();i++) 373 { 374 temp.erase((*j)[i]); 375 codimOne.insert(temp); 376 temp.insert((*j)[i]); 377 } 378 } 379 for(set<set<int> >::const_iterator i=codimOne.begin();i!=codimOne.end();i++) 380 { 381 list<int> additional; 382 for(set<IntegerVector>::const_iterator j=bases.begin();j!=bases.end();j++) 383 { 384 int theDifference; 385 if(differenceOne(*j,*i,theDifference)) 386 { 387 additional.push_back(theDifference); 388 } 389 } 390 if(additional.size()==2) 391 { 392 list<int> temp; 393 for(set<int>::const_iterator j=i->begin();j!=i->end();j++)temp.push_back(*j); 394 list<int>::const_iterator a=additional.begin(); 395 396 list<int> temp2=temp; 397 temp2.push_back(*a); 398 IntegerVector temp3(temp2.size()); 399 list<int>::const_iterator K=temp2.begin(); 400 for(int k=0;k<temp3.size();k++,K++)temp3[k]=*K; 401 if(determinantSign(subsetRows(Atransposed,temp3).getRows())>0) 402 { 403 temp.push_back(*a); 404 a++; 405 temp.push_back(*a); 406 } 407 else 408 { 409 a++; 410 temp.push_back(*a); 411 a--; 412 temp.push_back(*a); 413 } 414 ret.push_back(-determinantInequality(temp));//HERE//MAY4 415 } 416 } 417 } 418 #else 419 { 420 multimap<IntegerVector,pair<int,bool> > codimOne; 421 for(set<IntegerVector>::const_iterator j=bases.begin();j!=bases.end();j++) 422 { 423 bool orientation=0>determinantSign(subsetRows(Atransposed,*j).getRows()); 424 IntegerVector temp=j->subvector(1,j->size()); 425 for(int i=0;i<j->size();i++) 426 { 427 if(i)temp[i-1]=(*j)[i-1]; 428 IntegerVector temp2=temp; 429 int nswaps=mergeSort(temp2)+orientation; 430 codimOne.insert(pair<IntegerVector,pair<int,bool> >(temp2,pair<int,bool>((*j)[i],(nswaps+i+j->size())&1))); 431 } 432 } 433 for(multimap<IntegerVector,pair<int,bool> >::const_iterator i=codimOne.begin();i!=codimOne.end();i++) 434 { 435 multimap<IntegerVector,pair<int,bool> >::const_iterator next=i;next++; 436 if(next==codimOne.end())break; 437 if(i->first == next->first) 438 { 439 list<int> temp; 440 IntegerVector const &v=i->first; 441 for(int j=0;j<v.size();j++)temp.push_back(v[j]); 442 if(i->second.second) 443 { 444 temp.push_back(i->second.first); 445 temp.push_back(next->second.first); 446 } 447 else 448 { 449 temp.push_back(next->second.first); 450 temp.push_back(i->second.first); 451 } 452 ret.push_back(-determinantInequality(temp));//MAY4 453 i++; 454 } 455 } 456 } 457 #endif 458 return ret; 459 } inequalities()460 IntegerVectorList inequalities()const 461 { 462 int n=A.getWidth(); 463 IntegerVectorList ret; 464 for(set<IntegerVector>::const_iterator i=bases.begin();i!=bases.end();i++) 465 { 466 IntegerMatrix A2transposed(i->size()+1,A.getHeight()); 467 for(int j=0;j<i->size();j++) 468 { 469 A2transposed[j]=A.column((*i)[j]); 470 } 471 IntegerVector iComplement=complement(*i,n); 472 473 for(int k=0;k<iComplement.size();k++) 474 { 475 IntegerVector nyC; 476 477 A2transposed[i->size()]=A.column(iComplement[k]); 478 IntegerMatrix temp=A2transposed.transposed(); 479 nyC=vectorInKernel(temp); 480 481 482 IntegerVector c(n); 483 for(int j=0;j<i->size();j++) 484 c[(*i)[j]]=nyC[j]; 485 c[iComplement[k]]=nyC[i->size()]; 486 if(c[iComplement[k]]>0) 487 { 488 ret.push_back(c); 489 } 490 else 491 { 492 ret.push_back(-c); 493 } 494 } 495 } 496 return ret; 497 } totalVolume()498 int totalVolume()const 499 { 500 FieldElement s(Q); 501 for(set<IntegerVector>::const_iterator i=bases.begin();i!=bases.end();i++) 502 { 503 FieldElement vol=volume(*i,Atransposed); 504 s=s+vol; 505 } 506 cerr<<"retateawtat"<< s.toString()<<endl; 507 return toInteger(s); 508 } 509 /** 510 * Using the formula of Theorem 6.4 in [Dickenstein, Feichtner, Sturmfels] this computes this 511 * compute the exponent vector of the monomial in the resultant. Is this different from the 512 * vertex of the secondary polytope???!! 513 */ DFSResultantCoordinate()514 IntegerVector DFSResultantCoordinate()const 515 { 516 IntegerVector ret(A.getWidth()); 517 for(set<IntegerVector>::const_iterator i=bases.begin();i!=bases.end();i++) 518 { 519 FieldElement vol=volume(*i,Atransposed); 520 for(int j=0;j<i->size();j++)ret[(*i)[j]]+=toInteger(vol); 521 } 522 return ret; 523 } print(Printer & p)524 void print(Printer &p) 525 { 526 p.printString("Bases("); 527 p.printInteger(bases.size()); 528 p.printString(":\n"); 529 FieldElement s(Q); 530 for(set<IntegerVector>::const_iterator i=bases.begin();i!=bases.end();i++) 531 { 532 p.printVector(*i); 533 FieldElement vol=volume(*i,Atransposed); 534 p.printString(" Vol: "); 535 p.printFieldElement(vol); 536 p.printNewLine(); 537 s=s+vol; 538 } 539 p.printString("Total volume: "); 540 p.printFieldElement(s); 541 p.printNewLine(); 542 /* p.printString("Circuits (inequalities):\n"); 543 p.printVectorList(inequalities()); 544 p.printString("Facets:\n"); 545 p.printVectorList(facets()); 546 p.printString("Interior point:\n"); 547 p.printVector(interiorPoint()); 548 */ 549 } facets()550 IntegerVectorList facets()const 551 { 552 return fastNormals(inequalitiesFast()); 553 554 IntegerVectorList flipable; 555 IntegerVectorList normals=wallRemoveScaledInequalities(inequalitiesFast()); 556 // IntegerVectorList normals=wallRemoveScaledInequalities(inequalities()); 557 558 559 560 for(IntegerVectorList::iterator i=normals.begin();i!=normals.end();i++) 561 { 562 // if(!termOrder(*i,*i-*i)) 563 { 564 if(isFacet(normals,i)) 565 { 566 // if(wallContainsPositiveVector(*i)) 567 flipable.push_back(*i); 568 } 569 else 570 { 571 IntegerVectorList::iterator temp=i; 572 temp++; 573 normals.erase(i); 574 temp--; 575 i=temp; 576 } 577 } 578 } 579 /* AsciiPrinter Q(Stderr); 580 Q.printString("Bases:\n"); 581 for(set<IntegerVector>::const_iterator i=bases.begin();i!=bases.end();i++) 582 { 583 Q.printVector(*i); 584 Q.printNewLine(); 585 } 586 587 AsciiPrinter(Stderr).printVectorList(inequalities()); 588 AsciiPrinter(Stderr).printVectorList(flipable); 589 AsciiPrinter(Stderr).printVectorList(inequalitiesFast()); 590 log0 fprintf(stderr,"-----------------\n"); 591 */ 592 593 //log0 fprintf(stderr,"%i\n",0); 594 // log0 fprintf(stderr,"Number of facets:%i\n",flipable.size()); 595 596 return flipable; 597 } isEmpty()598 bool isEmpty()const 599 { 600 return bases.empty(); 601 } 602 /* void flip(IntegerVector const &normal) 603 { 604 AsciiPrinter P(Stderr); 605 gfan_log2 print(P); 606 //log0 P.printVector(normal); 607 int n=normal.size(); 608 // IntegerVectorList l=wallRemoveScaledInequalities(inequalities());// This is not needed - one circuit should be enough 609 IntegerVectorList l;l.push_back(normal);// Let's do this instead 610 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++) 611 if(dependent(*i,normal)) 612 { 613 gfan_log2 AsciiPrinter(Stderr).printVector(*i); 614 for(int k=0;k<normal.size();k++) 615 if((*i)[k]<0) 616 { 617 int s1=bases.size(); 618 IntegerVector i2=*i; 619 i2[k]=0; 620 bases.insert(i2.supportIndices()); 621 int s2=bases.size(); 622 assert(s2==s1+1); 623 } 624 else if((*i)[k]>0) 625 { 626 int s1=bases.size(); 627 IntegerVector i2=*i; 628 i2[k]=0; 629 bases.erase(i2.supportIndices()); 630 int s2=bases.size(); 631 assert(s2==s1-1); 632 } 633 } 634 }*/ toSet(IntegerVector const & v)635 static set<int> toSet(IntegerVector const &v) 636 { 637 set<int> ret; 638 for(int i=0;i<v.size();i++)ret.insert(v[i]); 639 return ret; 640 } toVector(set<int> const & s)641 static IntegerVector toVector(set<int> const &s) 642 { 643 IntegerVector ret(s.size()); 644 int I=0; 645 for(set<int>::const_iterator i=s.begin();i!=s.end();i++,I++)ret[I]=*i; 646 return ret; 647 } contains(IntegerVector const & v,set<int> const & s)648 bool contains(IntegerVector const &v, set<int> const &s) 649 { 650 for(set<int>::const_iterator i=s.begin();i!=s.end();i++) 651 { 652 bool found=false; 653 for(int j=0;j!=v.size();j++) 654 if(v[j]==*i) 655 { 656 found=true; 657 break; 658 } 659 if(!found)return false; 660 } 661 return true; 662 } uni(set<int> const & a,set<int> const & b)663 set<int> uni(set<int> const &a, set<int> const &b) 664 { 665 set<int> ret=a; 666 for(set<int>::const_iterator i=b.begin();i!=b.end();i++)ret.insert(*i); 667 return ret; 668 } printSet(set<int> const & s)669 void printSet(set<int> const &s) 670 { 671 cerr<<"{"; 672 for(set<int>::const_iterator i=s.begin();i!=s.end();i++) 673 { 674 if(i!=s.begin())cerr<<","; 675 cerr<<*i; 676 } 677 cerr<<"}"; 678 } printSetSet(set<set<int>> const & s)679 void printSetSet(set<set<int> > const &s) 680 { 681 cerr<<"{"<<endl; 682 for(set<set<int> >::const_iterator i=s.begin();i!=s.end();i++) 683 { 684 if(i!=s.begin())cerr<<endl<<","; 685 printSet(*i); 686 } 687 cerr<<"}"<<endl; 688 } 689 flipNew(IntegerVector const & normal)690 void flipNew(IntegerVector const &normal) 691 { 692 set<set<int> > toBeInserted; 693 set<set<int> > toBeErased; 694 int n=normal.size(); 695 set<int> nonZeroIndices; 696 for(int i=0;i<n;i++) 697 if(normal[i])nonZeroIndices.insert(i); 698 699 for(set<int>::const_iterator i=nonZeroIndices.begin();i!=nonZeroIndices.end();i++) 700 { 701 set<int> temp=nonZeroIndices; 702 temp.erase(*i); 703 if(normal[*i]>0) 704 toBeErased.insert(temp); 705 else 706 toBeInserted.insert(temp); 707 } 708 709 assert(toBeErased.size()); 710 711 /* cerr<<"To be inserted"<<endl; 712 printSetSet(toBeInserted); 713 714 cerr<<"ToBeErased"<<endl; 715 printSetSet(toBeErased); 716 */ 717 718 set<int> A=*toBeErased.begin(); 719 set<set<int> > link; 720 for(set<IntegerVector>::const_iterator i=bases.begin();i!=bases.end();i++) 721 if(contains(*i,A)) 722 link.insert(difference(*i,A)); 723 724 for(set<IntegerVector>::iterator i=bases.begin();i!=bases.end();) 725 { 726 bool e=false; 727 for(set<set<int> >::const_iterator j=toBeErased.begin();j!=toBeErased.end();j++) 728 if(contains(*i,*j)) 729 { 730 e=true; 731 break; 732 } 733 if(e) 734 { 735 set<IntegerVector>::iterator I=i; 736 I++; 737 bases.erase(i); 738 i=I; 739 } 740 else 741 i++; 742 } 743 744 for(set<set<int> >::iterator i=link.begin();i!=link.end();i++) 745 { 746 for(set<set<int> >::iterator j=toBeInserted.begin();j!=toBeInserted.end();j++) 747 { 748 bases.insert(toVector(uni(*i,*j))); 749 } 750 } 751 } usedRays()752 list<int> usedRays()const 753 { 754 list<int> ret; 755 for(int i=0;i<A.getWidth();i++) 756 { 757 bool contains=false; 758 for(set<IntegerVector>::const_iterator j=bases.begin();j!=bases.end();j++) 759 { 760 for(int k=0;k<j->size();k++) 761 { 762 if((*j)[k]==i) 763 { 764 contains=true; 765 goto leave; 766 } 767 } 768 } 769 leave: 770 if(contains)ret.push_back(i); 771 } 772 return ret; 773 } hirschScore()774 float hirschScore()const 775 { 776 int timesAttained; 777 int nVertices=bases.size(); 778 int nEdges=coDimensionOneTriangles().size(); 779 int diameter=edgeGraph().diameter(×Attained); 780 int dimension=A.getHeight(); 781 int nFacets=usedRays().size(); 782 783 if(diameter+dimension-nFacets>0) 784 { 785 cerr<<"Counter example found\n"; 786 // print(); 787 assert(0); 788 } 789 790 return diameter+dimension-nFacets+(nFacets>120?(-0.01*nFacets):0); 791 } changeToTriangulationInducedBy(TermOrder const & T)792 void changeToTriangulationInducedBy(TermOrder const &T) 793 { 794 //TODO: use ray shooting to find facet instead 795 // log0 debug<<"changing\n"; 796 while(1) 797 { 798 IntegerVectorList l=facets(); 799 bool found=false; 800 for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++) 801 { 802 if(!T(-*i,*i-*i)) 803 { 804 // debug<<*i; 805 flipNew(*i); 806 807 found=true; 808 break; 809 } 810 } 811 if(!found)break; 812 } 813 // log0 debug<<"done changing\n"; 814 } makeConfigurationFullDimensional(IntegerMatrix & A)815 static void makeConfigurationFullDimensional(IntegerMatrix &A) 816 { 817 // debug<<"C1\n"; 818 /* If the vector configuration A does not have full rank then 819 change coordinates. */ 820 if(::rank(A)!=A.getHeight()) 821 { 822 // debug<<"C2\n"; 823 FieldMatrix M=integerMatrixToFieldMatrix(A,Q); 824 M.reduce(false,true);//force integer operations - preserving volume 825 M.removeZeroRows(); 826 A=fieldMatrixToIntegerMatrix(M); 827 } 828 } triangulate()829 void triangulate() 830 { 831 /* Convert a Triangulation to a Triangulation2 */ 832 bases=set<IntegerVector>(); 833 list<Triangulation::Cone> T=Triangulation::triangulate(A.transposed()); 834 for(list<Triangulation::Cone>::const_iterator i=T.begin();i!=T.end();i++) 835 { 836 IntegerVector v=i->size(); 837 int J=0; 838 for(Triangulation::Cone::const_iterator j=i->begin();j!=i->end();j++,J++) 839 v[J]=*j; 840 bases.insert(v); 841 } 842 } 843 }; 844 845 #endif 846