tropicalCurves.cc
Go to the documentation of this file.
1 #include <gfanlib/gfanlib_matrix.h>
2 #include <gfanlib/gfanlib_zcone.h>
5 #include <std_wrapper.h>
6 #include <containsMonomial.h>
7 #include <initial.h>
8 #include <witness.h>
9 #include <tropicalStrategy.h>
11 #include <tropicalCurves.h>
12 #include <set>
13 #ifndef NDEBUG
14 #include <bbfan.h>
15 #endif
16 
17 /***
18  * Given two sets of cones A,B and a dimensional bound d,
19  * computes the intersections of all cones of A with all cones of B,
20  * and throws away those of lower dimension than d.
21  **/
23  const ZConesSortedByDimension &setB,
24  int d=0)
25 {
26  if (setA.empty())
27  return setB;
28  if (setB.empty())
29  return setA;
31  for (ZConesSortedByDimension::iterator coneOfA=setA.begin(); coneOfA!=setA.end(); coneOfA++)
32  {
33  for (ZConesSortedByDimension::iterator coneOfB=setB.begin(); coneOfB!=setB.end(); coneOfB++)
34  {
35  gfan::ZCone coneOfIntersection = gfan::intersection(*coneOfA,*coneOfB);
36  if (coneOfIntersection.dimension()>=d)
37  {
38  coneOfIntersection.canonicalize();
39  setAB.insert(coneOfIntersection);
40  }
41  }
42  }
43  return setAB;
44 }
45 
46 /***
47  * Given a ring r, weights u, w, and a matrix E, returns a copy of r whose ordering is,
48  * for any ideal homogeneous with respect to u, weighted with respect to u and
49  * whose tiebreaker is genericly weighted with respect to v and E in the following sense:
50  * the ordering "lies" on the affine space A running through v and spanned by the row vectors of E,
51  * and it lies in a Groebner cone of dimension at least rank(E)=dim(A).
52  **/
53 static ring genericlyWeightedOrdering(const ring r, const gfan::ZVector &u, const gfan::ZVector &w,
54  const gfan::ZMatrix &W, const tropicalStrategy* currentStrategy)
55 {
56  int n = rVar(r);
57  int h = W.getHeight();
58 
59  /* create a copy s of r and delete its ordering */
60  ring s = rCopy0(r);
61  omFree(s->order);
62  s->order = (int*) omAlloc0((h+4)*sizeof(int));
63  omFree(s->block0);
64  s->block0 = (int*) omAlloc0((h+4)*sizeof(int));
65  omFree(s->block1);
66  s->block1 = (int*) omAlloc0((h+4)*sizeof(int));
67  for (int j=0; s->wvhdl[j]; j++) omFree(s->wvhdl[j]);
68  omFree(s->wvhdl);
69  s->wvhdl = (int**) omAlloc0((h+4)*sizeof(int*));
70 
71  /* construct a new ordering as describe above */
72  bool overflow = false;
73  s->order[0] = ringorder_a;
74  s->block0[0] = 1;
75  s->block1[0] = n;
76  gfan::ZVector uAdjusted = currentStrategy->adjustWeightForHomogeneity(u);
77  s->wvhdl[0] = ZVectorToIntStar(uAdjusted,overflow);
78  s->order[1] = ringorder_a;
79  s->block0[1] = 1;
80  s->block1[1] = n;
81  gfan::ZVector wAdjusted = currentStrategy->adjustWeightUnderHomogeneity(w,uAdjusted);
82  s->wvhdl[1] = ZVectorToIntStar(wAdjusted,overflow);
83  for (int j=0; j<h-1; j++)
84  {
85  s->order[j+2] = ringorder_a;
86  s->block0[j+2] = 1;
87  s->block1[j+2] = n;
88  wAdjusted = currentStrategy->adjustWeightUnderHomogeneity(W[j],uAdjusted);
89  s->wvhdl[j+2] = ZVectorToIntStar(wAdjusted,overflow);
90  }
91  s->order[h+1] = ringorder_wp;
92  s->block0[h+1] = 1;
93  s->block1[h+1] = n;
94  wAdjusted = currentStrategy->adjustWeightUnderHomogeneity(W[h-1],uAdjusted);
95  s->wvhdl[h+1] = ZVectorToIntStar(wAdjusted,overflow);
96  s->order[h+2] = ringorder_C;
97 
98  if (overflow)
99  {
100  WerrorS("genericlyWeightedOrdering: overflow in weight vector");
101  throw 0; // weightOverflow;
102  }
103 
104  /* complete the ring and return it */
105  rComplete(s);
106  rTest(s);
107  return s;
108 }
109 
110 
111 /***
112  * Let I be an ideal. Given a weight vector u in the relative interior
113  * of a one-codimensional cone of the tropical variety of I and
114  * the initial ideal inI with respect to it, computes the star of the tropical variety in u.
115  **/
116 ZConesSortedByDimension tropicalStar(ideal inI, const ring r, const gfan::ZVector &u,
117  const tropicalStrategy* currentStrategy)
118 {
119  int k = idSize(inI);
120  int d = currentStrategy->getExpectedDimension();
121 
122  /* Compute the common refinement over all tropical varieties
123  * of the polynomials in the generating set */
124  ZConesSortedByDimension C = tropicalVarietySortedByDimension(inI->m[0],r,currentStrategy);
125  for (int i=1; i<k; i++)
126  C = intersect(C,tropicalVarietySortedByDimension(inI->m[i],r,currentStrategy),d);
127 
128  /* Cycle through all maximal cones of the refinement.
129  * Pick a monomial ordering corresponding to a generic weight vector in it
130  * and check if the initial ideal is monomial free, generic meaning
131  * that it lies in a maximal Groebner cone in the maximal cone of the refinement.
132  * If the initial ideal is not monomial free, compute a witness for the monomial
133  * and compute the common refinement with its tropical variety.
134  * If all initial ideals are monomial free, then we have our tropical curve */
135  // gfan::ZFan* zf = toFanStar(C);
136  // std::cout << zf->toString(2+4+8+128) << std::endl;
137  // delete zf;
138  for (std::set<gfan::ZCone>::iterator zc=C.begin(); zc!=C.end();)
139  {
140  gfan::ZVector w = zc->getRelativeInteriorPoint();
141  gfan::ZMatrix W = zc->generatorsOfSpan();
142  // std::cout << zc->extremeRays() << std::endl;
143 
144  ring s = genericlyWeightedOrdering(r,u,w,W,currentStrategy);
145  nMapFunc identity = n_SetMap(r->cf,s->cf);
146  ideal inIs = idInit(k);
147  for (int j=0; j<k; j++)
148  inIs->m[j] = p_PermPoly(inI->m[j],NULL,r,s,identity,NULL,0);
149 
150  ideal inIsSTD = gfanlib_kStd_wrapper(inIs,s,isHomog);
151  id_Delete(&inIs,s);
152  ideal ininIs = initial(inIsSTD,s,w,W);
153 
154  std::pair<poly,int> mons = currentStrategy->checkInitialIdealForMonomial(ininIs,s);
155 
156  if (mons.first!=NULL)
157  {
158  poly gs;
159  if (mons.second>=0)
160  // cheap way out, ininIsSTD already contains a monomial in its generators
161  gs = inIsSTD->m[mons.second];
162  else
163  // compute witness
164  gs = witness(mons.first,inIsSTD,ininIs,s);
165 
166  C = intersect(C,tropicalVarietySortedByDimension(gs,s,currentStrategy),d);
167  nMapFunc mMap = n_SetMap(s->cf,r->cf);
168  poly gr = p_PermPoly(gs,NULL,s,r,mMap,NULL,0);
169  idInsertPoly(inI,gr);
170  k++;
171 
172  if (mons.second<0)
173  {
174  // if necessary, cleanup mons and gs
175  p_Delete(&mons.first,s);
176  p_Delete(&gs,s);
177  }
178  // cleanup rest, reset zc
179  id_Delete(&inIsSTD,s);
180  id_Delete(&ininIs,s);
181  rDelete(s);
182  zc = C.begin();
183  }
184  else
185  {
186  // cleanup remaining data of first stage
187  id_Delete(&inIsSTD,s);
188  id_Delete(&ininIs,s);
189  rDelete(s);
190 
191  gfan::ZVector wNeg = currentStrategy->negateWeight(w);
192  if (zc->contains(wNeg))
193  {
194  s = genericlyWeightedOrdering(r,u,wNeg,W,currentStrategy);
195  identity = n_SetMap(r->cf,s->cf);
196  inIs = idInit(k);
197  for (int j=0; j<k; j++)
198  inIs->m[j] = p_PermPoly(inI->m[j],NULL,r,s,identity,NULL,0);
199 
200  inIsSTD = gfanlib_kStd_wrapper(inIs,s,isHomog);
201  id_Delete(&inIs,s);
202  ininIs = initial(inIsSTD,s,wNeg,W);
203 
204  mons = currentStrategy->checkInitialIdealForMonomial(ininIs,s);
205  if (mons.first!=NULL)
206  {
207  poly gs;
208  if (mons.second>=0)
209  // cheap way out, ininIsSTD already contains a monomial in its generators
210  gs = inIsSTD->m[mons.second];
211  else
212  // compute witness
213  gs = witness(mons.first,inIsSTD,ininIs,s);
214 
215  C = intersect(C,tropicalVarietySortedByDimension(gs,s,currentStrategy),d);
216  nMapFunc mMap = n_SetMap(s->cf,r->cf);
217  poly gr = p_PermPoly(gs,NULL,s,r,mMap,NULL,0);
218  idInsertPoly(inI,gr);
219  k++;
220 
221  if (mons.second<0)
222  {
223  // if necessary, cleanup mons and gs
224  p_Delete(&mons.first,s);
225  p_Delete(&gs,s);
226  }
227  // reset zc
228  zc = C.begin();
229  }
230  else
231  zc++;
232  // cleanup remaining data of second stage
233  id_Delete(&inIsSTD,s);
234  id_Delete(&ininIs,s);
235  rDelete(s);
236  }
237  else
238  zc++;
239  }
240  }
241  return C;
242 }
243 
244 
245 gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector &u, const tropicalStrategy* currentStrategy)
246 {
247  ZConesSortedByDimension C = tropicalStar(I,r,u,currentStrategy);
248  // gfan::ZFan* zf = toFanStar(C);
249  // std::cout << zf->toString();
250  // delete zf;
251  gfan::ZMatrix raysOfC(0,u.size());
252  if (!currentStrategy->restrictToLowerHalfSpace())
253  {
254  for (ZConesSortedByDimension::iterator zc=C.begin(); zc!=C.end(); zc++)
255  {
256  assume(zc->dimensionOfLinealitySpace()+1 >= zc->dimension());
257  if (zc->dimensionOfLinealitySpace()+1 >= zc->dimension())
258  raysOfC.appendRow(zc->getRelativeInteriorPoint());
259  else
260  {
261  gfan::ZVector interiorPoint = zc->getRelativeInteriorPoint();
262  if (!currentStrategy->homogeneitySpaceContains(interiorPoint))
263  {
264  raysOfC.appendRow(interiorPoint);
265  raysOfC.appendRow(currentStrategy->negateWeight(interiorPoint));
266  }
267  else
268  {
269  gfan::ZMatrix zm = zc->generatorsOfLinealitySpace();
270  for (int i=0; i<zm.getHeight(); i++)
271  {
272  gfan::ZVector point = zm[i];
273  if (currentStrategy->homogeneitySpaceContains(point))
274  {
275  raysOfC.appendRow(point);
276  raysOfC.appendRow(currentStrategy->negateWeight(point));
277  break;
278  }
279  }
280  }
281  }
282  }
283  }
284  else
285  {
286  for (ZConesSortedByDimension::iterator zc=C.begin(); zc!=C.end(); zc++)
287  {
288  assume(zc->dimensionOfLinealitySpace()+2 >= zc->dimension());
289  if (zc->dimensionOfLinealitySpace()+2 == zc->dimension())
290  raysOfC.appendRow(zc->getRelativeInteriorPoint());
291  else
292  {
293  gfan::ZVector interiorPoint = zc->getRelativeInteriorPoint();
294  if (!currentStrategy->homogeneitySpaceContains(interiorPoint))
295  {
296  raysOfC.appendRow(interiorPoint);
297  raysOfC.appendRow(currentStrategy->negateWeight(interiorPoint));
298  }
299  else
300  {
301  gfan::ZMatrix zm = zc->generatorsOfLinealitySpace();
302  for (int i=0; i<zm.getHeight(); i++)
303  {
304  gfan::ZVector point = zm[i];
305  if (currentStrategy->homogeneitySpaceContains(point))
306  {
307  raysOfC.appendRow(point);
308  raysOfC.appendRow(currentStrategy->negateWeight(point));
309  break;
310  }
311  }
312  }
313  }
314  }
315  }
316  return raysOfC;
317 }
318 
319 
320 /***
321  * Computes the tropical curve of an x-homogeneous ideal I
322  * which is weighted homogeneous with respect to weight w in ring r
323  **/
324 #ifndef NDEBUG
326 {
327  leftv u = args;
328  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
329  {
330  leftv v = u->next;
331  if ((v!=NULL) && (v->Typ()==BIGINTMAT_CMD))
332  {
333  omUpdateInfo();
334  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
335  ideal inI = (ideal) u->CopyD();
336  bigintmat* u = (bigintmat*) v->CopyD();
337  tropicalStrategy currentCase(inI,currRing);
338  gfan::ZVector* v = bigintmatToZVector(u);
339  ZConesSortedByDimension C = tropicalStar(inI,currRing,*v,&currentCase);
340  id_Delete(&inI,currRing);
341  delete u;
342  delete v;
343  res->rtyp = NONE;
344  res->data = NULL;
345  // res->rtyp = fanID;
346  // res->data = (char*) toFanStar(C);
347  return FALSE;
348  }
349  }
350  WerrorS("tropicalStarDebug: unexpected parameters");
351  return TRUE;
352 }
353 #endif
implementation of the class tropicalStrategy
const CanonicalForm int s
Definition: facAbsFact.cc:55
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
#define Print
Definition: emacs.cc:83
static ring genericlyWeightedOrdering(const ring r, const gfan::ZVector &u, const gfan::ZVector &w, const gfan::ZMatrix &W, const tropicalStrategy *currentStrategy)
#define FALSE
Definition: auxiliary.h:140
Matrices of numbers.
Definition: bigintmat.h:51
ZConesSortedByDimension tropicalVarietySortedByDimension(const poly g, const ring r, const tropicalStrategy *currentCase)
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:540
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
#define TRUE
Definition: auxiliary.h:144
int * ZVectorToIntStar(const gfan::ZVector &v, bool &overflow)
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
gfan::ZVector negateWeight(const gfan::ZVector &w) const
int Typ()
Definition: subexpr.cc:969
int getExpectedDimension() const
returns the expected Dimension of the polyhedral output
void * data
Definition: subexpr.h:89
poly res
Definition: myNF.cc:322
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
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
const ring r
Definition: syzextra.cc:208
std::set< gfan::ZCone, ZConeCompareDimensionFirst > ZConesSortedByDimension
BOOLEAN rComplete(ring r, int force)
this needs to be called whenever a new ring is created: new fields in ring are created (like VarOffse...
Definition: ring.cc:3436
int j
Definition: myNF.cc:70
ideal gfanlib_kStd_wrapper(ideal I, ring r, tHomog h=testHomog)
Definition: std_wrapper.cc:5
#define omFree(addr)
Definition: omAllocDecl.h:261
omInfo_t om_Info
Definition: omStats.c:13
#define assume(x)
Definition: mod2.h:405
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1318
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:72
#define rTest(r)
Definition: ring.h:781
BOOLEAN idInsertPoly(ideal h1, poly h2)
insert h2 into h1 (if h2 is not the zero polynomial) return TRUE iff h2 was indeed inserted ...
bool homogeneitySpaceContains(const gfan::ZVector &v) const
returns true, if v is contained in the homogeneity space; false otherwise
static ZConesSortedByDimension intersect(const ZConesSortedByDimension &setA, const ZConesSortedByDimension &setB, int d=0)
int i
Definition: cfEzgcd.cc:123
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition: coeffs.h:720
gfan::ZVector adjustWeightUnderHomogeneity(gfan::ZVector v, gfan::ZVector w) const
Given strictly positive weight w and weight v, returns a strictly positive weight u such that on an i...
leftv next
Definition: subexpr.h:87
bool restrictToLowerHalfSpace() const
returns true, if valuation non-trivial, false otherwise
ZConesSortedByDimension tropicalStar(ideal inI, const ring r, const gfan::ZVector &u, const tropicalStrategy *currentStrategy)
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:850
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
poly p_PermPoly(poly p, const int *perm, const ring oldRing, const ring dst, nMapFunc nMap, const int *par_perm, int OldPar)
Definition: p_polys.cc:3926
gfan::ZVector adjustWeightForHomogeneity(gfan::ZVector w) const
Given weight w, returns a strictly positive weight u such that an ideal satisfying the valuation-sepc...
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
#define NULL
Definition: omList.c:10
poly witness(const poly m, const ideal I, const ideal inI, const ring r)
Let w be the uppermost weight vector in the matrix defining the ordering on r.
Definition: witness.cc:34
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:448
void omUpdateInfo()
Definition: omStats.c:24
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 & w
Definition: facAbsFact.cc:55
int rtyp
Definition: subexpr.h:92
static int idSize(const ideal id)
Count the effective size of an ideal (without the trailing allocated zero-elements) ...
Definition: ideals.h:43
BOOLEAN tropicalStarDebug(leftv res, leftv args)
polyrec * poly
Definition: hilb.h:10
static Poly * h
Definition: janet.cc:978
int BOOLEAN
Definition: auxiliary.h:131
#define NONE
Definition: tok.h:173
gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector &u, const tropicalStrategy *currentStrategy)
void * CopyD(int t)
Definition: subexpr.cc:676
gfan::ZVector * bigintmatToZVector(const bigintmat &bim)
#define omAlloc0(size)
Definition: omAllocDecl.h:211