Context Navigation

source: git/kernel/gfan.cc @ 7d10b3

spielwiese

Last change on this file since 7d10b3 was 7d10b3, checked in by Martin Monerjan, 15 years ago
Added and debugged method gcone::normalize to make the rows of a mpq_t matrix into intvecs BUG: Singular crashes, while Singularg terminates normally... git-svn-id: file:///usr/local/Singular/svn/trunk@11822 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 47.7 KB

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

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

Download in other formats: