Context Navigation

source: git/kernel/gfan.cc @ 5cd07b

spielwiese

Last change on this file since 5cd07b was 5cd07b, checked in by Martin Monerjan, 15 years ago
Two new methods: int gcone::dotProduct(intvec a, intvec b) bool gcone::isParallel(intvec a, intvec b) gcone::flip now computes the initial form as an array of polynomials. How to get this into an ideal? git-svn-id: file:///usr/local/Singular/svn/trunk@11615 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 16.6 KB

Line
1	/*
2	Compute the Groebner fan of an ideal
3	$Author: monerjan $
4	$Date: 2009-04-03 13:49:16 $
5	$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.28 2009-04-03 13:49:16 monerjan Exp $
6	$Id: gfan.cc,v 1.28 2009-04-03 13:49:16 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	this->fNormal = ivCopy(iv);
71	//return;
72	}
73
74	/** Hopefully returns the facet normal */
75	intvec *getFacetNormal()
76	{
77	//this->fNormal->show(); cout << endl;
78	return this->fNormal;
79	}
80
81	/** Method to print the facet normal*/
82	void printNormal()
83	{
84	fNormal->show();
85	}
86
87	/** \brief The Groebner basis on the other side of a shared facet
88	*
89	* In order not to have to compute the flipped GB twice we store the basis we already get
90	* when identifying search facets. Thus in the next step of the reverse search we can
91	* just copy the old cone and update the facet and the gcBasis
92	*/
93	ideal flibGB; //The Groebner Basis on the other side, computed via gcone::flip
94
95	bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i]
96	bool isIncoming; //Is the facet incoming or outgoing?
97	facet *next; //Pointer to next facet
98	};
99
100	/**
101	*\brief Implements the cone structure
102	*
103	* A cone is represented by a linked list of facet normals
104	* @see facet
105	*/
106	/*class gcone
107	finally this should become s.th. like gconelib.{h,cc} to provide an API
108	*/
109	class gcone
110	{
111	private:
112	int numFacets; //#of facets of the cone
113
114	public:
115	/** \brief Default constructor.
116	*
117	* Initialises this->next=NULL and this->facetPtr=NULL
118	*/
119	gcone()
120	{
121	this->next=NULL;
122	this->facetPtr=NULL;
123	}
124	~gcone(); //destructor
125
126	/** Pointer to the first facet */
127	facet *facetPtr; //Will hold the adress of the first facet; set by gcone::getConeNormals
128	poly gcMarkedTerm; //marked terms of the cone's Groebner basis
129	int numVars; //#of variables in the ring
130
131	/** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
132	ideal gcBasis; //GB of the cone, set by gcone::getGB();
133	gcone next; //Pointer to previous* cone in search tree
134
135	/** \brief Compute the normals of the cone
136	*
137	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
138	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
139	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
140	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
141	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
142	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
143	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
144	*/
145	void getConeNormals(ideal I)
146	{
147	#ifdef gfan_DEBUG
148	cout << "* Computing Inequalities... *" << endl;
149	#endif
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	#ifdef gfan_DEBUG
167	cout << "The Groebner basis has " << lengthGB << " elements" << endl;
168	cout << "The current ring has " << numvar << " variables" << endl;
169	#endif
170	cols = numvar;
171
172	//Compute the # inequalities i.e. rows of the matrix
173	rows=0; //Initialization
174	for (int ii=0;ii<IDELEMS(I);ii++)
175	{
176	aktpoly=(poly)I->m[ii];
177	rows=rows+pLength(aktpoly)-1;
178	}
179	#ifdef gfan_DEBUG
180	cout << "rows=" << rows << endl;
181	cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
182	#endif
183	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
184	dd_set_global_constants();
185	ddrows=rows;
186	ddcols=cols;
187	dd_MatrixPtr ddineq; //Matrix to store the inequalities
188	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
189
190	// We loop through each g\in GB and compute the resulting inequalities
191	for (int i=0; i<IDELEMS(I); i++)
192	{
193	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
194	pCompCount=pLength(aktpoly); //How many terms does aktpoly consist of?
195	cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
196
197	int v=(int )omAlloc((numvar+1)*sizeof(int));
198	pGetExpV(aktpoly,v); //find the exp.vect in v[1],...,v[n], use pNext(p)
199
200	//Store leadexp for aktpoly
201	for (int kk=0;kk<numvar;kk++)
202	{
203	leadexp[kk]=v[kk+1];
204	//Since we need to know the difference of leadexp with the other expvects we do nothing here
205	//but compute the diff below
206	}
207
208
209	while (pNext(aktpoly)!=NULL) //move to next term until NULL
210	{
211	aktpoly=pNext(aktpoly);
212	pSetm(aktpoly); //doesn't seem to help anything
213	pGetExpV(aktpoly,v);
214	for (int kk=0;kk<numvar;kk++)
215	{
216	aktexp[kk]=v[kk+1];
217	//ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk]; //dito
218	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
219	}
220	aktmatrixrow=aktmatrixrow+1;
221	}//while
222
223	} //for
224
225	//Maybe add another row to contain the constraints of the standard simplex?
226
227	#ifdef gfan_DEBUG
228	cout << "The inequality matrix is" << endl;
229	dd_WriteMatrix(stdout, ddineq);
230	#endif
231
232	// The inequalities are now stored in ddineq
233	// Next we check for superflous rows
234	ddredrows = dd_RedundantRows(ddineq, &dderr);
235	if (dderr!=dd_NoError) // did an error occur?
236	{
237	dd_WriteErrorMessages(stderr,dderr); //if so tell us
238	} else
239	{
240	cout << "Redundant rows: ";
241	set_fwrite(stdout, ddredrows); //otherwise print the redundant rows
242	}//if dd_Error
243
244	//Remove reduntant rows here!
245	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
246	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
247	ddcols = ddineq->colsize;
248	#ifdef gfan_DEBUG
249	cout << "Having removed redundancies, the normals now read:" << endl;
250	dd_WriteMatrix(stdout,ddineq);
251	cout << "Rows = " << ddrows << endl;
252	cout << "Cols = " << ddcols << endl;
253	#endif
254	/Write the normals into class facet/
255	#ifdef gfan_DEBUG
256	cout << "Creating list of normals" << endl;
257	#endif
258	/The pointer fRoot should be the return value of this function*/
259	facet *fRoot = new facet(); //instantiate new facet
260	this->facetPtr = fRoot; //set variable facetPtr of class gcone to first facet
261	facet *fAct; //instantiate pointer to active facet
262	fAct = fRoot; //This does not seem to do the trick. fRoot and fAct have to point to the same adress!
263	std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
264	for (int kk = 0; kk<ddrows; kk++)
265	{
266	intvec *load = new intvec(numvar); //intvec to store a single facet normal that will then be stored via setFacetNormal
267	for (int jj = 1; jj <ddcols; jj++)
268	{
269	double *foo;
270	foo = (double*)ddineq->matrix[kk][jj]; //get entry from actual position
271	#ifdef gfan_DEBUG
272	std::cout << "fAct is " << *foo << " at " << fAct << std::endl;
273	#endif
274	(load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
275	//check for flipability here
276	if (jj<ddcols) //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
277	{
278	fAct->next = new facet(); //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor
279	}
280	}//for jj<ddcols
281	/Now load should be full and we can call setFacetNormal/
282	fAct->setFacetNormal(load);
283	//fAct->printNormal();
284	fAct=fAct->next; //this should definitely not be called in the above while-loop :D
285	delete load;
286	}
287	/*
288	Now we should have a linked list containing the facet normals of those facets that are
289	-irredundant
290	-flipable
291	Adressing is done via *facetPtr
292	*/
293
294	//Clean up but don't delete the return value! (Whatever it will turn out to be)
295	dd_FreeMatrix(ddineq);
296	set_free(ddredrows);
297	free(ddnewpos);
298	set_free(ddlinset);
299	dd_free_global_constants();
300
301	}//method getConeNormals(ideal I)
302
303	/*bool isParallel(int v[], intvec iv)
304	{
305	}*/
306
307	/** \brief Compute the Groebner Basis on the other side of a shared facet
308	*
309	* Implements algorithm 4.3.2 from Anders' thesis.
310	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
311	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
312	* is parallel to \f$ leadexp(g) \f$
313	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
314	* Other possibilities includes computing the rank of the matrix consisting of the vectors in question and
315	* computing an interior point of the facet and taking all terms having the same weight with respect
316	* to this interior point.
317	*\param ideal, facet
318	* Input: a marked,reduced Groebner basis and a facet
319	*/
320	void flip(ideal gb, facet *f) //Compute "the other side"
321	{
322
323	intvec *fNormal = new intvec(this->numVars); //facet normal, check for parallelity
324	fNormal = f->getFacetNormal(); //read this->fNormal;
325	#ifdef gfan_DEBUG
326	cout << "fNormal=";
327	fNormal->show();
328	cout << endl;
329	#endif
330	/1st step: Compute the initial ideal/
331	poly initialFormElement[IDELEMS(gb)]; //array of #polys in GB to store initial form
332	ideal initialForm;
333	poly aktpoly;
334	intvec *check = new intvec(this->numVars); //array to store the difference of LE and v
335
336	for (int ii=0;ii<IDELEMS(gb);ii++)
337	{
338	aktpoly = (poly)gb->m[ii];
339	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
340	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
341	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
342	//pGetExpV(aktpoly,ivLeadExpV);
343	initialFormElement[ii]=pHead(aktpoly);
344
345	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
346	{
347	aktpoly=pNext(aktpoly); //next term
348	pSetm(aktpoly);
349	pGetExpV(aktpoly,v);
350	/* Convert (int)v into (intvec)check */
351	for (int jj=0;jj<this->numVars;jj++)
352	{
353	//cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
354	//cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
355	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
356	}
357	cout << "check=";
358	check->show();
359	cout << endl;
360	if (isParallel(check,fNormal)) //pass *check when
361	{
362	cout << "Parallel vector found, adding to initialFormElement" << endl;
363	initialFormElement[ii] = pAdd(initialFormElement[ii],(poly)pHead(aktpoly));
364	}
365	}//while
366	cout << "Initial Form=";
367	pWrite(initialFormElement[ii]);
368	cout << "---" << endl;
369	/Now initialFormElement must be added to (ideal)initialForm /
370	}//for
371	}//void flip(ideal gb, facet *f)
372
373	/** \brief Compute a Groebner Basis
374	*
375	* Computes the Groebner basis and stores the result in gcone::gcBasis
376	*\param ideal
377	*\return void
378	*/
379	void getGB(ideal inputIdeal)
380	{
381	ideal gb;
382	gb=kStd(inputIdeal,NULL,testHomog,NULL);
383	idSkipZeroes(gb);
384	this->gcBasis=gb; //write the GB into gcBasis
385	}
386
387	ideal GenGrbWlk(ideal, ideal); //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
388
389
390	/** \brief Compute the dot product of two intvecs
391	*
392	* THIS IS WEIRD - BUT WORKS
393	*/
394	int dotProduct(intvec a, intvec b)
395	{
396	intvec iva=*a;
397	intvec ivb=*b;
398	int res=0;
399	for (int i=0;i<this->numVars;i++)
400	{
401	res = res+(iva[i]*ivb[i]);
402	}
403	return res;
404	}
405
406	/** \brief Check whether two intvecs are parallel
407	*
408	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
409	*/
410	bool isParallel(intvec a, intvec b)
411	{
412	//bool res=FALSE;
413	//intvec iva(this->numVars);
414	//intvec ivb(this->numVars);
415	int lhs,rhs;
416	//lhs=0;
417	//rhs=0;
418	//iva = a;
419	//ivb = b;
420	//lhs=dotProduct(&iva,&ivb)*dotProduct(&iva,&ivb);
421	//lhs=dotProduct(&iva,&iva)*dotProduct(&ivb,&ivb);
422	lhs=dotProduct(&a,&b)*dotProduct(&a,&b);
423	rhs=dotProduct(&a,&a)*dotProduct(&b,&b);
424	cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
425	if (lhs==rhs)
426	{
427	return TRUE;
428	}
429	else
430	{
431	return FALSE;
432	}
433	}//bool isParallel
434
435	/** \brief Check two intvecs for equality --- PROBABLY NOT NEEDED
436	*
437	* \f$ \alpha=\beta\Leftrightarrow\langle\alpha-\beta,\alpha-\beta\rangle=0 \f$
438	*/
439	bool isEqual(intvec a, intvec b)
440	{
441	intvec *ivdiff;
442	int res;
443
444	ivdiff=ivSub(&a,&b);
445	cout << "ivdiff=";
446	ivdiff->show();
447	cout << endl;
448	//res=dotProduct(ivdiff,ivdiff);
449	if (res==0)
450	{
451	return TRUE;
452	}
453	else
454	{
455	return FALSE;
456	}
457	}//bool isEqual
458	};//class gcone
459
460	ideal gfan(ideal inputIdeal)
461	{
462	int numvar = pVariables;
463
464	#ifdef gfan_DEBUG
465	cout << "Now in subroutine gfan" << endl;
466	#endif
467	ring inputRing=currRing; // The ring the user entered
468	ring rootRing; // The ring associated to the target ordering
469	ideal res;
470	//matrix ineq; //Matrix containing the boundary inequalities
471	facet *fRoot;
472
473	/*
474	1. Select target order, say dp.
475	2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
476	3. getConeNormals
477	*/
478
479
480	/* Construct a new ring which will serve as our root
481	Does not yet work as expected. Will work fine with order dp,Dp but otherwise hangs in getGB
482	*/
483	rootRing=rCopy0(currRing);
484	rootRing->order[0]=ringorder_dp;
485	rComplete(rootRing);
486	//rChangeCurrRing(rootRing);
487	ideal rootIdeal;
488	/* Fetch the inputIdeal into our rootRing */
489	map m=(map)idInit(IDELEMS(inputIdeal),0);
490	//rootIdeal=fast_map(inputIdeal,inputRing,(ideal)m, currRing);
491	#ifdef gfan_DEBUG
492	cout << "Root ideal is " << endl;
493	idPrint(rootIdeal);
494	cout << "The current ring is " << endl;
495	rWrite(rootRing);
496	cout << endl;
497	#endif
498
499	gcone *gcRoot = new gcone(); //Instantiate the sink
500	gcone *gcAct;
501	gcAct = gcRoot;
502	gcAct->numVars=pVariables;
503	gcAct->getGB(inputIdeal);
504	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
505	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
506	/*Now it is time to compute the search facets, respectively start the reverse search.
507	But since we are in the root all facets should be search facets. IS THIS TRUE?
508	MIND: AS OF NOW, THE LIST OF FACETS IS NOT PURGED OF NON-FLIPPAPLE FACETS
509	*/
510
511	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
512	The return type will then be of type LIST_CMD
513	Assume gfan has finished, thus we have enumerated all the cones
514	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
515	Groebner Basis and merge this somehow with LIST_CMD
516	=> Count the cones!
517	*/
518
519	res=gcAct->gcBasis;
520	//cout << fRoot << endl;
521	return res;
522	//return GBlist;
523	}
524	/*
525	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
526	*/
527	#endif

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

Download in other formats: