Changeset 472eeb in git for kernel/gfan.cc

Timestamp:

May 15, 2009, 3:54:46 PM (15 years ago)

Author:

Martin Monerjan

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

294312ffcc6c4d753fa40971861665f124f63851

Parents:

8f63e4c7ff0fb85b12623e48dd7beeb869ab2257

Message:

facets2Matrix


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

File:

• kernel/gfan.cc (modified) (15 diffs)

kernel/gfan.cc

@@ -2 +2 @@
 Compute the Groebner fan of an ideal
 $Author: monerjan $
-$Date: 2009-05-14 09:55:36 $
-$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.51 2009-05-14 09:55:36 monerjan Exp $
-$Id: gfan.cc,v 1.51 2009-05-14 09:55:36 monerjan Exp $
+$Date: 2009-05-15 13:54:46 $
+$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.52 2009-05-15 13:54:46 monerjan Exp $
+$Id: gfan.cc,v 1.52 2009-05-15 13:54:46 monerjan Exp $
 */
 
@@ -212 +212 @@
 class gcone
 {
-        private:
-                int numFacets;          //#of facets of the cone
+        private:                
                 ring rootRing;          //good to know this -> generic walk
                 ideal inputIdeal;       //the original
@@ -222 +221 @@
                 
         public: 
-                /** Pointer to the first facet */
+                /** \brief Pointer to the first facet */
                 facet *facetPtr;        //Will hold the adress of the first facet; set by gcone::getConeNormals
                 
                 /** # of variables in the ring */
                 int numVars;            //#of variables in the ring
+                
+                /** # of facets of the cone
+                * This value is set by gcone::getConeNormals
+                */
+                int numFacets;          //#of facets of the cone
                 
                 /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
@@ -276 +280 @@
                         
                         intvec *ivtmp = new intvec(this->numVars);                                              
-                        ivtmp = gc.facetPtr->getFacetNormal();
-                        //ivtmp->show();                        
+                        ivtmp = gc.facetPtr->getFacetNormal();                  
                         
                         ideal gb;
                         gb=gc.facetPtr->getFlipGB();                    
-                        this->gcBasis=gb;//gc.facetPtr->getFlipGB();    //this cone's GB is the flipped GB
-                        //idShow(this->gcBasis);
+                        this->gcBasis=gb;       //this cone's GB is the flipped GB                      
                         
                         /*Reverse direction of the facet normal to make it an inner normal*/                    
@@ -394 +396 @@
                                 aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
                                 pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
+#ifdef gfan_DEBUG
                                 cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
+#endif
                 
                                 int *v=(int *)omAlloc((numvar+1)*sizeof(int));
@@ -469 +473 @@
                                 for (int jj = 1; jj <ddcols; jj++)
                                 {
-#ifdef GMPRATIONAL
                                         double foo;
                                         foo = mpq_get_d(ddineq->matrix[kk][jj]);
@@ -475 +478 @@
                                         std::cout << "fAct is " << foo << " at " << fAct << std::endl;
 #endif
-                                        (*load)[jj-1] = (int)foo;               //store typecasted entry at pos jj-1 of load
-#endif                  
-
-                                        //(*load)[jj-1] = (int)foo;             //store typecasted entry at pos jj-1 of load
-                                        //check for flipability here
-//                                      if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
-//                                      {
-                                                //fAct->next = new facet();     //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                                
-//                                      }
+                                        (*load)[jj-1] = (int)foo;               //store typecasted entry at pos jj-1 of load            
                                 }//for (int jj = 1; jj <ddcols; jj++)
+                                
                                 /*Quick'n'dirty hack for flippability*/ 
                                 bool isFlippable=FALSE;                         
@@ -512 +508 @@
                                 delete load;
                         }//for (int kk = 0; kk<ddrows; kk++)
