Context Navigation

source: git/kernel/gfan.cc @ 224d153

spielwiese

Last change on this file since 224d153 was 224d153, checked in by Martin Monerjan, 14 years ago
More bugfixes in normalize Implemented list of (codim-2)-facets using the facet structure. NOTE: Need to save this when deleting cones in order to preserve the list of comparable facets git-svn-id: file:///usr/local/Singular/svn/trunk@11852 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 49.0 KB

Line
1	/*
2	Compute the Groebner fan of an ideal
3	$Author: monerjan $
4	$Date: 2009-05-29 07:52:12 $
5	$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.57 2009-05-29 07:52:12 monerjan Exp $
6	$Id: gfan.cc,v 1.57 2009-05-29 07:52:12 monerjan Exp $
7	*/
8
9	#include "mod2.h"
10
11	#ifdef HAVE_GFAN
12
13	#include "kstd1.h"
14	#include "kutil.h"
15	#include "intvec.h"
16	#include "polys.h"
17	#include "ideals.h"
18	#include "kmatrix.h"
19	#include "fast_maps.h" //Mapping of ideals
20	#include "maps.h"
21	#include "ring.h"
22	#include "prCopy.h"
23	#include <iostream>
24	#include <bitset>
25	#include <fstream> //read-write cones to files
26	#include <gmp.h>
27	//#include <gmpxx.h>
28
29	/DO NOT REMOVE THIS/
30	#ifndef GMPRATIONAL
31	#define GMPRATIONAL
32	#endif
33
34	//Hacks for different working places
35	#define ITWM
36
37	#ifdef UNI
38	#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
39	#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
40	#endif
41
42	#ifdef HOME
43	#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
44	#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
45	#endif
46
47	#ifdef ITWM
48	#include "/u/slg/monerjan/cddlib/include/setoper.h"
49	#include "/u/slg/monerjan/cddlib/include/cdd.h"
50	#include "/u/slg/monerjan/cddlib/include/cddmp.h"
51	#endif
52
53	#ifndef gfan_DEBUG
54	#define gfan_DEBUG
55	#endif
56
57	//#include gcone.h
58	using namespace std;
59
60	/**
61	*\brief Class facet
62	* Implements the facet structure as a linked list
63	*
64	*/
65	class facet
66	{
67	private:
68	/** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
69	intvec *fNormal;
70
71	/** \brief An interior point of the facet*/
72	intvec *interiorPoint;
73
74	/** \brief Universal Cone Number
75	* The number of the cone the facet belongs to
76	*/
77	int UCN;
78
79	/** \brief The codim of the facet
80	*/
81	int codim;
82
83	/** \brief The Groebner basis on the other side of a shared facet
84	*
85	* In order not to have to compute the flipped GB twice we store the basis we already get
86	* when identifying search facets. Thus in the next step of the reverse search we can
87	* just copy the old cone and update the facet and the gcBasis.
88	* facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
89	*/
90	ideal flipGB; //The Groebner Basis on the other side, computed via gcone::flip
91
92	public:
93	//bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i]
94	bool isIncoming; //Is the facet incoming or outgoing?
95	facet *next; //Pointer to next facet
96	facet *codim2Ptr; //Pointer to (codim-2)-facet. Bit of recursion here ;-)
97	//intvec codim2Normals =(intvec)omAlloc0(sizeof(intvec*)); //Integer matrix containing the (codim-2)-facets
98
99	/** The default constructor. Do I need a constructor of type facet(intvec)? */
100	facet()
101	{
102	// Pointer to next facet. */
103	/* Defaults to NULL. This way there is no need to check explicitly */
104	this->next=NULL;
105	this->UCN=NULL;
106	this->codim2Ptr=NULL;
107	this->codim=1; //default for (codim-1)-facets
108	}
109
110	/** \brief Constructor for facets of codim >= 2
111	*/
112	facet(int const &n)
113	{
114	this->next=NULL;
115	this->UCN=NULL;
116	this->codim2Ptr=NULL;
117	if(n>1)
118	{
119	this->codim=n;
120	}//NOTE Handle exception here!
121	}
122
123	/** The default destructor */
124	~facet(){;}
125
126	/** Stores the facet normal \param intvec*/
127	void setFacetNormal(intvec *iv)
128	{
129	this->fNormal = ivCopy(iv);
130	}
131
132	/** Hopefully returns the facet normal */
133	intvec *getFacetNormal()
134	{
135	return this->fNormal;
136	}
137
138	/** Method to print the facet normal*/
139	void printNormal()
140	{
141	fNormal->show();
142	}
143
144	/** Store the flipped GB*/
145	void setFlipGB(ideal I)
146	{
147	this->flipGB=I;
148	}
149
150	/** Return the flipped GB*/
151	ideal getFlipGB()
152	{
153	return this->flipGB;
154	}
155
156	/** Print the flipped GB*/
157	void printFlipGB()
158	{
159	idShow(this->flipGB);
160	}
161
162	void setUCN(int n)
163	{
164	this->UCN=n;
165	}
166
167	int getUCN()
168	{
169	return this->UCN;
170	}
171
172	void setInteriorPoint(intvec *iv)
173	{
174	this->interiorPoint = ivCopy(iv);
175	}
176
177	intvec *getInteriorPoint()
178	{
179	return this->interiorPoint;
180	}
181
182	/*bool isFlippable(intvec &load)
183	{
184	bool res=TRUE;
185	int jj;
186	for (int jj = 0; jj<load.length(); jj++)
187	{
188	intvec *ivCanonical = new intvec(load.length());
189	(*ivCanonical)[jj]=1;
190	if (ivMult(&load,ivCanonical)<0)
191	{
192	res=FALSE;
193	break;
194	}
195	}
196	return res;
197
198	/*while (dotProduct(load,ivCanonical)>=0)
199	{
200	if (jj!=this->numVars)
201	{
202	intvec *ivCanonical = new intvec(this->numVars);
203	(*ivCanonical)[jj]=1;
204	res=TRUE;
205	jj += 1;
206	}
207	}
208	if (jj==this->numVars)
209	{
210	delete ivCanonical;
211	return FALSE;
212	}
213	else
214	{
215	delete ivCanonical;
216	return TRUE;
217	}*/
218	//}//bool isFlippable(facet &f)
219
220
221	friend class gcone; //Bad style
222	};
223
224	/**
225	*\brief Implements the cone structure
226	*
227	* A cone is represented by a linked list of facet normals
228	* @see facet
229	*/
230	/*class gcone
231	finally this should become s.th. like gconelib.{h,cc} to provide an API
232	*/
233	class gcone
234	{
235	private:
236	ring rootRing; //good to know this -> generic walk
237	ideal inputIdeal; //the original
238	ring baseRing; //the basering of the cone
239	/* TODO in order to save memory use pointers to rootRing and inputIdeal instead */
240	intvec *ivIntPt; //an interior point of the cone
241	int UCN; //unique number of the cone
242
243	public:
244	/** \brief Pointer to the first facet */
245	facet *facetPtr; //Will hold the adress of the first facet; set by gcone::getConeNormals
246
247	/** # of variables in the ring */
248	int numVars; //#of variables in the ring
249
250	/** # of facets of the cone
251	* This value is set by gcone::getConeNormals
252	*/
253	int numFacets; //#of facets of the cone
254
255	/** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
256	ideal gcBasis; //GB of the cone, set by gcone::getGB();
257	gcone next; //Pointer to previous* cone in search tree
258
259	/** \brief Default constructor.
260	*
261	* Initialises this->next=NULL and this->facetPtr=NULL
262	*/
263	gcone()
264	{
265	this->next=NULL;
266	this->facetPtr=NULL; //maybe this->facetPtr = new facet();
267	this->baseRing=currRing;
268	this->UCN=1;
269	this->numFacets=0;
270	}
271
272	/** \brief Constructor with ring and ideal
273	*
274	* This constructor takes the root ring and the root ideal as parameters and stores
275	* them in the private members gcone::rootRing and gcone::inputIdeal
276	* Since knowledge of the root ring is only needed when using reverse search,
277	* this constructor is not needed when using the "second" method
278	*/
279	gcone(ring r, ideal I)
280	{
281	this->next=NULL;
282	this->facetPtr=NULL;
283	this->rootRing=r;
284	this->inputIdeal=I;
285	this->baseRing=currRing;
286	this->UCN=1;
287	this->numFacets=0;
288	}
289
290	/** \brief Copy constructor
291	*
292	* Copies one cone, sets this->gcBasis to the flipped GB and reverses the
293	* direction of the according facet normal
294	*/
295	//NOTE Prolly need to specify the facet to flip over
296	gcone(const gcone& gc, const facet &f)
297	{
298	this->next=NULL;
299	this->numVars=gc.numVars;
300	this->UCN=(gc.UCN)+1; //add 1 to the UCN of previous cone. This is NOT UNIQUE!
301	facet *fAct= new facet();
302	this->facetPtr=fAct;
303
304	intvec *ivtmp = new intvec(this->numVars);
305	ivtmp = gc.facetPtr->getFacetNormal();
306
307	ideal gb;
308	gb=gc.facetPtr->getFlipGB();
309	this->gcBasis=gb; //this cone's GB is the flipped GB
310
311	/Reverse direction of the facet normal to make it an inner normal/
312	for (int ii=0; ii<this->numVars;ii++)
313	{
314	(ivtmp)[ii]=-(ivtmp)[ii];
315	}
316
317	fAct->setFacetNormal(ivtmp);
318	delete ivtmp;
319	delete fAct;
320	}
321
322	/** \brief Default destructor
323	*/
324	~gcone()
325	{
326	//NOTE SAVE THE FACET STRUCTURE!!!
327	}
328
329
330	/** \brief Set the interior point of a cone */
331	void setIntPoint(intvec *iv)
332	{
333	this->ivIntPt=ivCopy(iv);
334	}
335
336	/** \brief Return the interior point */
337	intvec *getIntPoint()
338	{
339	return this->ivIntPt;
340	}
341
342	void showIntPoint()
343	{
344	ivIntPt->show();
345	}
346
347	void showFacets()
348	{
349	facet *f = new facet();
350	f = this->facetPtr;
351	intvec *iv = new intvec(this->numVars);
352
353	while (f->next!=NULL)
354	{
355	iv = f->getFacetNormal();
356	iv->show();
357	f=f->next;
358	}
359	//delete iv;
360	//delete f;
361	}
362
363	/** \brief Set gcone::numFacets */
364	void setNumFacets()
365	{
366	}
367
368	/** \brief Get gcone::numFacets */
369	int getNumFacets()
370	{
371	return this->numFacets;
372	}
373
374	int getUCN()
375	{
376	return this->UCN;
377	}
378
379	/** \brief Compute the normals of the cone
380	*
381	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
382	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
383	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
384	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
385	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
386	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
387	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
388	*
389	* Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
390	* an interior point of the cone.
391	*/
392	void getConeNormals(ideal const &I, bool compIntPoint=FALSE)
393	{
394	#ifdef gfan_DEBUG
395	std::cout << "* Computing Inequalities... *" << std::endl;
396	#endif
397	//All variables go here - except ineq matrix and *v, see below
398	int lengthGB=IDELEMS(I); // # of polys in the groebner basis
399	int pCompCount; // # of terms in a poly
400	poly aktpoly;
401	int numvar = pVariables; // # of variables in a polynomial (or ring?)
402	int leadexp[numvar]; // dirty hack of exp.vects
403	int aktexp[numvar];
404	int cols,rows; // will contain the dimensions of the ineq matrix - deprecated by
405	dd_rowrange ddrows;
406	dd_colrange ddcols;
407	dd_rowset ddredrows; // # of redundant rows in ddineq
408	dd_rowset ddlinset; // the opposite
409	dd_rowindex ddnewpos; // all to make dd_Canonicalize happy
410	dd_NumberType ddnumb=dd_Integer; //Number type
411	dd_ErrorType dderr=dd_NoError; //
412	// End of var declaration
413	#ifdef gfan_DEBUG
414	cout << "The Groebner basis has " << lengthGB << " elements" << endl;
415	cout << "The current ring has " << numvar << " variables" << endl;
416	#endif
417	cols = numvar;
418
419	//Compute the # inequalities i.e. rows of the matrix
420	rows=0; //Initialization
421	for (int ii=0;ii<IDELEMS(I);ii++)
422	{
423	aktpoly=(poly)I->m[ii];
424	rows=rows+pLength(aktpoly)-1;
425	}
426	#ifdef gfan_DEBUG
427	cout << "rows=" << rows << endl;
428	cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
429	#endif
430	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
431	dd_set_global_constants();
432	ddrows=rows;
433	ddcols=cols;
434	dd_MatrixPtr ddineq; //Matrix to store the inequalities
435	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
436
437	// We loop through each g\in GB and compute the resulting inequalities
438	for (int i=0; i<IDELEMS(I); i++)
439	{
440	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
441	pCompCount=pLength(aktpoly); //How many terms does aktpoly consist of?
442	#ifdef gfan_DEBUG
443	cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
444	#endif
445
446	int v=(int )omAlloc((numvar+1)*sizeof(int));
447	pGetExpV(aktpoly,v); //find the exp.vect in v[1],...,v[n], use pNext(p)
448
449	//Store leadexp for aktpoly
450	for (int kk=0;kk<numvar;kk++)
451	{
452	leadexp[kk]=v[kk+1];
453	//Since we need to know the difference of leadexp with the other expvects we do nothing here
454	//but compute the diff below
455	}
456
457
458	while (pNext(aktpoly)!=NULL) //move to next term until NULL
459	{
460	aktpoly=pNext(aktpoly);
461	pSetm(aktpoly); //doesn't seem to help anything
462	pGetExpV(aktpoly,v);
463	for (int kk=0;kk<numvar;kk++)
464	{
465	aktexp[kk]=v[kk+1];
466	//ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk]; //dito
467	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
468	}
469	aktmatrixrow=aktmatrixrow+1;
470	}//while
471
472	} //for
473
474	//Maybe add another row to contain the constraints of the standard simplex?
475
476	#ifdef gfan_DEBUG
477	cout << "The inequality matrix is" << endl;
478	dd_WriteMatrix(stdout, ddineq);
479	#endif
480
481	// The inequalities are now stored in ddineq
482	// Next we check for superflous rows
483	ddredrows = dd_RedundantRows(ddineq, &dderr);
484	if (dderr!=dd_NoError) // did an error occur?
485	{
486	dd_WriteErrorMessages(stderr,dderr); //if so tell us
487	} else
488	{
489	cout << "Redundant rows: ";
490	set_fwrite(stdout, ddredrows); //otherwise print the redundant rows
491	}//if dd_Error
492
493	//Remove reduntant rows here!
494	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
495	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
496	ddcols = ddineq->colsize;
497	#ifdef gfan_DEBUG
498	cout << "Having removed redundancies, the normals now read:" << endl;
499	dd_WriteMatrix(stdout,ddineq);
500	cout << "Rows = " << ddrows << endl;
501	cout << "Cols = " << ddcols << endl;
502	#endif
503
504	/Write the normals into class facet/
505	#ifdef gfan_DEBUG
506	cout << "Creating list of normals" << endl;
507	#endif
508	/The pointer fRoot should be the return value of this function*/
509	facet *fRoot = new facet(); //instantiate new facet
510	this->facetPtr = fRoot; //set variable facetPtr of class gcone to first facet
511	facet *fAct; //instantiate pointer to active facet
512	fAct = fRoot; //Seems to do the trick. fRoot and fAct have to point to the same adress!
513	//std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
514	for (int kk = 0; kk<ddrows; kk++)
515	{
516	intvec *load = new intvec(this->numVars); //intvec to store a single facet normal that will then be stored via setFacetNormal
517	for (int jj = 1; jj <ddcols; jj++)
518	{
519	double foo;
520	foo = mpq_get_d(ddineq->matrix[kk][jj]);
521	#ifdef gfan_DEBUG
522	std::cout << "fAct is " << foo << " at " << fAct << std::endl;
523	#endif
524	(*load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
525	}//for (int jj = 1; jj <ddcols; jj++)
526
527	/Quick'n'dirty hack for flippability/
528	bool isFlippable=FALSE;
529	for (int jj = 0; jj<load->length(); jj++)
530	{
531	intvec *ivCanonical = new intvec(load->length());
532	(*ivCanonical)[jj]=1;
533	cout << "dotProd=" << dotProduct(load,ivCanonical) << endl;
534	if (dotProduct(load,ivCanonical)<0)
535	//if (ivMult(load,ivCanonical)<0)
536	{
537	isFlippable=TRUE;
538	break; //URGHS
539	}
540	}
541	if (isFlippable==FALSE)
542	{
543	cout << "Ignoring facet";
544	load->show();
545	//fAct->next=NULL;
546	}
547	else
548	{ /Now load should be full and we can call setFacetNormal/
549	fAct->setFacetNormal(load);
550	fAct->next = new facet();
551	fAct=fAct->next; //this should definitely not be called in the above while-loop :D
552	this->numFacets++;
553	}//if (isFlippable==FALSE)
554	//delete load;
555	}//for (int kk = 0; kk<ddrows; kk++)
556
557	/*
558	Now we should have a linked list containing the facet normals of those facets that are
559	-irredundant
560	-flipable
561	Adressing is done via *facetPtr
562	*/
563
564	if (compIntPoint==TRUE)
565	{
566	intvec *iv = new intvec(this->numVars);
567	interiorPoint(ddineq, *iv); //NOTE ddineq contains non-flippable facets
568	this->setIntPoint(iv); //stores the interior point in gcone::ivIntPt
569	//delete iv;
570	}
571
572	//Compute the number of facets
573	// wrong because ddineq->rowsize contains only irredundant facets. There can still be
574	// non-flippable ones!
575	//this->numFacets=ddineq->rowsize;
576
577	//Clean up but don't delete the return value! (Whatever it will turn out to be)
578	//dd_FreeMatrix(ddineq);
579	//set_free(ddredrows);
580	//free(ddnewpos);
581	//set_free(ddlinset);
582	//NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
583	//THIS SUCKS
584	//dd_free_global_constants();
585
586	}//method getConeNormals(ideal I)
587
588
589	/** \brief Compute the Groebner Basis on the other side of a shared facet
590	*
591	* Implements algorithm 4.3.2 from Anders' thesis.
592	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
593	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
594	* is parallel to \f$ leadexp(g) \f$
595	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
596	* Other possibilities include computing the rank of the matrix consisting of the vectors in question and
597	* computing an interior point of the facet and taking all terms having the same weight with respect
598	* to this interior point.
599	*\param ideal, facet
600	* Input: a marked,reduced Groebner basis and a facet
601	*/
602	void flip(ideal gb, facet *f) //Compute "the other side"
603	{
604	intvec *fNormal = new intvec(this->numVars); //facet normal, check for parallelity
605	fNormal = f->getFacetNormal(); //read this->fNormal;
606	#ifdef gfan_DEBUG
607	std::cout << "===" << std::endl;
608	std::cout << "running gcone::flip" << std::endl;
609	// std::cout << "fNormal=";
610	// fNormal->show();
611	// std::cout << std::endl;
612	#endif
613	/1st step: Compute the initial ideal/
614	poly initialFormElement[IDELEMS(gb)]; //array of #polys in GB to store initial form
615	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
616	poly aktpoly;
617	intvec *check = new intvec(this->numVars); //array to store the difference of LE and v
618
619	for (int ii=0;ii<IDELEMS(gb);ii++)
620	{
621	aktpoly = (poly)gb->m[ii];
622	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
623	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
624	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
625	initialFormElement[ii]=pHead(aktpoly);
626
627	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
628	{
629	aktpoly=pNext(aktpoly); //next term
630	pSetm(aktpoly);
631	pGetExpV(aktpoly,v);
632	/* Convert (int)v into (intvec)check */
633	for (int jj=0;jj<this->numVars;jj++)
634	{
635	//cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
636	//cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
637	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
638	}
639	#ifdef gfan_DEBUG
640	// cout << "check=";
641	// check->show();
642	// cout << endl;
643	#endif
644	//TODO why not check, fNormal????
645	if (isParallel(check,fNormal)) //pass *check when
646	{
647	// cout << "Parallel vector found, adding to initialFormElement" << endl;
648	initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
649	}
650	}//while
651	#ifdef gfan_DEBUG
652	cout << "Initial Form=";
653	pWrite(initialFormElement[ii]);
654	cout << "---" << endl;
655	#endif
656	/Now initialFormElement must be added to (ideal)initialForm /
657	initialForm->m[ii]=initialFormElement[ii];
658	}//for
659	#ifdef gfan_DEBUG
660	cout << "Initial ideal is: " << endl;
661	idShow(initialForm);
662	//f->printFlipGB();
663	cout << "===" << endl;
664	#endif
665	//delete check;
666
667	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
668	/*Substep 2.1
669	compute $G_{-\alpha}(in_v(I))
670	see journal p. 66
671	*/
672	ring srcRing=currRing;
673
674	//intvec *negfNormal = new intvec(this->numVars);
675	//negfNormal=ivNeg(fNormal);
676	ring tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
677	rChangeCurrRing(tmpRing);
678
679	//rWrite(currRing); cout << endl;
680
681	ideal ina;
682	ina=idrCopyR(initialForm,srcRing);
683	#ifdef gfan_DEBUG
684	cout << "ina=";
685	idShow(ina); cout << endl;
686	#endif
687	ideal H;
688	H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
689	idSkipZeroes(H);
690	#ifdef gfan_DEBUG
691	// cout << "H="; idShow(H); cout << endl;
692	#endif
693	/*Substep 2.2
694	do the lifting and mark according to H
695	*/
696	rChangeCurrRing(srcRing);
697	ideal srcRing_H;
698	ideal srcRing_HH;
699	srcRing_H=idrCopyR(H,tmpRing);
700	#ifdef gfan_DEBUG
701	// cout << "srcRing_H = ";
702	// idShow(srcRing_H); cout << endl;
703	#endif
704	srcRing_HH=ffG(srcRing_H,this->gcBasis);
705	#ifdef gfan_DEBUG
706	// cout << "srcRing_HH = ";
707	// idShow(srcRing_HH); cout << endl;
708	#endif
709	/*Substep 2.2.1
710	Mark according to G_-\alpha
711	Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
712	we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
713	represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
714	Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
715	compute the difference accordingly
716	*/
717	dd_set_global_constants();
718	bool markingsAreCorrect=FALSE;
719	dd_MatrixPtr intPointMatrix;
720	int iPMatrixRows=0;
721	dd_rowrange aktrow=0;
722	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
723	{
724	poly aktpoly=(poly)srcRing_HH->m[ii];
725	iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
726	}
727	/* additionally one row for the standard-simplex and another for a row that becomes 0 during
728	construction of the differences
729	*/
730	intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
731	intPointMatrix->numbtype=dd_Integer; //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
732
733	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
734	{
735	markingsAreCorrect=FALSE; //crucial to initialise here
736	poly aktpoly=srcRing_HH->m[ii];
737	/Comparison of leading monomials is done via exponent vectors/
738	for (int jj=0;jj<IDELEMS(H);jj++)
739	{
740	int src_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
741	int dst_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
742	pGetExpV(aktpoly,src_ExpV);
743	rChangeCurrRing(tmpRing); //this ring change is crucial!
744	pGetExpV(pCopy(H->m[ii]),dst_ExpV);
745	rChangeCurrRing(srcRing);
746	bool expVAreEqual=TRUE;
747	for (int kk=1;kk<=this->numVars;kk++)
748	{
749	#ifdef gfan_DEBUG
750	//cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
751	#endif
752	if (src_ExpV[kk]!=dst_ExpV[kk])
753	{
754	expVAreEqual=FALSE;
755	}
756	}
757	//if (src_ExpV == dst_ExpV)
758	if (expVAreEqual==TRUE)
759	{
760	markingsAreCorrect=TRUE; //everything is fine
761	#ifdef gfan_DEBUG
762	// cout << "correct markings" << endl;
763	#endif
764	}//if (pHead(aktpoly)==pHead(H->m[jj])
765	delete src_ExpV;
766	delete dst_ExpV;
767	}//for (int jj=0;jj<IDELEMS(H);jj++)
768
769	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
770	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
771	if (markingsAreCorrect==TRUE)
772	{
773	pGetExpV(aktpoly,leadExpV);
774	}
775	else
776	{
777	rChangeCurrRing(tmpRing);
778	pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
779	rChangeCurrRing(srcRing);
780	}
781	/compute differences of the expvects/
782	while (pNext(aktpoly)!=NULL)
783	{
784	/*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
785	is not omitted when computing the differences*/
786	if(markingsAreCorrect==TRUE)
787	{
788	aktpoly=pNext(aktpoly);
789	pGetExpV(aktpoly,v);
790	}
791	else
792	{
793	pGetExpV(pHead(aktpoly),v);
794	markingsAreCorrect=TRUE;
795	}
796
797	for (int jj=0;jj<this->numVars;jj++)
798	{
799	/Store into ddMatrix/
800	dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
801	}
802	aktrow +=1;
803	}
804	delete v;
805	delete leadExpV;
806	}//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
807	/Now we add the constraint for the standard simplex/
808	/NOTE:Might actually work without the standard simplex/
809	dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
810	for (int jj=1;jj<=this->numVars;jj++)
811	{
812	dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
813	}
814	dd_WriteMatrix(stdout,intPointMatrix);
815	intvec *iv_weight = new intvec(this->numVars);
816	interiorPoint(intPointMatrix, *iv_weight); //iv_weight now contains the interior point
817	dd_FreeMatrix(intPointMatrix);
818	dd_free_global_constants();
819
820	/*Step 3
821	turn the minimal basis into a reduced one
822	*/
823	int i,j;
824	ring dstRing=rCopy0(srcRing);
825	i=rBlocks(srcRing);
826
827	dstRing->order=(int )omAlloc((i+1)sizeof(int));
828	for(j=i;j>0;j--)
829	{
830	dstRing->order[j]=srcRing->order[j-1];
831	dstRing->block0[j]=srcRing->block0[j-1];
832	dstRing->block1[j]=srcRing->block1[j-1];
833	if (srcRing->wvhdl[j-1] != NULL)
834	{
835	dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
836	}
837	}
838	dstRing->order[0]=ringorder_a;
839	dstRing->order[1]=ringorder_dp;
840	dstRing->order[2]=ringorder_C;
841	dstRing->wvhdl[0] =( int )omAlloc((iv_weight->length())sizeof(int));
842
843	for (int ii=0;ii<this->numVars;ii++)
844	{
845	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
846	}
847	rComplete(dstRing);
848
849	// NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
850	// Thus:
851	//ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
852	//cout << "PLING" << endl;
853	/*ring dstRing=rCopy0(srcRing);
854	for (int ii=0;ii<this->numVars;ii++)
855	{
856	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
857	}*/
858	rChangeCurrRing(dstRing);
859	#ifdef gfan_DEBUG
860	rWrite(dstRing); cout << endl;
861	#endif
862	ideal dstRing_I;
863	dstRing_I=idrCopyR(srcRing_HH,srcRing);
864	//validOpts<1>=TRUE;
865	#ifdef gfan_DEBUG
866	//idShow(dstRing_I);
867	#endif
868	BITSET save=test;
869	test\|=Sy_bit(OPT_REDSB);
870	test\|=Sy_bit(6); //OPT_DEBUG
871	dstRing_I=kStd(idrCopyR(this->inputIdeal,this->rootRing),NULL,testHomog,NULL);
872	kInterRed(dstRing_I);
873	idSkipZeroes(dstRing_I);
874	test=save;
875	/End of step 3 - reduction/
876
877	f->setFlipGB(dstRing_I);//store the flipped GB
878	#ifdef gfan_DEBUG
879	cout << "Flipped GB is: " << endl;
880	f->printFlipGB();
881	#endif
882	}//void flip(ideal gb, facet *f)
883
884	/** \brief Compute the remainder of a polynomial by a given ideal
885	*
886	* Compute \f$ f^{\mathcal{G}} \f$
887	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
888	* However, since we are only interested in the remainder, there is no need to
889	* compute the factors \f$ a_i \f$
890	*/
891	//NOTE: Should be replaced by kNF or kNF2
892	poly restOfDiv(poly const &f, ideal const &I)
893	{
894	cout << "Entering restOfDiv" << endl;
895	poly p=f;
896	//pWrite(p);
897	//poly r=kCreateZeroPoly(,currRing,currRing); //The 0-polynomial, hopefully
898	poly r=NULL; //The zero polynomial
899	int ii;
900	bool divOccured;
901
902	while (p!=NULL)
903	{
904	ii=1;
905	divOccured=FALSE;
906
907	while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
908	{
909	if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
910	{
911	poly step1,step2,step3;
912	//cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
913	step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
914	//cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
915	step2 = ppMult_qq(step1, I->m[ii-1]);
916	step3 = pSub(pCopy(p), step2);
917	//p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
918	//pSetm(p);
919	pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
920	p=step3;
921	//pWrite(p);
922	divOccured=TRUE;
923	}
924	else
925	{
926	ii += 1;
927	}//if (pLmDivisibleBy(I->m[ii],p,currRing))
928	}//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
929	if (divOccured==FALSE)
930	{
931	//cout << "TICK 5" << endl;
932	r=pAdd(pCopy(r),pHead(p));
933	pSetm(r);
934	pSort(r);
935	//cout << "r="; pWrite(r); cout << endl;
936
937	if (pLength(p)!=1)
938	{
939	p=pSub(pCopy(p),pHead(p)); //Here it may occur that p=0 instead of p=NULL
940	}
941	else
942	{
943	p=NULL; //Hack to correct this situation
944	}
945	//cout << "p="; pWrite(p);
946	}//if (divOccured==FALSE)
947	}//while (p!=0)
948	return r;
949	}//poly restOfDiv(poly const &f, ideal const &I)
950
951	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
952	*/
953	//NOTE: use kNF or kNF2 instead of restOfDivision
954	ideal ffG(ideal const &H, ideal const &G)
955	{
956	cout << "Entering ffG" << endl;
957	int size=IDELEMS(H);
958	ideal res=idInit(size,1);
959	poly temp1, temp2, temp3; //polys to temporarily store values for pSub
960	for (int ii=0;ii<size;ii++)
961	{
962	res->m[ii]=restOfDiv(H->m[ii],G);
963	//res->m[ii]=kNF(H->m[ii],G);
964	temp1=H->m[ii];
965	temp2=res->m[ii];
966	temp3=pSub(temp1, temp2);
967	res->m[ii]=temp3;
968	//res->m[ii]=pSub(temp1,temp2); //buggy
969	//pSort(res->m[ii]);
970	//pSetm(res->m[ii]);
971	//cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
972	}
973	return res;
974	}
975
976	/** \brief Compute a Groebner Basis
977	*
978	* Computes the Groebner basis and stores the result in gcone::gcBasis
979	*\param ideal
980	*\return void
981	*/
982	void getGB(ideal const &inputIdeal)
983	{
984	ideal gb;
985	gb=kStd(inputIdeal,NULL,testHomog,NULL);
986	idSkipZeroes(gb);
987	this->gcBasis=gb; //write the GB into gcBasis
988	}//void getGB
989
990	/** \brief The Generic Groebner Walk due to FJLT
991	* Needed for computing the search facet
992	*/
993	ideal GenGrbWlk(ideal, ideal)
994	{
995	}//GGW
996
997	/** \brief Compute the negative of a given intvec
998	*/
999	intvec *ivNeg(intvec const &iv)
1000	{
1001	intvec *res = new intvec(this->numVars);
1002	for(int ii=0;ii<this->numVars;ii++)
1003	{
1004	res[ii] = -(int)iv[ii];
1005	}
1006	return res;
1007	}
1008
1009
1010	/** \brief Compute the dot product of two intvecs
1011	*
1012	*/
1013	int dotProduct(intvec const &iva, intvec const &ivb)
1014	{
1015	//intvec iva=a;
1016	//intvec ivb=b;
1017	int res=0;
1018	for (int i=0;i<this->numVars;i++)
1019	{
1020	res = res+(iva[i]*ivb[i]);
1021	}
1022	return res;
1023	}//int dotProduct
1024
1025	/** \brief Check whether two intvecs are parallel
1026	*
1027	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
1028	*/
1029	bool isParallel(intvec const &a, intvec const &b)
1030	{
1031	int lhs,rhs;
1032	lhs=dotProduct(a,b)*dotProduct(a,b);
1033	rhs=dotProduct(a,a)*dotProduct(b,b);
1034	//cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
1035	if (lhs==rhs)
1036	{
1037	return TRUE;
1038	}
1039	else
1040	{
1041	return FALSE;
1042	}
1043	}//bool isParallel
1044
1045	/** \brief Compute an interior point of a given cone
1046	* Result will be written into intvec iv
1047	*/
1048	void interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
1049	{
1050	dd_LPPtr lp,lpInt;
1051	dd_ErrorType err=dd_NoError;
1052	dd_LPSolverType solver=dd_DualSimplex;
1053	dd_LPSolutionPtr lpSol=NULL;
1054	dd_rowset ddlinset,ddredrows; //needed for dd_FindRelativeInterior
1055	dd_rowindex ddnewpos;
1056	dd_NumberType numb;
1057	//M->representation=dd_Inequality;
1058	//M->objective-dd_LPMin; //Not sure whether this is needed
1059
1060	//NOTE: Make this n-dimensional!
1061	//dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1062
1063	//dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
1064	//if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1065	//cout << "Tick 2" << endl;
1066	//dd_WriteMatrix(stdout,M);
1067
1068	lp=dd_Matrix2LP(M, &err);
1069	if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}
1070	if (lp==NULL){cout << "LP is NULL" << endl;}
1071	#ifdef gfan_DEBUG
1072	// dd_WriteLP(stdout,lp);
1073	#endif
1074
1075	lpInt=dd_MakeLPforInteriorFinding(lp);
1076	if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
1077	#ifdef gfan_DEBUG
1078	// dd_WriteLP(stdout,lpInt);
1079	#endif
1080
1081	dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1082	if (err!=dd_NoError)
1083	{
1084	cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
1085	dd_WriteErrorMessages(stdout, err);
1086	}
1087
1088	//dd_LPSolve(lpInt,solver,&err); //This will not result in a point from the relative interior
1089	if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
1090	//cout << "Tick 5" << endl;
1091
1092	//lpSol=dd_CopyLPSolution(lpInt);
1093	if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}
1094	//cout << "Tick 6" << endl;
1095	#ifdef gfan_DEBUG
1096	cout << "Interior point: ";
1097	#endif
1098	for (int ii=1; ii<(lpSol->d)-1;ii++)
1099	{
1100	#ifdef gfan_DEBUG
1101	dd_WriteNumber(stdout,lpSol->sol[ii]);
1102	#endif
1103	/* NOTE This works only as long as gmp returns fractions with the same denominator*/
1104	(iv)[ii-1]=(int)mpz_get_d(mpq_numref(lpSol->sol[ii])); //looks evil, but does the trick
1105	}
1106	dd_FreeLPSolution(lpSol);
1107	dd_FreeLPData(lpInt);
1108	dd_FreeLPData(lp);
1109	set_free(ddlinset);
1110	set_free(ddredrows);
1111
1112	}//void interiorPoint(dd_MatrixPtr const &M)
1113
1114	/** \brief Copy a ring and add a weighted ordering in first place
1115	* Kudos to walkSupport.cc
1116	*/
1117	ring rCopyAndAddWeight(ring const &r, intvec const *ivw)
1118	{
1119	ring res=(ring)omAllocBin(ip_sring_bin);
1120	memcpy4(res,r,sizeof(ip_sring));
1121	res->VarOffset = NULL;
1122	res->ref=0;
1123
1124	if (r->algring!=NULL)
1125	r->algring->ref++;
1126	if (r->parameter!=NULL)
1127	{
1128	res->minpoly=nCopy(r->minpoly);
1129	int l=rPar(r);
1130	res->parameter=(char *)omAlloc(lsizeof(char_ptr));
1131	int i;
1132	for(i=0;i<rPar(r);i++)
1133	{
1134	res->parameter[i]=omStrDup(r->parameter[i]);
1135	}
1136	}
1137
1138	int i=rBlocks(r);
1139	int jj;
1140
1141	res->order =(int )omAlloc((i+1)sizeof(int));
1142	res->block0=(int )omAlloc((i+1)sizeof(int));
1143	res->block1=(int )omAlloc((i+1)sizeof(int));
1144	res->wvhdl =(int *)omAlloc((i+1)sizeof(int**));
1145	for(jj=0;jj<i;jj++)
1146	{
1147	if (r->wvhdl[jj] != NULL)
1148	{
1149	res->wvhdl[jj] = (int*) omMemDup(r->wvhdl[jj-1]);
1150	}
1151	else
1152	{
1153	res->wvhdl[jj+1]=NULL;
1154	}
1155	}
1156
1157	for (jj=0;jj<i;jj++)
1158	{
1159	res->order[jj+1]=r->order[jj];
1160	res->block0[jj+1]=r->block0[jj];
1161	res->block1[jj+1]=r->block1[jj];
1162	}
1163
1164	res->order[0]=ringorder_a;
1165	res->order[1]=ringorder_dp; //basically useless, since that should never be used
1166	int length=ivw->length();
1167	int A=(int )omAlloc(length*sizeof(int));
1168	for (jj=0;jj<length;jj++)
1169	{
1170	A[jj]=(*ivw)[jj];
1171	}
1172	res->wvhdl[0]=(int *)A;
1173	res->block0[0]=1;
1174	res->block1[0]=length;
1175
1176	res->names = (char *)omAlloc0(rVar(r) sizeof(char_ptr));
1177	for (i=rVar(res)-1;i>=0; i--)
1178	{
1179	res->names[i] = omStrDup(r->names[i]);
1180	}
1181	rComplete(res);
1182	return res;
1183	}//rCopyAndAdd
1184
1185	ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1186	{
1187	ring res=rCopy0(currRing);
1188	rComplete(res);
1189	rSetWeightVec(res,(int64*)ivw);
1190	//rChangeCurrRing(rnew);
1191	return res;
1192	}
1193
1194	/** \brief Checks whether a given facet is a search facet
1195	* Determines whether a given facet of a cone is the search facet of a neighbouring cone
1196	* This is done in the following way:
1197	* We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1198	* that is first crossed during the generic walk.
1199	* We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1200	* If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1201	*/
1202	bool isSearchFacet(gcone &gcTmp, facet *testfacet)
1203	{
1204	ring actRing=currRing;
1205	facet facetPtr=(facet)gcTmp.facetPtr;
1206	facet fMin=new facet(facetPtr); //Pointer to the "minimal" facet
1207	//facet *fMin = new facet(tmpcone.facetPtr);
1208	//fMin=tmpcone.facetPtr; //Initialise to first facet of tmpcone
1209	facet *fAct; //Ptr to alpha_i
1210	facet *fCmp; //Ptr to alpha_j
1211	fAct = fMin;
1212	fCmp = fMin->next;
1213
1214	rChangeCurrRing(this->rootRing); //because we compare the monomials in the rootring
1215	poly p=pInit();
1216	poly q=pInit();
1217	intvec *alpha_i = new intvec(this->numVars);
1218	intvec *alpha_j = new intvec(this->numVars);
1219	intvec *sigma = new intvec(this->numVars);
1220	sigma=gcTmp.getIntPoint();
1221
1222	int u=(int )omAlloc((this->numVars+1)*sizeof(int));
1223	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1224	u[0]=0; v[0]=0;
1225	int weight1,weight2;
1226	while(fAct->next->next!=NULL) //NOTE this is ugly. Can it be done without fCmp?
1227	{
1228	/* Get alpha_i and alpha_{i+1} */
1229	alpha_i=fAct->getFacetNormal();
1230	alpha_j=fCmp->getFacetNormal();
1231	#ifdef gfan_DEBUG
1232	alpha_i->show();
1233	alpha_j->show();
1234	#endif
1235	/Compute the dot product of sigma and alpha_{i,j}/
1236	weight1=dotProduct(sigma,alpha_i);
1237	weight2=dotProduct(sigma,alpha_j);
1238	#ifdef gfan_DEBUG
1239	cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
1240	#endif
1241	/Adjust alpha_i and alpha_i+1 accordingly/
1242	for(int ii=1;ii<=this->numVars;ii++)
1243	{
1244	u[ii]=weight1(alpha_i)[ii-1];
1245	v[ii]=weight2(alpha_j)[ii-1];
1246	}
1247
1248	/Now p_weight and q_weight need to be compared as exponent vectors/
1249	pSetCoeff0(p,nInit(1));
1250	pSetCoeff0(q,nInit(1));
1251	pSetExpV(p,u);
1252	pSetm(p);
1253	pSetExpV(q,v);
1254	pSetm(q);
1255	#ifdef gfan_DEBUG
1256	pWrite(p);pWrite(q);
1257	#endif
1258	/*We want to check whether x^p < x^q
1259	=> want to check for return value 1 */
1260	if (pLmCmp(p,q)==1) //i.e. x^q is smaller
1261	{
1262	fMin=fCmp;
1263	fAct=fMin;
1264	fCmp=fCmp->next;
1265	}
1266	else
1267	{
1268	//fAct=fAct->next;
1269	if(fCmp->next!=NULL)
1270	{
1271	fCmp=fCmp->next;
1272	}
1273	else
1274	{
1275	fAct=fAct->next;
1276	}
1277	}
1278	//fAct=fAct->next;
1279	}//while(fAct.facetPtr->next!=NULL)
1280	delete alpha_i,alpha_j,sigma;
1281
1282	/If testfacet was minimal then fMin should still point there /
1283
1284	//if(fMin->getFacetNormal()==ivNeg(testfacet.getFacetNormal()))
1285	#ifdef gfan_DEBUG
1286	cout << "Checking for parallelity" << endl <<" fMin is";
1287	fMin->printNormal();
1288	cout << "testfacet is ";
1289	testfacet->printNormal();
1290	cout << endl;
1291	#endif
1292	if (fMin==gcTmp.facetPtr)
1293	//if(areEqual(fMin->getFacetNormal(),ivNeg(testfacet.getFacetNormal())))
1294	//if (isParallel(fMin->getFacetNormal(),testfacet->getFacetNormal()))
1295	{
1296	cout << "Parallel" << endl;
1297	rChangeCurrRing(actRing);
1298	//delete alpha_min, test;
1299	return TRUE;
1300	}
1301	else
1302	{
1303	cout << "Not parallel" << endl;
1304	rChangeCurrRing(actRing);
1305	//delete alpha_min, test;
1306	return FALSE;
1307	}
1308	}//bool isSearchFacet
1309
1310	/** \brief Check for equality of two intvecs
1311	*/
1312	bool areEqual(intvec const &a, intvec const &b)
1313	{
1314	bool res=TRUE;
1315	for(int ii=0;ii<this->numVars;ii++)
1316	{
1317	if(a[ii]!=b[ii])
1318	{
1319	res=FALSE;
1320	break;
1321	}
1322	}
1323	return res;
1324	}
1325
1326	/** \brief The reverse search algorithm
1327	*/
1328	void reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
1329	{
1330	facet *fAct=new facet();
1331	fAct = gcAct->facetPtr;
1332
1333	while(fAct->next!=NULL) //NOTE NOT SURE WHETHER THIS IS RIGHT! Do I reach EVERY facet or only all but the last?
1334	{
1335	cout << "==========================================================================================="<< endl;
1336	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1337	gcone gcTmp = new gcone(gcAct);
1338	//idShow(gcTmp->gcBasis);
1339	gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1340	#ifdef gfan_DEBUG
1341	facet *f = new facet();
1342	f=gcTmp->facetPtr;
1343	while(f->next!=NULL)
1344	{
1345	f->printNormal();
1346	f=f->next;
1347	}
1348	#endif
1349	gcTmp->showIntPoint();
1350	/recursive part goes gere/
1351	if (isSearchFacet(gcTmp,(facet)gcAct->facetPtr))
1352	{
1353	gcAct->next=gcTmp;
1354	cout << "PING"<< endl;
1355	reverseSearch(gcTmp);
1356	}
1357	else
1358	{
1359	delete gcTmp;
1360	/NOTE remove fAct from linked list. It's no longer needed/
1361	}
1362	/recursion ends/
1363	fAct = fAct->next;
1364	}//while(fAct->next!=NULL)
1365	}//reverseSearch
1366
1367	/** \brief The new method of Markwig and Jensen
1368	* Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
1369	* In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
1370	* Keep a list of facets as a linked list containing an intvec and an integer matrix.
1371	* Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
1372	* soon as we compute this facet again. Comparison of facets is done by...
1373	*/
1374	void noRevS(gcone &gc, bool usingIntPoint=FALSE)
1375	{
1376	facet *fListPtr = new facet(); //The list containing ALL facets we come across
1377	facet *fAct = new facet();
1378	fAct = gc.facetPtr;
1379
1380	#ifdef gfan_DEBUG
1381	cout << "NoRevs" << endl;
1382	cout << "Facets are:" << endl;
1383	gc.showFacets();
1384	#endif
1385
1386	dd_set_global_constants();
1387	dd_rowrange ddrows;
1388	dd_colrange ddcols;
1389	ddrows=2; //Each facet is described by its normal
1390	ddcols=gc.numVars+1; // plus one for the standard simplex
1391	if(usingIntPoint){
1392	while(fAct->next!=NULL)
1393	{
1394	/Compute an interior point for each facet/
1395	dd_MatrixPtr ddineq;
1396	ddineq=dd_CreateMatrix(ddrows,ddcols);
1397	intvec *comp = new intvec(this->numVars);
1398	comp=fAct->getFacetNormal();
1399	for(int ii=0; ii<this->numVars; ii++)
1400	{
1401	dd_set_si(ddineq->matrix[0][ii+1],(*comp)[ii]);
1402	dd_set_si(ddineq->matrix[1][ii+1],1); //Assure we search in the pos. orthant
1403	}
1404	set_addelem(ddineq->linset,1); //We want equality in the first row
1405	//dd_WriteMatrix(stdout,ddineq);
1406	interiorPoint(ddineq,*comp);
1407	/**/
1408	#ifdef gfan_DEBUG
1409	cout << "IP is";
1410	comp->show(); cout << endl;
1411	#endif
1412	//Store the interior point and the UCN
1413	fListPtr->setInteriorPoint( comp );
1414	fListPtr->setUCN( gc.getUCN() );
1415
1416	dd_FreeMatrix(ddineq);
1417	fAct=fAct->next; //iterate
1418	}
1419	}
1420	else/Facets of facets/
1421	{
1422	dd_MatrixPtr ddineq,P,ddakt;
1423	dd_rowset impl_linset, redset;
1424	dd_ErrorType err;
1425	dd_rowindex newpos;
1426
1427	fAct = fListPtr;
1428
1429	#ifdef gfan_DEBUG
1430	cout << "before f2M" << endl;
1431	gc.showFacets();
1432	ddineq = facets2Matrix(gc);
1433	cout << "after f2M" << endl;
1434	gc.showFacets();
1435	#endif
1436
1437	/Now set appropriate linearity/
1438	dd_PolyhedraPtr ddpolyh;
1439	for (int ii=0; ii<this->numFacets; ii++)
1440	{
1441	cout << endl << "------------" << endl;
1442	ddakt = dd_CopyMatrix(ddineq);
1443	set_addelem(ddakt->linset,ii+1);
1444
1445	dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
1446
1447	//dd_WriteMatrix(stdout,ddakt);
1448	ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
1449	P=dd_CopyGenerators(ddpolyh);
1450	dd_WriteMatrix(stdout,P);
1451
1452	/* We loop through each row of P
1453	* normalize it by making all entries integer ones
1454	* and add the resulting vector to the int matrix facet::codim2Facets
1455	*/
1456	this->facetPtr->codim2Ptr = new facet(2); //construct a (codim-2)-facet
1457	facet *codim2Act;
1458	codim2Act = this->facetPtr->codim2Ptr;
1459	for (int jj=1;jj<=P->rowsize;jj++)
1460	{
1461	intvec *n = new intvec(this->numVars);
1462	normalize(P,jj,*n);
1463	codim2Act->setFacetNormal(n);
1464	codim2Act->next = new facet(2);
1465	codim2Act = codim2Act->next;
1466	n->show();
1467	delete n;
1468	}
1469
1470	dd_FreeMatrix(ddakt);
1471	dd_FreePolyhedra(ddpolyh);
1472	}
1473	}
1474
1475
1476	//NOTE Hm, comment in and get a crash for free...
1477	//dd_free_global_constants();
1478	gc.writeConeToFile(gc);
1479
1480	/*2nd step
1481	Choose a facet from fListPtr, flip it and forget the previous cone
1482	*/
1483	fAct = fListPtr;
1484	//gcone *gcTmp = new gcone(gc); //copy
1485	//gcTmp->flip(gcTmp->gcBasis,fAct);
1486	//NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
1487
1488	}//void noRevS(gcone &gc)
1489
1490
1491	/** \brief Make a set of rational vectors into integers
1492	*
1493	* We compute the lcm of the denominators and multiply with this to get integer values
1494	* \param dd_MatrixPtr,intvec
1495	*/
1496	void normalize(dd_MatrixPtr const &M, int const line, intvec &n)
1497	{
1498	mpz_t denom[this->numVars];
1499	for(int ii=0;ii<this->numVars;ii++)
1500	{
1501	mpz_init_set_str(denom[ii],"0",10);
1502	}
1503
1504	mpz_t kgV,tmp;
1505	mpz_init(kgV);
1506	mpz_init(tmp);
1507
1508	for (int ii=0;ii<(M->colsize)-1;ii++)
1509	{
1510	mpz_t z;
1511	mpz_init(z);
1512	mpq_get_den(z,M->matrix[line-1][ii+1]);
1513	mpz_set( denom[ii], z);
1514	mpz_clear(z);
1515	}
1516
1517	/Compute lcm of the denominators/
1518	mpz_set(tmp,denom[0]);
1519	for (int ii=0;ii<(M->colsize)-1;ii++)
1520	{
1521	mpz_lcm(kgV,tmp,denom[ii]);
1522	}
1523
1524	/Multiply the nominators by kgV/
1525	mpq_t qkgV,res;
1526	mpq_init(qkgV);
1527	mpq_set_str(qkgV,"1",10);
1528	mpq_canonicalize(qkgV);
1529
1530	mpq_init(res);
1531	mpq_set_str(res,"1",10);
1532	mpq_canonicalize(res);
1533
1534	mpq_set_num(qkgV,kgV);
1535
1536	// mpq_canonicalize(qkgV);
1537	for (int ii=0;ii<(M->colsize)-1;ii++)
1538	{
1539	mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
1540	n[ii]=(int)mpz_get_d(mpq_numref(res));
1541	}
1542	//mpz_clear(denom[this->numVars]);
1543	mpz_clear(kgV);
1544	mpq_clear(qkgV); mpq_clear(res);
1545
1546	}
1547
1548	/** \brief Construct a dd_MatrixPtr from a cone's list of facets
1549	*
1550	*/
1551	dd_MatrixPtr facets2Matrix(gcone const &gc)
1552	{
1553	facet *fAct = new facet();
1554	fAct = gc.facetPtr;
1555
1556	dd_MatrixPtr M;
1557	dd_rowrange ddrows;
1558	dd_colrange ddcols;
1559	ddcols=(this->numVars)+1;
1560	ddrows=this->numFacets;
1561	dd_NumberType numb = dd_Integer;
1562	M=dd_CreateMatrix(ddrows,ddcols);
1563
1564	//dd_set_global_constants();
1565
1566	intvec *comp = new intvec(this->numVars);
1567
1568	for (int jj=0; jj<M->rowsize; jj++)
1569	{
1570	comp = fAct->getFacetNormal();
1571	for (int ii=0; ii<this->numVars; ii++)
1572	{
1573	dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
1574	}
1575	if(fAct->next!=NULL)
1576	{
1577	fAct = fAct->next;
1578	}
1579	}//jj
1580
1581	//delete fAct;
1582	//delete comp;
1583	return M;
1584	}
1585
1586	/** \brief Write information about a cone into a file on disk
1587	*
1588	* This methods writes the information needed for the "second" method into a file.
1589	* The file's is divided in sections containing information on
1590	* 1) the ring
1591	* 2) the cone's Groebner Basis
1592	* 3) the cone's facets
1593	* Each line contains exactly one date
1594	* Each section starts with its name in CAPITALS
1595	*/
1596	void writeConeToFile(gcone const &gc, bool usingIntPoints=FALSE)
1597	{
1598	ofstream gcOutputFile("/tmp/cone1.gc");
1599	facet *fAct = new facet();
1600	fAct = gc.facetPtr;
1601
1602	if (!gcOutputFile)
1603	{
1604	cout << "Error opening file for writing in writeConeToFile" << endl;
1605	}
1606	else
1607	{ /*gcOutputFile << "UCN" << endl;
1608	gcOutputFile << this->UCN << endl;*/
1609	gcOutputFile << "RING" << endl;
1610	//Write this->gcBasis into file
1611	gcOutputFile << "GCBASIS" << endl;
1612	for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
1613	{
1614	gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
1615	}
1616
1617	gcOutputFile << "FACETS" << endl;
1618	while(fAct->next!=NULL)
1619	{
1620	intvec *iv = new intvec(gc.numVars);
1621	iv=fAct->getFacetNormal();
1622	for (int ii=0;ii<iv->length();ii++)
1623	{
1624	if (ii<iv->length()-1)
1625	{
1626	gcOutputFile << (*iv)[ii] << ",";
1627	}
1628	else
1629	{
1630	gcOutputFile << (*iv)[ii] << endl;
1631	}
1632	}
1633	fAct=fAct->next;
1634	//delete iv; iv=NULL;
1635	}
1636
1637	gcOutputFile.close();
1638	//delete fAct; fAct=NULL;
1639	}
1640
1641	}//writeConeToFile(gcone const &gc)
1642
1643	/** \brief Reads a cone from a file identified by its number
1644	*/
1645	void readConeFromFile(int gcNum)
1646	{
1647	}
1648
1649	friend class facet;
1650	};//class gcone
1651
1652	ideal gfan(ideal inputIdeal)
1653	{
1654	int numvar = pVariables;
1655
1656	enum searchMethod {
1657	reverseSearch,
1658	noRevS
1659	};
1660
1661	searchMethod method;
1662	method = noRevS;
1663	// method = reverseSearch;
1664
1665	#ifdef gfan_DEBUG
1666	cout << "Now in subroutine gfan" << endl;
1667	#endif
1668	ring inputRing=currRing; // The ring the user entered
1669	ring rootRing; // The ring associated to the target ordering
1670	ideal res;
1671	facet *fRoot;
1672
1673	if (method==reverseSearch)
1674	{
1675
1676	/* Construct a new ring which will serve as our root*/
1677	rootRing=rCopy0(currRing);
1678	rootRing->order[0]=ringorder_lp;
1679
1680	rComplete(rootRing);
1681	rChangeCurrRing(rootRing);
1682
1683	/* Fetch the inputIdeal into our rootRing */
1684	map theMap=(map)idMaxIdeal(1); //evil hack!
1685	theMap->preimage=NULL; //neccessary?
1686	ideal rootIdeal;
1687	rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
1688	#ifdef gfan_DEBUG
1689	cout << "Root ideal is " << endl;
1690	idShow(rootIdeal);
1691	cout << "The root ring is " << endl;
1692	rWrite(rootRing);
1693	cout << endl;
1694	#endif
1695
1696	//gcone *gcRoot = new gcone(); //Instantiate the sink
1697	gcone *gcRoot = new gcone(rootRing,rootIdeal);
1698	gcone *gcAct;
1699	gcAct = gcRoot;
1700	gcAct->numVars=pVariables;
1701	gcAct->getGB(rootIdeal); //sets gcone::gcBasis
1702	idShow(gcAct->gcBasis);
1703	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
1704	//gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1705	/*Now it is time to compute the search facets, respectively start the reverse search.
1706	But since we are in the root all facets should be search facets. IS THIS TRUE?
1707	NOTE: Check for flippability is not very sophisticated
1708	*/
1709	//gcAct->reverseSearch(gcAct);
1710	rChangeCurrRing(rootRing);
1711	res=gcRoot->gcBasis;
1712	}//if method==reverSearch
1713
1714	if(method==noRevS)
1715	{
1716	gcone *gcRoot = new gcone(currRing,inputIdeal);
1717	gcone *gcAct;
1718	gcAct = gcRoot;
1719	gcAct->numVars=pVariables;
1720	gcAct->getGB(inputIdeal);
1721	gcAct->getConeNormals(gcAct->gcBasis);
1722	gcAct->noRevS(*gcAct);
1723	res=gcAct->gcBasis;
1724	}
1725
1726	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
1727	The return type will then be of type LIST_CMD
1728	Assume gfan has finished, thus we have enumerated all the cones
1729	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
1730	Groebner Basis and merge this somehow with LIST_CMD
1731	=> Count the cones!
1732	*/
1733	//rChangeCurrRing(rootRing);
1734	//res=gcAct->gcBasis;
1735	//res=gcRoot->gcBasis;
1736	return res;
1737	//return GBlist;
1738	}
1739	/*
1740	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
1741	*/
1742	#endif

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

Download in other formats: