My Project
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 
13 /***
14  * Given two sets of cones A,B and a dimensional bound d,
15  * computes the intersections of all cones of A with all cones of B,
16  * and throws away those of lower dimension than d.
17  **/
19  const ZConesSortedByDimension &setB,
20  int d=0)
21 {
22  if (setA.empty())
23  return setB;
24  if (setB.empty())
25  return setA;
27  for (ZConesSortedByDimension::iterator coneOfA=setA.begin(); coneOfA!=setA.end(); coneOfA++)
28  {
29  for (ZConesSortedByDimension::iterator coneOfB=setB.begin(); coneOfB!=setB.end(); coneOfB++)
30  {
31  gfan::ZCone coneOfIntersection = gfan::intersection(*coneOfA,*coneOfB);
32  if (coneOfIntersection.dimension()>=d)
33  {
34  coneOfIntersection.canonicalize();
35  setAB.insert(coneOfIntersection);
36  }
37  }
38  }
39  return setAB;
40 }
41 
42 /***
43  * Given a ring r, weights u, w, and a matrix E, returns a copy of r whose ordering is,
44  * for any ideal homogeneous with respect to u, weighted with respect to u and
45  * whose tiebreaker is genericly weighted with respect to v and E in the following sense:
46  * the ordering "lies" on the affine space A running through v and spanned by the row vectors of E,
47  * and it lies in a Groebner cone of dimension at least rank(E)=dim(A).
48  **/
49 static ring genericlyWeightedOrdering(const ring r, const gfan::ZVector &u, const gfan::ZVector &w,
50  const gfan::ZMatrix &W, const tropicalStrategy* currentStrategy)
51 {
52  int n = rVar(r);
53  int h = W.getHeight();
54 
55  /* create a copy s of r and delete its ordering */
56  ring s = rCopy0(r,FALSE,FALSE);
57  s->order = (rRingOrder_t*) omAlloc0((h+4)*sizeof(rRingOrder_t));
58  s->block0 = (int*) omAlloc0((h+4)*sizeof(int));
59  s->block1 = (int*) omAlloc0((h+4)*sizeof(int));
60  s->wvhdl = (int**) omAlloc0((h+4)*sizeof(int*));
61 
62  /* construct a new ordering as describe above */
63  bool overflow = false;
64  s->order[0] = ringorder_a;
65  s->block0[0] = 1;
66  s->block1[0] = n;
67  gfan::ZVector uAdjusted = currentStrategy->adjustWeightForHomogeneity(u);
68  s->wvhdl[0] = ZVectorToIntStar(uAdjusted,overflow);
69  s->order[1] = ringorder_a;
70  s->block0[1] = 1;
71  s->block1[1] = n;
72  gfan::ZVector wAdjusted = currentStrategy->adjustWeightUnderHomogeneity(w,uAdjusted);
73  s->wvhdl[1] = ZVectorToIntStar(wAdjusted,overflow);
74  for (int j=0; j<h-1; j++)
75  {
76  s->order[j+2] = ringorder_a;
77  s->block0[j+2] = 1;
78  s->block1[j+2] = n;
79  wAdjusted = currentStrategy->adjustWeightUnderHomogeneity(W[j],uAdjusted);
80  s->wvhdl[j+2] = ZVectorToIntStar(wAdjusted,overflow);
81  }
82  s->order[h+1] = ringorder_wp;
83  s->block0[h+1] = 1;
84  s->block1[h+1] = n;
85  wAdjusted = currentStrategy->adjustWeightUnderHomogeneity(W[h-1],uAdjusted);
86  s->wvhdl[h+1] = ZVectorToIntStar(wAdjusted,overflow);
87  s->order[h+2] = ringorder_C;
88 
89  if (overflow)
90  {
91  WerrorS("genericlyWeightedOrdering: overflow in weight vector");
92  throw 0; // weightOverflow;
93  }
94 
95  /* complete the ring and return it */
96  rComplete(s);
97  rTest(s);
98  return s;
99 }
100 
101 
102 /***
103  * Let I be an ideal. Given a weight vector u in the relative interior
104  * of a one-codimensional cone of the tropical variety of I and
105  * the initial ideal inI with respect to it, computes the star of the tropical variety in u.
106  **/
107 ZConesSortedByDimension tropicalStar(ideal inI, const ring r, const gfan::ZVector &u,
108  const tropicalStrategy* currentStrategy)
109 {
110  int k = IDELEMS(inI);
111  int d = currentStrategy->getExpectedDimension();
112 
113  /* Compute the common refinement over all tropical varieties
114  * of the polynomials in the generating set */
115  ZConesSortedByDimension C = tropicalVarietySortedByDimension(inI->m[0],r,currentStrategy);
116  int PayneOsserman = rVar(r)-1;
117  for (int i=0; i<k; i++)
118  {
119  if(inI->m[i]!=NULL)
120  {
121  PayneOsserman--;
122  C = intersect(C,tropicalVarietySortedByDimension(inI->m[i],r,currentStrategy),si_max(PayneOsserman,d));
123  }
124  }
125 
126  /* Cycle through all maximal cones of the refinement.
127  * Pick a monomial ordering corresponding to a generic weight vector in it
128  * and check if the initial ideal is monomial free, generic meaning
129  * that it lies in a maximal Groebner cone in the maximal cone of the refinement.
130  * If the initial ideal is not monomial free, compute a witness for the monomial
131  * and compute the common refinement with its tropical variety.
132  * If all initial ideals are monomial free, then we have our tropical curve */
133  // gfan::ZFan* zf = toFanStar(C);
134  // std::cout << zf->toString(2+4+8+128) << std::endl;
135  // delete zf;
136  for (std::set<gfan::ZCone>::iterator zc=C.begin(); zc!=C.end();)
137  {
138  gfan::ZVector w = zc->getRelativeInteriorPoint();
139  gfan::ZMatrix W = zc->generatorsOfSpan();
140  // std::cout << zc->extremeRays() << std::endl;
141 
142  ring s = genericlyWeightedOrdering(r,u,w,W,currentStrategy);
143  nMapFunc identity = n_SetMap(r->cf,s->cf);
144  ideal inIs = idInit(k);
145  for (int j=0; j<k; j++)
146  inIs->m[j] = p_PermPoly(inI->m[j],NULL,r,s,identity,NULL,0);
147 
148  ideal inIsSTD = gfanlib_kStd_wrapper(inIs,s,isHomog);
149  id_Delete(&inIs,s);
150  ideal ininIs = initial(inIsSTD,s,w,W);
151 
152  std::pair<poly,int> mons = currentStrategy->checkInitialIdealForMonomial(ininIs,s);
153 
154  if (mons.first!=NULL)
155  {
156  poly gs;
157  if (mons.second>=0)
158  // cheap way out, ininIsSTD already contains a monomial in its generators
159  gs = inIsSTD->m[mons.second];
160  else
161  // compute witness
162  gs = witness(mons.first,inIsSTD,ininIs,s);
163 
164  C = intersect(C,tropicalVarietySortedByDimension(gs,s,currentStrategy),d);
165  nMapFunc mMap = n_SetMap(s->cf,r->cf);
166  poly gr = p_PermPoly(gs,NULL,s,r,mMap,NULL,0);
167  idInsertPoly(inI,gr);
168  k++;
169 
170  if (mons.second<0)
171  {
172  // if necessary, cleanup mons and gs
173  p_Delete(&mons.first,s);
174  p_Delete(&gs,s);
175  }
176  // cleanup rest, reset zc
177  id_Delete(&inIsSTD,s);
178  id_Delete(&ininIs,s);
179  rDelete(s);
180  zc = C.begin();
181  }
182  else
183  {
184  // cleanup remaining data of first stage
185  id_Delete(&inIsSTD,s);
186  id_Delete(&ininIs,s);
187  rDelete(s);
188 
189  gfan::ZVector wNeg = currentStrategy->negateWeight(w);
190  if (zc->contains(wNeg))
191  {
192  s = genericlyWeightedOrdering(r,u,wNeg,W,currentStrategy);
193  identity = n_SetMap(r->cf,s->cf);
194  inIs = idInit(k);
195  for (int j=0; j<k; j++)
196  inIs->m[j] = p_PermPoly(inI->m[j],NULL,r,s,identity,NULL,0);
197 
198  inIsSTD = gfanlib_kStd_wrapper(inIs,s,isHomog);
199  id_Delete(&inIs,s);
200  ininIs = initial(inIsSTD,s,wNeg,W);
201 
202  mons = currentStrategy->checkInitialIdealForMonomial(ininIs,s);
203  if (mons.first!=NULL)
204  {
205  poly gs;
206  if (mons.second>=0)
207  // cheap way out, ininIsSTD already contains a monomial in its generators
208  gs = inIsSTD->m[mons.second];
209  else
210  // compute witness
211  gs = witness(mons.first,inIsSTD,ininIs,s);
212 
213  C = intersect(C,tropicalVarietySortedByDimension(gs,s,currentStrategy),d);
214  nMapFunc mMap = n_SetMap(s->cf,r->cf);
215  poly gr = p_PermPoly(gs,NULL,s,r,mMap,NULL,0);
216  idInsertPoly(inI,gr);
217  k++;
218 
219  if (mons.second<0)
220  {
221  // if necessary, cleanup mons and gs
222  p_Delete(&mons.first,s);
223  p_Delete(&gs,s);
224  }
225  // reset zc
226  zc = C.begin();
227  }
228  else
229  zc++;
230  // cleanup remaining data of second stage
231  id_Delete(&inIsSTD,s);
232  id_Delete(&ininIs,s);
233  rDelete(s);
234  }
235  else
236  zc++;
237  }
238  }
239  return C;
240 }
241 
242 
243 gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector &u, const tropicalStrategy* currentStrategy)
244 {
245  ZConesSortedByDimension C = tropicalStar(I,r,u,currentStrategy);
246  // gfan::ZFan* zf = toFanStar(C);
247  // std::cout << zf->toString();
248  // delete zf;
249  gfan::ZMatrix raysOfC(0,u.size());
250  if (!currentStrategy->restrictToLowerHalfSpace())
251  {
252  for (ZConesSortedByDimension::iterator zc=C.begin(); zc!=C.end(); zc++)
253  {
254  assume(zc->dimensionOfLinealitySpace()+1 >= zc->dimension());
255  if (zc->dimensionOfLinealitySpace()+1 >= zc->dimension())
256  raysOfC.appendRow(zc->getRelativeInteriorPoint());
257  else
258  {
259  gfan::ZVector interiorPoint = zc->getRelativeInteriorPoint();
260  if (!currentStrategy->homogeneitySpaceContains(interiorPoint))
261  {
262  raysOfC.appendRow(interiorPoint);
263  raysOfC.appendRow(currentStrategy->negateWeight(interiorPoint));
264  }
265  else
266  {
267  gfan::ZMatrix zm = zc->generatorsOfLinealitySpace();
268  for (int i=0; i<zm.getHeight(); i++)
269  {
270  gfan::ZVector point = zm[i];
271  if (currentStrategy->homogeneitySpaceContains(point))
272  {
273  raysOfC.appendRow(point);
274  raysOfC.appendRow(currentStrategy->negateWeight(point));
275  break;
276  }
277  }
278  }
279  }
280  }
281  }
282  else
283  {
284  for (ZConesSortedByDimension::iterator zc=C.begin(); zc!=C.end(); zc++)
285  {
286  assume(zc->dimensionOfLinealitySpace()+2 >= zc->dimension());
287  if (zc->dimensionOfLinealitySpace()+2 == zc->dimension())
288  raysOfC.appendRow(zc->getRelativeInteriorPoint());
289  else
290  {
291  gfan::ZVector interiorPoint = zc->getRelativeInteriorPoint();
292  if (!currentStrategy->homogeneitySpaceContains(interiorPoint))
293  {
294  raysOfC.appendRow(interiorPoint);
295  raysOfC.appendRow(currentStrategy->negateWeight(interiorPoint));
296  }
297  else
298  {
299  gfan::ZMatrix zm = zc->generatorsOfLinealitySpace();
300  for (int i=0; i<zm.getHeight(); i++)
301  {
302  gfan::ZVector point = zm[i];
303  if (currentStrategy->homogeneitySpaceContains(point))
304  {
305  raysOfC.appendRow(point);
306  raysOfC.appendRow(currentStrategy->negateWeight(point));
307  break;
308  }
309  }
310  }
311  }
312  }
313  }
314  return raysOfC;
315 }
static int si_max(const int a, const int b)
Definition: auxiliary.h:124
#define FALSE
Definition: auxiliary.h:96
int * ZVectorToIntStar(const gfan::ZVector &v, bool &overflow)
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
bool homogeneitySpaceContains(const gfan::ZVector &v) const
returns true, if v is contained in the homogeneity space; false otherwise
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...
int getExpectedDimension() const
returns the expected Dimension of the polyhedral output
gfan::ZVector adjustWeightForHomogeneity(gfan::ZVector w) const
Given weight w, returns a strictly positive weight u such that an ideal satisfying the valuation-sepc...
gfan::ZVector negateWeight(const gfan::ZVector &w) 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...
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:723
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:74
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm & w
Definition: facAbsFact.cc:51
int j
Definition: facHensel.cc:110
void WerrorS(const char *s)
Definition: feFopen.cc:24
BOOLEAN idInsertPoly(ideal h1, poly h2)
insert h2 into h1 (if h2 is not the zero polynomial) return TRUE iff h2 was indeed inserted
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
STATIC_VAR Poly * h
Definition: janet.cc:971
#define assume(x)
Definition: mod2.h:387
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
poly p_PermPoly(poly p, const int *perm, const ring oldRing, const ring dst, nMapFunc nMap, const int *par_perm, int OldPar, BOOLEAN use_mult)
Definition: p_polys.cc:4156
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:861
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:3403
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1366
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:449
rRingOrder_t
order stuff
Definition: ring.h:68
@ ringorder_a
Definition: ring.h:70
@ ringorder_C
Definition: ring.h:73
@ ringorder_wp
Definition: ring.h:81
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:597
#define rTest(r)
Definition: ring.h:790
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
#define IDELEMS(i)
Definition: simpleideals.h:23
ideal gfanlib_kStd_wrapper(ideal I, ring r, tHomog h=testHomog)
Definition: std_wrapper.cc:6
@ isHomog
Definition: structs.h:42
static ZConesSortedByDimension intersect(const ZConesSortedByDimension &setA, const ZConesSortedByDimension &setB, int d=0)
ZConesSortedByDimension tropicalStar(ideal inI, const ring r, const gfan::ZVector &u, const tropicalStrategy *currentStrategy)
static ring genericlyWeightedOrdering(const ring r, const gfan::ZVector &u, const gfan::ZVector &w, const gfan::ZMatrix &W, const tropicalStrategy *currentStrategy)
gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector &u, const tropicalStrategy *currentStrategy)
implementation of the class tropicalStrategy
ZConesSortedByDimension tropicalVarietySortedByDimension(const poly g, const ring r, const tropicalStrategy *currentCase)
std::set< gfan::ZCone, ZConeCompareDimensionFirst > ZConesSortedByDimension
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