+                        
                         /*
                         Now we should have a linked list containing the facet normals of those facets that are
@@ -526 +523 @@
                                 delete iv;
                         }
+                        
+                        //Compute the number of facets
+                        this->numFacets=ddineq->rowsize;
                         
                         //Clean up but don't delete the return value! (Whatever it will turn out to be)                 
@@ -1142 +1142 @@
                 }
                 
-                /** 
+                /** \brief Checks whether a given facet is a search facet
                 * Determines whether a given facet of a cone is the search facet of a neighbouring cone
                 * This is done in the following way:
@@ -1319 +1319 @@
                 void noRevS(gcone &gc, bool usingIntPoint=FALSE)
                 {
-                        facet *fListPtr = new facet();                  
+                        facet *fListPtr = new facet();  //The list containing ALL facets we come across                 
                         facet *fAct = new facet();
                         fAct = gc.facetPtr;
@@ -1356 +1356 @@
                                         fAct=fAct->next;        //iterate
                                 }       
-                        }               
+                        }
+                        else/*Facets of facets*/
+                        {
+                                dd_MatrixPtr ddineq,P,ddakt;
+                                dd_rowset impl_linset, redset;
+                                dd_ErrorType err;
+                                dd_rowindex newpos;
+                                
+                                ddineq = facets2Matrix(gc);
+                                
+                                /*Now set appropriate linearity*/
+                                dd_PolyhedraPtr ddpolyh;
+                                for (int ii=0; ii<ddineq->rowsize; ii++)
+                                {       
+                                        cout << "------------" << endl;
+                                        ddakt = dd_CopyMatrix(ddineq);
+                                        set_addelem(ddakt->linset,ii+1);
+//                                      dd_WriteMatrix(stdout,ddakt);
+                                                                                
+                                        dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);                    
+                                        
+//                                      dd_WriteMatrix(stdout,ddakt);
+                                        ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
+                                        P=dd_CopyGenerators(ddpolyh);
+                                        dd_WriteMatrix(stdout,P);
+                                        
+                                        intvec *load = new intvec(this->numVars);
+                                        
+                                        dd_FreeMatrix(ddakt);
+                                        dd_FreePolyhedra(ddpolyh);
+                                }
+                                
+                        }                               
                         
                         //NOTE Hm, comment in and get a crash for free...
@@ -1367 +1399 @@
                         fAct = fListPtr;
                         gcone *gcTmp = new gcone(gc);   //copy
-                        gcTmp->flip(gcTmp->gcBasis,fAct);
+                        //gcTmp->flip(gcTmp->gcBasis,fAct);
                         //NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
                         
                 }//void noRevS(gcone &gc)
+                
+                /** \brief Construct a ddMatrix from a cone's list of facets
+                * 
+                */
+                dd_MatrixPtr facets2Matrix(gcone const &gc)
+                {
+                        facet *fAct = new facet();
+                        fAct = gc.facetPtr;                     
+                        dd_MatrixPtr M;
+                        dd_rowrange ddrows;
+                        dd_colrange ddcols;
+                        ddcols=(this->numVars)+1;
+                        ddrows=this->numFacets;
+                        M=dd_CreateMatrix(ddrows,ddcols);
+                        
+                        //dd_set_global_constants();                                            
+                                                
+                        intvec *comp = new intvec(this->numVars);                       
+                        
+                        for (int jj=0; jj<M->rowsize; jj++)
+                        {
+                                comp = fAct->getFacetNormal();
+                                for (int ii=0; ii<this->numVars; ii++)
+                                {
+                                        dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
+                                }
+                                fAct = fAct->next;                              
+                        }//jj
+                        delete comp;            
+                        
+                        return M;
+                }
                 
                 /** \brief Write information about a cone into a file on disk
@@ -1390 +1454 @@
                         }
                         else
-                        {                       
+                        {       /*gcOutputFile << "UCN" << endl;
+                                gcOutputFile << this->UCN << endl;*/
                                 gcOutputFile << "RING" << endl;                         
                                 //Write this->gcBasis into file
@@ -1402 +1467 @@
                                 fAct = gc.facetPtr;
                                 gcOutputFile << "FACETS" << endl;                                                               
-                                while(fAct->next!=NULL)
+                                while(fAct!=NULL)
                                 {       
                                         intvec *iv = new intvec(gc.numVars);