source: git/kernel/gfan.cc @ fc4f9a

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