source: git/Singular/dyn_modules/gfanlib/bbpolytope.cc @ fbaeb5

spielwiese
Last change on this file since fbaeb5 was fbaeb5, checked in by Yue <ren@…>, 8 years ago
fix: mixedVolume small mistakes
  • Property mode set to 100644
File size: 14.6 KB
Line 
1#include <kernel/mod2.h>
2
3#if HAVE_GFANLIB
4
5#include <Singular/ipid.h>
6#include <Singular/ipshell.h>
7#include <Singular/blackbox.h>
8#include <misc/intvec.h>
9#include <coeffs/bigintmat.h>
10
11#include <callgfanlib_conversion.h>
12#include <sstream>
13
14#include <gfanlib/gfanlib.h>
15#include <gfanlib/gfanlib_q.h>
16
17int polytopeID;
18
19std::string bbpolytopeToString(gfan::ZCone const &c)
20{
21  std::stringstream s;
22  gfan::ZMatrix i=c.getInequalities();
23  gfan::ZMatrix e=c.getEquations();
24  s<<"AMBIENT_DIM"<<std::endl;
25  s<<c.ambientDimension()-1<<std::endl;
26  s<<"INEQUALITIES"<<std::endl;
27  s<<toString(i)<<std::endl;
28  s<<"EQUATIONS"<<std::endl;
29  s<<toString(e)<<std::endl;
30  return s.str();
31}
32
33void *bbpolytope_Init(blackbox* /*b*/)
34{
35  return (void*)(new gfan::ZCone());
36}
37
38BOOLEAN bbpolytope_Assign(leftv l, leftv r)
39{
40  gfan::ZCone* newZc;
41  if (r==NULL)
42  {
43    if (l->Data()!=NULL)
44    {
45      gfan::ZCone* zd = (gfan::ZCone*)l->Data();
46      delete zd;
47    }
48    newZc = new gfan::ZCone();
49  }
50  else if (r->Typ()==l->Typ())
51  {
52    if (l->Data()!=NULL)
53    {
54      gfan::ZCone* zd = (gfan::ZCone*)l->Data();
55      delete zd;
56    }
57    gfan::ZCone* zc = (gfan::ZCone*)r->Data();
58    newZc = new gfan::ZCone(*zc);
59  }
60  // else if (r->Typ()==INT_CMD)  TODO:r->Typ()==BIGINTMAT_CMD
61  // {
62  //   int ambientDim = (int)(long)r->Data();
63  //   if (ambientDim < 0)
64  //   {
65  //     Werror("expected an int >= 0, but got %d", ambientDim);
66  //     return TRUE;
67  //   }
68  //   if (l->Data()!=NULL)
69  //   {
70  //     gfan::ZCone* zd = (gfan::ZCone*)l->Data();
71  //     delete zd;
72  //   }
73  //   newZc = new gfan::ZCone(ambientDim);
74  // }
75  else
76  {
77    Werror("assign Type(%d) = Type(%d) not implemented",l->Typ(),r->Typ());
78    return TRUE;
79  }
80
81  if (l->rtyp==IDHDL)
82  {
83    IDDATA((idhdl)l->data) = (char*) newZc;
84  }
85  else
86  {
87    l->data=(void *)newZc;
88  }
89  return FALSE;
90}
91
92char* bbpolytope_String(blackbox* /*b*/, void *d)
93{ if (d==NULL) return omStrDup("invalid object");
94   else
95   {
96     gfan::ZCone* zc = (gfan::ZCone*)d;
97     std::string s=bbpolytopeToString(*zc);
98     return omStrDup(s.c_str());
99   }
100}
101
102void bbpolytope_destroy(blackbox* /*b*/, void *d)
103{
104  if (d!=NULL)
105  {
106    gfan::ZCone* zc = (gfan::ZCone*) d;
107    delete zc;
108  }
109}
110
111void* bbpolytope_Copy(blackbox* /*b*/, void *d)
112{
113  gfan::ZCone* zc = (gfan::ZCone*)d;
114  gfan::ZCone* newZc = new gfan::ZCone(*zc);
115  return newZc;
116}
117
118static BOOLEAN ppCONERAYS1(leftv res, leftv v)
119{
120  /* method for generating a cone object from half-lines
121     (cone = convex hull of the half-lines; note: there may be
122     entire lines in the cone);
123     valid parametrizations: (bigintmat) */
124  bigintmat* rays = NULL;
125  if (v->Typ() == INTMAT_CMD)
126  {
127    intvec* rays0 = (intvec*) v->Data();
128    rays = iv2bim(rays0,coeffs_BIGINT);
129  }
130  else
131    rays = (bigintmat*) v->Data();
132
133  gfan::ZMatrix* zm = bigintmatToZMatrix(rays);
134  gfan::ZCone* zc = new gfan::ZCone();
135  *zc = gfan::ZCone::givenByRays(*zm, gfan::ZMatrix(0, zm->getWidth()));
136  res->rtyp = polytopeID;
137  res->data = (void*) zc;
138
139  delete zm;
140  if (v->Typ() == INTMAT_CMD)
141    delete rays;
142  return FALSE;
143}
144
145static BOOLEAN ppCONERAYS3(leftv res, leftv u, leftv v)
146{
147  /* method for generating a cone object from half-lines
148     (any point in the cone being the sum of a point
149     in the convex hull of the half-lines and a point in the span
150     of the lines), and an integer k;
151     valid parametrizations: (bigintmat, int);
152     Errors will be invoked in the following cases:
153     - k not 0 or 1;
154     if the k=1, then the extreme rays are known:
155     each half-line spans a (different) extreme ray */
156  bigintmat* rays = NULL;
157  if (u->Typ() == INTMAT_CMD)
158  {
159    intvec* rays0 = (intvec*) u->Data();
160    rays = iv2bim(rays0,coeffs_BIGINT);
161  }
162  else
163    rays = (bigintmat*) u->Data();
164  int k = (int)(long)v->Data();
165
166  if ((k < 0) || (k > 1))
167  {
168    WerrorS("expected int argument in [0..1]");
169    return TRUE;
170  }
171  k=k*2;
172  gfan::ZMatrix* zm = bigintmatToZMatrix(rays);
173  gfan::ZCone* zc = new gfan::ZCone();
174  *zc = gfan::ZCone::givenByRays(*zm,gfan::ZMatrix(0, zm->getWidth()));
175  //k should be passed on to zc; not available yet
176  res->rtyp = polytopeID;
177  res->data = (void*) zc;
178
179  delete zm;
180  if (v->Typ() == INTMAT_CMD)
181    delete rays;
182  return FALSE;
183}
184
185BOOLEAN polytopeViaVertices(leftv res, leftv args)
186{
187  gfan::initializeCddlibIfRequired();
188  leftv u = args;
189  if ((u != NULL) && ((u->Typ() == BIGINTMAT_CMD) || (u->Typ() == INTMAT_CMD)))
190  {
191    if (u->next == NULL) return ppCONERAYS1(res, u);
192    leftv v = u->next;
193    if ((v != NULL) && (v->Typ() == INT_CMD))
194    {
195      if (v->next == NULL) return ppCONERAYS3(res, u, v);
196    }
197  }
198  WerrorS("polytopeViaPoints: unexpected parameters");
199  return TRUE;
200}
201
202static BOOLEAN ppCONENORMALS1(leftv res, leftv v)
203{
204  /* method for generating a cone object from inequalities;
205     valid parametrizations: (bigintmat) */
206  bigintmat* ineq = NULL;
207  if (v->Typ() == INTMAT_CMD)
208  {
209    intvec* ineq0 = (intvec*) v->Data();
210    ineq = iv2bim(ineq0,coeffs_BIGINT);
211  }
212  else
213    ineq = (bigintmat*) v->Data();
214  gfan::ZMatrix* zm = bigintmatToZMatrix(ineq);
215  gfan::ZCone* zc = new gfan::ZCone(*zm, gfan::ZMatrix(0, zm->getWidth()));
216  delete zm;
217  if (v->Typ() == INTMAT_CMD)
218    delete ineq;
219  res->rtyp = polytopeID;
220  res->data = (void*) zc;
221  return FALSE;
222}
223
224static BOOLEAN ppCONENORMALS2(leftv res, leftv u, leftv v)
225{
226  /* method for generating a cone object from iequalities,
227     and equations (...)
228     valid parametrizations: (bigintmat, bigintmat)
229     Errors will be invoked in the following cases:
230     - u and v have different numbers of columns */
231  bigintmat* ineq = NULL; bigintmat* eq = NULL;
232  if (u->Typ() == INTMAT_CMD)
233  {
234    intvec* ineq0 = (intvec*) u->Data();
235    ineq = iv2bim(ineq0,coeffs_BIGINT);
236  }
237  else
238    ineq = (bigintmat*) u->Data();
239  if (v->Typ() == INTMAT_CMD)
240  {
241    intvec* eq0 = (intvec*) v->Data();
242    eq = iv2bim(eq0,coeffs_BIGINT);
243  }
244  else
245    eq = (bigintmat*) v->Data();
246
247  if (ineq->cols() != eq->cols())
248  {
249    Werror("expected same number of columns but got %d vs. %d",
250           ineq->cols(), eq->cols());
251    return TRUE;
252  }
253  gfan::ZMatrix* zm1 = bigintmatToZMatrix(ineq);
254  gfan::ZMatrix* zm2 = bigintmatToZMatrix(eq);
255  gfan::ZCone* zc = new gfan::ZCone(*zm1, *zm2);
256  delete zm1;
257  delete zm2;
258  if (u->Typ() == INTMAT_CMD)
259    delete ineq;
260  if (v->Typ() == INTMAT_CMD)
261    delete eq;
262
263  res->rtyp = polytopeID;
264  res->data = (void*) zc;
265  return FALSE;
266}
267
268static BOOLEAN ppCONENORMALS3(leftv res, leftv u, leftv v, leftv w)
269{
270  /* method for generating a cone object from inequalities, equations,
271     and an integer k;
272     valid parametrizations: (bigintmat, bigintmat, int);
273     Errors will be invoked in the following cases:
274     - u and v have different numbers of columns,
275     - k not in [0..3];
276     if the 2^0-bit of k is set, then ... */
277  bigintmat* ineq = NULL; bigintmat* eq = NULL;
278  if (u->Typ() == INTMAT_CMD)
279  {
280    intvec* ineq0 = (intvec*) u->Data();
281    ineq = iv2bim(ineq0,coeffs_BIGINT);
282  }
283  else
284    ineq = (bigintmat*) u->Data();
285  if (v->Typ() == INTMAT_CMD)
286  {
287    intvec* eq0 = (intvec*) v->Data();
288    eq = iv2bim(eq0,coeffs_BIGINT);
289  }
290  else
291    eq = (bigintmat*) v->Data();
292
293  if (ineq->cols() != eq->cols())
294  {
295    Werror("expected same number of columns but got %d vs. %d",
296           ineq->cols(), eq->cols());
297    return TRUE;
298  }
299  int k = (int)(long)w->Data();
300  if ((k < 0) || (k > 3))
301  {
302    WerrorS("expected int argument in [0..3]");
303    return TRUE;
304  }
305  gfan::ZMatrix* zm1 = bigintmatToZMatrix(ineq);
306  gfan::ZMatrix* zm2 = bigintmatToZMatrix(eq);
307  gfan::ZCone* zc = new gfan::ZCone(*zm1, *zm2, k);
308  delete zm1;
309  delete zm2;
310  if (u->Typ() == INTMAT_CMD)
311    delete ineq;
312  if (v->Typ() == INTMAT_CMD)
313    delete eq;
314
315  res->rtyp = polytopeID;
316  res->data = (void*) zc;
317  return FALSE;
318}
319
320BOOLEAN polytopeViaNormals(leftv res, leftv args)
321{
322  gfan::initializeCddlibIfRequired();
323  leftv u = args;
324  if ((u != NULL) && ((u->Typ() == BIGINTMAT_CMD) || (u->Typ() == INTMAT_CMD)))
325  {
326    if (u->next == NULL) return ppCONENORMALS1(res, u);
327  }
328  leftv v = u->next;
329  if ((v != NULL) && ((v->Typ() == BIGINTMAT_CMD) || (v->Typ() == INTMAT_CMD)))
330  {
331    if (v->next == NULL) return ppCONENORMALS2(res, u, v);
332  }
333  leftv w = v->next;
334  if ((w != NULL) && (w->Typ() == INT_CMD))
335  {
336    if (w->next == NULL) return ppCONENORMALS3(res, u, v, w);
337  }
338  WerrorS("polytopeViaInequalities: unexpected parameters");
339  return TRUE;
340}
341
342BOOLEAN vertices(leftv res, leftv args)
343{
344  gfan::initializeCddlibIfRequired();
345  leftv u = args;
346  if ((u != NULL) && (u->Typ() == polytopeID))
347    {
348      gfan::ZCone* zc = (gfan::ZCone*)u->Data();
349      gfan::ZMatrix zmat = zc->extremeRays();
350      res->rtyp = BIGINTMAT_CMD;
351      res->data = (void*)zMatrixToBigintmat(zmat);
352      return FALSE;
353    }
354  WerrorS("vertices: unexpected parameters");
355  return TRUE;
356}
357
358int getAmbientDimension(gfan::ZCone* zc) // zc is meant to represent a polytope here
359{                                        // hence ambientDimension-1
360  return zc->ambientDimension()-1;
361}
362
363int getCodimension(gfan::ZCone *zc)
364{
365  return zc->codimension();
366}
367
368int getDimension(gfan::ZCone* zc)
369{
370  return zc->dimension()-1;
371}
372
373gfan::ZVector intStar2ZVectorWithLeadingOne(const int d, const int* i)
374{
375  gfan::ZVector zv(d+1);
376  zv[0]=1;
377  for(int j=1; j<=d; j++)
378  {
379    zv[j]=i[j];
380  }
381  return zv;
382}
383
384gfan::ZCone newtonPolytope(poly p, ring r)
385{
386  int N = rVar(r);
387  gfan::ZMatrix zm(0,N+1);
388  int *leadexpv = (int*)omAlloc((N+1)*sizeof(int));
389  while (p!=NULL)
390  {
391    p_GetExpV(p,leadexpv,r);
392    gfan::ZVector zv = intStar2ZVectorWithLeadingOne(N, leadexpv);
393    zm.appendRow(zv);
394    pIter(p);
395  }
396  omFreeSize(leadexpv,(N+1)*sizeof(int));
397  gfan::ZCone Delta = gfan::ZCone::givenByRays(zm,gfan::ZMatrix(0, zm.getWidth()));
398  return Delta;
399}
400
401BOOLEAN newtonPolytope(leftv res, leftv args)
402{
403  gfan::initializeCddlibIfRequired();
404  leftv u = args;
405  if ((u != NULL) && (u->Typ() == POLY_CMD))
406  {
407    poly p = (poly)u->Data();
408    res->rtyp = polytopeID;
409    res->data = (void*) new gfan::ZCone(newtonPolytope(p,currRing));
410    return FALSE;
411  }
412  WerrorS("newtonPolytope: unexpected parameters");
413  return TRUE;
414}
415
416BOOLEAN scalePolytope(leftv res, leftv args)
417{
418  gfan::initializeCddlibIfRequired();
419  leftv u = args;
420  if ((u != NULL) && (u->Typ() == INT_CMD))
421  {
422    leftv v = u->next;
423    if ((v != NULL) && (v->Typ() == polytopeID))
424    {
425      int s = (int)(long) u->Data();
426      gfan::ZCone* zp = (gfan::ZCone*) v->Data();
427      gfan::ZMatrix zm = zp->extremeRays();
428      for (int i=0; i<zm.getHeight(); i++)
429        for (int j=1; j<zm.getWidth(); j++)
430          zm[i][j]*=s;
431      gfan::ZCone* zq = new gfan::ZCone();
432      *zq = gfan::ZCone::givenByRays(zm,gfan::ZMatrix(0, zm.getWidth()));
433      res->rtyp = polytopeID;
434      res->data = (void*) zq;
435      return FALSE;
436    }
437  }
438  WerrorS("scalePolytope: unexpected parameters");
439  return TRUE;
440}
441
442BOOLEAN dualPolytope(leftv res, leftv args)
443{
444  gfan::initializeCddlibIfRequired();
445  leftv u = args;
446  if ((u != NULL) && (u->Typ() == polytopeID))
447  {
448    gfan::ZCone* zp = (gfan::ZCone*) u->Data();
449    gfan::ZCone* zq = new gfan::ZCone(zp->dualCone());
450    res->rtyp = polytopeID;
451    res->data = (void*) zq;
452    return FALSE;
453  }
454  WerrorS("dualPolytope: unexpected parameters");
455  return TRUE;
456}
457
458BOOLEAN mixedVolume(leftv res, leftv args)
459{
460  gfan::initializeCddlibIfRequired();
461  leftv u = args;
462  if ((u != NULL) && (u->Typ() == LIST_CMD))
463  {
464    lists l = (lists) u->Data();
465    int k = lSize(l)+1;
466    std::vector<gfan::IntMatrix> P(k);
467    for (int i=0; i<k; i++)
468    {
469      if (l->m[i].Typ() == polytopeID)
470      {
471        gfan::ZCone* p = (gfan::ZCone*) l->m[i].Data();
472        gfan::ZMatrix pv = p->extremeRays();
473        int r = pv.getHeight();
474        int c = pv.getWidth();
475        gfan::IntMatrix pw(r,c-1);
476        for (int n=0; n<r; n++)
477          for (int m=1; m<c; m++)
478            pw[n][m-1] = pv[n][m].toInt();
479        P[i]=pw.transposed();
480      } else if (l->m[i].Typ() == POLY_CMD)
481      {
482        poly p = (poly) l->m[i].Data();
483        int N = rVar(currRing);
484        gfan::IntMatrix pw(0,N);
485        int *leadexpv = (int*)omAlloc((N+1)*sizeof(int));
486        while (p!=NULL)
487        {
488          p_GetExpV(p,leadexpv,currRing);
489          gfan::IntVector zv(N);
490          for (int i=0; i<N; i++)
491            zv[i] = leadexpv[i+1];
492          pw.appendRow(zv);
493          pIter(p);
494        }
495        P[i]=pw.transposed();
496        omFreeSize(leadexpv,(N+1)*sizeof(int));
497      }
498      else
499      {
500        WerrorS("mixedVolume: entries of unsupported type in list");
501        return TRUE;
502      }
503    }
504    gfan::Integer mv = gfan::mixedVolume(P);
505
506    res->rtyp = BIGINT_CMD;
507    res->data = (void*) integerToNumber(mv);
508    return FALSE;
509  }
510  WerrorS("mixedVolume: unexpected parameters");
511  return TRUE;
512}
513
514
515
516void bbpolytope_setup(SModulFunctions* p)
517{
518  blackbox *b=(blackbox*)omAlloc0(sizeof(blackbox));
519  // all undefined entries will be set to default in setBlackboxStuff
520  // the default Print is quite usefule,
521  // all other are simply error messages
522  b->blackbox_destroy=bbpolytope_destroy;
523  b->blackbox_String=bbpolytope_String;
524  //b->blackbox_Print=blackbox_default_Print;
525  b->blackbox_Init=bbpolytope_Init;
526  b->blackbox_Copy=bbpolytope_Copy;
527  b->blackbox_Assign=bbpolytope_Assign;
528  p->iiAddCproc("","polytopeViaPoints",FALSE,polytopeViaVertices);
529  p->iiAddCproc("","polytopeViaInequalities",FALSE,polytopeViaNormals);
530  p->iiAddCproc("","vertices",FALSE,vertices);
531  p->iiAddCproc("","newtonPolytope",FALSE,newtonPolytope);
532  p->iiAddCproc("","scalePolytope",FALSE,scalePolytope);
533  p->iiAddCproc("","dualPolytope",FALSE,dualPolytope);
534  p->iiAddCproc("","mixedVolume",FALSE,mixedVolume);
535  /********************************************************/
536  /* the following functions are implemented in bbcone.cc */
537  // iiAddCproc("","getAmbientDimension",FALSE,getAmbientDimension);
538  // iiAddCproc("","getCodimension",FALSE,getAmbientDimension);
539  // iiAddCproc("","getDimension",FALSE,getDimension);
540  /********************************************************/
541  /* the following functions are identical to those in bbcone.cc */
542  // iiAddCproc("","facets",FALSE,facets);
543  // iiAddCproc("","setLinearForms",FALSE,setLinearForms);
544  // iiAddCproc("","getLinearForms",FALSE,getLinearForms);
545  // iiAddCproc("","setMultiplicity",FALSE,setMultiplicity);
546  // iiAddCproc("","getMultiplicity",FALSE,getMultiplicity);
547  // iiAddCproc("","hasFace",FALSE,hasFace);
548  /***************************************************************/
549  // iiAddCproc("","getEquations",FALSE,getEquations);
550  // iiAddCproc("","getInequalities",FALSE,getInequalities);
551  polytopeID=setBlackboxStuff(b,"polytope");
552  //Print("created type %d (polytope)\n",polytopeID);
553}
554
555#endif
Note: See TracBrowser for help on using the repository browser.