1 #include "parser.h" 2 #include "printer.h" 3 #include "polynomial.h" 4 #include "division.h" 5 #include "buchberger.h" 6 #include "wallideal.h" 7 #include "lp.h" 8 #include "reversesearch.h" 9 #include "breadthfirstsearch.h" 10 #include "termorder.h" 11 #include "ep_standard.h" 12 #include "ep_xfig.h" 13 #include "gfanapplication.h" 14 #include "timer.h" 15 #include "log.h" 16 #include "matrix.h" 17 #include "lll.h" 18 #include "polyhedralfan.h" 19 #include "linalg.h" 20 #include "determinant.h" 21 #include "triangulation.h" 22 #include "intsinpolytope.h" 23 #include "graph.h" 24 25 #include "triangulation2.h" 26 27 #include "traverser_secondaryfan.h" 28 #include "symmetrictraversal.h" 29 30 #include <iostream> 31 #include <algorithm> 32 33 #define NOOUTPUT 0 34 35 class SecondaryFanApplication : public GFanApplication 36 { 37 SimpleOption hirschOption; 38 SimpleOption searchOption; 39 IntegerOption scaleOption; 40 StringOption optionRestrictingFan; 41 SimpleOption symmetryOption; 42 SimpleOption optionIgnoreCones; 43 InterruptOption optionInterrupt; 44 public: helpText()45 const char *helpText() 46 { 47 return "This program computes the secondary fan of a vector configuration. The configuration is given as an ordered list of vectors. In order to compute the secondary fan of a point configuration an additional coordinate of ones must be added. For example {(1,0),(1,1),(1,2),(1,3)}.\n"; 48 } SecondaryFanApplication()49 SecondaryFanApplication(): 50 searchOption("--unimodular","Use heuristics to search for unimodular triangulation rather than computing the complete secondary fan"), 51 scaleOption("--scale","Assuming that the first coordinate of each vector is 1, this option will take the polytope in the 1 plane and scale it. The point configuration will be all lattice points in that scaled polytope. The polytope must have maximal dimension. When this option is used the vector configuration must have full rank. This option may be removed in the future."), 52 symmetryOption("--symmetry","Tells the program to read in generators for a group of symmetries (subgroup of $S_n$) after having read in the vector configuration. The program checks that the configuration stays fixed when permuting the variables with respect to elements in the group. The output is grouped according to the symmetry.\n"), 53 optionRestrictingFan("--restrictingfan","Specify the name of a file containing a polyhedral fan in Polymake format. The computation of the Secondary fan will be restricted to this fan. If the --symmetry option is used then this restricting fan must be invariant under the symmetry and the orbits in the file must be with respect to the specified group of symmetries. The orbits of maximal cones of the file are then read in rather than the maximal cones.\n",0), 54 optionIgnoreCones("--nocones","Tells the program not to output the CONES and MAXIMAL_CONES sections, but still output CONES_COMPRESSED and MAXIMAL_CONES_COMPRESSED if --symmetry is used."), 55 hirschOption("--hirsch","") 56 { 57 hirschOption.hide(); 58 registerOptions(); 59 } 60 name()61 const char *name() 62 { 63 return "_secondaryfan"; 64 } 65 enumerate(Triangulation2 const & t)66 PolyhedralFan enumerate(Triangulation2 const &t) 67 { 68 PolyhedralFan ret(t.getN()); 69 list<Triangulation2> active; 70 active.push_back(t); 71 IntegerVectorList interiorPoints; 72 interiorPoints.push_back(t.interiorPoint()); 73 while(!active.empty()) 74 { 75 Triangulation2 a=active.front(); 76 PolyhedralCone C=a.secondaryCone(); 77 78 79 // if(active.size()>100)break;//SLETMIGGGGG 80 //log0 fprintf(stderr,"a\n"); 81 /* { 82 PolyhedralCone C2=C; 83 C2.canonicalize(); 84 }*/ 85 C.canonicalize(); 86 //log0 fprintf(stderr,"b\n"); 87 ret.insert(C); 88 AsciiPrinter P(Stderr); 89 // C.print(&P); 90 active.pop_front(); 91 // fprintf(stderr,"pop\n"); 92 IntegerVectorList flips=a.facets(); 93 for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++) 94 { 95 { 96 IntegerVectorList t=C.getEquations(); 97 t.push_back(*i); 98 PolyhedralCone CF(C.getHalfSpaces(),t); 99 CF.findFacets(); 100 // CF.canonicalize(); 101 } 102 103 if(!i->isNonNegative()) //is this the right condition or should i be negated? 104 // if(!(-*i).isNonNegative()) //is this the right condition or should i be negated? 105 { 106 Triangulation2 b=a; 107 log1 AsciiPrinter(Stderr)<<*i; 108 /*fprintf(stderr,"Before:"); 109 b.print(P);*/ 110 // b.flip(*i); 111 b.flipNew(-*i); 112 /*fprintf(stderr,"After:"); 113 b.print(P);*/ 114 if(!b.isEmpty()) 115 { 116 IntegerVectorList inequalities=b.inequalities(); 117 bool isKnown=false; 118 for(IntegerVectorList::const_iterator j=interiorPoints.begin();j!=interiorPoints.end();j++) 119 { 120 bool match=true; 121 for(IntegerVectorList::const_iterator k=inequalities.begin();k!=inequalities.end();k++) 122 { 123 if(dotLong(-*k,*j)<0) 124 { 125 match=false; 126 break; 127 } 128 } 129 if(match)isKnown=true; 130 } 131 if(!isKnown) 132 { 133 active.push_back(b); 134 interiorPoints.push_back(b.interiorPoint()); 135 } 136 } 137 } 138 } 139 } 140 return ret; 141 } interactive(Triangulation2 const & t)142 PolyhedralFan interactive(Triangulation2 const &t) 143 { 144 Triangulation2 a=t; 145 146 while(1) 147 { 148 fprintf(stdout,"Triangles in current triangulation:%i\n",(int)a.bases.size()); 149 PolyhedralCone C=a.secondaryCone(); 150 151 C.canonicalize(); 152 153 154 AsciiPrinter Pstd(Stderr); 155 IntegerVectorList flips=a.facets(); 156 int I=0; 157 for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++,I++) 158 { 159 fprintf(stdout,"%i:\n",I); 160 Pstd.printVector(*i); 161 162 163 if(!i->isNonNegative()) 164 { 165 Triangulation2 b=a; 166 //fprintf(stderr,"Before:"); 167 //b.print(P); 168 // b.flip(*i); 169 /* b.flipNew(*i); 170 //fprintf(stderr,"After:"); 171 //b.print(P); 172 fprintf(stdout,"Triangles in new triangulation:%i\n",b.bases.size()); 173 */ 174 } 175 176 fprintf(stdout,"\n"); 177 } 178 int s; 179 //cin >> s; 180 int err=fscanf(stdin,"%i",&s); 181 assert(err==1); 182 if((s>=0)&&(s<I)) 183 { 184 IntegerVectorList::const_iterator i=flips.begin(); 185 for(int I=0;I<s;I++)i++; 186 a.flipNew(*i); 187 } 188 else 189 { 190 a.print(Pstd); 191 } 192 } 193 PolyhedralFan ret(0); 194 195 return ret; 196 } automatic(Triangulation2 const & t,int abortAtSize)197 PolyhedralFan automatic(Triangulation2 const &t, int abortAtSize) 198 { 199 Triangulation2 a=t; 200 201 cout << "Looking for triangulation with "<<abortAtSize<<" simplices."<<endl; 202 203 while(1) 204 { 205 fprintf(stdout,"Triangles in current triangulation:%i\n",(int)a.bases.size()); 206 // PolyhedralCone C=a.secondaryCone(); 207 208 // C.canonicalize(); 209 210 211 AsciiPrinter Pstd(Stderr); 212 IntegerVectorList flips=a.facets(); 213 int I=0; 214 cerr << "Number of facets of secondary cone: " << flips.size() <<endl; 215 for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++,I++) 216 { 217 if(!i->isNonNegative()) 218 { 219 Triangulation2 b=a; 220 b.flipNew(-*i); 221 fprintf(stdout,"Triangles in new triangulation:%i\n",(int)b.bases.size()); 222 223 if(b.bases.size()==abortAtSize) 224 { 225 b.print(Pstd); 226 227 exit(0); 228 } 229 230 if((b.bases.size()>a.bases.size())||((rand()&127)==0)) 231 { 232 a=b; 233 break; 234 } 235 } 236 } 237 } 238 PolyhedralFan ret(0); 239 240 return ret; 241 } automaticHirsch(Triangulation2 const & t)242 PolyhedralFan automaticHirsch(Triangulation2 const &t) 243 { 244 Triangulation2 a=t; 245 while(1) 246 { 247 int nVertices=a.bases.size(); 248 int nEdges=a.coDimensionOneTriangles().size(); 249 int diameter=a.edgeGraph().diameter(); 250 int dimension=a.getD(); 251 int nFacets=a.usedRays().size(); 252 fprintf(stdout,"NVER: %i NEDGES: %i DIAMETER:%i DIMENSION:%i NFACETS:%i\n",nVertices,nEdges,diameter,dimension,nFacets); 253 254 AsciiPrinter Pstd(Stderr); 255 IntegerVectorList flips=a.facets(); 256 int I=0; 257 float currentScore=a.hirschScore(); 258 cerr << "Current score: " << currentScore <<endl; 259 for(IntegerVectorList::const_iterator i=flips.begin();i!=flips.end();i++,I++) 260 { 261 if(!i->isNonNegative()) 262 { 263 Triangulation2 b=a; 264 b.flipNew(-*i); 265 float bScore=b.hirschScore(); 266 fprintf(stdout,"New score:%f\n",bScore); 267 268 if((bScore>currentScore)||((rand()&31)==0)) 269 { 270 a=b; 271 break; 272 } 273 } 274 } 275 } 276 PolyhedralFan ret(0); 277 278 return ret; 279 } main()280 int main() 281 { 282 283 284 // debug<<(int)sizeof(SymmetricComplex::Cone); 285 // debug<<(int)sizeof(IntegerVector); 286 // debug<<(int)sizeof(vector<int>); 287 IntegerMatrix A=rowsToIntegerMatrix(FileParser(Stdin).parseIntegerVectorList()).transposed(); 288 289 int n=A.getWidth(); 290 291 SymmetryGroup s(n); 292 if(symmetryOption.getValue()) 293 { 294 IntegerVectorList generators=FileParser(Stdin).parseIntegerVectorList(); 295 for(IntegerVectorList::const_iterator i=generators.begin();i!=generators.end();i++) 296 { 297 assert(i->size()==n); 298 FieldMatrix M1=integerMatrixToFieldMatrix(A,Q); 299 FieldMatrix M2=integerMatrixToFieldMatrix(rowsToIntegerMatrix(SymmetryGroup::permuteIntegerVectorList(A.getRows(),*i)),Q); 300 M1.reduce(); 301 M1.REformToRREform(true); 302 M2.reduce(); 303 M2.REformToRREform(true); 304 305 if(!(M1==M2)) 306 { 307 AsciiPrinter(Stderr) << "Permutation "<< *i << 308 " does not keep the configuration fixed.\n"; 309 assert(0); 310 } 311 } 312 s.computeClosure(generators); 313 s.createTrie(); 314 } 315 316 317 if(scaleOption.getValue()) 318 { 319 if(::rank(A)!=A.getHeight()) 320 { 321 cerr << "The vector configuration must have full rank in order to use the scale option.\n"; 322 assert(0); 323 } 324 325 int s=scaleOption.getValue(); 326 327 cout << "Input configuration:" << endl; 328 AsciiPrinter(Stdout)<<A.transposed().getRows(); 329 330 IntegerVectorList empty; 331 PolyhedralCone dual(A.transposed().getRows(),empty,n); 332 dual.canonicalize(); 333 assert(dual.dimensionOfLinealitySpace()==0); 334 IntegerVectorList inequalities=dual.extremeRays(); 335 336 IntegerMatrix M1=rowsToIntegerMatrix(inequalities); 337 IntegerMatrix M2=-1*M1.submatrix(0,1,M1.getHeight(),M1.getWidth()); 338 IntegerVector rightHandSide=s*M1.transposed().submatrix(0,0,1,M1.getWidth())[0].toVector(); 339 IntegerVector v=s*A.submatrix(1,0,A.getHeight(),1).transposed()[0].toVector(); 340 341 IntegerVectorList l=intsInPolytopeGivenIneqAndPt(M2,rightHandSide,v); 342 343 IntegerVectorList lT=rowsToIntegerMatrix(l).transposed().getRows(); 344 lT.push_front(IntegerVector::allOnes(l.size())); 345 l=rowsToIntegerMatrix(lT).transposed().getRows(); 346 347 cout << "New configuration:" << endl; 348 AsciiPrinter(Stdout)<<l; 349 A=rowsToIntegerMatrix(l).transposed(); 350 } 351 352 /* If the vector configuration A does not have full rank then 353 change coordinates. */ 354 Triangulation2::makeConfigurationFullDimensional(A); 355 356 Triangulation2 t(A); 357 t.triangulate(); 358 359 if(searchOption.getValue()) 360 { 361 PolyhedralFan f=automatic(t,t.totalVolume()); 362 } 363 else if(hirschOption.getValue()) 364 { 365 PolyhedralFan f=automaticHirsch(t); 366 } 367 else 368 { 369 #if !NOOUTPUT 370 SymmetricTargetFanBuilder target(n,s); 371 #else 372 SymmetricTargetNull target(n,s); 373 #endif 374 375 if(!optionRestrictingFan.getValue()) 376 { 377 SecondaryFanTraverser traverser(t); 378 SymmetricTargetCounterInterrupted target2(target,optionInterrupt.getValue()); 379 symmetricTraverse(traverser,target2,&s); 380 } 381 else 382 { 383 PolyhedralFan f1=PolyhedralFan::readFan(optionRestrictingFan.getValue(),true,0,0,/*optionSymmetry.getValue()?&s:0*/0,false/*true*/); 384 385 for(PolyhedralFan::coneIterator i=f1.conesBegin();i!=f1.conesEnd();i++) 386 { 387 static int a; 388 log2 cerr<<"Processing Cone "<<a++<<" which has dimension "<<i->dimension()<<endl; 389 SecondaryFanTraverser traverser(triangulationWithFullDimensionalIntersection(t,*i),*i); 390 SymmetricTargetCounterInterrupted target2(target,optionInterrupt.getValue()); 391 symmetricTraverse(traverser,target2,&s); 392 } 393 } 394 395 #if !NOOUTPUT 396 target.getFanRef().printWithIndices(&pout, 397 (symmetryOption.getValue()?FPF_group|FPF_conesCompressed:0)| 398 (optionIgnoreCones.getValue()?0:FPF_conesExpanded)| 399 FPF_maximalCones|FPF_cones, 400 &s); 401 #endif 402 /* target.getFanRef().printWithIndices(&p, 403 FPF_default| 404 (symmetryOption.getValue()?FPF_group|FPF_conesCompressed:0), 405 &s); 406 */ 407 /* PolyhedralFan f=enumerate(t).negated();//Changing sign 408 409 f.printWithIndices(&p, 410 FPF_default| 411 (symmetryOption.getValue()?FPF_group|FPF_conesCompressed:0), 412 &s); 413 */ 414 } 415 return 0; 416 } 417 }; 418 419 static SecondaryFanApplication theApplication; 420 421