Changeset 008e81 in git


Ignore:
Timestamp:
Jun 12, 2009, 6:32:00 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '00e2e9c41af3fde1273eb3633f4c0c7c3db2579d')
Children:
cf14b9d47b88351f19f0fd851564bb4cb6d6a586
Parents:
351764983bd71ea77010993df1b11657f12eda1d
Message:
introduced method gcone::areEqual to check == of facets. Probably still buggy.
New copy constructor
Modified ivNeg
More on noRevs, where the code for intPoint was removed
Some cleanup


git-svn-id: file:///usr/local/Singular/svn/trunk@11897 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r3517649 r008e81  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-06-09 15:23:08 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.58 2009-06-09 15:23:08 monerjan Exp $
    6 $Id: gfan.cc,v 1.58 2009-06-09 15:23:08 monerjan Exp $
     4$Date: 2009-06-12 16:32:00 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.59 2009-06-12 16:32:00 monerjan Exp $
     6$Id: gfan.cc,v 1.59 2009-06-12 16:32:00 monerjan Exp $
    77*/
    88
    … …  
    9595                facet *next;            //Pointer to next facet
    9696                facet *codim2Ptr;       //Pointer to (codim-2)-facet. Bit of recursion here ;-)
    97                 //intvec **codim2Normals =(intvec**)omAlloc0(sizeof(intvec*));  //Integer matrix containing the (codim-2)-facets
     97                int numCodim2Facets;
     98                ring flipRing;          //the ring on the other side of the facet
    9899                                               
    99100                /** The default constructor. Do I need a constructor of type facet(intvec)? */
    … …  
    106107                        this->codim2Ptr=NULL;
    107108                        this->codim=1;          //default for (codim-1)-facets
     109                        this->numCodim2Facets=0;
    108110                }
    109111               
    … …  
    119121                                this->codim=n;
    120122                        }//NOTE Handle exception here!
     123                        this->numCodim2Facets=0;
    121124                }
    122125               
    123126                /** The default destructor */
    124127                ~facet(){;}
     128               
     129                /** \brief Comparison of facets*/
     130                bool areEqual(facet &f, facet &g)
     131                {
     132                        bool res = TRUE;
     133                        intvec *fNormal = new intvec(pVariables);
     134                        intvec *gNormal = new intvec(pVariables);
     135                        fNormal = f.getFacetNormal();
     136                        gNormal = g.getFacetNormal();
     137                        if((fNormal == gNormal))//||(gcone::isParallel(fNormal,gNormal)))
     138                        {
     139                                if(f.numCodim2Facets==g.numCodim2Facets)
     140                                {
     141                                        facet *f2Act;
     142                                        facet *g2Act;
     143                                        f2Act = f.codim2Ptr;
     144                                        g2Act = g.codim2Ptr;
     145                                        intvec *f2N = new intvec(pVariables);
     146                                        intvec *g2N = new intvec(pVariables);
     147                                        while(f2Act->next!=NULL && g2Act->next!=NULL)
     148                                        {
     149                                                for(int ii=0;ii<f2N->length();ii++)
     150                                                {
     151                                                        if(f2Act->getFacetNormal() != g2Act->getFacetNormal())
     152                                                        {
     153                                                                res=FALSE;                                                             
     154                                                        }
     155                                                        if (res==FALSE)
     156                                                                break;
     157                                                }
     158                                                if(res==FALSE)
     159                                                        break;
     160                                               
     161                                                f2Act = f2Act->next;
     162                                                g2Act = g2Act->next;
     163                                        }//while
     164                                }//if           
     165                        }//if
     166                        else
     167                        {
     168                                res = FALSE;
     169                        }
     170                        return res;
     171                }
    125172               
    126173                /** Stores the facet normal \param intvec*/
    … …  
    295342                */                     
    296343                //NOTE Prolly need to specify the facet to flip over
    297                 gcone(const gcone& gc, const facet &f)         
     344//              gcone(const gcone& gc, const facet &f)         
     345//              {
     346//                      this->next=NULL;
     347//                      this->numVars=gc.numVars;
     348//                      this->UCN=(gc.UCN)+1;   //add 1 to the UCN of previous cone. This is NOT UNIQUE!
     349//                      facet *fAct= new facet();                       
     350//                      this->facetPtr=fAct;
     351//                     
     352//                      intvec *ivtmp = new intvec(this->numVars);
     353//                      //NOTE ivtmp = f->getFacetNormal();                                             
     354//                      ivtmp = gc.facetPtr->getFacetNormal();                 
     355//                     
     356//                      ideal gb;
     357//                      //NOTE gb=f->getFlipGB();
     358//                      gb=gc.facetPtr->getFlipGB();                   
     359//                      this->gcBasis=gb;       //this cone's GB is the flipped GB                     
     360//                     
     361//                      /*Reverse direction of the facet normal to make it an inner normal*/                   
     362//                      for (int ii=0; ii<this->numVars;ii++)
     363//                      {
     364//                              (*ivtmp)[ii]=-(*ivtmp)[ii];                             
     365//                      }
     366//                     
     367//                      fAct->setFacetNormal(ivtmp);
     368//                      delete ivtmp;
     369//                      delete fAct;                                           
     370//              }
     371               
     372                gcone(const gcone& gc, const facet &f)
    298373                {
    299374                        this->next=NULL;
    300375                        this->numVars=gc.numVars;
    301                         this->UCN=(gc.UCN)+1;   //add 1 to the UCN of previous cone. This is NOT UNIQUE!
    302                         facet *fAct= new facet();                       
    303                         this->facetPtr=fAct;
    304                        
    305                         intvec *ivtmp = new intvec(this->numVars);
    306                         //NOTE ivtmp = f->getFacetNormal();                                             
    307                         ivtmp = gc.facetPtr->getFacetNormal();                 
    308                        
    309                         ideal gb;
    310                         //NOTE gb=f->getFlipGB();
    311                         gb=gc.facetPtr->getFlipGB();                   
    312                         this->gcBasis=gb;       //this cone's GB is the flipped GB                     
    313                        
    314                         /*Reverse direction of the facet normal to make it an inner normal*/                   
    315                         for (int ii=0; ii<this->numVars;ii++)
    316                         {
    317                                 (*ivtmp)[ii]=-(*ivtmp)[ii];                             
    318                         }
    319                        
    320                         fAct->setFacetNormal(ivtmp);
    321                         delete ivtmp;
    322                         delete fAct;                                           
     376                        this->UCN=(gc.UCN)+1;
     377                        this->facetPtr=NULL;
     378                        this->gcBasis=gc.gcBasis;
     379                        this->baseRing=f.flipRing;
     380                        rChangeCurrRing(this->baseRing);
    323381                }
    324382               
    … …  
    328386                {
    329387                        //NOTE SAVE THE FACET STRUCTURE!!!
    330                 }
    331                        
     388                }                       
    332389               
    333390                /** \brief Set the interior point of a cone */
    … …  
    602659                void getCodim2Normals(gcone const &gc)
    603660                {
    604                         this->facetPtr->codim2Ptr = new facet(2);       //instantiate a (codim-2)-facet
     661                        this->facetPtr->codim2Ptr = new facet(2);       //instantiate a (codim-2)-facet                 
    605662                        facet *codim2Act;
    606663                        codim2Act = this->facetPtr->codim2Ptr;
    … …  
    638695                                        makeInt(P,jj,*n);                                               
    639696                                        codim2Act->setFacetNormal(n);
     697                                        this->facetPtr->numCodim2Facets++;      //Honor the creation of a codim-2-facet
    640698                                        codim2Act->next = new facet(2);
    641699                                        codim2Act = codim2Act->next;                                           
    … …  
    937995                       
    938996                        f->setFlipGB(dstRing_I);//store the flipped GB
     997                        f->flipRing=dstRing;    //store the ring on the other side
    939998#ifdef gfan_DEBUG
    940999                        cout << "Flipped GB is: " << endl;
    9411000                        f->printFlipGB();
    9421001#endif                 
     1002                        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
    9431003                }//void flip(ideal gb, facet *f)
    9441004                               
    … …  
    10581118                /** \brief Compute the negative of a given intvec
    10591119                */             
    1060                 intvec *ivNeg(intvec const &iv)
    1061                 {
    1062                         intvec *res = new intvec(this->numVars);
    1063                         for(int ii=0;ii<this->numVars;ii++)
    1064                         {
    1065                                 res[ii] = -(int)iv[ii];
    1066                         }
     1120                intvec *ivNeg(const intvec *iv)
     1121                {
     1122                        intvec *res = new intvec(iv->length());
     1123                        res=ivCopy(iv);
     1124                        *res *= (int)-1;
     1125                        //for(int ii=0;ii<this->numVars;ii++)
     1126                        //{                     
     1127                                //(res)[ii] = (*res[ii])*(int)(-1);
     1128                        //}
     1129                        res->show();                   
    10671130                        return res;
    10681131                }
    … …  
    14341497                */
    14351498                void noRevS(gcone &gc, bool usingIntPoint=FALSE)
    1436                 {
    1437                         facet *fListPtr = new facet();  //The list containing ALL facets we come across                 
    1438                         facet *fAct;// = new facet();
    1439                         fAct = gc.facetPtr;
    1440                         fListPtr = gc.facetPtr;
     1499                {       
     1500                        facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
     1501                        facet *SearchListAct;
     1502                        SearchListAct = SearchListRoot;
     1503                       
     1504                        facet *fAct;
     1505                        fAct = gc.facetPtr;                     
     1506                       
    14411507#ifdef gfan_DEBUG
    14421508                        cout << "NoRevs" << endl;
    … …  
    14441510                        gc.showFacets();
    14451511#endif                 
    1446                         //dd_set_global_constants();
    1447                         dd_rowrange ddrows;
    1448                         dd_colrange ddcols;
    1449                         ddrows=2;       //Each facet is described by its normal
    1450                         ddcols=gc.numVars+1;    // plus one for the standard simplex
    1451                         if(usingIntPoint){
    1452                                 while(fAct->next!=NULL)
    1453                                 {
    1454                                         /*Compute an interior point for each facet*/                           
    1455                                         dd_MatrixPtr ddineq;   
    1456                                         ddineq=dd_CreateMatrix(ddrows,ddcols);
    1457                                         intvec *comp = new intvec(this->numVars);
    1458                                         comp=fAct->getFacetNormal();                           
    1459                                         for(int ii=0; ii<this->numVars; ii++)                                   
    1460                                         {                                                                               
    1461                                                 dd_set_si(ddineq->matrix[0][ii+1],(*comp)[ii]);
    1462                                                 dd_set_si(ddineq->matrix[1][ii+1],1);   //Assure we search in the pos. orthant         
    1463                                         }
    1464                                         set_addelem(ddineq->linset,1);  //We want equality in the first row
    1465                                         //dd_WriteMatrix(stdout,ddineq);
    1466                                         interiorPoint(ddineq,*comp);                           
    1467                                         /**/
    1468 #ifdef gfan_DEBUG
    1469                                         cout << "IP is";
    1470                                         comp->show(); cout << endl;
    1471 #endif
    1472                                         //Store the interior point and the UCN
    1473                                         fListPtr->setInteriorPoint( comp );                             
    1474                                         fListPtr->setUCN( gc.getUCN() );       
    1475                                                                
    1476                                         dd_FreeMatrix(ddineq);
    1477                                         fAct=fAct->next;        //iterate
    1478                                 }       
    1479                         }
    1480                         else/*Facets of facets*/
    1481                         {
     1512                       
     1513                        this->getCodim2Normals(gc);
    14821514                               
    1483                                 this->getCodim2Normals(gc);
    1484                                
    1485                                 //Compute unique representation of codim-2-facets
    1486                                 this->normalize();
    1487                                
    1488                                 gc.writeConeToFile(gc);
    1489                                
    1490                                 /*2nd step
    1491                                 Choose a facet from fListPtr, flip it and forget the previous cone
    1492                                 We always choose the first facet from fListPtr as facet to be flipped
    1493                                 */
    1494                                
    1495                                
    1496                                
    1497                         }//else usingIntPoint                   
    1498                        
    1499                        
    1500        
    1501                        
     1515                        //Compute unique representation of codim-2-facets
     1516                        this->normalize();
     1517                       
     1518                        /*Make a copy of the facet list of first cone
     1519                        Since the operations getCodim2Normals and normalize affect the facets
     1520                        we must not memcpy them before these ops!
     1521                        */
     1522                        while(fAct->next!=NULL)
     1523                        {                               
     1524                                memcpy(SearchListAct,fAct,sizeof(facet));
     1525                                SearchListAct->next = new facet();
     1526                                SearchListAct = SearchListAct->next;
     1527                                fAct = fAct->next;
     1528                        }
     1529                        SearchListAct = SearchListRoot; //Set to beginning of list
     1530                       
     1531                        fAct = gc.facetPtr;
     1532                        if(areEqual(fAct->getFacetNormal(),fAct->next->getFacetNormal()))
     1533                        {
     1534                                cout <<"HI" << endl;
     1535                        }
     1536                       
     1537                        gc.writeConeToFile(gc);
     1538                       
     1539                        /*2nd step
     1540                        Choose a facet from fListPtr, flip it and forget the previous cone
     1541                        We always choose the first facet from fListPtr as facet to be flipped
     1542                        */
     1543                        while(SearchListAct->next!=NULL)
     1544                        {
     1545                                gc.flip(gc.gcBasis,SearchListAct);
     1546                                gcone *gcTmp = new gcone::gcone(gc,*SearchListAct);
     1547                       
     1548                        //gcTmp->getConeNormals();
     1549                        //add new facets
     1550                        //fListPtr = fListPtr->next;
     1551                                SearchListAct = SearchListAct->next;
     1552                        }
     1553               
    15021554                        //NOTE Hm, comment in and get a crash for free...
    15031555                        //dd_free_global_constants();                           
    … …  
    16421694                        M=dd_CreateMatrix(ddrows,ddcols);                       
    16431695                       
    1644                         //dd_set_global_constants();                                           
    1645                                                
    1646                         //intvec *comp = new intvec(this->numVars);
    1647                        
    1648                         /*for (int jj=0; jj<M->rowsize; jj++)
    1649                         {
    1650                                 intvec *comp = new intvec(this->numVars);
    1651                                 comp = fAct->getFacetNormal();
    1652                                 for (int ii=0; ii<this->numVars; ii++)
    1653                                 {
    1654                                         dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
    1655                                 }
    1656                                 if(fAct->next!=NULL)
    1657                                 {
    1658                                         fAct = fAct->next;
    1659                                 }
    1660                                 delete comp;
    1661                         }//jj
    1662                         */
    16631696                        int jj=0;
    16641697                        while(fAct->next!=NULL)
    … …  
    16931726                        fAct = gc.facetPtr;
    16941727                       
     1728                        char *ringString = rString(currRing);
     1729                       
    16951730                        if (!gcOutputFile)
    16961731                        {
    … …  
    16981733                        }
    16991734                        else
    1700                         {       /*gcOutputFile << "UCN" << endl;
    1701                                 gcOutputFile << this->UCN << endl;*/
    1702                                 gcOutputFile << "RING" << endl;                         
     1735                        {       gcOutputFile << "UCN" << endl;
     1736                                gcOutputFile << this->UCN << endl;
     1737                                gcOutputFile << "RING" << endl;
     1738                                gcOutputFile << ringString << endl;                     
    17031739                                //Write this->gcBasis into file
    17041740                                gcOutputFile << "GCBASIS" << endl;                             
Note: See TracChangeset for help on using the changeset viewer.