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