Changeset f26ea0 in git


Ignore:
Timestamp:
Jun 9, 2009, 5:23:08 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')
Children:
3231f3dc7a66f525efae0e24ca084a7703a2bb7f
Parents:
44ab6bc29303a0aa0c6a2c8d1708d78ac34d6f28
Message:
Maintenance and bugfixes
method gcone::showFacets() now takes the codim as an optional parameter
New method gcone::getCodim2Normals
Made UCN global and modified constructor of gcone accordingly


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r44ab6b rf26ea0  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-05-29 07:52:12 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.57 2009-05-29 07:52:12 monerjan Exp $
    6 $Id: gfan.cc,v 1.57 2009-05-29 07:52:12 monerjan Exp $
     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 $
    77*/
    88
    … …  
    239239                /* TODO in order to save memory use pointers to rootRing and inputIdeal instead */
    240240                intvec *ivIntPt;        //an interior point of the cone
    241                 int UCN;                //unique number of the cone
     241                static int UCN;         //unique number of the cone
    242242               
    243243        public:
    … …  
    264264                {
    265265                        this->next=NULL;
    266                         this->facetPtr=NULL;    //maybe this->facetPtr = new facet();
     266                        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();                   
    267267                        this->baseRing=currRing;
    268                         this->UCN=1;
     268                        this->UCN++;
    269269                        this->numFacets=0;
    270270                }
    … …  
    280280                {
    281281                        this->next=NULL;
    282                         this->facetPtr=NULL;
     282                        this->facetPtr=NULL;                   
    283283                        this->rootRing=r;
    284284                        this->inputIdeal=I;
    285285                        this->baseRing=currRing;
    286                         this->UCN=1;
     286                        this->UCN++;
    287287                        this->numFacets=0;
    288288                }
    … …  
    290290                /** \brief Copy constructor
    291291                *
    292                 * Copies one cone, sets this->gcBasis to the flipped GB and reverses the
     292                * Copies a cone, sets this->gcBasis to the flipped GB and reverses the
    293293                * direction of the according facet normal
     294                * Call this only after a successfull call to gcone::flip which sets facet::flipGB
    294295                */                     
    295296                //NOTE Prolly need to specify the facet to flip over
    … …  
    302303                        this->facetPtr=fAct;
    303304                       
    304                         intvec *ivtmp = new intvec(this->numVars);                                             
     305                        intvec *ivtmp = new intvec(this->numVars);
     306                        //NOTE ivtmp = f->getFacetNormal();                                             
    305307                        ivtmp = gc.facetPtr->getFacetNormal();                 
    306308                       
    307309                        ideal gb;
     310                        //NOTE gb=f->getFlipGB();
    308311                        gb=gc.facetPtr->getFlipGB();                   
    309312                        this->gcBasis=gb;       //this cone's GB is the flipped GB                     
    … …  
    345348                }
    346349               
    347                 void showFacets()
    348                 {
    349                         facet *f = new facet();
    350                         f = this->facetPtr;
     350                void showFacets(short codim=1)
     351                {
     352                        facet *f;
     353                        switch(codim)
     354                        {
     355                                case 1:
     356                                        f = this->facetPtr;
     357                                        break;
     358                                case 2:
     359                                        f = this->facetPtr->codim2Ptr;
     360                                        break;
     361                        }                       
     362                       
    351363                        intvec *iv = new intvec(this->numVars);
    352364                       
    … …  
    357369                                f=f->next;
    358370                        }
    359                         //delete iv;
    360                         //delete f;
    361                 }               
     371                        //delete iv;                   
     372                }
    362373               
    363374                /** \brief Set gcone::numFacets */
    … …  
    586597                }//method getConeNormals(ideal I)       
    587598               
     599                /** \brief Compute the (codim-2)-facets of a given cone
     600                * This method is used during noRevS
     601                */
     602                void getCodim2Normals(gcone const &gc)
     603                {
     604                        this->facetPtr->codim2Ptr = new facet(2);       //instantiate a (codim-2)-facet
     605                        facet *codim2Act;
     606                        codim2Act = this->facetPtr->codim2Ptr;
     607                       
     608                        dd_set_global_constants();
     609                       
     610                        dd_MatrixPtr ddineq,P,ddakt;
     611                        dd_rowset impl_linset, redset;
     612                        dd_ErrorType err;
     613                        dd_rowindex newpos;             
     614
     615                        ddineq = facets2Matrix(gc);     //get a matrix representation of the cone
     616                               
     617                        /*Now set appropriate linearity*/
     618                        dd_PolyhedraPtr ddpolyh;
     619                        for (int ii=0; ii<this->numFacets; ii++)
     620                        {                               
     621                                ddakt = dd_CopyMatrix(ddineq);
     622                                set_addelem(ddakt->linset,ii+1);
     623                                                                               
     624                                dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);                   
     625                                       
     626                                        //dd_WriteMatrix(stdout,ddakt);
     627                                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
     628                                P=dd_CopyGenerators(ddpolyh);
     629                                dd_WriteMatrix(stdout,P);
     630                                       
     631                                /* We loop through each row of P
     632                                * normalize it by making all entries integer ones
     633                                * and add the resulting vector to the int matrix facet::codim2Facets
     634                                */
     635                                for (int jj=1;jj<=P->rowsize;jj++)
     636                                {
     637                                        intvec *n = new intvec(this->numVars);
     638                                        makeInt(P,jj,*n);                                               
     639                                        codim2Act->setFacetNormal(n);
     640                                        codim2Act->next = new facet(2);
     641                                        codim2Act = codim2Act->next;                                           
     642                                        delete n;                                                                       
     643                                }                                                                               
     644                                       
     645                                dd_FreeMatrix(ddakt);
     646                                dd_FreePolyhedra(ddpolyh);
     647                        }
     648                }
    588649               
    589650                /** \brief Compute the Groebner Basis on the other side of a shared facet
    … …  
    13751436                {
    13761437                        facet *fListPtr = new facet();  //The list containing ALL facets we come across                 
    1377                         facet *fAct = new facet();
     1438                        facet *fAct;// = new facet();
    13781439                        fAct = gc.facetPtr;
    1379                        
     1440                        fListPtr = gc.facetPtr;
    13801441#ifdef gfan_DEBUG
    13811442                        cout << "NoRevs" << endl;
    13821443                        cout << "Facets are:" << endl;
    13831444                        gc.showFacets();
    1384 #endif
    1385                        
    1386                         dd_set_global_constants();
     1445#endif                 
     1446                        //dd_set_global_constants();
    13871447                        dd_rowrange ddrows;
    13881448                        dd_colrange ddcols;
    … …  
    14201480                        else/*Facets of facets*/
    14211481                        {
    1422                                 dd_MatrixPtr ddineq,P,ddakt;
    1423                                 dd_rowset impl_linset, redset;
    1424                                 dd_ErrorType err;
    1425                                 dd_rowindex newpos;                     
    14261482                               
    1427                                 fAct = fListPtr;
    1428 
    1429 #ifdef gfan_DEBUG
    1430                                 cout << "before f2M" << endl;
    1431                                 gc.showFacets();
    1432                                 ddineq = facets2Matrix(gc);
    1433                                 cout << "after f2M" << endl;
    1434                                 gc.showFacets();
    1435 #endif
     1483                                this->getCodim2Normals(gc);
    14361484                               
    1437                                 /*Now set appropriate linearity*/
    1438                                 dd_PolyhedraPtr ddpolyh;
    1439                                 for (int ii=0; ii<this->numFacets; ii++)
    1440                                 {       
    1441                                         cout << endl << "------------" << endl;
    1442                                         ddakt = dd_CopyMatrix(ddineq);
    1443                                         set_addelem(ddakt->linset,ii+1);
    1444                                                                                
    1445                                         dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);                   
    1446                                        
    1447                                         //dd_WriteMatrix(stdout,ddakt);
    1448                                         ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
    1449                                         P=dd_CopyGenerators(ddpolyh);
    1450                                         dd_WriteMatrix(stdout,P);
    1451                                        
    1452                                         /* We loop through each row of P
    1453                                         * normalize it by making all entries integer ones
    1454                                         * and add the resulting vector to the int matrix facet::codim2Facets
    1455                                         */
    1456                                         this->facetPtr->codim2Ptr = new facet(2);       //construct a (codim-2)-facet
    1457                                         facet *codim2Act;
    1458                                         codim2Act = this->facetPtr->codim2Ptr;
    1459                                         for (int jj=1;jj<=P->rowsize;jj++)
    1460                                         {
    1461                                                 intvec *n = new intvec(this->numVars);
    1462                                                 normalize(P,jj,*n);                                             
    1463                                                 codim2Act->setFacetNormal(n);
    1464                                                 codim2Act->next = new facet(2);
    1465                                                 codim2Act = codim2Act->next;
    1466                                                 n->show();
    1467                                                 delete n;                                                                       
    1468                                         }                                       
    1469                                        
    1470                                         dd_FreeMatrix(ddakt);
    1471                                         dd_FreePolyhedra(ddpolyh);
    1472                                 }                               
    1473                         }                       
     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                       
    14741500       
    14751501                       
    14761502                        //NOTE Hm, comment in and get a crash for free...
    14771503                        //dd_free_global_constants();                           
    1478                         gc.writeConeToFile(gc);
    1479                        
    1480                         /*2nd step
    1481                         Choose a facet from fListPtr, flip it and forget the previous cone
    1482                         */
    1483                         fAct = fListPtr;
     1504                        //gc.writeConeToFile(gc);
     1505                       
     1506                       
     1507                        //fAct = fListPtr;
    14841508                        //gcone *gcTmp = new gcone(gc); //copy
    14851509                        //gcTmp->flip(gcTmp->gcBasis,fAct);
    … …  
    14911515                /** \brief Make a set of rational vectors into integers
    14921516                *
    1493                 * We compute the lcm of the denominators and multiply with this to get integer values
     1517                * We compute the lcm of the denominators and multiply with this to get integer values.         
    14941518                * \param dd_MatrixPtr,intvec
    14951519                */
    1496                 void normalize(dd_MatrixPtr const &M, int const line, intvec &n)
     1520                void makeInt(dd_MatrixPtr const &M, int const line, intvec &n)
    14971521                {                       
    14981522                        mpz_t denom[this->numVars];
    … …  
    15451569                       
    15461570                }
     1571                /**
     1572                 * We compute the gcd of the components of the codim-2-facets and
     1573                 * multiply the each codim-2facet by it. Thus we get a normalized representation of each
     1574                 * (codim-2)-facet normal.
     1575                */
     1576                void normalize()
     1577                {
     1578                        int ggT=1;
     1579                        facet *codim2Act;
     1580                        codim2Act = this->facetPtr->codim2Ptr;
     1581                        intvec *n = new intvec(this->numVars);
     1582                       
     1583                        while(codim2Act->next!=NULL)
     1584                        {                               
     1585                                n=codim2Act->getFacetNormal();
     1586                                for(int ii=0;ii<this->numVars;ii++)
     1587                                {
     1588                                        ggT = intgcd(ggT,(*n)[ii]);
     1589                                }
     1590                                codim2Act = codim2Act->next;                           
     1591                        }
     1592                        //delete n;
     1593                        codim2Act = this->facetPtr->codim2Ptr;  //reset to start of linked list
     1594                        while(codim2Act->next!=NULL)
     1595                        {
     1596                                //intvec *n = new intvec(this->numVars);
     1597                                n=codim2Act->getFacetNormal();
     1598                                intvec *new_n = new intvec(this->numVars);
     1599                                for(int ii=0;ii<this->numVars;ii++)
     1600                                {
     1601                                        (*new_n)[ii] = (*n)[ii]*ggT;
     1602                                }
     1603                                codim2Act->setFacetNormal(new_n);
     1604                                codim2Act = codim2Act->next;
     1605                                //delete n;
     1606                                //delete new_n;
     1607                        }                       
     1608                }
     1609               
     1610                /** \brief Compute the gcd of two ints
     1611                */
     1612                int intgcd(int a, int b)
     1613                {
     1614                        int r, p=a, q=b;
     1615                        if(p < 0)
     1616                                p = -p;
     1617                        if(q < 0)
     1618                                q = -q;
     1619                        while(q != 0)
     1620                        {
     1621                                r = p % q;
     1622                                p = q;
     1623                                q = r;
     1624                        }
     1625                        return p;
     1626                }               
    15471627               
    15481628                /** \brief Construct a dd_MatrixPtr from a cone's list of facets
    … …  
    15511631                dd_MatrixPtr facets2Matrix(gcone const &gc)
    15521632                {
    1553                         facet *fAct = new facet();
     1633                        facet *fAct;
    15541634                        fAct = gc.facetPtr;
    15551635       
    … …  
    15641644                        //dd_set_global_constants();                                           
    15651645                                               
    1566                         intvec *comp = new intvec(this->numVars);                       
    1567                        
    1568                         for (int jj=0; jj<M->rowsize; jj++)
    1569                         {                               
     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);
    15701651                                comp = fAct->getFacetNormal();
    15711652                                for (int ii=0; ii<this->numVars; ii++)
    … …  
    15771658                                        fAct = fAct->next;
    15781659                                }
     1660                                delete comp;
    15791661                        }//jj
    1580                        
    1581                         //delete fAct;
    1582                         //delete comp;                 
     1662                        */
     1663                        int jj=0;
     1664                        while(fAct->next!=NULL)
     1665                        {
     1666                                intvec *comp;
     1667                                comp = fAct->getFacetNormal();
     1668                                for(int ii=0;ii<this->numVars;ii++)
     1669                                {
     1670                                        dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
     1671                                }
     1672                                jj++;
     1673                                fAct=fAct->next;                               
     1674                        }                       
     1675                       
    15831676                        return M;
    15841677                }
    … …  
    16491742        friend class facet;
    16501743};//class gcone
     1744
     1745int gcone::UCN=0;
    16511746
    16521747ideal gfan(ideal inputIdeal)
Note: See TracChangeset for help on using the changeset viewer.