Changeset bb503c7 in git for kernel/gfan.cc


Ignore:
Timestamp:
Oct 20, 2009, 5:14:02 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
a388ae5c0bb2238bcecc75e032fb91a5a79c3b00
Parents:
d2cd5158d8448e8c4da85855996e71834afaec9d
Message:
Cleanup
idNorm in gcone::getGB


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    rd2cd515 rbb503c7  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-10-20 10:39:18 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.98 2009-10-20 10:39:18 monerjan Exp $
    6 $Id: gfan.cc,v 1.98 2009-10-20 10:39:18 monerjan Exp $
     4$Date: 2009-10-20 15:14:02 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.99 2009-10-20 15:14:02 monerjan Exp $
     6$Id: gfan.cc,v 1.99 2009-10-20 15:14:02 monerjan Exp $
    77*/
    88
    … …  
    7575*/
    7676               
    77 /** The default constructor. Do I need a constructor of type facet(intvec)? */
     77/** \brief The default constructor for facets
     78* Do I need a constructor of type facet(intvec)?
     79*/
    7880facet::facet()                 
    7981{
    … …  
    266268        cout << "Facet normal = (";
    267269        fNormal->show(1,1);
    268         cout << ")"<<endl;
    269         /*for(int ii=0;ii<pVariables;ii++)
    270         {
    271                 cout << (*fNormal)[ii] << ",";
    272         }
    273         if(this->isFlippable==TRUE)
    274                 cout << ")" << endl;
    275         else
    276                 cout << ")*" << endl;*/ //This case should never happen in SLA!
     270        cout << ")"<<endl;     
    277271        cout << "-----------------------" << endl;
    278272        cout << "Codim2 facets:" << endl;
    … …  
    282276                cout << "(";
    283277                f2Normal->show(1,0);
    284 //              for(int jj=0;jj<pVariables;jj++)
    285 //              {
    286 //                              cout << (*f2Normal)[jj] << ",";
    287 //              }
    288278                cout << ")" << endl;
    289279                codim2Act = codim2Act->next;
    … …  
    444434        facet *codim2Act;
    445435        codim2Act = fAct->codim2Ptr;
    446         intvec *fNormal;// = new intvec(this->numVars);
     436        intvec *fNormal;
    447437        intvec *f2Normal;
    448438        cout << endl;
    … …  
    520510        std::cout << "*** Computing Inequalities... ***" << std::endl;
    521511#endif         
    522                         //All variables go here - except ineq matrix and *v, see below
     512        //All variables go here - except ineq matrix and *v, see below
    523513        int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
    524514        int pCompCount;                 // # of terms in a poly
    … …  
    535525        dd_NumberType ddnumb=dd_Integer; //Number type
    536526        dd_ErrorType dderr=dd_NoError;  //
    537                         // End of var declaration
     527        // End of var declaration
    538528#ifdef gfan_DEBUG
    539529//                      cout << "The Groebner basis has " << lengthGB << " elements" << endl;
    … …  
    560550        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
    561551               
    562                         // We loop through each g\in GB and compute the resulting inequalities
     552        // We loop through each g\in GB and compute the resulting inequalities
    563553        for (int i=0; i<IDELEMS(I); i++)
    564554        {
    … …  
    572562                pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
    573563                               
    574                                 //Store leadexp for aktpoly
     564                //Store leadexp for aktpoly
    575565                for (int kk=0;kk<numvar;kk++)
    576566                {
    577567                        leadexp[kk]=v[kk+1];
    578                                         //Since we need to know the difference of leadexp with the other expvects we do nothing here
    579                                         //but compute the diff below
     568                        //Since we need to know the difference of leadexp with the other expvects we do nothing here
     569                        //but compute the diff below
    580570                }
    581571               
    … …  
    604594#endif
    605595
    606                         // The inequalities are now stored in ddineq
    607                         // Next we check for superflous rows
     596        // The inequalities are now stored in ddineq
     597        // Next we check for superflous rows
    608598        ddredrows = dd_RedundantRows(ddineq, &dderr);
    609599        if (dderr!=dd_NoError)                  // did an error occur?
    … …  
    640630                //fRoot->setUCN(this->getUCN());
    641631                //this->facetPtr = fRoot;               //set variable facetPtr of class gcone to first facet
    642         facet *fAct;                    //pointer to active facet
     632        facet *fAct;    //pointer to active facet
    643633                        //fAct = fRoot;                 //Seems to do the trick. fRoot and fAct have to point to the same adress!
    644634                        //std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
    … …  
    715705                        fAct->setUCN(this->getUCN());                                   
    716706                }//if (isFlippable==FALSE)
    717                                 //delete load;                         
     707                //delete load;                         
    718708        }//for (int kk = 0; kk<ddrows; kk++)
    719709                       
    720                         //In cases like I=<x-1,y-1> there are only non-flippable facets...
     710        //In cases like I=<x-1,y-1> there are only non-flippable facets...
    721711        if(numNonFlip==this->numFacets)
    722712        {                                       
    … …  
    747737        }
    748738                       
    749                 //Compute the number of facets
    750                 // wrong because ddineq->rowsize contains only irredundant facets. There can still be
    751                 // non-flippable ones!
    752                 //this->numFacets=ddineq->rowsize;
    753                
    754                 //Clean up but don't delete the return value! (Whatever it will turn out to be)                 
    755                 dd_FreeMatrix(ddineq);
    756                 set_free(ddredrows);
    757                 //free(ddnewpos);
    758                 //set_free(ddlinset);
    759                 //NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
    760                 //THIS SUCKS
    761                 //dd_free_global_constants();
     739        //Compute the number of facets
     740        // wrong because ddineq->rowsize contains only irredundant facets. There can still be
     741        // non-flippable ones!
     742        //this->numFacets=ddineq->rowsize;
     743       
     744        //Clean up but don't delete the return value! (Whatever it will turn out to be)                 
     745        dd_FreeMatrix(ddineq);
     746        set_free(ddredrows);
     747        //free(ddnewpos);
     748        //set_free(ddlinset);
     749        //NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
     750        //THIS SUCKS
     751        //dd_free_global_constants();
    762752
    763753}//method getConeNormals(ideal I)       
    … …  
    770760void gcone::getCodim2Normals(gcone const &gc)
    771761{
    772                         //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet
     762        //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet
    773763        facet *fAct;
    774764        fAct = this->facetPtr;         
    775765        facet *codim2Act;
    776                         //codim2Act = this->facetPtr->codim2Ptr;
     766        //codim2Act = this->facetPtr->codim2Ptr;
    777767                       
    778768        dd_set_global_constants();
    … …  
    783773        dd_rowindex newpos;             
    784774
    785                         //ddineq = facets2Matrix(gc);   //get a matrix representation of the cone
     775        //ddineq = facets2Matrix(gc);   //get a matrix representation of the cone
    786776        ddineq = dd_CopyMatrix(gc.ddFacets);
    787777                               
    … …  
    923913//                                      cout << endl;
    924914#endif
    925                                         //TODO why not *check, *fNormal????
    926915                        if (isParallel(*check,*fNormal)) //pass *check when
    927916                        {
    928 //                                              cout << "Parallel vector found, adding to initialFormElement" << endl;                 
     917//                              cout << "Parallel vector found, adding to initialFormElement" << endl;                 
    929918                                initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
    930919                        }                                               
    … …  
    939928        }//for                 
    940929#ifdef gfan_DEBUG
    941 /*                      cout << "Initial ideal is: " << endl;
     930/*      cout << "Initial ideal is: " << endl;
    942931        idShow(initialForm);
    943932        //f->printFlipGB();*/
    944 //                      cout << "===" << endl;
    945 #endif
    946                         //delete check;
     933//      cout << "===" << endl;
     934#endif
     935        //delete check;
    947936                       
    948937        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
    949                         /*Substep 2.1
     938        /*Substep 2.1
    950939        compute $G_{-\alpha}(in_v(I))
    951940        see journal p. 66
    952941        NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
    953942        srcRing already has a weighted ordering
    954                         */
     943        */
    955944        ring srcRing=currRing;
    956945        ring tmpRing;
    … …  
    976965        rChangeCurrRing(tmpRing);
    977966                       
    978                         //rWrite(currRing); cout << endl;
     967        //rWrite(currRing); cout << endl;
    979968                       
    980969        ideal ina;                     
    … …  
    991980//                      cout << "H="; idShow(H); cout << endl;
    992981#endif
    993                         /*Substep 2.2
     982        /*Substep 2.2
    994983        do the lifting and mark according to H
    995                         */
     984        */
    996985        rChangeCurrRing(srcRing);
    997986        ideal srcRing_H;
    … …  
    1007996//                      idShow(srcRing_HH); cout << endl;
    1008997#endif
    1009                         /*Substep 2.2.1
    1010         Mark according to G_-\alpha
    1011         Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
    1012         we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
    1013         represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
    1014         Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
    1015         compute the difference accordingly
    1016                         */
     998        /*Substep 2.2.1
     999         * Mark according to G_-\alpha
     1000         * Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
     1001         * we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
     1002         * represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
     1003         * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
     1004         * compute the difference accordingly
     1005        */
    10171006        dd_set_global_constants();
    10181007        bool markingsAreCorrect=FALSE;
    … …  
    10251014                iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
    10261015        }
    1027                         /* additionally one row for the standard-simplex and another for a row that becomes 0 during
    1028         construction of the differences
    1029                         */
     1016        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
     1017         * construction of the differences
     1018         */
    10301019        intPointMatrix = dd_CreateMatrix(iPMatrixRows+10,this->numVars+1); //iPMatrixRows+10;
    10311020        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
    … …  
    10821071                while (pNext(aktpoly)!=NULL)
    10831072                {
    1084                                         /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
     1073                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
    10851074                        is not omitted when computing the differences*/
    10861075                        if(markingsAreCorrect==TRUE)
    … …  
    11121101                dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
    11131102        }
    1114                         //Let's make sure we compute interior points from the positive orthant
     1103        //Let's make sure we compute interior points from the positive orthant
    11151104        dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
    11161105        int jj=1;
    … …  
    11871176//#ifdef gfan_DEBUG
    11881177        cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
    1189         f->printFlipGB();
    1190 //      this->idDebugPrint(dstRing_I);
     1178//      f->printFlipGB();
     1179        this->idDebugPrint(dstRing_I);
    11911180        cout << endl;
    11921181//#endif                       
    … …  
    12991288        ideal gb;
    13001289        gb=kStd(inputIdeal,NULL,testHomog,NULL);
     1290        idNorm(gb);
    13011291        idSkipZeroes(gb);
    13021292        this->gcBasis=gb;       //write the GB into gcBasis
    … …  
    24032393        string UCNstr = ss.str();               
    24042394                       
    2405         char prefix[]="/tmp/cone";
    2406         char *UCNstring;
    2407         strcpy(UCNstring,UCNstr.c_str());
    2408         char suffix[]=".sg";
    2409         char *filename=strcat(prefix,UCNstring);
    2410         strcat(filename,suffix);
     2395        string prefix="/tmp/cone";
     2396        string suffix=".sg";
     2397        string filename;
     2398        filename.append(prefix);
     2399        filename.append(UCNstr);
     2400        filename.append(suffix);
    24112401                                       
    2412         ofstream gcOutputFile(filename);
    2413         facet *fAct;// = new facet(); //NOTE Why new?
     2402        ofstream gcOutputFile(filename.c_str());
     2403        facet *fAct;
    24142404        fAct = gc.facetPtr;                     
    24152405        facet *f2Act;
    … …  
    24382428                                       
    24392429                gcOutputFile << "FACETS" << endl;                                                               
    2440                 //while(fAct->next!=NULL)
     2430               
    24412431                while(fAct!=NULL)
    24422432                {       
    2443                         intvec *iv;// = new intvec(gc.numVars);
     2433                        intvec *iv;
    24442434                        intvec *iv2;
    24452435                        iv=fAct->getFacetNormal();
    … …  
    24742464                        gcOutputFile << endl;
    24752465                        fAct=fAct->next;
    2476 //                      delete iv; iv=NULL;
    24772466                }                       
    2478                                
    24792467        gcOutputFile.close();
    2480 //      delete fAct; fAct=NULL;
    24812468        }
    24822469                       
    … …  
    25042491        bool hasNegCoeff = FALSE;       //or a negative one?
    25052492        size_t found;                   //used for find_first_(not)_of
    2506                
    2507 //      char prefix[]="/tmp/cone";
     2493
    25082494        string prefix="/tmp/cone";
    2509 //      const char *UCNstring = UCNstr.c_str();
    2510         //strcpy(UCNstring,UCNstr.c_str());
    2511        
    2512 //      char suffix[]=".sg";
    25132495        string suffix=".sg";
    2514 //      char *filename=strcat(prefix,UCNstring);
    25152496        string filename;
    25162497        filename.append(prefix);
    25172498        filename.append(UCNstr);
    25182499        filename.append(suffix);
    2519 //      strcat(filename,suffix);
    25202500                                       
    25212501        ifstream gcInputFile(filename.c_str(), ifstream::in);
    … …  
    25392519                        getline(gcInputFile, line);
    25402520                        strGcBasisLength = line;
    2541 //                       >> gcBasisLength;
    25422521                        gcBasisLength=atoi(strGcBasisLength.c_str());
    25432522                }
Note: See TracChangeset for help on using the changeset viewer.