Changeset f5b077 in git for kernel/gfan.cc


Ignore:
Timestamp:
Apr 6, 2009, 4:57:18 PM (15 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
1288efb1bc6f003f059484c99a625cf47c1ed08b
Parents:
f3c31a15740169790d1f07f8cb43045dd6f97991
Message:
facet::flibGB is now private; set & print method exist
gcone::flip now computes the initial ideal of the respective facet.
Started work on the lifting step. Not sure how to set wvhdl though.
Made some methods call-by-reference
Switching of rings works for the root ring


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    rf3c31a rf5b077  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-04-03 13:49:16 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.28 2009-04-03 13:49:16 monerjan Exp $
    6 $Id: gfan.cc,v 1.28 2009-04-03 13:49:16 monerjan Exp $
     4$Date: 2009-04-06 14:57:18 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.29 2009-04-06 14:57:18 monerjan Exp $
     6$Id: gfan.cc,v 1.29 2009-04-06 14:57:18 monerjan Exp $
    77*/
    88
    … …  
    5252{
    5353        private:
    54                 /** inner normal, describing the facet uniquely up to isomorphism */
    55                 intvec *fNormal;               
     54                /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
     55                intvec *fNormal;
     56               
     57                /** \brief The Groebner basis on the other side of a shared facet
     58                 *
     59                 * In order not to have to compute the flipped GB twice we store the basis we already get
     60                 * when identifying search facets. Thus in the next step of the reverse search we can
     61                 * just copy the old cone and update the facet and the gcBasis.
     62                 * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
     63                 */
     64                ideal flipGB;           //The Groebner Basis on the other side, computed via gcone::flip
     65
     66                       
    5667        public:
    5768                /** The default constructor. Do I need a constructor of type facet(intvec)? */
    … …  
    8596                }
    8697               
    87                 /** \brief The Groebner basis on the other side of a shared facet
    88                 *
    89                 * In order not to have to compute the flipped GB twice we store the basis we already get
    90                 * when identifying search facets. Thus in the next step of the reverse search we can
    91                 * just copy the old cone and update the facet and the gcBasis
    92                 */
    93                 ideal flibGB;           //The Groebner Basis on the other side, computed via gcone::flip
     98                /** Store the flipped GB*/
     99                void setFlipGB(ideal I)
     100                {
     101                        this->flipGB=I;
     102                }
     103               
     104                /** Print the flipped GB*/
     105                void printFlipGB()
     106                {
     107                        idShow(this->flipGB);
     108                }
    94109               
    95110                bool isFlippable;       //flippable facet? Want to have cone->isflippable.facet[i]
    … …  
    143158                * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
    144159                */
    145                 void getConeNormals(ideal I)
     160                void getConeNormals(const ideal &I)
    146161                {
    147162#ifdef gfan_DEBUG
    … …  
    330345                        /*1st step: Compute the initial ideal*/
    331346                        poly initialFormElement[IDELEMS(gb)];   //array of #polys in GB to store initial form
    332                         ideal initialForm;
     347                        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
    333348                        poly aktpoly;
    334349                        intvec *check = new intvec(this->numVars);      //array to store the difference of LE and v
    … …  
    340355                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
    341356                                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
    342                                 //pGetExpV(aktpoly,ivLeadExpV);
    343357                                initialFormElement[ii]=pHead(aktpoly);
    344358                               
    … …  
    355369                                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
    356370                                        }
     371#ifdef gfan_DEBUG
    357372                                        cout << "check=";                       
    358373                                        check->show();
    359374                                        cout << endl;
     375#endif
    360376                                        if (isParallel(check,fNormal)) //pass *check when
    361377                                        {
    … …  
    367383                                pWrite(initialFormElement[ii]);
    368384                                cout << "---" << endl;
    369                                 /*Now initialFormElement must be added to (ideal)initialForm */                                                 
     385                                /*Now initialFormElement must be added to (ideal)initialForm */
     386                                //f->flibGB->m[ii]=(poly)initialFormElement[ii];
     387                                //(poly)initialForm->m[ii]=pAdd(initialForm[ii],initialFormElement[ii]);
     388                                initialForm->m[ii]=initialFormElement[ii];
    370389                        }//for
     390                        f->setFlipGB(initialForm);                     
     391#ifdef gfan_DEBUG
     392                        cout << "Initial ideal is: " << endl;
     393                        //idShow(initialForm);
     394                        f->printFlipGB();
     395                        cout << "===" << endl;
     396#endif
     397                        delete check;
     398                       
     399                        /*2nd step: lift initial ideal to a GB of the neighbouring cone*/
     400                        ring tmpring=rCopy0(currRing);
     401                        tmpring->order[0]=ringorder_a;
     402                        tmpring->order[1]=ringorder_dp;
     403                        tmpring->order[2]=ringorder_C;
     404                        //rWrite(tmpring);
     405                        tmpring->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
     406                       
     407                        for (int ii=0;ii<this->numVars;ii++)
     408                        {
     409                                cout << "ping" << endl;
     410                                tmpring->wvhdl[0][ii]=(*fNormal)[ii];   //What exactly am I doing here?
     411                                cout << "pong" << endl;
     412                                cout << *tmpring->wvhdl[ii] << endl;
     413                        }
     414                        //tmpring->wvhdl=(int**)(fNormal);
     415                        rComplete(tmpring);
     416                        rWrite(tmpring);
     417                        /*setring(tmpring);
     418                        ideal ina=initialForm;
     419                        ideal H;
     420                        H=kstd(ina,NULL,testHomog,NULL);
     421                        idShow(H);*/
     422                       
    371423                }//void flip(ideal gb, facet *f)
    372424                               
    … …  
    377429                *\return void
    378430                */
    379                 void getGB(ideal inputIdeal)
     431                void getGB(const ideal &inputIdeal)
    380432                {
    381433                        ideal gb;
    … …  
    410462                bool isParallel(intvec *a, intvec *b)
    411463                {                       
    412                         //bool res=FALSE;
    413                         //intvec iva(this->numVars);
    414                         //intvec ivb(this->numVars);
    415464                        int lhs,rhs;
    416                         //lhs=0;
    417                         //rhs=0;
    418                         //iva = a;
    419                         //ivb = b;
    420                         //lhs=dotProduct(&iva,&ivb)*dotProduct(&iva,&ivb);
    421                         //lhs=dotProduct(&iva,&iva)*dotProduct(&ivb,&ivb);
    422465                        lhs=dotProduct(&a,&b)*dotProduct(&a,&b);
    423466                        rhs=dotProduct(&a,&a)*dotProduct(&b,&b);
    … …  
    476519        3. getConeNormals
    477520        */
    478 
    479521       
    480522        /* Construct a new ring which will serve as our root
    … …  
    482524        */
    483525        rootRing=rCopy0(currRing);
    484         rootRing->order[0]=ringorder_dp;
     526        rootRing->order[0]=ringorder_lp;
    485527        rComplete(rootRing);
    486         //rChangeCurrRing(rootRing);
     528        rChangeCurrRing(rootRing);
     529       
     530        /* Fetch the inputIdeal into our rootRing */
     531        map theMap=(map)idMaxIdeal(1);  //evil hack!
     532        //idShow(idMaxIdeal(1));
     533        /*for (int ii=0;ii<pVariables;ii++)
     534        {
     535                theMap->m[ii]=inputIdeal->m[ii];
     536        }*/
     537        theMap->preimage=NULL;
    487538        ideal rootIdeal;
    488         /* Fetch the inputIdeal into our rootRing */
    489         map m=(map)idInit(IDELEMS(inputIdeal),0);
    490         //rootIdeal=fast_map(inputIdeal,inputRing,(ideal)m, currRing);
     539        rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
    491540#ifdef gfan_DEBUG
    492541        cout << "Root ideal is " << endl;
    493         idPrint(rootIdeal);
    494         cout << "The current ring is " << endl;
     542        idShow(rootIdeal);
     543        cout << "The root ring is " << endl;
    495544        rWrite(rootRing);
    496545        cout << endl;
    … …  
    501550        gcAct = gcRoot;
    502551        gcAct->numVars=pVariables;
    503         gcAct->getGB(inputIdeal);
     552        gcAct->getGB(rootIdeal);        //sets gcone::gcBasis
     553        idShow(gcAct->gcBasis);
    504554        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
    505555        gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
Note: See TracChangeset for help on using the changeset viewer.