[5443c1] | 1 | #include <stddef.h> |
---|
[5cea7a] | 2 | #include "gfanlib_traversal.h" |
---|
| 3 | #include "gfanlib_symmetry.h" |
---|
| 4 | |
---|
| 5 | #include <map> |
---|
| 6 | #include <algorithm> |
---|
| 7 | #include <iostream> |
---|
| 8 | |
---|
| 9 | //#include "log.h" |
---|
| 10 | |
---|
| 11 | using namespace std; |
---|
| 12 | using namespace gfan; |
---|
| 13 | |
---|
| 14 | |
---|
| 15 | static list<ZVector> rowsToList(ZMatrix const &m) |
---|
[15813d] | 16 | { |
---|
| 17 | list<ZVector> ret; |
---|
| 18 | for(int i=0;i<m.getHeight();i++) |
---|
| 19 | ret.push_back(m[i]); |
---|
| 20 | return ret; |
---|
| 21 | } |
---|
[5cea7a] | 22 | |
---|
| 23 | bool FanTraverser::hasNoState()const |
---|
| 24 | { |
---|
| 25 | return false; |
---|
| 26 | } |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | /** |
---|
| 30 | * FanBuilder |
---|
| 31 | */ |
---|
| 32 | |
---|
| 33 | FanBuilder::FanBuilder(int n, SymmetryGroup const &sym): |
---|
[15813d] | 34 | coneCollection(n) |
---|
[5cea7a] | 35 | { |
---|
| 36 | } |
---|
| 37 | |
---|
| 38 | |
---|
| 39 | bool FanBuilder::process(FanTraverser &traverser) |
---|
| 40 | { |
---|
[15813d] | 41 | ZCone cone2=traverser.refToPolyhedralCone(); |
---|
| 42 | cone2.canonicalize(); |
---|
| 43 | coneCollection.insert(cone2); |
---|
| 44 | return true; |
---|
[5cea7a] | 45 | } |
---|
| 46 | |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | |
---|
| 50 | /** |
---|
| 51 | * Classes |
---|
| 52 | */ |
---|
| 53 | |
---|
| 54 | /** The hypergraph of ridges and facets can be considered as a usual |
---|
| 55 | bipartite graph where the right nodes are the ridges and the left |
---|
| 56 | nodes are the facets. We wish to make a traversal of this |
---|
| 57 | bipartite graph keeping track of the boundary edges of the |
---|
| 58 | traversed set. The ConeOrbit object represents the orbit of a |
---|
| 59 | ridge. The edges of the ridge are listed but only those which |
---|
| 60 | belong to the boundary of the set of ridges seen so far. When a |
---|
| 61 | ridge is discovered the ConeOrbit object will be created with all |
---|
| 62 | its edges present (except the one it was reached by). As progress |
---|
| 63 | in the computation is made these edges will be deleted. |
---|
| 64 | */ |
---|
| 65 | |
---|
| 66 | |
---|
| 67 | class Boundary |
---|
| 68 | { |
---|
| 69 | typedef pair<ZVector,ZVector> EFirst; |
---|
| 70 | struct ESecond{// L2 maybe zero, in that case i1==i2 |
---|
| 71 | list<ZVector>* L1; |
---|
| 72 | list<ZVector>::iterator i1; |
---|
| 73 | list<ZVector> *L2; |
---|
| 74 | list<ZVector>::iterator i2; |
---|
| 75 | ESecond():L1(0),L2(0){}; |
---|
| 76 | |
---|
| 77 | ESecond(list<ZVector>* L1_,list<ZVector>::iterator i1_,list<ZVector>* L2_,list<ZVector>::iterator i2_): |
---|
| 78 | L1(L1_), |
---|
| 79 | i1(i1_), |
---|
| 80 | L2(L2_), |
---|
| 81 | i2(i2_) |
---|
| 82 | { |
---|
| 83 | } |
---|
| 84 | }; |
---|
| 85 | SymmetryGroup const &sym; |
---|
| 86 | map<EFirst,ESecond > theSet; |
---|
| 87 | int theSetSize; |
---|
| 88 | public: |
---|
| 89 | Boundary(SymmetryGroup const &sym_): |
---|
| 90 | sym(sym_), |
---|
| 91 | theSetSize(0) |
---|
| 92 | { |
---|
| 93 | } |
---|
| 94 | int size()const |
---|
| 95 | { |
---|
| 96 | return theSetSize; |
---|
| 97 | } |
---|
| 98 | pair<ZVector,ZVector> normalForm(ZVector const &ridge, ZVector const &ray)const |
---|
| 99 | { |
---|
| 100 | pair<ZVector,ZVector> ret; |
---|
| 101 | Permutation perm(ridge.size()); |
---|
| 102 | ret.first=sym.orbitRepresentative(ridge,&perm); |
---|
| 103 | ret.second=sym.orbitRepresentativeFixing(perm.apply(ray),ret.first); |
---|
| 104 | return ret; |
---|
| 105 | } |
---|
| 106 | bool containsFlip(ZVector const &ridge, ZVector const &ray, list<ZVector> *storedInList_, list<ZVector>::iterator listIterator_, list<ZVector> *storedInList2_, list<ZVector>::iterator listIterator2_) |
---|
| 107 | { |
---|
| 108 | assert(ridge.size()==ray.size()); |
---|
| 109 | EFirst p=normalForm(ridge,ray); |
---|
| 110 | if(theSet.count(p)==1) |
---|
| 111 | { |
---|
| 112 | theSet[p].L1->erase(theSet[p].i1); |
---|
| 113 | if(theSet[p].L2)theSet[p].L2->erase(theSet[p].i2); |
---|
| 114 | theSet.erase(p); |
---|
| 115 | theSetSize--; |
---|
| 116 | return true; |
---|
| 117 | } |
---|
| 118 | theSet[p]=ESecond(storedInList_,listIterator_,storedInList2_,listIterator2_); |
---|
| 119 | theSetSize++; |
---|
| 120 | return false; |
---|
| 121 | } |
---|
| 122 | /** |
---|
| 123 | * This routine remove rays from rays, such that only one ridge-ray pair is left for each orbit. |
---|
| 124 | * The routine allows an additional list of vectors with the same number of elements as rays to be passed. |
---|
| 125 | * The routine will remove those vectors from this set which correspond to rays removed from rays. |
---|
| 126 | * |
---|
| 127 | * To do this it must know the symmetry group. |
---|
| 128 | */ |
---|
| 129 | void removeDuplicates(ZVector const &ridge, list<ZVector> &rays, list<ZVector> *normals=0)const |
---|
| 130 | { |
---|
[15813d] | 131 | list<ZVector> ret; |
---|
| 132 | list<ZVector> normalsRet; |
---|
| 133 | set<ZVector> representatives; |
---|
| 134 | list<ZVector>::const_iterator I; |
---|
[5cea7a] | 135 | if(normals)I=normals->begin(); |
---|
| 136 | for(list<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++) |
---|
| 137 | { |
---|
[15813d] | 138 | ZVector rep=sym.orbitRepresentativeFixing(*i,ridge); |
---|
| 139 | if(representatives.count(rep)==0) |
---|
| 140 | { |
---|
| 141 | representatives.insert(rep); |
---|
| 142 | ret.push_back(*i); |
---|
| 143 | if(normals)normalsRet.push_back(*I); |
---|
| 144 | } |
---|
| 145 | if(normals)I++; |
---|
[5cea7a] | 146 | } |
---|
| 147 | rays=ret; |
---|
| 148 | if(normals)*normals=normalsRet; |
---|
| 149 | } |
---|
| 150 | void print()const |
---|
| 151 | { |
---|
| 152 | cerr<< "Boundary" <<endl; |
---|
| 153 | for(map<EFirst,ESecond>::const_iterator i=theSet.begin();i!=theSet.end();i++) |
---|
| 154 | { |
---|
[15813d] | 155 | cerr << i->first.first << i->first.second; |
---|
| 156 | cerr << endl; |
---|
[5cea7a] | 157 | } |
---|
| 158 | cerr<<endl<<endl; |
---|
| 159 | } |
---|
| 160 | }; |
---|
| 161 | |
---|
| 162 | /** |
---|
| 163 | Rewrite these comments. |
---|
| 164 | |
---|
| 165 | During traversal the path from the current facet to the starting |
---|
| 166 | facet is stored on a stack. The elements of the stack are objects |
---|
| 167 | of the class pathStep. The top element on the stack tells through |
---|
| 168 | which ridge the current facet was reached. This is done by the |
---|
| 169 | value parent ridge which is the unique ray of the ridge. In order |
---|
| 170 | not to recompute the ridge the path facet contains rays of the link |
---|
| 171 | of the ridge represented by their unique vector. - or rather only |
---|
| 172 | the rays that are yet to be processed are stored in ridgeRays. In |
---|
| 173 | order to trace the path back the unique point of the ray from which |
---|
| 174 | the ridge was reached is stored in parentRay. |
---|
| 175 | */ |
---|
| 176 | struct pathStepRidge |
---|
| 177 | { |
---|
| 178 | ZVector parentRidge; |
---|
| 179 | list<ZVector> rays; |
---|
| 180 | ZVector parentRay; |
---|
| 181 | }; |
---|
| 182 | |
---|
| 183 | struct pathStepFacet |
---|
| 184 | { |
---|
[15813d] | 185 | list<ZVector> ridges; |
---|
| 186 | list<ZVector> ridgesRayUniqueVector;//stores the ray of the link that we came from |
---|
[5cea7a] | 187 | }; |
---|
| 188 | |
---|
| 189 | /** |
---|
| 190 | We need to simulate two mutually recursive functions. An actual |
---|
| 191 | implementation of these two functions would probably not work since |
---|
| 192 | the recursion depth could easily be 10000. |
---|
| 193 | |
---|
| 194 | Here follows a sketch of the simulation |
---|
| 195 | |
---|
| 196 | lav kegle |
---|
| 197 | find ridges |
---|
| 198 | skriv ned i objekt |
---|
| 199 | put paa stakken |
---|
| 200 | |
---|
| 201 | L1: |
---|
| 202 | if ridges in top element |
---|
| 203 | compute tropical curve |
---|
| 204 | construct stak object with rays; set parrentRidge,ridgeRays |
---|
| 205 | push ridge |
---|
| 206 | else |
---|
| 207 | pop facet |
---|
| 208 | if empty break; |
---|
| 209 | |
---|
| 210 | goto L2 |
---|
| 211 | |
---|
| 212 | L2: |
---|
| 213 | if(ridgeRays not empty) |
---|
| 214 | change CONE |
---|
| 215 | <---entry point |
---|
| 216 | push facet |
---|
| 217 | else |
---|
| 218 | pop ridge |
---|
| 219 | change CONE |
---|
| 220 | |
---|
| 221 | goto L1 |
---|
| 222 | |
---|
| 223 | |
---|
| 224 | The strategy for marking is as follows Before a vertex is pushed the |
---|
| 225 | edges that needs to be taken are written in its data. A edge is only |
---|
| 226 | written if its orbit has not been marked. Each time an edge is written |
---|
| 227 | it is simultaneously marked. |
---|
| 228 | |
---|
| 229 | */ |
---|
| 230 | |
---|
| 231 | |
---|
| 232 | |
---|
| 233 | static void printStack(list<pathStepFacet> const &facetStack, list<pathStepRidge> const &ridgeStack) |
---|
| 234 | { |
---|
| 235 | list<pathStepFacet>::const_iterator i=facetStack.begin(); |
---|
| 236 | list<pathStepRidge>::const_iterator j=ridgeStack.begin(); |
---|
| 237 | cerr<<"STACK:"<<endl; |
---|
| 238 | if(facetStack.size()>ridgeStack.size())goto entry; |
---|
| 239 | do |
---|
| 240 | { |
---|
| 241 | cerr<<"RIDGE:"<<endl; |
---|
| 242 | // cerr<<j->parentRidge<<j->rays<<j->parentRay; |
---|
| 243 | cerr<<endl; |
---|
| 244 | |
---|
| 245 | j++; |
---|
| 246 | entry: |
---|
| 247 | cerr<<"FACET:"<<endl; |
---|
| 248 | // cerr<<i->ridges; |
---|
| 249 | cerr<<endl; |
---|
| 250 | i++; |
---|
| 251 | } |
---|
| 252 | while(i!=facetStack.end()); |
---|
| 253 | |
---|
| 254 | int a; |
---|
| 255 | //cin >> a; |
---|
| 256 | } |
---|
| 257 | |
---|
| 258 | |
---|
| 259 | |
---|
| 260 | void traverse(FanTraverser &traverser, Target &target, SymmetryGroup const *symmetryGroup) |
---|
| 261 | { |
---|
| 262 | //TO DO: at the moment a conetraverser can only report that it has no state if it is traversing a complete fan. |
---|
| 263 | //This is because symmetricTraverse needs go BACK to compute the links of previously seen facets. |
---|
| 264 | //Alternative the links should be computed and stored the first time a facet is seen. |
---|
| 265 | //Or the conetraverse should be given more info about the ridge to make computations quicker. |
---|
[15813d] | 266 | int lastNumberOfEdges=0; |
---|
| 267 | float averageEdge=0; |
---|
| 268 | int n=traverser.refToPolyhedralCone().ambientDimension();//symmetryGroup->sizeOfBaseSet(); |
---|
| 269 | SymmetryGroup localSymmetryGroup(n); |
---|
| 270 | if(!symmetryGroup)symmetryGroup=&localSymmetryGroup; |
---|
[5cea7a] | 271 | |
---|
[15813d] | 272 | ZMatrix linealitySpaceGenerators=traverser.refToPolyhedralCone().generatorsOfLinealitySpace(); |
---|
[5cea7a] | 273 | |
---|
[15813d] | 274 | int d=traverser.refToPolyhedralCone().dimension(); |
---|
[5cea7a] | 275 | |
---|
[15813d] | 276 | Boundary boundary(*symmetryGroup); |
---|
| 277 | list<pathStepFacet> facetStack; |
---|
| 278 | list<pathStepRidge> ridgeStack; |
---|
[5cea7a] | 279 | |
---|
[15813d] | 280 | int numberOfCompletedFacets=0; |
---|
| 281 | int numberOfCompletedRidges=0; |
---|
| 282 | int stackSize=0; |
---|
[5cea7a] | 283 | |
---|
[15813d] | 284 | ZVector facetUniqueVector; |
---|
| 285 | goto entry; |
---|
| 286 | while(1) |
---|
| 287 | { |
---|
| 288 | L1: |
---|
| 289 | // printStack(facetStack,ridgeStack); |
---|
| 290 | //if we have more ProcessRidge calls to do |
---|
| 291 | if(!facetStack.front().ridges.empty()) |
---|
| 292 | { |
---|
| 293 | //ProcessRidge "called" |
---|
| 294 | pathStepRidge top; |
---|
[5cea7a] | 295 | |
---|
| 296 | |
---|
| 297 | if(traverser.hasNoState()) |
---|
| 298 | top.parentRay=facetStack.front().ridgesRayUniqueVector.front(); |
---|
| 299 | else |
---|
| 300 | { |
---|
| 301 | ZCone link=traverser.refToPolyhedralCone().link(facetStack.front().ridges.front()); |
---|
| 302 | link.canonicalize(); |
---|
| 303 | top.parentRay=link.getUniquePoint(); |
---|
| 304 | } |
---|
| 305 | |
---|
[15813d] | 306 | top.parentRidge=facetStack.front().ridges.front(); |
---|
| 307 | // AsciiPrinter(Stderr)<<top.parentRay<<"--------------------------------++++\n"; |
---|
| 308 | list<ZVector> rays; |
---|
[5cea7a] | 309 | if(traverser.hasNoState()) |
---|
| 310 | { |
---|
| 311 | rays.push_back(top.parentRay); |
---|
| 312 | rays.push_back(-top.parentRay); |
---|
| 313 | } |
---|
| 314 | else |
---|
| 315 | rays=traverser.link(facetStack.front().ridges.front()); |
---|
| 316 | |
---|
[15813d] | 317 | assert(!rays.empty()); |
---|
| 318 | boundary.removeDuplicates(top.parentRidge,rays); |
---|
| 319 | ridgeStack.push_front(top);stackSize++; |
---|
| 320 | ZVector ridgeRidgeRidge=facetStack.front().ridges.front(); |
---|
| 321 | for(list<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++) |
---|
| 322 | { |
---|
| 323 | ridgeStack.front().rays.push_front(*i); |
---|
| 324 | if(boundary.containsFlip(ridgeRidgeRidge,*i,&ridgeStack.front().rays,ridgeStack.front().rays.begin(),0,ridgeStack.front().rays.begin())) |
---|
| 325 | ridgeStack.front().rays.pop_front(); |
---|
| 326 | } |
---|
| 327 | // "state saved" ready to do calls to ProcessFacet |
---|
| 328 | numberOfCompletedRidges++; |
---|
| 329 | } |
---|
| 330 | else |
---|
| 331 | { |
---|
| 332 | // No more calls to do - we now return from ProcessFacet |
---|
| 333 | //THIS IS THE PLACE TO CHANGE CONE BACK |
---|
| 334 | facetStack.pop_front();stackSize--; |
---|
| 335 | if(facetStack.empty())break; |
---|
| 336 | // log1 cerr<<"BACK"<<endl; |
---|
| 337 | if(!traverser.hasNoState())traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().parentRay); |
---|
| 338 | } |
---|
| 339 | L2: |
---|
| 340 | // printStack(facetStack,ridgeStack); |
---|
| 341 | //check if ProcessRidge needs to make more ProcessFacet calls |
---|
| 342 | if(!ridgeStack.front().rays.empty()) |
---|
| 343 | { |
---|
| 344 | // log1 cerr<<"FORWARD"<<endl; |
---|
| 345 | traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().rays.front()); |
---|
| 346 | entry: |
---|
| 347 | //ProcessFacet() |
---|
| 348 | averageEdge=0.99*averageEdge+0.01*(boundary.size()-lastNumberOfEdges); |
---|
| 349 | // log1 fprintf(Stderr,"\n-------------------------------------\n"); |
---|
| 350 | // log1 fprintf(Stderr,"Boundary edges in bipartite graph: %i, Completed ridges: %i, Completed facets: %i, Recursion depth:%i Average new edge/facet:%0.2f\n",boundary.size(),numberOfCompletedRidges,numberOfCompletedFacets,stackSize,averageEdge); |
---|
| 351 | // log1 fprintf(Stderr,"-------------------------------------\n"); |
---|
| 352 | lastNumberOfEdges=boundary.size(); |
---|
| 353 | |
---|
| 354 | // target.process(traverser);//Postponed until extrem rays have been computed |
---|
| 355 | ZMatrix extremeRays=traverser.refToPolyhedralCone().extremeRays(&linealitySpaceGenerators); |
---|
[5cea7a] | 356 | target.process(traverser); |
---|
| 357 | |
---|
[15813d] | 358 | // IntegerVectorList inequalities=traverser.refToPolyhedralCone().getHalfSpaces(); |
---|
[5cea7a] | 359 | ZMatrix equations=traverser.refToPolyhedralCone().getEquations(); |
---|
[15813d] | 360 | // facetUniqueVector=traverser.refToPolyhedralCone().getUniquePoint(); |
---|
| 361 | facetUniqueVector=traverser.refToPolyhedralCone().getUniquePointFromExtremeRays(extremeRays); |
---|
| 362 | list<ZVector> facetNormals=rowsToList(traverser.refToPolyhedralCone().getFacets()); |
---|
| 363 | |
---|
| 364 | pathStepFacet stepFacet; |
---|
| 365 | list<ZVector> ridges; |
---|
| 366 | |
---|
| 367 | for(list<ZVector>::iterator i=facetNormals.begin();i!=facetNormals.end();i++) |
---|
| 368 | { |
---|
| 369 | ZVector v(n); |
---|
| 370 | // for(IntegerVectorList::const_iterator j=extremeRays.begin();j!=extremeRays.end();j++) |
---|
| 371 | for(int j=0;j<extremeRays.getHeight();j++) |
---|
| 372 | if(dot(*i,extremeRays[j]).isZero())v+=extremeRays[j]; |
---|
| 373 | ridges.push_back(v); |
---|
| 374 | } |
---|
| 375 | |
---|
| 376 | ZVector temp(n); |
---|
| 377 | // boundary.removeDuplicates(temp,ridges);//use facetUniqueVector instead |
---|
| 378 | boundary.removeDuplicates(facetUniqueVector,ridges,&facetNormals);//use facetUniqueVector instead |
---|
| 379 | |
---|
| 380 | facetStack.push_front(stepFacet);stackSize++; |
---|
| 381 | list<ZVector>::const_iterator I=facetNormals.begin(); |
---|
| 382 | for(list<ZVector>::const_iterator i=ridges.begin();i!=ridges.end();i++,I++) |
---|
| 383 | { |
---|
| 384 | ZVector rayUniqueVector; |
---|
| 385 | |
---|
| 386 | if(d==n) |
---|
| 387 | { |
---|
| 388 | rayUniqueVector =I->normalized(); |
---|
| 389 | // if(dotLong(rayUniqueVector,*I) |
---|
| 390 | } |
---|
| 391 | else |
---|
| 392 | { |
---|
| 393 | ZCone rayCone=traverser.refToPolyhedralCone().link(*i); |
---|
| 394 | rayCone.canonicalize(); |
---|
| 395 | rayUniqueVector=rayCone.getUniquePoint(); |
---|
| 396 | // debug<<traverser.refToPolyhedralCone(); |
---|
| 397 | // debug<<rayCone; |
---|
| 398 | } |
---|
| 399 | facetStack.front().ridges.push_front(*i); |
---|
| 400 | if(traverser.hasNoState())facetStack.front().ridgesRayUniqueVector.push_front(rayUniqueVector); |
---|
| 401 | |
---|
| 402 | |
---|
| 403 | if(!traverser.hasNoState()) |
---|
| 404 | { |
---|
| 405 | if(boundary.containsFlip(*i,rayUniqueVector,&facetStack.front().ridges,facetStack.front().ridges.begin(),0,facetStack.front().ridges.begin())) |
---|
[5cea7a] | 406 | { |
---|
| 407 | facetStack.front().ridges.pop_front(); |
---|
| 408 | } |
---|
[15813d] | 409 | } |
---|
| 410 | else |
---|
| 411 | { |
---|
| 412 | if(boundary.containsFlip(*i,rayUniqueVector,&facetStack.front().ridges,facetStack.front().ridges.begin(),&facetStack.front().ridgesRayUniqueVector,facetStack.front().ridgesRayUniqueVector.begin())) |
---|
| 413 | { |
---|
| 414 | facetStack.front().ridges.pop_front(); |
---|
| 415 | facetStack.front().ridgesRayUniqueVector.pop_front(); |
---|
| 416 | } |
---|
| 417 | } |
---|
| 418 | } |
---|
| 419 | //"State pushed" ready to call ProcessRidge |
---|
| 420 | |
---|
| 421 | numberOfCompletedFacets++; |
---|
| 422 | } |
---|
| 423 | else |
---|
| 424 | { |
---|
| 425 | //ProcessRidge is done making its calls to ProcessFacet so we can return from ProcessRidge |
---|
| 426 | // cerr<<"BACK"<<endl; |
---|
| 427 | // traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().parentRay); |
---|
| 428 | ridgeStack.pop_front();stackSize--; |
---|
| 429 | } |
---|
| 430 | }//goto L1 |
---|
| 431 | // log1 fprintf(Stderr,"\n-------------------------------------\n"); |
---|
| 432 | // log1 fprintf(Stderr,"Boundary edges in bipartite graph: %i, Completed ridges: %i, Completed facets: %i, Recursion depth:%i\n",boundary.size(),numberOfCompletedRidges,numberOfCompletedFacets,stackSize); |
---|
| 433 | // log1 fprintf(Stderr,"-------------------------------------\n"); |
---|
[5cea7a] | 434 | } |
---|