1 | #include <stddef.h> |
---|
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) |
---|
16 | { |
---|
17 | list<ZVector> ret; |
---|
18 | for(int i=0;i<m.getHeight();i++) |
---|
19 | ret.push_back(m[i]); |
---|
20 | return ret; |
---|
21 | } |
---|
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): |
---|
34 | coneCollection(n) |
---|
35 | { |
---|
36 | } |
---|
37 | |
---|
38 | |
---|
39 | bool FanBuilder::process(FanTraverser &traverser) |
---|
40 | { |
---|
41 | ZCone cone2=traverser.refToPolyhedralCone(); |
---|
42 | cone2.canonicalize(); |
---|
43 | coneCollection.insert(cone2); |
---|
44 | return true; |
---|
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 | { |
---|
131 | list<ZVector> ret; |
---|
132 | list<ZVector> normalsRet; |
---|
133 | set<ZVector> representatives; |
---|
134 | list<ZVector>::const_iterator I; |
---|
135 | if(normals)I=normals->begin(); |
---|
136 | for(list<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++) |
---|
137 | { |
---|
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++; |
---|
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 | { |
---|
155 | cerr << i->first.first << i->first.second; |
---|
156 | cerr << endl; |
---|
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 | { |
---|
185 | list<ZVector> ridges; |
---|
186 | list<ZVector> ridgesRayUniqueVector;//stores the ray of the link that we came from |
---|
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. |
---|
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; |
---|
271 | |
---|
272 | ZMatrix linealitySpaceGenerators=traverser.refToPolyhedralCone().generatorsOfLinealitySpace(); |
---|
273 | |
---|
274 | int d=traverser.refToPolyhedralCone().dimension(); |
---|
275 | |
---|
276 | Boundary boundary(*symmetryGroup); |
---|
277 | list<pathStepFacet> facetStack; |
---|
278 | list<pathStepRidge> ridgeStack; |
---|
279 | |
---|
280 | int numberOfCompletedFacets=0; |
---|
281 | int numberOfCompletedRidges=0; |
---|
282 | int stackSize=0; |
---|
283 | |
---|
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; |
---|
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 | |
---|
306 | top.parentRidge=facetStack.front().ridges.front(); |
---|
307 | // AsciiPrinter(Stderr)<<top.parentRay<<"--------------------------------++++\n"; |
---|
308 | list<ZVector> rays; |
---|
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 | |
---|
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); |
---|
356 | target.process(traverser); |
---|
357 | |
---|
358 | // IntegerVectorList inequalities=traverser.refToPolyhedralCone().getHalfSpaces(); |
---|
359 | ZMatrix equations=traverser.refToPolyhedralCone().getEquations(); |
---|
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())) |
---|
406 | { |
---|
407 | facetStack.front().ridges.pop_front(); |
---|
408 | } |
---|
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"); |
---|
434 | } |
---|