source: git/gfanlib/gfanlib_traversal.cpp @ 47a774

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