Changeset 5684df in git


Ignore:
Timestamp:
May 4, 2009, 5:14:38 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
b2a6e3e0d28bf571be753914eaccf9f5c629080f
Parents:
e2efbe929b401e118eadf323f44b5bd894b457b0
Message:
found bug in flippability test


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    re2efbe9 r5684df  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-04-30 09:32:01 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.43 2009-04-30 09:32:01 monerjan Exp $
    6 $Id: gfan.cc,v 1.43 2009-04-30 09:32:01 monerjan Exp $
     4$Date: 2009-05-04 15:14:38 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.44 2009-05-04 15:14:38 monerjan Exp $
     6$Id: gfan.cc,v 1.44 2009-05-04 15:14:38 monerjan Exp $
    77*/
    88
    … …  
    210210               
    211211                /** \brief Default destructor */
    212                 ~gcone();               //destructor
     212                ~gcone(){;}             //destructor
    213213               
    214214                /** Pointer to the first facet */
    … …  
    403403                                }//for (int jj = 1; jj <ddcols; jj++)
    404404                                /*Quick'n'dirty hack for flippability*/
    405                                 bool isFlippable;                       
     405                                bool isFlippable;
     406                                //NOTE BUG HERE
     407                                /* The flippability check is wrong!
     408                                (1,-4) will pass, but (-1,7) will not.
     409                                REWRITE THIS
     410                                */                     
    406411                                for (int jj = 0; jj<this->numVars; jj++)
    407412                                {                                       
    … …  
    514519                                        cout << endl;
    515520#endif
    516                                         if (isParallel(check,fNormal)) //pass *check when
     521                                        //TODO why not *check, *fNormal????
     522                                        if (isParallel(*check,*fNormal)) //pass *check when
    517523                                        {
    518524                                                cout << "Parallel vector found, adding to initialFormElement" << endl;                 
    … …  
    10221028                /**
    10231029                * Determines whether a given facet of a cone is the search facet of a neighbouring cone
    1024                 *
     1030                * This is done in the following way:
     1031                * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
     1032                * that is first crossed during the generic walk.
     1033                * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
     1034                * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
    10251035                */
    1026                 bool isSearchFacet(gcone &tmpcone, facet *testfacet)
    1027                 {                       
     1036                bool isSearchFacet(gcone &tmpcone, facet &testfacet)
     1037                {                              
    10281038                        ring actRing=currRing;
    1029                        
    1030                         facet *fMin=new facet();        //Pointer to the "minimal" facet
     1039                        facet *facetPtr=(facet*)tmpcone.facetPtr;                       
     1040                        facet *fMin=new facet(*facetPtr);       //Pointer to the "minimal" facet
     1041                        //facet *fMin = new facet(tmpcone.facetPtr);
     1042                        //fMin=tmpcone.facetPtr;                //Initialise to first facet of tmpcone
    10311043                        facet *fAct;
    1032                         fMin = testfacet;
     1044                        //fMin = testfacet;
     1045                        //fMin->next=testfacet->next;
    10331046                        fAct = fMin;
    10341047                       
    1035                         rChangeCurrRing(this->rootRing);
     1048                        cout << endl<< fMin->next << endl;
     1049                        cout << tmpcone.facetPtr->next << endl;
     1050                        tmpcone.facetPtr->printNormal();
     1051                        fAct=tmpcone.facetPtr->next;
     1052                        fAct->printNormal();
     1053                        /*fMin->printNormal();
     1054                        fMin->next->printNormal();*/
     1055                       
     1056                        rChangeCurrRing(this->rootRing);        //because we compare the monomials in the rootring                     
    10361057                        poly p=NULL;
    10371058                        poly q=NULL;
    … …  
    10401061                        intvec *sigma = new intvec(this->numVars);
    10411062                        sigma=tmpcone.getIntPoint();
     1063                        cout << "pling" << endl;
    10421064                        int *u=(int *)omAlloc((this->numVars+1)*sizeof(int));
    10431065                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    10441066                        int weight1,weight2;
    1045                         while(tmpcone.facetPtr->next!=NULL)
     1067                        while(fAct->next!=NULL) //ersetzen durch fAct
    10461068                        {
    10471069                                /* Get alpha_i and alpha_{i+1} */
    10481070                                p_weight=fMin->getFacetNormal();
    10491071                                q_weight=fMin->next->getFacetNormal();
    1050                                
     1072                                p_weight->show();
     1073                                q_weight->show();
    10511074                                /*Compute the dot product of sigma and alpha_{i,j}*/
    10521075                                weight1=dotProduct(sigma,p_weight);
    10531076                                weight2=dotProduct(sigma,q_weight);
    1054                                
     1077                                cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
    10551078                                /*Adjust alpha_i and alpha_i+1 accordingly*/
    10561079                                for(int ii=1;ii<=this->numVars;ii++)
    … …  
    10611084                                        v[ii]=weight2*(*q_weight)[ii];
    10621085                                }
     1086                                cout << "PLONK" << endl;
    10631087                               
    10641088                                /*Now p_weight and q_weight need to be compared as exponent vectors*/
    10651089                                pSetExpV(p,u);
    10661090                                pSetExpV(q,v);
     1091                                cout << "AARGH" << endl;
    10671092                                /*We want to check whether x^p < x^q
    10681093                                => want to check for return value 1 */
    … …  
    10761101                                        fAct=fAct->next;
    10771102                                }
     1103                                fAct=fAct->next;
    10781104                        }//while(tmpcone.facetPtr->next!=NULL)
    1079                         rChangeCurrRing(actRing);
     1105                       
     1106                        /*If testfacet was minimal then fMin should still point there */
     1107                        //if (isParallel(fMin->getFacetNormal,testfacet.getFacetNormal))
     1108                        if (fMin==tmpcone.facetPtr)
     1109                        {
     1110                                rChangeCurrRing(actRing);
     1111                                return TRUE;
     1112                        }
     1113                        else
     1114                        {
     1115                                rChangeCurrRing(actRing);
     1116                                return FALSE;
     1117                        }
    10801118                }//bool isSearchFacet
    10811119               
    1082                 void reverseSearch(gcone gcAct) //no const possible here since we call gcAct->flip
     1120                void reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
    10831121                {
    10841122                        facet *fAct=new facet();
    1085                         fAct = gcAct.facetPtr;                 
    1086                        
    1087                         while(fAct->next!=NULL)
    1088                         {
    1089                                 gcAct.flip(gcAct.gcBasis,gcAct.facetPtr);
    1090                                 gcone *gcTmp = new gcone(gcAct);
     1123                        fAct = gcAct->facetPtr;                 
     1124                       
     1125                        while(fAct->next!=NULL)  //NOTE NOT SURE WHETHER THIS IS RIGHT! Do I reach EVERY facet or only all but the last?
     1126                        {
     1127                                cout << "==========================================================================================="<< endl;
     1128                                gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
     1129                                gcone *gcTmp = new gcone(*gcAct);
    10911130                                idShow(gcTmp->gcBasis);
    10921131                                gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
     1132#ifdef gfan_DEBUG
     1133                                facet *f = new facet();
     1134                                f=gcTmp->facetPtr;
     1135                                while(f->next!=NULL)
     1136                                {
     1137                                        f->printNormal();
     1138                                        f=f->next;                                     
     1139                                }
     1140#endif
    10931141                                gcTmp->showIntPoint();
    10941142                                /*recursive part goes gere*/
    1095                                 if (isSearchFacet(gcTmp,gcAct.facetPtr))
    1096                                 {
    1097                                         gcAct.next=gcTmp;
    1098                                         reverseSearch(*gcTmp);
     1143                                if (isSearchFacet(*gcTmp,(facet&)gcAct->facetPtr))
     1144                                {
     1145                                        gcAct->next=gcTmp;
     1146                                        cout << "PING"<< endl;
     1147                                        reverseSearch(gcTmp);
    10991148                                }
    11001149                                else
    … …  
    11131162        int numvar = pVariables;
    11141163       
    1115         #ifdef gfan_DEBUG
     1164#ifdef gfan_DEBUG
    11161165        cout << "Now in subroutine gfan" << endl;
    1117         #endif
     1166#endif
    11181167        ring inputRing=currRing;        // The ring the user entered
    11191168        ring rootRing;                  // The ring associated to the target ordering
    … …  
    11791228                fAct = fAct->next;             
    11801229        }*/
    1181         gcAct->reverseSearch(*gcAct);
     1230        gcAct->reverseSearch(gcAct);
    11821231        /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
    11831232        The return type will then be of type LIST_CMD
Note: See TracChangeset for help on using the changeset viewer.