source: git/Singular/dyn_modules/gfanlib/tropical.cc @ 76d8bb

spielwiese
Last change on this file since 76d8bb was 76d8bb, checked in by Yue <ren@…>, 8 years ago
fix: maximalGroebnerCone for ideals with generator 0 or generator monomial
  • Property mode set to 100644
File size: 13.8 KB
Line 
1#include <polys/monomials/p_polys.h>
2#include <coeffs/coeffs.h>
3
4#include <callgfanlib_conversion.h>
5#include <bbcone.h>
6#include <ppinitialReduction.h>
7#include <containsMonomial.h>
8#include <initial.h>
9#include <witness.h>
10#include <tropicalCurves.h>
11#include <tropicalStrategy.h>
12#include <startingCone.h>
13#include <groebnerFan.h>
14#include <groebnerComplex.h>
15#include <tropicalVariety.h>
16
17int tropicalVerboseLevel = 0;
18
19gfan::ZCone homogeneitySpace(ideal I, ring r)
20{
21  int n = rVar(r);
22  poly g;
23  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
24  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
25  gfan::ZVector leadexpw = gfan::ZVector(n);
26  gfan::ZVector tailexpw = gfan::ZVector(n);
27  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
28  for (int i=0; i<IDELEMS(I); i++)
29  {
30    g = (poly) I->m[i];
31    if (g)
32    {
33      p_GetExpV(g,leadexpv,r);
34      leadexpw = intStar2ZVector(n,leadexpv);
35      pIter(g);
36      while (g)
37      {
38        p_GetExpV(g,tailexpv,r);
39        tailexpw = intStar2ZVector(n,tailexpv);
40        equations.appendRow(leadexpw-tailexpw);
41        pIter(g);
42      }
43    }
44  }
45  omFreeSize(leadexpv,(n+1)*sizeof(int));
46  omFreeSize(tailexpv,(n+1)*sizeof(int));
47  return gfan::ZCone(gfan::ZMatrix(0, equations.getWidth()),equations);
48}
49
50
51BOOLEAN homogeneitySpace(leftv res, leftv args)
52{
53  leftv u = args;
54  if ((u != NULL) && (u->Typ() == POLY_CMD))
55  {
56    leftv v = u->next;
57    if (v == NULL)
58    {
59      poly g = (poly) u->Data();
60      ideal I = idInit(1);
61      I->m[0] = g;
62      res->rtyp = coneID;
63      res->data = (void*) new gfan::ZCone(homogeneitySpace(I,currRing));
64      I->m[0] = NULL;
65      id_Delete(&I,currRing);
66      return FALSE;
67    }
68  }
69  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
70  {
71    leftv v = u->next;
72    if (v == NULL)
73    {
74      ideal I = (ideal) u->Data();
75      res->rtyp = coneID;
76      res->data = (void*) new gfan::ZCone(homogeneitySpace(I,currRing));
77      return FALSE;
78    }
79  }
80  WerrorS("homogeneitySpace: unexpected parameters");
81  return TRUE;
82}
83
84
85gfan::ZCone lowerHomogeneitySpace(ideal I, ring r)
86{
87  int n = rVar(r);
88  poly g;
89  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
90  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
91  gfan::ZVector leadexpw = gfan::ZVector(n);
92  gfan::ZVector tailexpw = gfan::ZVector(n);
93  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
94  for (int i=0; i<IDELEMS(I); i++)
95  {
96    g = (poly) I->m[i];
97    if (g)
98    {
99      p_GetExpV(g,leadexpv,r);
100      leadexpw = intStar2ZVector(n,leadexpv);
101      pIter(g);
102      while (g)
103      {
104        p_GetExpV(g,tailexpv,r);
105        tailexpw = intStar2ZVector(n,tailexpv);
106        equations.appendRow(leadexpw-tailexpw);
107        pIter(g);
108      }
109    }
110  }
111  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
112  gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
113  lowerHalfSpaceCondition[0] = -1;
114  inequalities.appendRow(lowerHalfSpaceCondition);
115
116  omFreeSize(leadexpv,(n+1)*sizeof(int));
117  omFreeSize(tailexpv,(n+1)*sizeof(int));
118  return gfan::ZCone(inequalities,equations);
119}
120
121
122BOOLEAN lowerHomogeneitySpace(leftv res, leftv args)
123{
124  leftv u = args;
125  if ((u != NULL) && (u->Typ() == POLY_CMD))
126  {
127    leftv v = u->next;
128    if (v == NULL)
129    {
130      poly g = (poly) u->Data();
131      ideal I = idInit(1);
132      I->m[0] = g;
133      res->rtyp = coneID;
134      res->data = (void*) new gfan::ZCone(lowerHomogeneitySpace(I,currRing));
135      I->m[0] = NULL;
136      id_Delete(&I,currRing);
137      return FALSE;
138    }
139  }
140  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
141  {
142    leftv v = u->next;
143    if (v == NULL)
144    {
145      ideal I = (ideal) u->Data();
146      res->rtyp = coneID;
147      res->data = (void*) new gfan::ZCone(lowerHomogeneitySpace(I,currRing));
148      return FALSE;
149    }
150  }
151  WerrorS("lowerHomogeneitySpace: unexpected parameters");
152  return TRUE;
153}
154
155
156gfan::ZCone groebnerCone(const ideal I, const ring r, const gfan::ZVector &w)
157{
158  int n = rVar(r);
159  poly g = NULL;
160  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
161  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
162  gfan::ZVector leadexpw = gfan::ZVector(n);
163  gfan::ZVector tailexpw = gfan::ZVector(n);
164
165  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
166  for (int i=0; i<IDELEMS(I); i++)
167  {
168    g = (poly) I->m[i];
169    if (g!=NULL)
170    {
171      p_GetExpV(g,leadexpv,currRing);
172      leadexpw = intStar2ZVector(n, leadexpv);
173      pIter(g);
174      while (g!=NULL)
175      {
176        pGetExpV(g,tailexpv);
177        tailexpw = intStar2ZVector(n, tailexpv);
178        inequalities.appendRow(leadexpw-tailexpw);
179        pIter(g);
180      }
181    }
182  }
183
184  ideal inI = initial(I,currRing,w);
185  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
186  for (int i=0; i<IDELEMS(I); i++)
187  {
188    g = (poly) inI->m[i];
189    if (g!=NULL)
190    {
191      p_GetExpV(g,leadexpv,currRing);
192      leadexpw = intStar2ZVector(n, leadexpv);
193      pIter(g);
194      while (g!=NULL)
195      {
196        pGetExpV(g,tailexpv);
197        tailexpw = intStar2ZVector(n, tailexpv);
198        equations.appendRow(leadexpw-tailexpw);
199        pIter(g);
200      }
201    }
202  }
203
204  omFreeSize(leadexpv,(n+1)*sizeof(int));
205  omFreeSize(tailexpv,(n+1)*sizeof(int));
206  id_Delete(&inI,currRing);
207  return gfan::ZCone(inequalities,equations);
208}
209
210
211BOOLEAN groebnerCone(leftv res, leftv args)
212{
213  leftv u = args;
214  if ((u != NULL) && (u->Typ() == POLY_CMD))
215  {
216    leftv v = u->next;
217    if ((v !=NULL) && ((v->Typ() == BIGINTMAT_CMD) || (v->Typ() == INTVEC_CMD)))
218    {
219      try
220      {
221        poly g = (poly) u->Data();
222        ideal I = idInit(1);
223        I->m[0] = g;
224        gfan::ZVector* weightVector;
225        if (v->Typ() == INTVEC_CMD)
226        {
227          intvec* w0 = (intvec*) v->Data();
228          bigintmat* w1 = iv2bim(w0,coeffs_BIGINT);
229          w1->inpTranspose();
230          weightVector = bigintmatToZVector(*w1);
231          delete w1;
232        }
233        else
234        {
235          bigintmat* w1 = (bigintmat*) v->Data();
236          weightVector = bigintmatToZVector(*w1);
237        }
238        res->rtyp = coneID;
239        res->data = (void*) new gfan::ZCone(groebnerCone(I,currRing,*weightVector));
240        delete weightVector;
241        I->m[0] = NULL;
242        id_Delete(&I,currRing);
243        return FALSE;
244      }
245      catch (const std::exception& ex)
246      {
247        Werror("ERROR: %s",ex.what());
248        return TRUE;
249      }
250    }
251  }
252  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
253  {
254    leftv v = u->next;
255    if ((v !=NULL) && ((v->Typ() == BIGINTMAT_CMD) || (v->Typ() == INTVEC_CMD)))
256    {
257      try
258      {
259        ideal I = (ideal) u->Data();
260        gfan::ZVector* weightVector;
261        if (v->Typ() == INTVEC_CMD)
262        {
263          intvec* w0 = (intvec*) v->Data();
264          bigintmat* w1 = iv2bim(w0,coeffs_BIGINT);
265          w1->inpTranspose();
266          weightVector = bigintmatToZVector(*w1);
267          delete w1;
268        }
269        else
270        {
271          bigintmat* w1 = (bigintmat*) v->Data();
272          weightVector = bigintmatToZVector(*w1);
273        }
274        res->rtyp = coneID;
275        res->data = (void*) new gfan::ZCone(groebnerCone(I,currRing,*weightVector));
276        delete weightVector;
277        return FALSE;
278      }
279      catch (const std::exception& ex)
280      {
281        Werror("ERROR: %s",ex.what());
282        return TRUE;
283      }
284    }
285  }
286  WerrorS("groebnerCone: unexpected parameters");
287  return TRUE;
288}
289
290
291gfan::ZCone maximalGroebnerCone(const ideal &I, const ring &r)
292{
293  int n = rVar(r);
294  poly g = NULL;
295  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
296  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
297  gfan::ZVector leadexpw = gfan::ZVector(n);
298  gfan::ZVector tailexpw = gfan::ZVector(n);
299  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
300  for (int i=0; i<IDELEMS(I); i++)
301  {
302    g = (poly) I->m[i];
303    if (g != NULL && pNext(g) != NULL)
304    {
305      pGetExpV(g,leadexpv);
306      leadexpw = intStar2ZVector(n, leadexpv);
307      pIter(g);
308      while (g != NULL)
309      {
310        pGetExpV(g,tailexpv);
311        tailexpw = intStar2ZVector(n, tailexpv);
312        inequalities.appendRow(leadexpw-tailexpw);
313        pIter(g);
314      }
315    }
316  }
317  omFreeSize(leadexpv,(n+1)*sizeof(int));
318  omFreeSize(tailexpv,(n+1)*sizeof(int));
319  return gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
320}
321
322
323BOOLEAN maximalGroebnerCone(leftv res, leftv args)
324{
325  leftv u = args;
326  if ((u != NULL) && (u->Typ() == POLY_CMD))
327  {
328    leftv v = u->next;
329    if (v == NULL)
330    {
331      try
332      {
333        poly g = (poly) u->Data();
334        ideal I = idInit(1);
335        I->m[0] = g;
336        res->rtyp = coneID;
337        res->data = (void*) new gfan::ZCone(maximalGroebnerCone(I,currRing));
338        I->m[0] = NULL;
339        id_Delete(&I,currRing);
340        return FALSE;
341      }
342      catch (const std::exception& ex)
343      {
344        Werror("ERROR: %s",ex.what());
345        return TRUE;
346      }
347    }
348  }
349  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
350  {
351    leftv v = u->next;
352    if (v == NULL)
353    {
354      try
355      {
356        ideal I = (ideal) u->Data();
357        res->rtyp = coneID;
358        res->data = (void*) new gfan::ZCone(maximalGroebnerCone(I,currRing));
359        return FALSE;
360      }
361      catch (const std::exception& ex)
362      {
363        Werror("ERROR: %s",ex.what());
364        return TRUE;
365      }
366    }
367  }
368  WerrorS("maximalGroebnerCone: unexpected parameters");
369  return TRUE;
370}
371
372
373BOOLEAN initial(leftv res, leftv args)
374{
375  leftv u = args;
376  if ((u != NULL) && (u->Typ() == POLY_CMD))
377  {
378    leftv v = u->next;
379    if ((v !=NULL) && ((v->Typ() == BIGINTMAT_CMD) || (v->Typ() == INTVEC_CMD)))
380    {
381      poly p = (poly) u->Data();
382      gfan::ZVector* weightVector;
383      if (v->Typ() == INTVEC_CMD)
384      {
385        intvec* w0 = (intvec*) v->Data();
386        bigintmat* w1 = iv2bim(w0,coeffs_BIGINT);
387        w1->inpTranspose();
388        weightVector = bigintmatToZVector(*w1);
389        delete w1;
390      }
391      else
392      {
393        bigintmat* w1 = (bigintmat*) v->Data();
394        weightVector = bigintmatToZVector(*w1);
395      }
396      res->rtyp = POLY_CMD;
397      res->data = (void*) initial(p, currRing, *weightVector);
398      delete weightVector;
399      return FALSE;
400    }
401  }
402  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
403  {
404    leftv v = u->next;
405    if ((v !=NULL) && ((v->Typ() == BIGINTMAT_CMD) || (v->Typ() == INTVEC_CMD)))
406    {
407      try
408      {
409        ideal I = (ideal) u->Data();
410        gfan::ZVector* weightVector;
411        if (v->Typ() == INTVEC_CMD)
412        {
413          intvec* w0 = (intvec*) v->Data();
414          bigintmat* w1 = iv2bim(w0,coeffs_BIGINT);
415          w1->inpTranspose();
416          weightVector = bigintmatToZVector(*w1);
417          delete w1;
418        }
419        else
420        {
421          bigintmat* w1 = (bigintmat*) v->Data();
422          weightVector = bigintmatToZVector(*w1);
423        }
424        res->rtyp = IDEAL_CMD;
425        res->data = (void*) initial(I, currRing, *weightVector);
426        delete weightVector;
427        return FALSE;
428      }
429      catch (const std::exception& ex)
430      {
431        Werror("ERROR: %s",ex.what());
432        return TRUE;
433      }
434    }
435  }
436  WerrorS("initial: unexpected parameters");
437  return TRUE;
438}
439
440
441void tropical_setup(SModulFunctions* p)
442{
443  p->iiAddCproc("","groebnerCone",FALSE,groebnerCone);
444  p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
445  p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
446  p->iiAddCproc("","lowerHomogeneitySpace",FALSE,lowerHomogeneitySpace);
447  p->iiAddCproc("","initial",FALSE,initial);
448  p->iiAddCproc("","tropicalVariety",FALSE,tropicalVariety);
449  p->iiAddCproc("","groebnerFan",FALSE,groebnerFan);
450  p->iiAddCproc("","groebnerComplex",FALSE,groebnerComplex);
451#ifndef NDEBUG
452  // p->iiAddCproc("","tropicalNeighbours",FALSE,tropicalNeighbours);
453  // p->iiAddCproc("","initial0",FALSE,initial0);
454  p->iiAddCproc("","pReduceDebug",FALSE,pReduceDebug);
455  // p->iiAddCproc("","ppreduceInitially0",FALSE,ppreduceInitially0);
456  // p->iiAddCproc("","ppreduceInitially1",FALSE,ppreduceInitially1);
457  // p->iiAddCproc("","ppreduceInitially2",FALSE,ppreduceInitially2);
458  p->iiAddCproc("","ptNormalize",FALSE,ptNormalize);
459  p->iiAddCproc("","ppreduceInitially3",FALSE,ppreduceInitially3);
460  // p->iiAddCproc("","ppreduceInitially4",FALSE,ppreduceInitially4);
461  // p->iiAddCproc("","ttpReduce",FALSE,ttpReduce);
462  // p->iiAddCproc("","ttreduceInitially0",FALSE,ttreduceInitially0);
463  // p->iiAddCproc("","ttreduceInitially1",FALSE,ttreduceInitially1);
464  // p->iiAddCproc("","ttreduceInitially2",FALSE,ttreduceInitially2);
465  // p->iiAddCproc("","ttreduceInitially3",FALSE,ttreduceInitially3);
466  // p->iiAddCproc("","ttreduceInitially4",FALSE,ttreduceInitially4);
467  // p->iiAddCproc("","checkForMonomial",FALSE,checkForMonomial);
468  // p->iiAddCproc("","dwr0",FALSE,dwr0);
469  // p->iiAddCproc("","witness0",FALSE,witness0);
470  // p->iiAddCproc("","tropicalVariety00",FALSE,tropicalVariety00);
471  // p->iiAddCproc("","tropicalVariety01",FALSE,tropicalVariety01);
472  // p->iiAddCproc("","tropicalCurve0",FALSE,tropicalCurve0);
473  // p->iiAddCproc("","tropicalCurve1",FALSE,tropicalCurve1);
474  p->iiAddCproc("","reduceInitiallyDebug",FALSE,reduceInitiallyDebug);
475  p->iiAddCproc("","computeWitnessDebug",FALSE,computeWitnessDebug);
476  p->iiAddCproc("","computeFlipDebug",FALSE,computeFlipDebug);
477  p->iiAddCproc("","flipConeDebug",FALSE,flipConeDebug);
478  // p->iiAddCproc("","groebnerNeighboursDebug",FALSE,groebnerNeighboursDebug);
479  // p->iiAddCproc("","tropicalNeighboursDebug",FALSE,tropicalNeighboursDebug);
480  p->iiAddCproc("","tropicalStarDebug",FALSE,tropicalStarDebug);
481  p->iiAddCproc("","tropicalStartingPoint",FALSE,tropicalStartingPoint);
482  p->iiAddCproc("","positiveTropicalStartingPoint",FALSE,positiveTropicalStartingPoint);
483  p->iiAddCproc("","nonNegativeTropicalStartingPoint",FALSE,nonNegativeTropicalStartingPoint);
484  p->iiAddCproc("","negativeTropicalStartingPoint",FALSE,negativeTropicalStartingPoint);
485  p->iiAddCproc("","nonPositiveTropicalStartingPoint",FALSE,nonPositiveTropicalStartingPoint);
486  p->iiAddCproc("","tropicalStartingCone",FALSE,tropicalStartingCone);
487#endif //NDEBUG
488  // p->iiAddCproc("","ppreduceInitially",FALSE,ppreduceInitially);
489  // p->iiAddCproc("","ttreduceInitially",FALSE,ttreduceInitially);
490}
Note: See TracBrowser for help on using the repository browser.