Context Navigation

source: git/kernel/gfan.cc @ 2b7bb5d

fieker-DuValspielwiese

Last change on this file since 2b7bb5d was 2b7bb5d, checked in by Martin Monerjan, 15 years ago
Continued work on gcone::flip - HIGHLY EXPERIMENTAL CODE! Added public variable gcone::numVars to store #of variables Began work on bool gcone::isParallel git-svn-id: file:///usr/local/Singular/svn/trunk@11609 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 14.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: