My Project
Loading...
Searching...
No Matches
groebnerCone.cc
Go to the documentation of this file.
1
3#include "kernel/ideals.h"
4#include "Singular/ipid.h"
5
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"
17#include "containsMonomial.h"
18#include "tropicalCurves.h"
19#include "bbcone.h"
20#include "tropicalDebug.h"
21
22#ifndef NDEBUG
23bool 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
72 polyhedralCone(gfan::ZCone(0)),
73 interiorPoint(gfan::ZVector(0)),
75{
76}
77
78groebnerCone::groebnerCone(const ideal I, const ring r, const tropicalStrategy& currentCase):
81 currentStrategy(&currentCase)
82{
84 if (r) polynomialRing = rCopy(r);
85 if (I)
86 {
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
131groebnerCone::groebnerCone(const ideal I, const ring r, const gfan::ZVector& w, const tropicalStrategy& currentCase):
134 currentStrategy(&currentCase)
135{
137 if (r) polynomialRing = rCopy(r);
138 if (I)
139 {
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 {
155 gfan::ZVector leadexpv = intStar2ZVector(n,expv);
156 long d = wDeg(g,polynomialRing,w);
157 for (pIter(g); g; pIter(g))
158 {
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
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 **/
189groebnerCone::groebnerCone(const ideal I, const ring r, const gfan::ZVector& u, const gfan::ZVector& w, const tropicalStrategy& currentCase):
192 currentStrategy(&currentCase)
193{
196 if (r) polynomialRing = rCopy(r);
197 if (I)
198 {
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 {
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 {
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
240 polyhedralCone.canonicalize();
241 interiorPoint = polyhedralCone.getRelativeInteriorPoint();
243}
244
245
246groebnerCone::groebnerCone(const ideal I, const ideal inI, const ring r, const tropicalStrategy& currentCase):
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));
292 if (currentStrategy->restrictToLowerHalfSpace())
293 {
294 gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
295 lowerHalfSpaceCondition[0] = -1;
296 inequalities.appendRow(lowerHalfSpaceCondition);
297 }
298
300 polyhedralCone.canonicalize();
301 interiorPoint = polyhedralCone.getRelativeInteriorPoint();
303}
304
318
327
340
341/**
342 * Returns true if Groebner cone contains w, false otherwise
343 */
344bool 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 **/
354gfan::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 */
383groebnerCone 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];
416 if (currentStrategy->restrictToLowerHalfSpace())
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
428bool 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
461gfan::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}
std::pair< gfan::ZMatrix, gfan::ZMatrix > interiorPointsAndNormalsOfFacets(const gfan::ZCone zc, const std::set< gfan::ZVector > &exceptThesePoints, const bool onlyLowerHalfSpace)
Definition bbcone.cc:1853
gfan::ZMatrix interiorPointsOfFacets(const gfan::ZCone &zc, const std::set< gfan::ZVector > &exceptThese)
Definition bbcone.cc:1799
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
gfan::ZVector intStar2ZVector(const int d, const int *i)
gfan::ZVector expvToZVector(const int n, const int *expv)
int i
Definition cfEzgcd.cc:132
g
Definition cfModGcd.cc:4098
const tropicalStrategy * currentStrategy
const tropicalStrategy * getTropicalStrategy() const
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.
groebnerCone & operator=(const groebnerCone &sigma)
bool contains(const gfan::ZVector &w) const
Returns true if Groebner cone contains w, false otherwise.
gfan::ZVector interiorPoint
gfan::ZCone getPolyhedralCone() const
bool checkFlipConeInput(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const
Debug tools.
ideal polynomialIdeal
ideal to which this Groebner cone belongs to
bool pointsOutwards(const gfan::ZVector w) const
Return 1 if w points is in the dual of the polyhedral cone, 0 otherwise.
ideal getPolynomialIdeal() const
gfan::ZCone polyhedralCone
ring getPolynomialRing() const
gfan::ZVector getInteriorPoint() const
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
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
const CanonicalForm int s
Definition facAbsFact.cc:51
const CanonicalForm & w
Definition facAbsFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
int j
Definition facHensel.cc:110
gfan::ZFan * toFanStar(groebnerCones setOfCones)
implementation of the class groebnerCone
std::set< groebnerCone, groebnerCone_compare > groebnerCones
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:389
#define pIter(p)
Definition monomials.h:37
#define omFreeSize(addr, size)
#define omAlloc(size)
#define NULL
Definition omList.c:12
static void p_Delete(poly *p, const ring r)
Definition p_polys.h:903
static void p_GetExpV(poly p, int *ev, const ring r)
Definition p_polys.h:1536
VAR 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:454
static int sign(int x)
Definition ring.cc:3503
ring rCopy(ring r)
Definition ring.cc:1736
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition ring.h:598
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
#define IDELEMS(i)
#define R
Definition sirandom.c:27
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