# source:git/gfanlib/gfanlib_traversal.cpp@5443c1

jengelh-datetimespielwiese
Last change on this file since 5443c1 was 5443c1, checked in by Hans Schoenemann <hannes@…>, 4 years ago
porting: gcc-5 (Ubuntu 16.04.5 LTS)
• Property mode set to `100644`
File size: 15.6 KB
Line
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
11using namespace std;
12using namespace gfan;
13
14
15static 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
23bool FanTraverser::hasNoState()const
24{
25  return false;
26}
27
28
29/**
30 * FanBuilder
31 */
32
33FanBuilder::FanBuilder(int n, SymmetryGroup const &sym):
34        coneCollection(n)
35{
36}
37
38
39bool 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
67class 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;
88public:
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/**
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 */
176struct pathStepRidge
177{
178  ZVector parentRidge;
179  list<ZVector> rays;
180  ZVector parentRay;
181};
182
183struct 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
196lav kegle
197find ridges
198skriv ned i objekt
199put paa stakken
200
201L1:
202if ridges in top element
203  compute tropical curve
204  construct stak object with rays; set parrentRidge,ridgeRays
205  push ridge
206else
207  pop facet
208  if empty break;
209
210goto L2
211
212L2:
213if(ridgeRays not empty)
214   change CONE
215   <---entry point
216   push facet
217else
218  pop ridge
219  change CONE
220
221goto L1
222
223
224The strategy for marking is as follows Before a vertex is pushed the
225edges that needs to be taken are written in its data. A edge is only
226written if its orbit has not been marked. Each time an edge is written
227it is simultaneously marked.
228
229*/
230
231
232
233static 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
260void 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.
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                    {
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
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);
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                          {
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}
Note: See TracBrowser for help on using the repository browser.