Context Navigation

source: git/kernel/gfan.cc @ 472eeb

spielwiese

Last change on this file since 472eeb was 472eeb, checked in by Martin Monerjan, 15 years ago
facets2Matrix git-svn-id: file:///usr/local/Singular/svn/trunk@11807 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 45.1 KB

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

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

Download in other formats: