source: git/kernel/gfan.cc @ 19567b

spielwiese
Last change on this file since 19567b was 90adba8, checked in by Martin Monerjan, 15 years ago
extra.cc: prepared for LIST_CMD gfan.cc: removed functions getGB and *getConeNormals. These are now methods of class gcone. Construction of rings works partially. When entering a ring with ordering lp the program hangs after rWrite(rootRing) line 347. Prepared for return type LIST_CMD git-svn-id: file:///usr/local/Singular/svn/trunk@11601 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 11.7 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
4$Date: 2009-03-30 14:07:21 $
5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.24 2009-03-30 14:07:21 monerjan Exp $
6$Id: gfan.cc,v 1.24 2009-03-30 14:07:21 monerjan Exp $
7*/
8
9#include "mod2.h"
10
11#ifdef HAVE_GFAN
12
13#include "kstd1.h"
14#include "intvec.h"
15#include "polys.h"
16#include "ideals.h"
17#include "kmatrix.h"
18#include "fast_maps.h"  //Mapping of ideals
19#include "iostream.h"   //deprecated
20
21//Hacks for different working places
22#define ITWM
23
24#ifdef UNI
25#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
26#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
27#endif
28
29#ifdef HOME
30#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
31#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
32#endif
33
34#ifdef ITWM
35#include "/u/slg/monerjan/cddlib/include/setoper.h"
36#include "/u/slg/monerjan/cddlib/include/cdd.h"
37#endif
38
39#ifndef gfan_DEBUG
40#define gfan_DEBUG
41#endif
42
43//#include gcone.h
44
45/**
46*\brief Class facet
47*       Implements the facet structure as a linked list
48*
49*/
50class facet
51{
52        private:
53                /** inner normal, describing the facet uniquely up to isomorphism */
54                intvec *fNormal;               
55        public:
56                /** The default constructor. Do I need a constructor of type facet(intvec)? */
57                facet()                 
58                {
59                        // Pointer to next facet.  */
60                        /* Defaults to NULL. This way there is no need to check explicitly */
61                        this->next=NULL; 
62                }
63               
64                /** The default destructor */
65                ~facet(){;}
66               
67                /** Stores the facet normal \param intvec*/
68                void setFacetNormal(intvec *iv){
69                        fNormal = iv;
70                        //return;
71                }
72               
73                /** Method to print the facet normal*/
74                void printNormal()
75                {
76                        fNormal->show();
77                }
78               
79                bool isFlippable;       //flippable facet? Want to have cone->isflippable.facet[i]
80                bool isIncoming;        //Is the facet incoming or outgoing?
81                facet *next;            //Pointer to next facet
82};
83
84/**
85*\brief Implements the cone structure
86*
87* A cone is represented by a linked list of facet normals
88* @see facet
89*/
90/*class gcone
91finally this should become s.th. like gconelib.{h,cc} to provide an API
92*/
93class gcone
94{
95        private:
96                int numFacets;          //#of facets of the cone
97
98        public:
99                /** \brief Default constructor.
100                *
101                * Initialises this->next=NULL and this->facetPtr=NULL
102                */
103                gcone()
104                {
105                        this->next=NULL;
106                        this->facetPtr=NULL;
107                }
108                ~gcone();               //destructor
109                /** Pointer to the first facet */
110                facet *facetPtr;        //Will hold the adress of the first facet
111                poly gcMarkedTerm;      //marked terms of the cone's Groebner basis
112                ideal gcBasis;          //GB of the cone
113                gcone *next;            //Pointer to *previous* cone in search tree
114                /** \brief Compute the normals of the cone
115                *
116                * This method computes a representation of the cone in terms of facet normals. It takes an ideal
117                * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
118                * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
119                * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
120                * each row to represent an inequality of type const+x1+...+xn <= 0
121                * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
122                */
123                void getConeNormals(ideal I)
124                {
125#ifdef gfan_DEBUG
126                        cout << "*** Computing Inequalities... ***" << endl;
127#endif         
128                        //All variables go here - except ineq matrix and *v, see below
129                        int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
130                        int pCompCount;                 // # of terms in a poly
131                        poly aktpoly;
132                        int numvar = pVariables;        // # of variables in a polynomial (or ring?)
133                        int leadexp[numvar];            // dirty hack of exp.vects
134                        int aktexp[numvar];
135                        int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
136                        dd_rowrange ddrows;
137                        dd_colrange ddcols;
138                        dd_rowset ddredrows;            // # of redundant rows in ddineq
139                        dd_rowset ddlinset;             // the opposite
140                        dd_rowindex ddnewpos;           // all to make dd_Canonicalize happy
141                        dd_NumberType ddnumb=dd_Real;   //Number type
142                        dd_ErrorType dderr=dd_NoError;  //
143                        // End of var declaration
144#ifdef gfan_DEBUG
145                        cout << "The Groebner basis has " << lengthGB << " elements" << endl;
146                        cout << "The current ring has " << numvar << " variables" << endl;
147#endif
148                        cols = numvar;
149               
150                        //Compute the # inequalities i.e. rows of the matrix
151                        rows=0; //Initialization
152                        for (int ii=0;ii<IDELEMS(I);ii++)
153                        {
154                                aktpoly=(poly)I->m[ii];
155                                rows=rows+pLength(aktpoly)-1;
156                        }
157#ifdef gfan_DEBUG
158                        cout << "rows=" << rows << endl;
159                        cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
160#endif
161                        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
162                        dd_set_global_constants();
163                        ddrows=rows;
164                        ddcols=cols;
165                        dd_MatrixPtr ddineq;            //Matrix to store the inequalities
166                        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
167               
168                        // We loop through each g\in GB and compute the resulting inequalities
169                        for (int i=0; i<IDELEMS(I); i++)
170                        {
171                                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
172                                pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
173                                cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
174               
175                                int *v=(int *)omAlloc((numvar+1)*sizeof(int));
176                                pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
177                               
178                                //Store leadexp for aktpoly
179                                for (int kk=0;kk<numvar;kk++)
180                                {
181                                        leadexp[kk]=v[kk+1];
182                                        //Since we need to know the difference of leadexp with the other expvects we do nothing here
183                                        //but compute the diff below
184                                }
185               
186                               
187                                while (pNext(aktpoly)!=NULL) //move to next term until NULL
188                                {
189                                        aktpoly=pNext(aktpoly);
190                                        pSetm(aktpoly);         //doesn't seem to help anything
191                                        pGetExpV(aktpoly,v);
192                                        for (int kk=0;kk<numvar;kk++)
193                                        {
194                                                aktexp[kk]=v[kk+1];
195                                                //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
196                                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
197                                        }
198                                        aktmatrixrow=aktmatrixrow+1;
199                                }//while
200               
201                        } //for
202               
203                        //Maybe add another row to contain the constraints of the standard simplex?
204
205#ifdef gfan_DEBUG
206                        cout << "The inequality matrix is" << endl;
207                        dd_WriteMatrix(stdout, ddineq);
208#endif
209
210                        // The inequalities are now stored in ddineq
211                        // Next we check for superflous rows
212                        ddredrows = dd_RedundantRows(ddineq, &dderr);
213                        if (dderr!=dd_NoError)                  // did an error occur?
214                        {
215                                dd_WriteErrorMessages(stderr,dderr);    //if so tell us
216                        } else
217                        {
218                                cout << "Redundant rows: ";
219                                set_fwrite(stdout, ddredrows);          //otherwise print the redundant rows
220                        }//if dd_Error
221               
222                        //Remove reduntant rows here!
223                        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
224                        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
225                        ddcols = ddineq->colsize;
226#ifdef gfan_DEBUG
227                        cout << "Having removed redundancies, the normals now read:" << endl;
228                        dd_WriteMatrix(stdout,ddineq);
229                        cout << "Rows = " << ddrows << endl;
230                        cout << "Cols = " << ddcols << endl;
231#endif
232                        /*Write the normals into class facet*/
233#ifdef gfan_DEBUG
234                        cout << "Creating list of normals" << endl;
235#endif
236                        /*The pointer *fRoot should be the return value of this function*/
237                        facet *fRoot = new facet();             //instantiate new facet with intvec with numvar rows, one column and initial values all 0
238                        facetPtr = fRoot;               //set variable facetPtr of class gcone to first facet
239                        facet *fAct;                    //instantiate pointer to active facet
240                        fAct = fRoot;           //This does not seem to do the trick. fRoot and fAct have to point to the same adress!
241                        std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
242                        for (int kk = 0; kk<ddrows; kk++)
243                        {
244                                intvec *load = new intvec(numvar);      //intvec to store a single facet normal that will then be stored via setFacetNormal
245                                for (int jj = 1; jj <ddcols; jj++)
246                                {
247                                        double *foo;
248                                        foo = (double*)ddineq->matrix[kk][jj];  //get entry from actual position
249#ifdef gfan_DEBUG
250                                        std::cout << "fAct is " << *foo << " at " << fAct << std::endl;
251#endif 
252                                        (*load)[jj-1] = (int)*foo;              //store typecasted entry at pos jj-1 of load
253                                        //check for flipability here
254                                        if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
255                                        {
256                                                fAct->next = new facet();       //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor
257                                                fAct = fAct->next;              //scary :)
258                                        }
259                                }//for jj<ddcols
260                                /*Now load should be full and we can call setFacetNormal*/
261                                fAct->setFacetNormal(load);
262                                fAct->printNormal();
263                        }
264                        /*
265                        Now we should have a linked list containing the facet normals of those facets that are
266                        -irredundant
267                        -flipable
268                        Adressing is done via *facetPtr
269                        */
270                       
271                        //Clean up but don't delete the return value! (Whatever it will turn out to be)
272                        dd_FreeMatrix(ddineq);
273                        set_free(ddredrows);
274                        free(ddnewpos);
275                        set_free(ddlinset);
276                        dd_free_global_constants();
277
278                }//method getConeNormals(ideal I)       
279               
280                void flip();            //Compute "the other side"
281                               
282                /** \brief Compute a Groebner Basis
283                *
284                * Computes the Groebner basis and stores the result in this->gcBasis
285                *\param ideal
286                *\return void
287                */
288                void getGB(ideal inputIdeal)
289                {
290                        ideal gb;
291                        gb=kStd(inputIdeal,NULL,testHomog,NULL);
292                        idSkipZeroes(gb);
293                        this->gcBasis=gb;       //write the GB into gcBasis
294                }
295               
296                ideal GenGrbWlk(ideal, ideal);  //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
297               
298
299};//class gcone
300
301/*
302function getGB incorporated into class gcone with rev 1.24
303*/
304
305//DEPRECATED since rev 1.24 with existence of gcone::getConeNormals(ideal I);
306//Kept for unknown reasons ;)
307facet *getConeNormals(ideal I)
308{
309        return NULL;
310}
311
312ideal gfan(ideal inputIdeal)
313{
314        int numvar = pVariables; 
315       
316        #ifdef gfan_DEBUG
317        cout << "Now in subroutine gfan" << endl;
318        #endif
319        ring inputRing=currRing;        // The ring the user entered
320        ring rootRing;                  // The ring associated to the target ordering
321        ideal res;
322        //matrix ineq;                  //Matrix containing the boundary inequalities
323        facet *fRoot;
324       
325        /*
326        1. Select target order, say dp.
327        2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
328        3. getConeNormals
329        */
330
331       
332        /* Construct a new ring which will serve as our root
333        Does not yet work as expected. Will work fine with order dp,Dp but otherwise hangs in getGB
334        */
335        rootRing=rCopy0(currRing);
336        rootRing->order[0]=ringorder_dp;
337        rComplete(rootRing);
338        rChangeCurrRing(rootRing);
339        ideal rootIdeal;
340        /* Fetch the inputIdeal into our rootRing */
341        map m=(map)idInit(IDELEMS(inputIdeal),0);
342        rootIdeal=fast_map(inputIdeal,inputRing,(ideal)m, currRing);
343#ifdef gfan_DEBUG
344        cout << "Root ideal is " << endl;
345        idPrint(rootIdeal);
346        cout << "The current ring is " << endl;
347        rWrite(rootRing);
348        cout << endl;
349#endif 
350       
351        gcone *gcRoot = new gcone();    //Instantiate the sink
352        gcone *gcAct;
353        gcAct = gcRoot;
354        gcAct->getGB(inputIdeal);
355        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
356        /*Now it is time to compute the search facets, respectively start the reverse search.
357        But since we are in the root all facets should be search facets.
358        MIND: AS OF NOW, THE LIST OF FACETS IS NOT PURGED OF NON-FLIPPAPLE FACETS
359        */
360       
361        /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
362        The return type will then be of type LIST_CMD
363        Assume gfan has finished, thus we have enumerated all the cones
364        Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
365        Groebner Basis and merge this somehow with LIST_CMD
366        => Count the cones!
367        */
368       
369        res=gcAct->gcBasis;
370        //cout << fRoot << endl;
371        return res;
372        //return GBlist;
373}
374/*
375Since gfan.cc is #included from extra.cc there must not be a int main(){} here
376*/
377#endif
Note: See TracBrowser for help on using the repository browser.