source: git/dyn_modules/callpolymake/polymake_conversion.cc @ 62de8de

fieker-DuValspielwiese
Last change on this file since 62de8de was 62de8de, checked in by Yue Ren <ren@…>, 11 years ago
chg: updated callpolymake to newest version
  • Property mode set to 100644
File size: 13.9 KB
Line 
1#include <gmpxx.h>
2
3#include <polymake/Main.h>
4#include <polymake/Matrix.h>
5#include <polymake/Rational.h>
6#include <polymake/Integer.h>
7#include <polymake/Set.h>
8#include <polymake/common/lattice_tools.h>
9#include <polymake/IncidenceMatrix.h>
10
11#include <gfanlib/gfanlib.h>
12#include <gfanlib/gfanlib_q.h>
13
14#include <kernel/mod2.h>
15#include <libpolys/misc/intvec.h>
16#include <libpolys/coeffs/numbers.h>
17#include <libpolys/coeffs/bigintmat.h>
18#include <Singular/lists.h>
19#include <Singular/ipid.h> // for bigints,
20// is there really nothing better than this?
21
22/* Functions for converting Integers, Rationals and their Matrices
23   in between C++, gfan, polymake and singular */
24
25/* gfan -> polymake */
26
27polymake::Integer GfInteger2PmInteger (const gfan::Integer& gi)
28{
29  mpz_t cache; mpz_init(cache);
30  gi.setGmp(cache);
31  polymake::Integer pi(cache);
32  return pi;
33}
34
35polymake::Rational GfRational2PmRational (const gfan::Rational& gr)
36{
37  mpq_t cache; mpq_init(cache);
38  gr.setGmp(cache);
39  polymake::Rational pr(cache);
40  return pr;
41}
42
43polymake::Vector<polymake::Integer> Intvec2PmVectorInteger (const intvec* iv)
44{
45  polymake::Vector<polymake::Integer> vi(iv->length());
46  for(int i=1; i<=iv->length(); i++)
47  {
48    vi[i-1]=(*iv)[i-1];
49  }
50  return vi;
51}
52
53polymake::Matrix<polymake::Integer> GfZMatrix2PmMatrixInteger (const gfan::ZMatrix* zm)
54{
55  int rows=zm->getHeight();
56  int cols=zm->getWidth();
57  polymake::Matrix<polymake::Integer> mi(rows,cols);
58  for(int r=1; r<=rows; r++)
59    for(int c=1; c<=cols; c++)
60      mi(r-1,c-1) = GfInteger2PmInteger((*zm)[r-1][c-1]);
61  return mi;
62}
63
64polymake::Matrix<polymake::Rational> GfQMatrix2PmMatrixRational (const gfan::QMatrix* qm)
65{
66  int rows=qm->getHeight();
67  int cols=qm->getWidth();
68  polymake::Matrix<polymake::Rational> mr(rows,cols);
69  for(int r=1; r<=rows; r++)
70    for(int c=1; c<=cols; c++)
71      mr(r-1,c-1) = GfRational2PmRational((*qm)[r-1][c-1]);
72  return mr;
73}
74
75/* gfan <- polymake */
76
77gfan::Integer PmInteger2GfInteger (const polymake::Integer& pi)
78{
79  mpz_class cache(pi.get_rep());
80  gfan::Integer gi(cache.get_mpz_t());
81  return gi;
82}
83
84gfan::Rational PmRational2GfRational (const polymake::Rational& pr)
85{
86  mpq_class cache(pr.get_rep());
87  gfan::Rational gr(cache.get_mpq_t());
88  return gr;
89}
90
91gfan::ZMatrix PmMatrixInteger2GfZMatrix (const polymake::Matrix<polymake::Integer>* mi)
92{
93  int rows=mi->rows();
94  int cols=mi->cols();
95  gfan::ZMatrix zm(rows,cols);
96  for(int r=1; r<=rows; r++)
97    for(int c=1; c<=cols; c++)
98      zm[r-1][c-1] = PmInteger2GfInteger((*mi)(r-1,c-1));
99  return zm;
100}
101
102gfan::QMatrix PmMatrixRational2GfQMatrix (const polymake::Matrix<polymake::Rational>* mr)
103{
104  int rows=mr->rows();
105  int cols=mr->cols();
106  gfan::QMatrix qm(rows,cols);
107  for(int r=1; r<=rows; r++)
108    for(int c=1; c<=cols; c++)
109      qm[r-1][c-1] = PmRational2GfRational((*mr)(r-1,c-1));
110  return qm;
111}
112
113/* polymake -> singular */
114
115int PmInteger2Int(const polymake::Integer& pi, bool &ok)
116{
117  int i=0;
118  try
119  {
120    i = pi.to_int();
121  }
122  catch (const std::exception& ex)
123  {
124    WerrorS("ERROR: "); WerrorS(ex.what()); WerrorS("\n");
125    ok = false;
126  }
127  return i;
128}
129
130number PmInteger2Number (const polymake::Integer& pi)
131{
132  mpz_class cache(pi.get_rep());
133  long m = 268435456;
134  if(mpz_cmp_si(cache.get_mpz_t(),m))
135  {
136    int temp = (int) mpz_get_si(cache.get_mpz_t());
137    return n_Init(temp,coeffs_BIGINT);
138  }
139  else
140    return n_InitMPZ(cache.get_mpz_t(),coeffs_BIGINT);
141}
142
143intvec* PmVectorInteger2Intvec (const polymake::Vector<polymake::Integer>* vi, bool &ok)
144{
145  intvec* iv = new intvec(vi->size());
146  for(int i=1; i<=vi->size(); i++)
147  {
148    (*iv)[i-1] = PmInteger2Int((*vi)[i-1],ok);
149  }
150  return iv;
151}
152
153intvec* PmMatrixInteger2Intvec (polymake::Matrix<polymake::Integer>* mi, bool &ok)
154{
155  int rows = mi->rows();
156  int cols = mi->cols();
157  intvec* iv = new intvec(rows,cols,0);
158  const polymake::Integer* pi = concat_rows(*mi).begin();
159  for (int r = 1; r <= rows; r++)
160    for (int c = 1; c <= cols; c++)
161    {
162      IMATELEM(*iv,r,c) = PmInteger2Int(*pi, ok);
163      pi++;
164    }
165  return iv;
166}
167
168bigintmat* PmMatrixInteger2Bigintmat (polymake::Matrix<polymake::Integer>* mi)
169{
170  int rows = mi->rows();
171  int cols = mi->cols();
172  bigintmat* bim= new bigintmat(rows,cols,coeffs_BIGINT);
173  const polymake::Integer* pi = concat_rows(*mi).begin();
174  for (int r = 1; r <= rows; r++)
175    for (int c = 1; c <= cols; c++)
176    {
177      number temp = PmInteger2Number(*pi);
178      bim->set(r,c,temp);
179      n_Delete(&temp,coeffs_BIGINT);
180      pi++;
181    }
182  return bim;
183}
184
185lists PmIncidenceMatrix2ListOfIntvecs (polymake::IncidenceMatrix<polymake::NonSymmetric>* icmat)
186{
187  int rows = icmat->rows();
188  int cols = icmat->cols();
189  lists L = (lists)omAllocBin(slists_bin);
190  L->Init(rows);
191
192  for (int r = 0; r < rows; r++)
193  {
194    intvec* iv = new intvec(cols); int i=0;
195    for (int c = 0; c < cols; c++)
196    {
197      if ((*icmat).row(r).exists(c))
198      { (*iv)[i]=c; i++; }
199    }
200    iv->resize(i);
201    L->m[r].rtyp = INTVEC_CMD;
202    L->m[r].data = (void*) iv;
203  }
204
205  return L;
206}
207
208lists PmAdjacencyMatrix2ListOfEdges (polymake::IncidenceMatrix<polymake::NonSymmetric>* icmat)
209{
210  int rows = icmat->rows();
211  int cols = icmat->cols();
212
213  // counting number of edges
214  int i=0; int r, c;
215  for (r=0; r<rows; r++)
216  {
217    for (c=0; c<cols; c++)
218    {
219      if ((*icmat).row(r).exists(c) && r<c)
220        i++;
221    }
222  }
223
224  lists L = (lists)omAllocBin(slists_bin);
225  L->Init(i);
226
227  i=0;
228  for (r=0; r<rows; r++)
229  {
230    for (c=0; c<cols; c++)
231    {
232      if ((*icmat).row(r).exists(c) && r<c)
233      {
234        intvec* iv = new intvec(2);
235        (*iv)[0]=r; (*iv)[1]=c;
236        L->m[i].rtyp = INTVEC_CMD;
237        L->m[i].data = (void*) iv;
238        i++;
239      }
240    }
241  }
242
243  return L;
244}
245
246intvec* PmSetInteger2Intvec (polymake::Set<polymake::Integer>* si, bool &b)
247{
248  polymake::Vector<polymake::Integer> vi(*si);
249  return PmVectorInteger2Intvec(&vi,b);
250}
251
252/* polymake <- singular */
253
254polymake::Matrix<polymake::Integer> Intvec2PmMatrixInteger (const intvec* im)
255{
256  int rows=im->rows();
257  int cols=im->cols();
258  polymake::Matrix<polymake::Integer> mi(rows,cols);
259  for(int r=0; r<rows; r++)
260    for(int c=0; c<cols; c++)
261      mi(r,c) = polymake::Integer(IMATELEM(*im, r+1, c+1));
262  return mi;
263}
264
265/* Functions for converting cones and fans in between gfan and polymake,
266   Singular shares the same cones and fans with gfan */
267
268gfan::ZCone* PmCone2ZCone (polymake::perl::Object* pc)
269{
270  if (pc->isa("Cone"))
271  {
272    polymake::Integer ambientdim1 = pc->give("CONE_AMBIENT_DIM");
273    bool ok=true; int ambientdim2 = PmInteger2Int(ambientdim1, ok);
274    if (!ok)
275    {
276      WerrorS("PmCone2ZCone: overflow while converting polymake::Integer to int");
277    }
278    polymake::Matrix<polymake::Rational> ineqrational = pc->give("FACETS");
279    polymake::Matrix<polymake::Rational> eqrational = pc->give("LINEAR_SPAN");
280    // polymake::Matrix<polymake::Rational> exraysrational = pc->give("RAYS");
281    // polymake::Matrix<polymake::Rational> linrational = pc->give("LINEALITY_SPACE");
282
283    gfan::ZMatrix zv, zw, zx, zy, zz;
284    // the following branching statements are to cover cases in which polymake returns empty matrices
285    // by convention, gfanlib ignores empty matrices, hence zero matrices of right dimensions have to be supplied
286    if (ineqrational.cols()!=0)
287    {
288      polymake::Matrix<polymake::Integer> ineqinteger = polymake::common::primitive(ineqrational);
289      zv = PmMatrixInteger2GfZMatrix(&ineqinteger);
290    }
291    else
292      zv = gfan::ZMatrix(0, ambientdim2);
293    if (eqrational.cols()!=0)
294    {
295      polymake::Matrix<polymake::Integer> eqinteger = polymake::common::primitive(eqrational);
296      zw = PmMatrixInteger2GfZMatrix(&eqinteger);
297    }
298    else
299      zw = gfan::ZMatrix(0, ambientdim2);
300    // if (exraysrational.cols()!=0)
301    // {
302    //   polymake::Matrix<polymake::Integer> exraysinteger = polymake::common::primitive(exraysrational);
303    //   zx = PmMatrixInteger2GfZMatrix(&exraysinteger);
304    // }
305    // else
306    //   zx = gfan::ZMatrix(0, ambientdim2);
307    // if (linrational.cols()!=0)
308    // {
309    //   polymake::Matrix<polymake::Integer> lininteger = polymake::common::primitive(linrational);
310    //   zy = PmMatrixInteger2GfZMatrix(&lininteger);
311    // }
312    // else
313    //   zy = gfan::ZMatrix(0, ambientdim2);
314
315    // gfan::ZCone* zc = new gfan::ZCone(zv,zw,zx,zy,zz,3);
316    gfan::ZCone* zc = new gfan::ZCone(zv,zw,3);
317    return zc;
318  }
319  WerrorS("PmCone2ZCone: unexpected parameters");
320  return NULL;
321}
322
323gfan::ZCone* PmPolytope2ZPolytope (polymake::perl::Object* pp)
324{
325  if (pp->isa("Polytope<Rational>"))
326  {
327    polymake::Integer ambientdim1 = pp->give("CONE_AMBIENT_DIM");
328    bool ok=true; int ambientdim2 = PmInteger2Int(ambientdim1, ok);
329    if (!ok)
330    {
331      WerrorS("overflow while converting polymake::Integer to int");
332    }
333    polymake::Matrix<polymake::Rational> ineqrational = pp->give("FACETS");
334    polymake::Matrix<polymake::Rational> eqrational = pp->give("AFFINE_HULL");
335    // polymake::Matrix<polymake::Rational> vertrational = pp->give("VERTICES");
336    // polymake::Matrix<polymake::Rational> linrational = pp->give("LINEALITY_SPACE");
337
338    gfan::ZMatrix zv, zw;
339    // the following branching statements are to cover the cases when polymake returns empty matrices
340    // by convention, gfanlib ignores empty matrices, hence zero matrices of right dimensions have to be supplied
341    if (ineqrational.cols()!=0)
342    {
343      polymake::Matrix<polymake::Integer> ineqinteger = polymake::common::primitive(ineqrational);
344      zv = PmMatrixInteger2GfZMatrix(&ineqinteger);
345    }
346    else
347      zv = gfan::ZMatrix(0, ambientdim2);
348
349    if (eqrational.cols()!=0)
350    {
351      polymake::Matrix<polymake::Integer> eqinteger = polymake::common::primitive(eqrational);
352      zw = PmMatrixInteger2GfZMatrix(&eqinteger);
353    }
354    else
355      zw = gfan::ZMatrix(0, ambientdim2);
356
357    // if (vertrational.cols()!=0)
358    // {
359    //   polymake::Matrix<polymake::Integer> vertinteger = polymake::common::primitive(vertrational);
360    //   zx = PmMatrixInteger2GfZMatrix(&vertinteger);
361    // }
362    // else
363    //   zx = gfan::ZMatrix(0, ambientdim2);
364    // if (linrational.cols()!=0)
365    //   {
366    //     polymake::Matrix<polymake::Integer> lininteger = polymake::common::primitive(linrational);
367    //     zy = PmMatrixInteger2GfZMatrix(&lininteger);
368    //   }
369    // else
370    //   zy = gfan::ZMatrix(0, ambientdim2);
371
372    // gfan::ZCone* zp = new gfan::ZCone(zv,zw,zx,zy,zz,3);
373    gfan::ZCone* zp = new gfan::ZCone(zv,zw,3);
374
375    return zp;
376  }
377  WerrorS("PmPolytope2ZPolytope: unexpected parameters");
378  return NULL;
379}
380
381gfan::ZFan* PmFan2ZFan (polymake::perl::Object* pf)
382{
383  if (pf->isa("PolyhedralFan"))
384  {
385    int d = (int) pf->give("FAN_AMBIENT_DIM");
386    gfan::ZFan* zf = new gfan::ZFan(d);
387
388    int n = pf->give("N_MAXIMAL_CONES");
389    for (int i=0; i<n; i++)
390      {
391        polymake::perl::Object pmcone=pf->CallPolymakeMethod("cone",i);
392        gfan::ZCone* zc=PmCone2ZCone(&pmcone);
393        zf->insert(*zc);
394      }
395    return zf;
396  }
397  WerrorS("PmFan2ZFan: unexpected parameters");
398  return NULL;
399}
400
401polymake::perl::Object* ZCone2PmCone (gfan::ZCone* zc)
402{
403  polymake::perl::Object* gc = new polymake::perl::Object("Cone<Rational>");
404
405  gfan::ZMatrix inequalities = zc->getInequalities();
406  gc->take("FACETS") << GfZMatrix2PmMatrixInteger(&inequalities);
407
408  gfan::ZMatrix equations = zc->getEquations();
409  gc->take("LINEAR_SPAN") << GfZMatrix2PmMatrixInteger(&equations);
410
411  // if(zc->areExtremeRaysKnown())
412  //   {
413  //     gfan::ZMatrix extremeRays = zc->extremeRays();
414  //     gc->take("RAYS") << GfZMatrix2PmMatrixInteger(&extremeRays);
415  //   }
416
417  // if(zc->areGeneratorsOfLinealitySpaceKnown())
418  //   {
419  //     gfan::ZMatrix lineality = zc->generatorsOfLinealitySpace();
420  //     gc->take("LINEALITY_SPACE") << GfZMatrix2PmMatrixInteger(&lineality);
421  //   }
422
423  return gc;
424}
425
426polymake::perl::Object* ZPolytope2PmPolytope (gfan::ZCone* zc)
427{
428  polymake::perl::Object* pp = new polymake::perl::Object("Polytope<Rational>");
429
430  gfan::ZMatrix inequalities = zc->getInequalities();
431  pp->take("FACETS") << GfZMatrix2PmMatrixInteger(&inequalities);
432
433  gfan::ZMatrix equations = zc->getEquations();
434  pp->take("LINEAR_SPAN") << GfZMatrix2PmMatrixInteger(&equations);
435
436  // if(zc->areExtremeRaysKnown())
437  //   {
438  //     gfan::ZMatrix vertices = zc->extremeRays();
439  //     pp->take("VERTICES") << GfZMatrix2PmMatrixInteger(&vertices);
440  //   }
441
442  return pp;
443}
444
445polymake::Matrix<polymake::Integer> raysOf(gfan::ZFan* zf)
446{
447  int d = zf->getAmbientDimension();
448  int n = zf->numberOfConesOfDimension(1,0,0);
449  gfan::ZMatrix zm(n,d);
450
451  for (int i=0; i<n; i++)
452    {
453      gfan::ZCone zc = zf->getCone(1,i,0,0);
454      gfan::ZMatrix ray = zc.extremeRays();
455      for (int j=0; j<d; j++)
456        {
457          zm[i][j]=ray[0][j];
458        }
459    }
460
461  return GfZMatrix2PmMatrixInteger(&zm);
462}
463
464int numberOfRaysOf(gfan::ZFan* zf)
465{
466  int n = zf->numberOfConesOfDimension(1,0,0);
467  return n;
468}
469
470int numberOfMaximalConesOf(gfan::ZFan* zf)
471{
472  int d = zf->getAmbientDimension();
473  int n = 0;
474
475  for (int i=0; i<=d; i++)
476    {
477      n = n + zf->numberOfConesOfDimension(i,0,1);
478    }
479
480  return n;
481}
482
483polymake::Array<polymake::Set<int> > conesOf(gfan::ZFan* zf)
484{
485  int r = numberOfMaximalConesOf(zf);
486
487  polymake::Matrix<polymake::Integer> pm=raysOf(zf);
488  polymake::Array<polymake::Set<int> > L(r);
489
490  int ii = 0;
491  for (int d=1; d<=zf->getAmbientDimension(); d++)
492    {
493      for (int i=0; i<zf->numberOfConesOfDimension(d,0,1); i++)
494        {
495          gfan::IntVector v = zf->getConeIndices(d,i,0,1);
496          polymake::Set<int> s;
497          for (int j=0; j<(int)v.size(); j++)
498            {
499              s = s+v[j];
500            }
501          L[ii] = s;
502          ii = ii + 1;
503        }
504    }
505  return L;
506}
507
508polymake::perl::Object* ZFan2PmFan (gfan::ZFan* zf)
509{
510  polymake::perl::Object* pf = new polymake::perl::Object("PolyhedralFan");
511
512  polymake::Matrix<polymake::Integer> zm = raysOf(zf);
513  pf->take("RAYS") << zm;  // using rays here instead of INPUT_RAYS prevents redundant computations
514
515  polymake::Array<polymake::Set<int> > ar = conesOf(zf);
516  pf->take("MAXIMAL_CONES") << ar;
517
518  return pf;
519}
Note: See TracBrowser for help on using the repository browser.