Changeset 15813d in git for gfanlib/gfanlib_traversal.cpp


Ignore:
Timestamp:
Aug 8, 2016, 2:16:11 PM (7 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
Children:
7cc3fdae6091a6fab65626857c1ed2c59f027dcd
Parents:
47a774fa4174b578843d43227fcc351165473747
Message:
format
File:
1 edited

Legend:

Unmodified
Added
Removed

  • gfanlib/gfanlib_traversal.cpp

    r47a774 r15813d  
    1313
    1414static 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                 }
     15                {
     16                        list<ZVector> ret;
     17                        for(int i=0;i<m.getHeight();i++)
     18                                ret.push_back(m[i]);
     19                        return ret;
     20                }
    2121
    2222bool FanTraverser::hasNoState()const
    … …  
    3131
    3232FanBuilder::FanBuilder(int n, SymmetryGroup const &sym):
    33         coneCollection(n)
     33        coneCollection(n)
    3434{
    3535}
    … …  
    3838bool FanBuilder::process(FanTraverser &traverser)
    3939{
    40         ZCone cone2=traverser.refToPolyhedralCone();
    41         cone2.canonicalize();
    42         coneCollection.insert(cone2);
    43         return true;
     40        ZCone cone2=traverser.refToPolyhedralCone();
     41        cone2.canonicalize();
     42        coneCollection.insert(cone2);
     43        return true;
    4444}
    4545
    … …  
    128128  void removeDuplicates(ZVector const &ridge, list<ZVector> &rays, list<ZVector> *normals=0)const
    129129  {
    130           list<ZVector> ret;
    131           list<ZVector> normalsRet;
    132           set<ZVector> representatives;
    133           list<ZVector>::const_iterator I;
     130          list<ZVector> ret;
     131          list<ZVector> normalsRet;
     132          set<ZVector> representatives;
     133          list<ZVector>::const_iterator I;
    134134    if(normals)I=normals->begin();
    135135    for(list<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++)
    136136      {
    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++;
     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++;
    145145      }
    146146    rays=ret;
    … …  
    152152    for(map<EFirst,ESecond>::const_iterator i=theSet.begin();i!=theSet.end();i++)
    153153      {
    154         cerr << i->first.first << i->first.second;
    155         cerr << endl;
     154            cerr << i->first.first << i->first.second;
     155            cerr << endl;
    156156      }
    157157    cerr<<endl<<endl;
    … …  
    182182struct pathStepFacet
    183183{
    184         list<ZVector> ridges;
    185         list<ZVector> ridgesRayUniqueVector;//stores the ray of the link that we came from
     184        list<ZVector> ridges;
     185        list<ZVector> ridgesRayUniqueVector;//stores the ray of the link that we came from
    186186};
    187187
    … …  
    263263  //Alternative the links should be computed and stored the first time a facet is seen.
    264264  //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;
     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;
    294294
    295295
    … …  
    303303                    }
    304304
    305                   top.parentRidge=facetStack.front().ridges.front();
    306 //                AsciiPrinter(Stderr)<<top.parentRay<<"--------------------------------++++\n";
    307                   list<ZVector> rays;
     305                  top.parentRidge=facetStack.front().ridges.front();
     306//                  AsciiPrinter(Stderr)<<top.parentRay<<"--------------------------------++++\n";
     307                  list<ZVector> rays;
    308308                  if(traverser.hasNoState())
    309309                    {
    … …  
    314314                    rays=traverser.link(facetStack.front().ridges.front());
    315315
    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);
     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);
    355355                  target.process(traverser);
    356356
    357 //                IntegerVectorList inequalities=traverser.refToPolyhedralCone().getHalfSpaces();
     357//                  IntegerVectorList inequalities=traverser.refToPolyhedralCone().getHalfSpaces();
    358358                  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()))
     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()))
    405405                        {
    406406                          facetStack.front().ridges.pop_front();
    407407                        }
    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");
     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");
    433433}
Note: See TracChangeset for help on using the changeset viewer.