Changeset e2e3ad in git for kernel/gfan.cc


Ignore:
Timestamp:
Nov 20, 2009, 6:01:20 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', 'e7cc1ebecb61be8b9ca6c18016352af89940b21a')
Children:
1d77dffe7911bafa212ddec6401c158b5357fe6c
Parents:
68b081826cda1d5edc4a5c9adb3f9ead2a95c338
Message:
bugfix in facet::getFacetNormal
gcone::rootRing removed
Introduced int gcone::pred to store predecessor
Cleanup



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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r68b081 re2e3ad  
    205205intvec *facet::getFacetNormal()
    206206{                               
    207         return this->fNormal;
     207        return ivCopy(this->fNormal);
    208208}
    209209               
    … …  
    321321 * This constructor takes the root ring and the root ideal as parameters and stores
    322322 * them in the private members gcone::rootRing and gcone::inputIdeal
     323 * This constructor is only called once in the computation of the Gröbner fan,
     324 * namely for the very first cone. Therefore pred is set to 1.
     325 * Might set it to this->UCN though...
    323326 * Since knowledge of the root ring is only needed when using reverse search,
    324327 * this constructor is not needed when using the "second" method
    … …  
    329332        this->prev=NULL;
    330333        this->facetPtr=NULL;                   
    331         this->rootRing=r;
     334//      this->rootRing=r;
    332335        this->inputIdeal=I;
    333336        this->baseRing=currRing;
    334337        this->counter++;
    335         this->UCN=this->counter;                       
     338        this->UCN=this->counter;
     339        this->pred=1;
    336340        this->numFacets=0;
    337341        this->ivIntPt=NULL;
    … …  
    350354        this->numVars=gc.numVars;                                               
    351355        this->counter++;
    352         this->UCN=this->counter;                       
     356        this->UCN=this->counter;
     357        this->pred=gc.UCN;                     
    353358        this->facetPtr=NULL;
    354359        this->gcBasis=idCopy(f.flipGB);
    … …  
    357362        this->numFacets=0;
    358363        this->ivIntPt=NULL;
    359         this->rootRing=NULL;
     364//      this->rootRing=NULL;
    360365        //rComplete(this->baseRing);
    361366        //rChangeCurrRing(this->baseRing);
    … …  
    370375        if(this->inputIdeal!=NULL)
    371376                idDelete((ideal *)&this->inputIdeal);
    372         if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
    373                 rDelete(this->rootRing);
     377//      if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
     378//              rDelete(this->rootRing);
    374379        //rDelete(this->baseRing);
    375380        facet *fAct;
    … …  
    532537                return -1;
    533538}
     539
     540int gcone::getPredUCN()
     541{
     542        return this->pred;
     543}
     544
    534545/** \brief Compute the normals of the cone
    535546 *
    … …  
    547558void gcone::getConeNormals(ideal const &I, bool compIntPoint)
    548559{
    549 #ifdef gfan_DEBUG
    550         std::cout << "*** Computing Inequalities... ***" << std::endl;
    551 #endif         
    552         //All variables go here - except ineq matrix and *v, see below
    553         int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
    554         int pCompCount;                 // # of terms in a poly
    555560        poly aktpoly;
    556         int numvar = pVariables;        // # of variables in currRing
    557 //      int leadexp[numvar];            // dirty hack of exp.vects
    558 //      int aktexp[numvar];
    559         int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
     561        int rows;                       // will contain the dimensions of the ineq matrix - deprecated by
    560562        dd_rowrange ddrows;
    561563        dd_colrange ddcols;
    … …  
    564566        dd_rowindex ddnewpos;             // all to make dd_Canonicalize happy
    565567        dd_NumberType ddnumb=dd_Integer; //Number type
    566         dd_ErrorType dderr=dd_NoError;  //
    567         // End of var declaration
    568 #ifdef gfan_DEBUG
    569 //      cout << "The Groebner basis has " << lengthGB << " elements" << endl;
    570 //      cout << "The current ring has " << numvar << " variables" << endl;
    571 #endif
    572         cols = numvar;
    573                
     568        dd_ErrorType dderr=dd_NoError;         
    574569        //Compute the # inequalities i.e. rows of the matrix
    575570        rows=0; //Initialization
    … …  
    579574                rows=rows+pLength(aktpoly)-1;
    580575        }
    581 #ifdef gfan_DEBUG
    582 //      cout << "rows=" << rows << endl;
    583 //      cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
    584 #endif
     576
    585577        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
    586578        dd_set_global_constants();
    587579        ddrows=rows;
    588         ddcols=cols;
     580        ddcols=this->numVars;
    589581        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                     
    590582        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
    … …  
    594586        {
    595587                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
    596                 pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
    597 #ifdef gfan_DEBUG
    598 //                              cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
    599 #endif
    600                
    601 //              int *v=(int *)omAlloc((numvar+1)*sizeof(int));
    602 //              pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
    603 //                             
    604 //              //Store leadexp for aktpoly
    605 //              for (int kk=0;kk<numvar;kk++)
    606 //              {
    607 //                      leadexp[kk]=v[kk+1];
    608 //                      //Since we need to know the difference of leadexp with the other expvects we do nothing here
    609 //                      //but compute the diff below
    610 //              }
    611 //             
    612 //                             
    613 //              while (pNext(aktpoly)!=NULL) //move to next term until NULL
    614 //              {
    615 //                      aktpoly=pNext(aktpoly);
    616 //                      pSetm(aktpoly);         //doesn't seem to help anything
    617 //                      pGetExpV(aktpoly,v);
    618 //                     
    619 //                      for (int kk=0;kk<numvar;kk++)
    620 //                      {
    621 //                              aktexp[kk]=v[kk+1];
    622 //                                              //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
    623 //                              dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
    624 //                      }
    625 //                      aktmatrixrow=aktmatrixrow+1;
    626 //              }//while
    627 //              omFree(v);
    628                
    629                 //simpler version of the above
    630                 int *leadexpv=(int*)omAlloc((numvar+1)*sizeof(int));
    631                 int *tailexpv=(int*)omAlloc((numvar+1)*sizeof(int));
     588                //simpler version of storing expvect diffs
     589                int *leadexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
     590                int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
    632591                pGetExpV(aktpoly,leadexpv);
    633592                while(pNext(aktpoly)!=NULL)
    … …  
    635594                        aktpoly=pNext(aktpoly);
    636595                        pGetExpV(aktpoly,tailexpv);
    637                         for(int kk=1;kk<=numvar;kk++)
     596                        for(int kk=1;kk<=this->numVars;kk++)
    638597                        {
    639598                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
    … …  
    642601                }
    643602                omFree(tailexpv);
    644                 omFree(leadexpv);
    645                
     603                omFree(leadexpv);       
    646604        } //for
    647                
    648605#if true
    649606        /*Let's make the preprocessing here. This could already be done in the above for-loop,
    … …  
    666623                        tmp=mpq_get_d(ddineq->matrix[ii][jj]);
    667624                        (*gamma)[jj-1]=(int)tmp;
    668 //                      (*gamma)[jj]=(ddineq->matrix[ii][jj]);
    669625                }
    670626                computeInv((ideal&)I,initialForm,*gamma);
    … …  
    751707        //Maybe add another row to contain the constraints of the standard simplex?
    752708
    753 #ifdef gfan_DEBUG
    754 //      cout << "The inequality matrix is" << endl;
    755 //      dd_WriteMatrix(stdout, ddineq);
    756 #endif
    757 
    758         // The inequalities are now stored in ddineq
    759         // Next we check for superflous rows
    760 //      time_t canonicalizeTic, canonicalizeTac;
    761 //      time(&canonicalizeTic);
    762 //      ddredrows = dd_RedundantRows(ddineq, &dderr);
    763 //      if (dderr!=dd_NoError)                  // did an error occur?
    764 //      {
    765 //              dd_WriteErrorMessages(stderr,dderr);    //if so tell us
    766 //      }
    767 #ifdef gfan_DEBUG
    768 //      else
    769 //      {
    770 //                              cout << "Redundant rows: ";
    771 //                              set_fwrite(stdout, ddredrows);          //otherwise print the redundant rows
    772 //      }//if dd_Error
    773 #endif
    774         //Remove reduntant rows here!
    775         //Necessary check here! C.f. FJT p18
    776                
     709
     710        //Remove redundant rows here!
     711        //Necessary check here! C.f. FJT p18           
    777712        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
    778 //      time(&canonicalizeTac);
    779 //      cout << "dd_MatrixCanonicalize time: " << difftime(canonicalizeTac,canonicalizeTic) << "sec" << endl;
     713        //time(&canonicalizeTac);
     714        //cout << "dd_MatrixCanonicalize time: " << difftime(canonicalizeTac,canonicalizeTic) << "sec" << endl;
    780715        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
    781716        ddcols = ddineq->colsize;
    … …  
    783718        //ddCreateMatrix(ddrows,ddcols+1);
    784719        ddFacets = dd_CopyMatrix(ddineq);
    785 #ifdef gfan_DEBUG
    786 //      cout << "Having removed redundancies, the normals now read:" << endl;
    787 //      dd_WriteMatrix(stdout,ddineq);
    788 //      cout << "Rows = " << ddrows << endl;
    789 //      cout << "Cols = " << ddcols << endl;
    790 #endif
    791720                       
    792721        /*Write the normals into class facet*/
    793 #ifdef gfan_DEBUG
    794 //      cout << "Creating list of normals" << endl;
    795 #endif
    796722        /*The pointer *fRoot should be the return value of this function*/
    797                 //facet *fRoot = new facet();   //instantiate new facet
    798                 //fRoot->setUCN(this->getUCN());
    799                 //this->facetPtr = fRoot;               //set variable facetPtr of class gcone to first facet
     723        //facet *fRoot = new facet();   //instantiate new facet
     724        //fRoot->setUCN(this->getUCN());
     725        //this->facetPtr = fRoot;               //set variable facetPtr of class gcone to first facet
    800726        facet *fAct;    //pointer to active facet
    801727                        //fAct = fRoot;                 //Seems to do the trick. fRoot and fAct have to point to the same adress!
    … …  
    809735                        double foo;
    810736                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
    811 #ifdef gfan_DEBUG
    812 //                                      std::cout << "fAct is " << foo << " at " << fAct << std::endl;
    813 #endif
    814                         (*load)[jj-1] = (int)foo;               //store typecasted entry at pos jj-1 of load           
     737                        (*load)[jj-1] = (int)foo;       //store typecasted entry at pos jj-1 of load           
    815738                }//for (int jj = 1; jj <ddcols; jj++)
    816739                               
    817740                /*Quick'n'dirty hack for flippability*/
    818                 bool isFlip=FALSE;
    819                                                                
     741                bool isFlip=FALSE;                                                             
    820742                for (int jj = 0; jj<load->length(); jj++)
    821743                {
    … …  
    829751                                break;  //URGHS
    830752                        }
    831                 }       
     753                        delete ivCanonical;
     754                }/*End of check for flippability*/
    832755                if (isFlip==FALSE)
    833756                {
    … …  
    873796                        fAct->setUCN(this->getUCN());                                   
    874797                }//if (isFlippable==FALSE)
    875                 //delete load;                         
     798                delete load;                           
    876799        }//for (int kk = 0; kk<ddrows; kk++)
    877800                       
    … …  
    879802        if(numNonFlip==this->numFacets)
    880803        {                                       
    881                 cerr << "Only non-flippable facets. Terminating..." << endl;
     804                WerrorS ("Only non-flippable facets. Terminating...\n");
    882805        }
    883806                       
    … …  
    902825                interiorPoint(ddineq, *iv);     //NOTE ddineq contains non-flippable facets
    903826                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
    904                 //delete iv;
     827                delete iv;
    905828        }
    906829                       
    … …  
    913836        dd_FreeMatrix(ddineq);
    914837        set_free(ddredrows);
     838        set_free(ddlinset);
    915839        //free(ddnewpos);
    916840        //set_free(ddlinset);
    … …  
    952876                set_addelem(ddakt->linset,ii+1);                               
    953877                dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);                   
    954                                        
    955 #ifdef gfan_DEBUG
    956 //              cout << "Codim2 matrix"<<endl;
    957 //              dd_WriteMatrix(stdout,ddakt);
    958 #endif
    959878                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
    960879                P=dd_CopyGenerators(ddpolyh);                           
    961 #ifdef gfan_DEBUG
    962 //              cout << "Codim2 facet:" << endl;
    963 //                      dd_WriteMatrix(stdout,P);
    964 //              cout << endl;
    965 #endif
    966                                        
    967                 /* We loop through each row of P
    968                 * normalize it by making all entries integer ones
    969                 * and add the resulting vector to the int matrix facet::codim2Facets
    970                 */
     880                /* We loop through each row of P normalize it by making all
     881                * entries integer ones and add the resulting vector to the
     882                * int matrix facet::codim2Facets */
    971883                for (int jj=1;jj<=P->rowsize;jj++)
    972884                {                                       
    … …  
    1000912                        dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
    1001913                }
    1002                 intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
    1003                                 //dd_WriteMatrix(stdout,intPointMatrix);
     914                intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);             
    1004915                interiorPoint(intPointMatrix,*iv_intPoint);
    1005916                for(int ll=0;ll<this->numVars;ll++)
    … …  
    1014925                fAct = fAct->next;     
    1015926                dd_FreeMatrix(ddakt);
    1016 //              dd_FreeMatrix(ddineq);
     927//              dd_FreeMatrix(ddineq);
     928                dd_FreeMatrix(shiftMatrix);
     929                dd_FreeMatrix(intPointMatrix);
    1017930                dd_FreePolyhedra(ddpolyh);
    1018931                delete iv_intPoint;
    1019         }//while
     932        }//for
     933        dd_FreeMatrix(ddineq);
    1020934}
    1021935               
    … …  
    1037951        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
    1038952        fNormal = f->getFacetNormal();  //read this->fNormal;
    1039 //#ifdef gfan_DEBUG
    1040                         //std::cout << "===" << std::endl;
     953
    1041954        std::cout << "running gcone::flip" << std::endl;
    1042955        std::cout << "flipping UCN " << this->getUCN() << endl;
    1043 //                      for(int ii=0;ii<IDELEMS(gb);ii++)       //not very handy with large examples
    1044 //                      {
    1045 //                              pWrite((poly)gb->m[ii]);
    1046 //                      }
    1047956        cout << "over facet (";
    1048957        fNormal->show(1,0);
    1049958        cout << ") with UCN " << f->getUCN();
    1050959        std::cout << std::endl;
    1051 //#endif                               
     960
    1052961        /*1st step: Compute the initial ideal*/
    1053962        /*poly initialFormElement[IDELEMS(gb)];*/       //array of #polys in GB to store initial form
    1054963        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
    1055         poly aktpoly;
    1056         /*intvec *check = new intvec(this->numVars);*/  //array to store the difference of LE and v
     964        poly aktpoly;   
    1057965        computeInv(gb,initialForm,*fNormal);
    1058966        //NOTE The code below went into gcone::computeInv
    … …  
    1062970        //f->printFlipGB();*/
    1063971//      cout << "===" << endl;
    1064 #endif
    1065         //delete check;
    1066                        
     972#endif                 
    1067973        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
    1068974        /*Substep 2.1
    … …  
    1093999                rComplete(tmpRing);     
    10941000        }
     1001        delete fNormal; //NOTE: Why does this crash?
    10951002        rChangeCurrRing(tmpRing);
    10961003                       
    … …  
    12661173        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
    12671174        // Thus:
    1268         //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
    1269         //cout << "PLING" << endl;
    1270         /*ring dstRing=rCopy0(srcRing);
    1271         for (int ii=0;ii<this->numVars;ii++)
    1272         {                               
    1273         dstRing->wvhdl[0][ii]=(*iv_weight)[ii];                         
    1274         }*/
     1175        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
    12751176        ring dstRing=rCopy0(tmpRing);
    12761177        int length=iv_weight->length();
    … …  
    12851186        //rDelete(tmpRing);
    12861187        delete iv_weight;
    1287 //#ifdef gfan_DEBUG
    1288         rWrite(dstRing); cout << endl;
    1289 //#endif
     1188
     1189        //rWrite(dstRing); cout << endl;
     1190
    12901191        ideal dstRing_I;                       
    12911192        dstRing_I=idrCopyR(srcRing_HH,srcRing);
    1292                         //dstRing_I=idrCopyR(inputIdeal,srcRing);
    1293                         //validOpts<1>=TRUE;
     1193        //dstRing_I=idrCopyR(inputIdeal,srcRing);
     1194        //validOpts<1>=TRUE;
    12941195#ifdef gfan_DEBUG
    1295                         //idShow(dstRing_I);
     1196        //idShow(dstRing_I);
    12961197#endif                 
    12971198        BITSET save=test;
    … …  
    13161217//#ifdef gfan_DEBUG
    13171218        cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
    1318 //      f->printFlipGB();
    1319         this->idDebugPrint(dstRing_I);
    1320         cout << endl;
     1219//      this->idDebugPrint(dstRing_I);
     1220//      cout << endl;
    13211221//#endif                       
    13221222        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
    … …  
    14381338intvec *gcone::ivNeg(const intvec *iv)
    14391339{
     1340        //NOTE: Can't this be done without new?
    14401341        intvec *res = new intvec(iv->length());
    14411342        res=ivCopy(iv);
    … …  
    17901691          We always choose the first facet from fListPtr as facet to be flipped
    17911692        */                     
    1792         time_t tic, tac;
     1693//      time_t tic, tac;
    17931694        while((SearchListAct!=NULL))// && counter<10)
    17941695        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
    … …  
    18041705                        rComplete(rTmp);                       
    18051706                        rChangeCurrRing(rTmp);
    1806                         gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);
     1707                        gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);//copy constructor!
    18071708                        /* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
    18081709                         * Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
    … …  
    26792580                res->m[ii].rtyp=LIST_CMD;
    26802581                lists l=(lists)omAllocBin(slists_bin);
    2681                 l->Init(5);
     2582                l->Init(6);
    26822583                l->m[0].rtyp=INT_CMD;
    26832584                l->m[0].data=(void*)gcAct->getUCN();
    … …  
    27082609                l->m[4].rtyp=RING_CMD;
    27092610                l->m[4].data=(void*)(gcAct->getBaseRing());             
    2710                
     2611                l->m[5].rtyp=INT_CMD;
     2612                l->m[5].data=(void*)gcAct->getPredUCN();
    27112613                res->m[ii].data=(void*)l;
    27122614                gcAct = gcAct->next;
Note: See TracChangeset for help on using the changeset viewer.