Changeset 0ad81f in git for kernel/gfan.cc


Ignore:
Timestamp:
Apr 7, 2009, 11:44:20 AM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
Children:
2ae96e40fc5453bcb155aec76d376d79dd549cbe
Parents:
f32177e9dea6c9833a7f2f8986a2cf6273f30fab
Message:
Ring construction & change works for the lifting step in gcone::flip;


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    rf32177 r0ad81f  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    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 $
     4$Date: 2009-04-07 09:44:20 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.30 2009-04-07 09:44:20 monerjan Exp $
     6$Id: gfan.cc,v 1.30 2009-04-07 09:44:20 monerjan Exp $
    77*/
    88
    … …  
    267267                        cout << "Cols = " << ddcols << endl;
    268268#endif
     269                       
    269270                        /*Write the normals into class facet*/
    270271#ifdef gfan_DEBUG
    … …  
    275276                        this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
    276277                        facet *fAct;                    //instantiate pointer to active facet
    277                         fAct = fRoot;           //This does not seem to do the trick. fRoot and fAct have to point to the same adress!
     278                        fAct = fRoot;                   //Seems to do the trick. fRoot and fAct have to point to the same adress!
    278279                        std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
    279280                        for (int kk = 0; kk<ddrows; kk++)
    … …  
    291292                                        if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
    292293                                        {
    293                                                 fAct->next = new facet();       //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                               
     294                                                //fAct->next = new facet();     //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                               
    294295                                        }
    295                                 }//for jj<ddcols
    296                                 /*Now load should be full and we can call setFacetNormal*/
    297                                 fAct->setFacetNormal(load);
    298                                 //fAct->printNormal();
    299                                 fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
     296                                }//for (int jj = 1; jj <ddcols; jj++)
     297                                /*Quick'n'dirty hack for flippability*/
     298                                bool isFlippable;                       
     299                                for (int jj = 0; jj<this->numVars; jj++)
     300                                {                                       
     301                                        intvec *ivCanonical = new intvec(this->numVars);
     302                                        (*ivCanonical)[jj]=1;                                   
     303                                        if (dotProduct(load,ivCanonical)>=0)
     304                                        {
     305                                                isFlippable=FALSE;                                             
     306                                        }
     307                                        else
     308                                        {
     309                                                isFlippable=TRUE;
     310                                        }                                       
     311                                        delete ivCanonical;
     312                                }//for (int jj = 0; jj<this->numVars; jj++)
     313                                if (isFlippable==FALSE)
     314                                {
     315                                        cout << "Ignoring facet";
     316                                        load->show();
     317                                        //fAct->next=NULL;
     318                                }
     319                                else
     320                                {       /*Now load should be full and we can call setFacetNormal*/
     321                                        fAct->setFacetNormal(load);
     322                                        fAct->next = new facet();
     323                                        //fAct->printNormal();
     324                                        fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
     325                                }//if (isFlippable==FALSE)
    300326                                delete load;
    301                         }
     327                        }//for (int kk = 0; kk<ddrows; kk++)
    302328                        /*
    303329                        Now we should have a linked list containing the facet normals of those facets that are
    … …  
    316342                }//method getConeNormals(ideal I)       
    317343               
    318                 /*bool isParallel(int v[], intvec iv)
    319                 {
    320 }*/
    321344               
    322345                /** \brief Compute the Groebner Basis on the other side of a shared facet
    … …  
    335358                void flip(ideal gb, facet *f)           //Compute "the other side"
    336359                {                       
    337                        
    338360                        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
    339361                        fNormal = f->getFacetNormal();  //read this->fNormal;
    340362#ifdef gfan_DEBUG
     363                        cout << "===" << endl;
     364                        cout << "running gcone::flip" << endl;
    341365                        cout << "fNormal=";
    342366                        fNormal->show();
    … …  
    384408                                cout << "---" << endl;
    385409                                /*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]);
    388410                                initialForm->m[ii]=initialFormElement[ii];
    389411                        }//for
    … …  
    391413#ifdef gfan_DEBUG
    392414                        cout << "Initial ideal is: " << endl;
    393                         //idShow(initialForm);
     415                        idShow(initialForm);
    394416                        f->printFlipGB();
    395417                        cout << "===" << endl;
    … …  
    397419                        delete check;
    398420                       
    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
     421                        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
     422                        /*Substep 2.1
     423                        compute $G_{-\alpha}(in_v(I))
     424                        see journal p. 66
     425                        */
     426                        ring srcRing=currRing;
     427                        ring dstRing=rCopy0(srcRing);
     428                        dstRing->order[0]=ringorder_a;
     429                        //tmpring->order[1]=ringorder_dp;
     430                        //tmpring->order[2]=ringorder_C;
     431                        dstRing->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
    406432                       
    407433                        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;
     434                        {                               
     435                                dstRing->wvhdl[0][ii]=-(*fNormal)[ii];  //What exactly am I doing here?
     436                                //cout << tmpring->wvhdl[0][ii] << endl;
    413437                        }
    414                         //tmpring->wvhdl=(int**)(fNormal);
    415                         rComplete(tmpring);
    416                         rWrite(tmpring);
    417                         /*setring(tmpring);
    418                         ideal ina=initialForm;
     438                        rComplete(dstRing);
     439                        rChangeCurrRing(dstRing);
     440                        map theMap=(map)idMaxIdeal(1);
     441                        rWrite(currRing); cout << endl;
     442                        ideal ina;
     443                        ina=fast_map(initialForm,srcRing,(ideal)theMap,dstRing);                       
     444                        cout << "ina=";
     445                        idShow(ina); cout << endl;
    419446                        ideal H;
    420                         H=kstd(ina,NULL,testHomog,NULL);
    421                         idShow(H);*/
     447                        H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
     448                        idSkipZeroes(H);
     449                        cout << "H="; idShow(H); cout << endl;
     450                        /*Substep 2.2
     451                        do the lifting
     452                        */
     453                        rChangeCurrRing(srcRing);
     454                        ideal srcRing_H;
     455                        ideal srcRing_HH;
     456                        //map theMap = (map)idMaxIdeal(1);
     457                        srcRing_H=fast_map(H,dstRing,(ideal)theMap,srcRing);
     458                        srcRing_HH=srcRing_H-stdred(srcRing_H,this->gcBasis);
     459                        /*Substep 2.3
     460                        turn the minimal basis into a reduced one
     461                        */
    422462                       
    423463                }//void flip(ideal gb, facet *f)
    … …  
    429469                *\return void
    430470                */
    431                 void getGB(const ideal &inputIdeal)
     471                void getGB(ideal const &inputIdeal)
    432472                {
    433473                        ideal gb;
    … …  
    435475                        idSkipZeroes(gb);
    436476                        this->gcBasis=gb;       //write the GB into gcBasis
    437                 }
     477                }//void getGB
    438478               
    439479                ideal GenGrbWlk(ideal, ideal);  //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
    … …  
    442482                /** \brief Compute the dot product of two intvecs
    443483                *
    444                 * THIS IS WEIRD - BUT WORKS
    445484                */
    446                 int dotProduct(intvec **a, intvec **b)                         
    447                 {
    448                         intvec iva=*a;
    449                         intvec ivb=*b;
     485                int dotProduct(intvec const &iva, intvec const &ivb)                           
     486                {
     487                        //intvec iva=a;
     488                        //intvec ivb=b;
    450489                        int res=0;
    451490                        for (int i=0;i<this->numVars;i++)
    … …  
    454493                        }
    455494                        return res;
    456                 }
     495                }//int dotProduct
    457496
    458497                /** \brief Check whether two intvecs are parallel
    … …  
    460499                * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
    461500                */
    462                 bool isParallel(intvec *a, intvec *b)
     501                bool isParallel(intvec const &a, intvec const &b)
    463502                {                       
    464503                        int lhs,rhs;
    465                         lhs=dotProduct(&a,&b)*dotProduct(&a,&b);
    466                         rhs=dotProduct(&a,&a)*dotProduct(&b,&b);
     504                        lhs=dotProduct(a,b)*dotProduct(a,b);
     505                        rhs=dotProduct(a,a)*dotProduct(b,b);
    467506                        cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
    468507                        if (lhs==rhs)
    … …  
    475514                        }
    476515                }//bool isParallel
    477                
    478                 /** \brief Check two intvecs for equality --- PROBABLY NOT NEEDED
    479                 *
    480                 * \f$ \alpha=\beta\Leftrightarrow\langle\alpha-\beta,\alpha-\beta\rangle=0 \f$
    481                 */
    482                 bool isEqual(intvec a, intvec b)
    483                 {
    484                         intvec *ivdiff;
    485                         int res;
    486                        
    487                         ivdiff=ivSub(&a,&b);
    488                         cout << "ivdiff=";
    489                         ivdiff->show();
    490                         cout << endl;
    491                         //res=dotProduct(ivdiff,ivdiff);
    492                         if (res==0)
    493                         {
    494                                 return TRUE;
    495                         }
    496                         else
    497                         {
    498                                 return FALSE;
    499                         }                       
    500                 }//bool isEqual
    501516};//class gcone
    502517
    … …  
    566581        => Count the cones!
    567582        */
    568        
     583        rChangeCurrRing(inputRing);
    569584        res=gcAct->gcBasis;
    570585        //cout << fRoot << endl;
Note: See TracChangeset for help on using the changeset viewer.