My Project  debian-1:4.1.1-p2+ds-4build4
groebnerCone.cc
Go to the documentation of this file.
1 
3 #include "kernel/ideals.h"
4 #include "Singular/ipid.h"
5 
7 #include "polys/monomials/ring.h"
8 #include "polys/prCopy.h"
9 
10 #include "gfanlib/gfanlib.h"
11 #include "gfanlib/gfanlib_matrix.h"
12 
13 #include "initial.h"
14 #include "tropicalStrategy.h"
15 #include "groebnerCone.h"
16 #include "callgfanlib_conversion.h"
17 #include "containsMonomial.h"
18 #include "tropicalCurves.h"
19 #include "bbcone.h"
20 #include "tropicalDebug.h"
21 
22 #ifndef NDEBUG
23 bool groebnerCone::checkFlipConeInput(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const
24 {
25  /* check first whether interiorPoint lies on the boundary of the cone */
26  if (!polyhedralCone.contains(interiorPoint))
27  {
28  std::cout << "ERROR: interiorPoint is not contained in the Groebner cone!" << std::endl
29  << "cone: " << std::endl
31  << "interiorPoint:" << std::endl
32  << interiorPoint << std::endl;
33  return false;
34  }
35  if (polyhedralCone.containsRelatively(interiorPoint))
36  {
37  std::cout << "ERROR: interiorPoint is contained in the interior of the maximal Groebner cone!" << std::endl
38  << "cone: " << std::endl
40  << "interiorPoint:" << std::endl
41  << interiorPoint << std::endl;
42  return false;
43  }
44  gfan::ZCone hopefullyAFacet = polyhedralCone.faceContaining(interiorPoint);
45  if (hopefullyAFacet.dimension()!=(polyhedralCone.dimension()-1))
46  {
47  std::cout << "ERROR: interiorPoint is not contained in the interior of a facet!" << std::endl
48  << "cone: " << std::endl
50  << "interiorPoint:" << std::endl
51  << interiorPoint << std::endl;
52  return false;
53  }
54  /* check whether facet normal points outwards */
55  gfan::ZCone dual = polyhedralCone.dualCone();
56  if(dual.containsRelatively(facetNormal))
57  {
58  std::cout << "ERROR: facetNormal is not pointing outwards!" << std::endl
59  << "cone: " << std::endl
61  << "facetNormal:" << std::endl
62  << facetNormal << std::endl;
63  return false;
64  }
65  return true;
66 }
67 #endif //NDEBUG
68 
70  polynomialIdeal(NULL),
71  polynomialRing(NULL),
72  polyhedralCone(gfan::ZCone(0)),
73  interiorPoint(gfan::ZVector(0)),
74  currentStrategy(NULL)
75 {
76 }
77 
78 groebnerCone::groebnerCone(const ideal I, const ring r, const tropicalStrategy& currentCase):
79  polynomialIdeal(NULL),
80  polynomialRing(NULL),
81  currentStrategy(&currentCase)
82 {
84  if (r) polynomialRing = rCopy(r);
85  if (I)
86  {
87  polynomialIdeal = id_Copy(I,r);
90  }
91 
92  int n = rVar(polynomialRing);
93  poly g = NULL;
94  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
95  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
96  gfan::ZVector leadexpw = gfan::ZVector(n);
97  gfan::ZVector tailexpw = gfan::ZVector(n);
98  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
99  for (int i=0; i<IDELEMS(polynomialIdeal); i++)
100  {
101  g = polynomialIdeal->m[i];
102  if (g)
103  {
104  p_GetExpV(g,leadexpv,r);
105  leadexpw = expvToZVector(n, leadexpv);
106  pIter(g);
107  while (g)
108  {
109  p_GetExpV(g,tailexpv,r);
110  tailexpw = expvToZVector(n, tailexpv);
111  inequalities.appendRow(leadexpw-tailexpw);
112  pIter(g);
113  }
114  }
115  }
116  omFreeSize(leadexpv,(n+1)*sizeof(int));
117  omFreeSize(tailexpv,(n+1)*sizeof(int));
118  // if (currentStrategy->restrictToLowerHalfSpace())
119  // {
120  // gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
121  // lowerHalfSpaceCondition[0] = -1;
122  // inequalities.appendRow(lowerHalfSpaceCondition);
123  // }
124 
125  polyhedralCone = gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
126  polyhedralCone.canonicalize();
127  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
129 }
130 
131 groebnerCone::groebnerCone(const ideal I, const ring r, const gfan::ZVector& w, const tropicalStrategy& currentCase):
132  polynomialIdeal(NULL),
133  polynomialRing(NULL),
134  currentStrategy(&currentCase)
135 {
137  if (r) polynomialRing = rCopy(r);
138  if (I)
139  {
140  polynomialIdeal = id_Copy(I,r);
143  }
144 
145  int n = rVar(polynomialRing);
146  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
147  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
148  int* expv = (int*) omAlloc((n+1)*sizeof(int));
149  for (int i=0; i<IDELEMS(polynomialIdeal); i++)
150  {
151  poly g = polynomialIdeal->m[i];
152  if (g)
153  {
154  p_GetExpV(g,expv,polynomialRing);
155  gfan::ZVector leadexpv = intStar2ZVector(n,expv);
156  long d = wDeg(g,polynomialRing,w);
157  for (pIter(g); g; pIter(g))
158  {
159  p_GetExpV(g,expv,polynomialRing);
160  gfan::ZVector tailexpv = intStar2ZVector(n,expv);
161  if (wDeg(g,polynomialRing,w)==d)
162  equations.appendRow(leadexpv-tailexpv);
163  else
164  {
166  inequalities.appendRow(leadexpv-tailexpv);
167  }
168  }
169  }
170  }
171  omFreeSize(expv,(n+1)*sizeof(int));
172  // if (currentStrategy->restrictToLowerHalfSpace())
173  // {
174  // gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
175  // lowerHalfSpaceCondition[0] = -1;
176  // inequalities.appendRow(lowerHalfSpaceCondition);
177  // }
178 
179  polyhedralCone = gfan::ZCone(inequalities,equations);
180  polyhedralCone.canonicalize();
181  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
183 }
184 
185 /***
186  * Computes the groebner cone of I around u+e*w for e>0 sufficiently small.
187  * Assumes that this cone is a face of the maximal Groenbner cone given by the ordering of r.
188  **/
189 groebnerCone::groebnerCone(const ideal I, const ring r, const gfan::ZVector& u, const gfan::ZVector& w, const tropicalStrategy& currentCase):
190  polynomialIdeal(NULL),
191  polynomialRing(NULL),
192  currentStrategy(&currentCase)
193 {
194  assume(checkWeightVector(I,r,u));
196  if (r) polynomialRing = rCopy(r);
197  if (I)
198  {
199  polynomialIdeal = id_Copy(I,r);
202  }
203 
204  int n = rVar(polynomialRing);
205  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
206  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
207  int* expv = (int*) omAlloc((n+1)*sizeof(int));
208  for (int i=0; i<IDELEMS(polynomialIdeal); i++)
209  {
210  poly g = polynomialIdeal->m[i];
211  if (g)
212  {
213  p_GetExpV(g,expv,polynomialRing);
214  gfan::ZVector leadexpv = intStar2ZVector(n,expv);
215  long d1 = wDeg(g,polynomialRing,u);
216  long d2 = wDeg(g,polynomialRing,w);
217  for (pIter(g); g; pIter(g))
218  {
219  p_GetExpV(g,expv,polynomialRing);
220  gfan::ZVector tailexpv = intStar2ZVector(n,expv);
221  if (wDeg(g,polynomialRing,u)==d1 && wDeg(g,polynomialRing,w)==d2)
222  equations.appendRow(leadexpv-tailexpv);
223  else
224  {
226  inequalities.appendRow(leadexpv-tailexpv);
227  }
228  }
229  }
230  }
231  omFreeSize(expv,(n+1)*sizeof(int));
232  // if (currentStrategy->restrictToLowerHalfSpace())
233  // {
234  // gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
235  // lowerHalfSpaceCondition[0] = -1;
236  // inequalities.appendRow(lowerHalfSpaceCondition);
237  // }
238 
239  polyhedralCone = gfan::ZCone(inequalities,equations);
240  polyhedralCone.canonicalize();
241  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
243 }
244 
245 
246 groebnerCone::groebnerCone(const ideal I, const ideal inI, const ring r, const tropicalStrategy& currentCase):
247  polynomialIdeal(id_Copy(I,r)),
248  polynomialRing(rCopy(r)),
249  currentStrategy(&currentCase)
250 {
253 
256 
257  int n = rVar(r);
258  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
259  int* expv = (int*) omAlloc((n+1)*sizeof(int));
260  for (int i=0; i<IDELEMS(inI); i++)
261  {
262  poly g = inI->m[i];
263  if (g)
264  {
265  p_GetExpV(g,expv,r);
266  gfan::ZVector leadexpv = intStar2ZVector(n,expv);
267  for (pIter(g); g; pIter(g))
268  {
269  p_GetExpV(g,expv,r);
270  gfan::ZVector tailexpv = intStar2ZVector(n,expv);
271  equations.appendRow(leadexpv-tailexpv);
272  }
273  }
274  }
275  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
276  for (int i=0; i<IDELEMS(polynomialIdeal); i++)
277  {
278  poly g = polynomialIdeal->m[i];
279  if (g)
280  {
281  p_GetExpV(g,expv,r);
282  gfan::ZVector leadexpv = intStar2ZVector(n,expv);
283  for (pIter(g); g; pIter(g))
284  {
285  p_GetExpV(g,expv,r);
286  gfan::ZVector tailexpv = intStar2ZVector(n,expv);
287  inequalities.appendRow(leadexpv-tailexpv);
288  }
289  }
290  }
291  omFreeSize(expv,(n+1)*sizeof(int));
293  {
294  gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
295  lowerHalfSpaceCondition[0] = -1;
296  inequalities.appendRow(lowerHalfSpaceCondition);
297  }
298 
299  polyhedralCone = gfan::ZCone(inequalities,equations);
300  polyhedralCone.canonicalize();
301  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
303 }
304 
306  polynomialIdeal(NULL),
307  polynomialRing(NULL),
308  polyhedralCone(gfan::ZCone(sigma.getPolyhedralCone())),
309  interiorPoint(gfan::ZVector(sigma.getInteriorPoint())),
310  currentStrategy(sigma.getTropicalStrategy())
311 {
317 }
318 
320 {
326 }
327 
329 {
338  return *this;
339 }
340 
341 /**
342  * Returns true if Groebner cone contains w, false otherwise
343  */
344 bool groebnerCone::contains(const gfan::ZVector &w) const
345 {
346  return polyhedralCone.contains(w);
347 }
348 
349 
350 /***
351  * Returns a point in the tropical variety, if the groebnerCone contains one.
352  * Returns an empty vector otherwise.
353  **/
354 gfan::ZVector groebnerCone::tropicalPoint() const
355 {
357  ideal I = polynomialIdeal;
358  ring r = polynomialRing;
359 
360  gfan::ZCone coneToCheck = polyhedralCone;
361  gfan::ZMatrix R = coneToCheck.extremeRays();
362  for (int i=0; i<R.getHeight(); i++)
363  {
364  assume(!currentStrategy->restrictToLowerHalfSpace() || R[i][0].sign()<=0);
365  if (currentStrategy->restrictToLowerHalfSpace() && R[i][0].sign()==0)
366  continue;
367  std::pair<poly,int> s = currentStrategy->checkInitialIdealForMonomial(I,r,R[i]);
368  if (s.first==NULL)
369  {
370  if (s.second<0)
371  // if monomial was initialized, delete it
372  p_Delete(&s.first,r);
373  return R[i];
374  }
375  }
376  return gfan::ZVector();
377 }
378 
379 /**
380  * Given an interior point on the facet and the outer normal factor on the facet,
381  * returns the adjacent groebnerCone sharing that facet
382  */
383 groebnerCone groebnerCone::flipCone(const gfan::ZVector &interiorPoint, const gfan::ZVector &facetNormal) const
384 {
385  assume(this->checkFlipConeInput(interiorPoint,facetNormal));
387  /* Note: the polynomial ring created will have a weighted ordering with respect to interiorPoint
388  * and with a weighted ordering with respect to facetNormal as tiebreaker.
389  * Hence it is sufficient to compute the initial form with respect to facetNormal,
390  * to obtain an initial form with respect to interiorPoint+e*facetNormal,
391  * for e>0 sufficiently small */
392  std::pair<ideal,ring> flipped = currentStrategy->computeFlip(polynomialIdeal,polynomialRing,interiorPoint,facetNormal);
393  assume(checkPolynomialInput(flipped.first,flipped.second));
394  groebnerCone flippedCone(flipped.first, flipped.second, interiorPoint, facetNormal, *currentStrategy);
395  id_Delete(&flipped.first,flipped.second);
396  rDelete(flipped.second);
397  return flippedCone;
398 }
399 
400 
401 /***
402  * Returns a complete list of neighboring Groebner cones.
403  **/
405 {
406  std::pair<gfan::ZMatrix, gfan::ZMatrix> facetsData = interiorPointsAndNormalsOfFacets(polyhedralCone);
407 
408  gfan::ZMatrix interiorPoints = facetsData.first;
409  gfan::ZMatrix facetNormals = facetsData.second;
410 
411  groebnerCones neighbours;
412  for (int i=0; i<interiorPoints.getHeight(); i++)
413  {
414  gfan::ZVector w = interiorPoints[i];
415  gfan::ZVector v = facetNormals[i];
417  {
418  assume(w[0].sign()<=0);
419  if (w[0].sign()==0 && v[0].sign()>0)
420  continue;
421  }
422  neighbours.insert(flipCone(w,v));
423  }
424  return neighbours;
425 }
426 
427 
428 bool groebnerCone::pointsOutwards(const gfan::ZVector w) const
429 {
430  gfan::ZCone dual = polyhedralCone.dualCone();
431  return (!dual.contains(w));
432 }
433 
434 
435 /***
436  * Returns a complete list of neighboring Groebner cones in the tropical variety.
437  **/
439 {
440  gfan::ZMatrix interiorPoints = interiorPointsOfFacets(polyhedralCone);
441  groebnerCones neighbours;
442  for (int i=0; i<interiorPoints.getHeight(); i++)
443  {
444  if (!(currentStrategy->restrictToLowerHalfSpace() && interiorPoints[i][0].sign()==0))
445  {
446  ideal initialIdeal = initial(polynomialIdeal,polynomialRing,interiorPoints[i]);
447  gfan::ZMatrix ray = raysOfTropicalStar(initialIdeal,polynomialRing,interiorPoints[i],currentStrategy);
448  for (int j=0; j<ray.getHeight(); j++)
449  if (pointsOutwards(ray[j]))
450  {
451  groebnerCone neighbour = flipCone(interiorPoints[i],ray[j]);
452  neighbours.insert(neighbour);
453  }
454  id_Delete(&initialIdeal,polynomialRing);
455  }
456  }
457  return neighbours;
458 }
459 
460 
461 gfan::ZFan* toFanStar(groebnerCones setOfCones)
462 {
463  if (setOfCones.size() > 0)
464  {
465  groebnerCones::iterator sigma = setOfCones.begin();
466  gfan::ZFan* zf = new gfan::ZFan(sigma->getPolyhedralCone().ambientDimension());
467  for (; sigma!=setOfCones.end(); sigma++)
468  {
469  gfan::ZCone zc = sigma->getPolyhedralCone();
470  // assume(isCompatible(zf,&zc));
471  zf->insert(zc);
472  }
473  return zf;
474  }
475  else
476  return new gfan::ZFan(gfan::ZFan(currRing->N));
477 }
gfan::ZMatrix interiorPointsOfFacets(const gfan::ZCone &zc, const std::set< gfan::ZVector > &exceptThese)
Definition: bbcone.cc:1848
BOOLEAN equations(leftv res, leftv args)
Definition: bbcone.cc:577
BOOLEAN inequalities(leftv res, leftv args)
Definition: bbcone.cc:560
std::string toString(const gfan::ZCone *const c)
Definition: bbcone.cc:27
std::pair< gfan::ZMatrix, gfan::ZMatrix > interiorPointsAndNormalsOfFacets(const gfan::ZCone zc, const std::set< gfan::ZVector > &exceptThesePoints, const bool onlyLowerHalfSpace)
Definition: bbcone.cc:1902
gfan::ZVector intStar2ZVector(const int d, const int *i)
gfan::ZVector expvToZVector(const int n, const int *expv)
int i
Definition: cfEzgcd.cc:125
g
Definition: cfModGcd.cc:4031
const tropicalStrategy * currentStrategy
Definition: groebnerCone.h:41
gfan::ZVector tropicalPoint() const
Returns a point in the tropical variety, if the groebnerCone contains one.
groebnerCones tropicalNeighbours() const
Returns a complete list of neighboring Groebner cones in the tropical variety.
groebnerCones groebnerNeighbours() const
Returns a complete list of neighboring Groebner cones.
bool pointsOutwards(const gfan::ZVector) const
groebnerCone & operator=(const groebnerCone &sigma)
bool contains(const gfan::ZVector &w) const
Returns true if Groebner cone contains w, false otherwise.
gfan::ZVector interiorPoint
Definition: groebnerCone.h:40
gfan::ZCone getPolyhedralCone() const
Definition: groebnerCone.h:64
bool checkFlipConeInput(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const
Debug tools.
Definition: groebnerCone.cc:23
ideal polynomialIdeal
ideal to which this Groebner cone belongs to
Definition: groebnerCone.h:34
ideal getPolynomialIdeal() const
Definition: groebnerCone.h:62
gfan::ZCone polyhedralCone
Definition: groebnerCone.h:39
const tropicalStrategy * getTropicalStrategy() const
Definition: groebnerCone.h:66
ring getPolynomialRing() const
Definition: groebnerCone.h:63
gfan::ZVector getInteriorPoint() const
Definition: groebnerCone.h:65
groebnerCone flipCone(const gfan::ZVector &interiorPoint, const gfan::ZVector &facetNormal) const
Given an interior point on the facet and the outer normal factor on the facet, returns the adjacent g...
ring polynomialRing
ring in which the ideal exists
Definition: groebnerCone.h:38
std::pair< ideal, ring > computeFlip(const ideal Ir, const ring r, const gfan::ZVector &interiorPoint, const gfan::ZVector &facetNormal) const
given an interior point of a groebner cone computes the groebner cone adjacent to it
bool reduce(ideal I, const ring r) const
reduces the generators of an ideal I so that the inequalities and equations of the Groebner cone can ...
void pReduce(ideal I, const ring r) const
bool restrictToLowerHalfSpace() const
returns true, if valuation non-trivial, false otherwise
std::pair< poly, int > checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector &w=0) const
If given w, assuming w is in the Groebner cone of the ordering on r and I is a standard basis with re...
const CanonicalForm int s
Definition: facAbsFact.cc:55
const CanonicalForm & w
Definition: facAbsFact.cc:55
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
int j
Definition: facHensel.cc:105
gfan::ZFan * toFanStar(groebnerCones setOfCones)
implementation of the class groebnerCone
std::set< groebnerCone, groebnerCone_compare > groebnerCones
Definition: groebnerCone.h:23
ideal id_Copy(ideal h1, const ring r)
copy an ideal
long wDeg(const poly p, const ring r, const gfan::ZVector &w)
various functions to compute the initial form of polynomials and ideals
Definition: initial.cc:6
poly initial(const poly p, const ring r, const gfan::ZVector &w)
Returns the initial form of p with respect to w.
Definition: initial.cc:30
#define assume(x)
Definition: mod2.h:390
#define pIter(p)
Definition: monomials.h:44
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define NULL
Definition: omList.c:10
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:857
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1466
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:439
static int sign(int x)
Definition: ring.cc:3328
ring rCopy(ring r)
Definition: ring.cc:1605
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:583
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
#define IDELEMS(i)
Definition: simpleideals.h:24
#define R
Definition: sirandom.c:26
gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector &u, const tropicalStrategy *currentStrategy)
bool checkWeightVector(const ideal I, const ring r, const gfan::ZVector &weightVector, bool checkBorder)
bool checkPolyhedralInput(const gfan::ZCone zc, const gfan::ZVector p)
bool checkOrderingAndCone(const ring r, const gfan::ZCone zc)
bool checkPolynomialInput(const ideal I, const ring r)
implementation of the class tropicalStrategy