source: git/kernel/gfan.cc @ aba53c

spielwiese
Last change on this file since aba53c was aba53c, checked in by Martin Monerjan, 15 years ago
Started to get normals from ddMatrix into intvec Esoteric things with pointers git-svn-id: file:///usr/local/Singular/svn/trunk@11586 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.9 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
4$Date: 2009-03-24 17:37:40 $
5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.19 2009-03-24 17:37:40 monerjan Exp $
6$Id: gfan.cc,v 1.19 2009-03-24 17:37:40 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 HOME
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
42class facet
43{
44        private:
45                intvec fNormal;         //inner normal, describing the facet uniquely
46        public:
47                facet()                 //default constructor. Do I need a constructor of type facet(intvec) ?
48                {
49                        //fNormal = FN;
50                        this->next=NULL;        //By default each facet is last
51                }
52               
53                ~facet(){;}             //destructor
54                void setFacetNormal(intvec iv){
55                        fNormal = iv;
56                        return;
57                }
58               
59                bool isFlippable;       //flippable facet? Want to have cone->isflippable.facet[i]
60                bool isIncoming;        //Is the facet incoming or outgoing?
61                facet *next;            //Pointer to next facet
62};
63
64/*class gcone
65finally this should become s.th. like gconelib.{h,cc} to provide an API
66*/
67class gcone
68{
69        private:
70                int numFacets;          //#of facets of the cone
71
72        public:
73                gcone(int);             //constructor with dimension
74                ~gcone();               //destructor
75                poly gcMarkedTerm;      //marked terms of the cone's Groebner basis
76                ideal gcBasis;          //GB of the cone
77                gcone *next;            //Pointer to *previous* cone in search tree
78       
79                void flip();            //Compute "the other side"
80                void remRedFacets();    //Remove redundant facets of the cone NOT NEEDED since this is already done by cddlib while compunting the normals
81               
82                ideal GenGrbWlk(ideal, ideal);  //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
83               
84
85};//class gcone
86
87ideal getGB(ideal inputIdeal)
88{
89        #ifdef gfan_DEBUG
90        cout << "Computing a groebner basis..." << endl;
91        #endif
92
93        ideal gb;
94        gb=kStd(inputIdeal,NULL,testHomog,NULL); //Possible to call without testHomog/isHomog?
95        idSkipZeroes(gb); //Get rid of zero entries
96
97        return gb;
98}
99
100/****** getConeNormals computes the inequalities ***/
101/*INPUT_TYPE: ideal                             */
102/*RETURN_TYPE: pointer to first facet           */
103/************************************************/
104facet getConeNormals(ideal I)
105{
106        #ifdef gfan_DEBUG
107        cout << "*** Computing Inequalities... ***" << endl;
108        #endif
109
110        //All variables go here - except ineq matrix and *v, see below
111        int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
112        int pCompCount;                 // # of terms in a poly
113        poly aktpoly;
114        int numvar = pVariables;        // # of variables in a polynomial (or ring?)
115        int leadexp[numvar];            // dirty hack of exp.vects
116        int aktexp[numvar];
117        int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
118        dd_rowrange ddrows;
119        dd_colrange ddcols;
120        dd_rowset ddredrows;            // # of redundant rows in ddineq
121        dd_rowset ddlinset;             // the opposite
122        dd_rowindex ddnewpos;           // all to make dd_Canonicalize happy
123        dd_NumberType ddnumb=dd_Real;   //Number type
124        dd_ErrorType dderr=dd_NoError;  //
125        // End of var declaration
126
127        cout << "The Groebner basis has " << lengthGB << " elements" << endl;
128        cout << "The current ring has " << numvar << " variables" << endl;
129        cols = numvar;
130
131        //Compute the # inequalities i.e. rows of the matrix
132        rows=0; //Initialization
133        for (int ii=0;ii<IDELEMS(I);ii++)
134        {
135                aktpoly=(poly)I->m[ii];
136                rows=rows+pLength(aktpoly)-1;
137        }
138        cout << "rows=" << rows << endl;
139        cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
140
141        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
142        dd_set_global_constants();
143        ddrows=rows;
144        ddcols=cols;
145        dd_MatrixPtr ddineq;            //Matrix to store the inequalities
146        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
147
148        // We loop through each g\in GB and compute the resulting inequalities
149        for (int i=0; i<IDELEMS(I); i++)
150        {
151                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
152                pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
153                cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
154
155                int *v=(int *)omAlloc((numvar+1)*sizeof(int));
156                pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
157               
158                //Store leadexp for aktpoly
159                for (int kk=0;kk<numvar;kk++)
160                {
161                        leadexp[kk]=v[kk+1];
162                        //printf("Leadexpcomp=%i\n",leadexp[kk]);
163                        //Since we need to know the difference of leadexp with the other expvects we do nothing here
164                        //but compute the diff below
165                }
166
167               
168                while (pNext(aktpoly)!=NULL) //move to next term until NULL
169                {
170                        aktpoly=pNext(aktpoly);
171                        pSetm(aktpoly);         //doesn't seem to help anything
172                        pGetExpV(aktpoly,v);
173                        for (int kk=0;kk<numvar;kk++)
174                        {
175//The ordering somehow gets lost here but this is acceptable, since we are only interested in the inequalities
176                                aktexp[kk]=v[kk+1];
177                                //printf("aktexpcomp=%i\n",aktexp[kk]);
178                                //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
179                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
180                        }
181                        aktmatrixrow=aktmatrixrow+1;
182                }//while
183
184        } //for
185
186        //Maybe add another row to contain the constraints of the standard simplex?
187
188        #ifdef gfan_DEBUG
189        cout << "The inequality matrix is" << endl;
190        dd_WriteMatrix(stdout, ddineq);
191        #endif
192
193        // The inequalities are now stored in ddineq
194        // Next we check for superflous rows
195        ddredrows = dd_RedundantRows(ddineq, &dderr);
196        if (dderr!=dd_NoError)                  // did an error occur?
197        {
198                 dd_WriteErrorMessages(stderr,dderr);   //if so tell us
199        } else
200        {
201                cout << "Redundant rows: ";
202                set_fwrite(stdout, ddredrows);          //otherwise print the redundant rows
203        }//if dd_Error
204
205        //Remove reduntant rows here!
206        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
207        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
208        ddcols = ddineq->colsize;
209        #ifdef gfan_DEBUG
210        cout << "Having removed redundancies, the normals now read:" << endl;
211        dd_WriteMatrix(stdout,ddineq);
212        cout << "Rows = " << ddrows << endl;
213        cout << "Cols = " << ddcols << endl;
214        #endif
215
216        /*dd_PolyhedraPtr ddpolyh;
217        dd_MatrixPtr G;
218        ddpolyh=dd_DDMatrix2Poly(ddineq, &dderr);
219        G=dd_CopyGenerators(ddpolyh);
220        printf("\nSpanning vectors = rows:\n");
221        dd_WriteMatrix(stdout, G);
222        */
223       
224       
225        /*Write the normals into class facet
226                How do I get the #rows in ddineq? \exists s.th. like dd_NumRows? dd_get_si(ddineq->matrix[i][j]) ?
227        Strange enough: facet *f=new facet(intvec(2,3)) does work for the constructor facet(intvec FN){fNormal = FN;}
228        */
229        #ifdef gfan_DEBUG
230        cout << "Creating list of normals" << endl;
231        #endif
232        /*The pointer *fRoot should be the return value of this function*/
233        facet *fRoot = new facet();             //instantiate new facet with intvec with numvar rows, one column and initial values all 0
234        //facet *fRoot;
235        facet *fAct;                    //instantiate pointer to active facet
236        fAct = (facet*)&fRoot; 
237        //fAct = fRoot;                 //Let fAct point to fRoot
238        for (int kk = 0; kk<ddrows; kk++)
239        {
240                intvec load;    //intvec to sto
241                for (int jj = 1; jj <ddcols; jj++)
242                {
243                        double *foo;
244                        foo = (double*)ddineq->matrix[kk][jj];  //get entry from actual position
245#ifdef gfan_DEBUG
246                        std::cout << "fAct is " << *foo << " at " << fAct << endl;
247#endif
248                        /*next two lines commented out. How to store values into intvec? */
249                        //load[jj] = (int)*foo;                 //store typecasted entry at pos jj of load
250                        //fAct->setFacetNormal(load);
251                        if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
252                        {
253                                fAct->next = new facet();       //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor
254                                fAct = fAct->next;              //scary :)
255                        }
256                }
257        }
258        /*
259        Now we should have a concatenated list containing the facet normals which we can adress via *f
260        */
261       
262        //ddineq->representation=dd_Inequality;         //We want our LP to be Ax>=0
263        //Clean up but don't delete the return value! (Whatever it will turn out to be)
264        dd_FreeMatrix(ddineq);
265        //set_free(ddrows);
266        //set_free(ddcols);
267        set_free(ddredrows);
268        //set_free(ddnumb);
269        //set_free(dderr);
270        free(ddnewpos);
271        //set_free(ddlinset);
272        dd_free_global_constants();
273
274
275        return *fRoot;
276}
277
278ideal gfan(ideal inputIdeal)
279{
280        #ifdef gfan_DEBUG
281        cout << "Now in subroutine gfan" << endl;
282        #endif
283        ideal res;
284        matrix ineq; //Matrix containing the boundary inequalities
285        /*
286        1. Select target order
287        2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
288        3. getConeNormals
289        */
290        res=getGB(inputIdeal);
291        getConeNormals(res);
292        return res;
293}
294#endif
Note: See TracBrowser for help on using the repository browser.