source: git/kernel/gfan.cc @ 53e33d9

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