volume.c

   1 #include <assert.h>
   2 #include <barvinok/polylib.h>
   3 #include <barvinok/barvinok.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/util.h>
   6 #include "reduce_domain.h"
   7 #include "param_util.h"
   8 #include "volume.h"
   9 #include "scale.h"
  10
  11 #define ALLOC(type) (type*)malloc(sizeof(type))
  12 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  13
  14 /* Computes an evalue representation of a coordinate
  15  * in a ParamVertices.
  16  */
  17 static evalue *vertex2evalue(Value *vertex, int nparam)
  18 {
  19     return affine2evalue(vertex, vertex[nparam+1], nparam);
  20 }
  21
  22 static void matrix_print(evalue ***matrix, int dim, int *cols,
  23                          const char * const *param_names)
  24 {
  25     int i, j;
  26
  27     for (i = 0; i < dim; ++i)
  28         for (j = 0; j < dim; ++j) {
  29             int k = cols ? cols[j] : j;
  30             fprintf(stderr, "%d %d: ", i, j);
  31             print_evalue(stderr, matrix[i][k], param_names);
  32             fprintf(stderr, "\n");
  33         }
  34 }
  35
  36 /* Compute determinant using Laplace's formula.
  37  * In particular, the determinant is expanded along the last row.
  38  * The cols array is a list of columns that remain in the currect submatrix.
  39  */
  40 static evalue *determinant_cols(evalue ***matrix, int dim, int *cols)
  41 {
  42     evalue *det, *tmp;
  43     evalue mone;
  44     int j;
  45     int *newcols;
  46
  47     if (dim == 1) {
  48         det = ALLOC(evalue);
  49         value_init(det->d);
  50         evalue_copy(det, matrix[0][cols[0]]);
  51         return det;
  52     }
  53
  54     value_init(mone.d);
  55     evalue_set_si(&mone, -1, 1);
  56     det = NULL;
  57     newcols = ALLOCN(int, dim-1);
  58     for (j = 1; j < dim; ++j)
  59         newcols[j-1] = cols[j];
  60     for (j = 0; j < dim; ++j) {
  61         if (j != 0)
  62             newcols[j-1] = cols[j-1];
  63         tmp = determinant_cols(matrix, dim-1, newcols);
  64         emul(matrix[dim-1][cols[j]], tmp);
  65         if ((dim+j)%2 == 0)
  66             emul(&mone, tmp);
  67         if (!det)
  68             det = tmp;
  69         else {
  70             eadd(tmp, det);
  71             free_evalue_refs(tmp);
  72             free(tmp);
  73         }
  74     }
  75     free(newcols);
  76     free_evalue_refs(&mone);
  77
  78     return det;
  79 }
  80
  81 static evalue *determinant(evalue ***matrix, int dim)
  82 {
  83     int i;
  84     int *cols = ALLOCN(int, dim);
  85     evalue *det;
  86
  87     for (i = 0; i < dim; ++i)
  88         cols[i] = i;
  89
  90     det = determinant_cols(matrix, dim, cols);
  91
  92     free(cols);
  93
  94     return det;
  95 }
  96
  97 /* Compute the facet of P that saturates constraint c.
  98  */
  99 static Polyhedron *facet(Polyhedron *P, int c, unsigned MaxRays)
 100 {
 101     Polyhedron *F;
 102     Vector *row = Vector_Alloc(1+P->Dimension+1);
 103     Vector_Copy(P->Constraint[c]+1, row->p+1, P->Dimension+1);
 104     F = AddConstraints(row->p, 1, P, MaxRays);
 105     Vector_Free(row);
 106     return F;
 107 }
 108
 109 /* Compute a dummy Param_Domain that contains all vertices of Param_Domain D
 110  * (which contains the vertices of P) that lie on the facet obtained by
 111  * saturating constraint c of P
 112  */
 113 static Param_Domain *face_vertices(Param_Polyhedron *PP, Param_Domain *D,
 114                                    Polyhedron *P, int c)
 115 {
 116     int nv;
 117     Param_Vertices *V;
 118     unsigned nparam = PP->V->Vertex->NbColumns-2;
 119     Vector *row = Vector_Alloc(1+nparam+1);
 120     Param_Domain *FD = ALLOC(Param_Domain);
 121     FD->Domain = 0;
 122     FD->next = 0;
 123
 124     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 125     FD->F = ALLOCN(unsigned, nv);
 126     memset(FD->F, 0, nv * sizeof(unsigned));
 127
 128     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _i, _ix, _bx internal counters */
 129         int n;
 130         Param_Inner_Product(P->Constraint[c], V->Vertex, row->p);
 131         if (First_Non_Zero(row->p+1, nparam+1) == -1)
 132             FD->F[_ix] |= _bx;
 133     END_FORALL_PVertex_in_ParamPolyhedron;
 134
 135     Vector_Free(row);
 136
 137     return FD;
 138 }
 139
 140 /* Substitute parameters by the corresponding element in subs
 141  */
 142 static evalue *evalue_substitute_new(evalue *e, evalue **subs)
 143 {
 144     evalue *res = NULL;
 145     evalue *c;
 146     int i;
 147
 148     if (value_notzero_p(e->d)) {
 149         res = ALLOC(evalue);
 150         value_init(res->d);
 151         evalue_copy(res, e);
 152         return res;
 153     }
 154     assert(e->x.p->type == polynomial);
 155
 156     res = evalue_zero();
 157     for (i = e->x.p->size-1; i > 0; --i) {
 158         c = evalue_substitute_new(&e->x.p->arr[i], subs);
 159         eadd(c, res);
 160         free_evalue_refs(c);
 161         free(c);
 162         emul(subs[e->x.p->pos-1], res);
 163     }
 164     c = evalue_substitute_new(&e->x.p->arr[0], subs);
 165     eadd(c, res);
 166     free_evalue_refs(c);
 167     free(c);
 168
 169     return res;
 170 }
 171
 172 struct parameter_point {
 173     Vector *coord;
 174     evalue **e;
 175 };
 176
 177 struct parameter_point *parameter_point_new(unsigned nparam)
 178 {
 179     struct parameter_point *point = ALLOC(struct parameter_point);
 180     point->coord = Vector_Alloc(nparam+1);
 181     point->e = NULL;
 182     return point;
 183 }
 184
 185 evalue **parameter_point_evalue(struct parameter_point *point)
 186 {
 187     int j;
 188     unsigned nparam = point->coord->Size-1;
 189
 190     if (point->e)
 191         return point->e;
 192
 193     point->e = ALLOCN(evalue *, nparam);
 194     for (j = 0; j < nparam; ++j) {
 195         point->e[j] = ALLOC(evalue);
 196         value_init(point->e[j]->d);
 197         evalue_set(point->e[j], point->coord->p[j], point->coord->p[nparam]);
 198     }
 199
 200     return point->e;
 201 }
 202
 203 void parameter_point_free(struct parameter_point *point)
 204 {
 205     int i;
 206     unsigned nparam = point->coord->Size-1;
 207
 208     Vector_Free(point->coord);
 209
 210     if (point->e) {
 211         for (i = 0; i < nparam; ++i) {
 212             free_evalue_refs(point->e[i]);
 213             free(point->e[i]);
 214         }
 215         free(point->e);
 216     }
 217     free(point);
 218 }
 219
 220 /* Computes point in pameter space where polyhedron is non-empty.
 221  */
 222 static struct parameter_point *non_empty_point(Param_Domain *D)
 223 {
 224     unsigned nparam = D->Domain->Dimension;
 225     struct parameter_point *point;
 226     Vector *v;
 227
 228     v = inner_point(D->Domain);
 229     point = parameter_point_new(nparam);
 230     Vector_Copy(v->p+1, point->coord->p, nparam+1);
 231     Vector_Free(v);
 232
 233     return point;
 234 }
 235
 236 static Matrix *barycenter(Param_Polyhedron *PP, Param_Domain *D)
 237 {
 238     int nbV;
 239     Matrix *center = NULL;
 240     Value denom;
 241     Value fc, fv;
 242     unsigned nparam;
 243     int i;
 244     Param_Vertices *V;
 245
 246     value_init(fc);
 247     value_init(fv);
 248     nbV = 0;
 249     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 250         ++nbV;
 251         if (!center) {
 252             center = Matrix_Copy(V->Vertex);
 253             nparam = center->NbColumns - 2;
 254         } else {
 255             for (i = 0; i < center->NbRows; ++i) {
 256                 value_assign(fc, center->p[i][nparam+1]);
 257                 value_lcm(fc, V->Vertex->p[i][nparam+1],
 258                             &center->p[i][nparam+1]);
 259                 value_division(fc, center->p[i][nparam+1], fc);
 260                 value_division(fv, center->p[i][nparam+1],
 261                                 V->Vertex->p[i][nparam+1]);
 262                 Vector_Combine(center->p[i], V->Vertex->p[i], center->p[i],
 263                                fc, fv, nparam+1);
 264             }
 265         }
 266     END_FORALL_PVertex_in_ParamPolyhedron;
 267     value_clear(fc);
 268     value_clear(fv);
 269
 270     value_init(denom);
 271     value_set_si(denom, nbV);
 272     for (i = 0; i < center->NbRows; ++i) {
 273         value_multiply(center->p[i][nparam+1], center->p[i][nparam+1], denom);
 274         Vector_Normalize(center->p[i], nparam+2);
 275     }
 276     value_clear(denom);
 277
 278     return center;
 279 }
 280
 281 static Matrix *triangulation_vertex(Param_Polyhedron *PP, Param_Domain *D,
 282                                     Polyhedron *F)
 283 {
 284     Param_Vertices *V;
 285
 286     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 287         return V->Vertex;
 288     END_FORALL_PVertex_in_ParamPolyhedron;
 289
 290     assert(0);
 291     return NULL;
 292 }
 293
 294 /* Compute dim! times the volume of polyhedron F in Param_Domain D.
 295  * If F is a simplex, then the volume is computed of a recursive pyramid
 296  * over F with the points already in matrix.
 297  * Otherwise, the barycenter of F is added to matrix and the function
 298  * is called recursively on the facets of F.
 299  *
 300  * The first row of matrix contain the _negative_ of the first point.
 301  * The remaining rows of matrix contain the distance of the corresponding
 302  * point to the first point.
 303  */
 304 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 305                                 unsigned dim, evalue ***matrix,
 306                                 struct parameter_point *point,
 307                                 int row, Polyhedron *F,
 308                                 struct barvinok_options *options);
 309
 310 static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
 311                                   unsigned dim, evalue ***matrix,
 312                                   struct parameter_point *point,
 313                                   int row, Polyhedron *F,
 314                                   struct barvinok_options *options)
 315 {
 316     int j;
 317     evalue *tmp;
 318     evalue *vol;
 319     evalue mone;
 320     Matrix *center;
 321     unsigned cut_MaxRays = options->MaxRays;
 322     unsigned nparam = PP->V->Vertex->NbColumns-2;
 323     Vector *v = NULL;
 324
 325     POL_UNSET(cut_MaxRays, POL_INTEGER);
 326
 327     value_init(mone.d);
 328     evalue_set_si(&mone, -1, 1);
 329
 330     if (options->volume_triangulate == BV_VOL_BARYCENTER)
 331         center = barycenter(PP, D);
 332     else
 333         center = triangulation_vertex(PP, D, F);
 334     for (j = 0; j < dim; ++j)
 335         matrix[row][j] = vertex2evalue(center->p[j], center->NbColumns - 2);
 336     if (options->volume_triangulate == BV_VOL_BARYCENTER)
 337         Matrix_Free(center);
 338     else
 339         v = Vector_Alloc(1+nparam+1);
 340
 341     if (row == 0) {
 342         for (j = 0; j < dim; ++j)
 343             emul(&mone, matrix[row][j]);
 344     } else {
 345         for (j = 0; j < dim; ++j)
 346             eadd(matrix[0][j], matrix[row][j]);
 347     }
 348
 349     vol = NULL;
 350     POL_ENSURE_FACETS(F);
 351     for (j = F->NbEq; j < F->NbConstraints; ++j) {
 352         Polyhedron *FF;
 353         Param_Domain *FD;
 354         if (First_Non_Zero(F->Constraint[j]+1, dim) == -1)
 355             continue;
 356         if (options->volume_triangulate != BV_VOL_BARYCENTER) {
 357             Param_Inner_Product(F->Constraint[j], center, v->p);
 358             if (First_Non_Zero(v->p+1, nparam+1) == -1)
 359                 continue;
 360         }
 361         FF = facet(F, j, options->MaxRays);
 362         FD = face_vertices(PP, D, F, j);
 363         tmp = volume_in_domain(PP, FD, dim, matrix, point,
 364                                row+1, FF, options);
 365         if (!vol)
 366             vol = tmp;
 367         else {
 368             eadd(tmp, vol);
 369             free_evalue_refs(tmp);
 370             free(tmp);
 371         }
 372         Polyhedron_Free(FF);
 373         Param_Domain_Free(FD);
 374     }
 375
 376     if (options->volume_triangulate != BV_VOL_BARYCENTER)
 377         Vector_Free(v);
 378
 379     for (j = 0; j < dim; ++j) {
 380         free_evalue_refs(matrix[row][j]);
 381         free(matrix[row][j]);
 382     }
 383
 384     free_evalue_refs(&mone);
 385     return vol;
 386 }
 387
 388 static evalue *volume_simplex(Param_Polyhedron *PP, Param_Domain *D,
 389                                 unsigned dim, evalue ***matrix,
 390                                 struct parameter_point *point,
 391                                 int row, struct barvinok_options *options)
 392 {
 393     evalue mone;
 394     Param_Vertices *V;
 395     evalue *vol, *val;
 396     int i, j;
 397
 398     options->stats->volume_simplices++;
 399
 400     value_init(mone.d);
 401     evalue_set_si(&mone, -1, 1);
 402
 403     i = row;
 404     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal counters */
 405         for (j = 0; j < dim; ++j) {
 406             matrix[i][j] = vertex2evalue(V->Vertex->p[j],
 407                                            V->Vertex->NbColumns - 2);
 408             if (i == 0)
 409                 emul(&mone, matrix[i][j]);
 410             else
 411                 eadd(matrix[0][j], matrix[i][j]);
 412         }
 413         ++i;
 414     END_FORALL_PVertex_in_ParamPolyhedron;
 415
 416     vol = determinant(matrix+1, dim);
 417
 418     val = evalue_substitute_new(vol, parameter_point_evalue(point));
 419
 420     assert(value_notzero_p(val->d));
 421     assert(value_notzero_p(val->x.n));
 422     if (value_neg_p(val->x.n))
 423         emul(&mone, vol);
 424
 425     free_evalue_refs(val);
 426     free(val);
 427
 428     for (i = row; i < dim+1; ++i) {
 429         for (j = 0; j < dim; ++j) {
 430             free_evalue_refs(matrix[i][j]);
 431             free(matrix[i][j]);
 432         }
 433     }
 434
 435     free_evalue_refs(&mone);
 436
 437     return vol;
 438 }
 439
 440 static evalue *volume_triangulate_lift(Param_Polyhedron *PP, Param_Domain *D,
 441                                         unsigned dim, evalue ***matrix,
 442                                         struct parameter_point *point,
 443                                         int row, struct barvinok_options *options)
 444 {
 445     const static int MAX_TRY=10;
 446     Param_Vertices *V;
 447     int nbV, nv;
 448     int i;
 449     int t = 0;
 450     Matrix *FixedRays, *M;
 451     Polyhedron *L;
 452     Param_Domain SD;
 453     Value tmp;
 454     evalue *s, *vol;
 455
 456     SD.Domain = 0;
 457     SD.next = 0;
 458     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 459     SD.F = ALLOCN(unsigned, nv);
 460
 461     FixedRays = Matrix_Alloc(PP->nbV+1, 1+dim+2);
 462     nbV = 0;
 463     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 464         unsigned nparam = V->Vertex->NbColumns - 2;
 465         Param_Vertex_Common_Denominator(V);
 466         for (i = 0; i < V->Vertex->NbRows; ++i) {
 467             value_multiply(FixedRays->p[nbV][1+i], V->Vertex->p[i][nparam],
 468                            point->coord->p[nparam]);
 469             Inner_Product(V->Vertex->p[i], point->coord->p, nparam,
 470                           &FixedRays->p[nbV][1+dim]);
 471             value_addto(FixedRays->p[nbV][1+i], FixedRays->p[nbV][1+i],
 472                         FixedRays->p[nbV][1+dim]);
 473         }
 474         value_multiply(FixedRays->p[nbV][1+dim+1], V->Vertex->p[0][nparam+1],
 475                        point->coord->p[nparam]);
 476         value_set_si(FixedRays->p[nbV][0], 1);
 477         ++nbV;
 478     END_FORALL_PVertex_in_ParamPolyhedron;
 479     value_set_si(FixedRays->p[nbV][0], 1);
 480     value_set_si(FixedRays->p[nbV][1+dim], 1);
 481     FixedRays->NbRows = nbV+1;
 482
 483     value_init(tmp);
 484     if (0) {
 485 try_again:
 486         /* Usually vol should still be NULL */
 487         if (vol) {
 488             free_evalue_refs(vol);
 489             free(vol);
 490         }
 491         Polyhedron_Free(L);
 492         ++t;
 493     }
 494     assert(t <= MAX_TRY);
 495     vol = NULL;
 496
 497     for (i = 0; i < nbV; ++i)
 498         value_set_si(FixedRays->p[i][1+dim], random_int((t+1)*dim*nbV)+1);
 499
 500     M = Matrix_Copy(FixedRays);
 501     L = Rays2Polyhedron(M, options->MaxRays);
 502     Matrix_Free(M);
 503
 504     POL_ENSURE_FACETS(L);
 505     for (i = 0; i < L->NbConstraints; ++i) {
 506         int r;
 507         /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
 508         if (value_negz_p(L->Constraint[i][1+dim]))
 509             continue;
 510
 511         memset(SD.F, 0, nv * sizeof(unsigned));
 512         nbV = 0;
 513         r = 0;
 514         FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal */
 515             Inner_Product(FixedRays->p[nbV]+1, L->Constraint[i]+1, dim+2, &tmp);
 516             if (value_zero_p(tmp)) {
 517                 if (r > dim-row)
 518                     goto try_again;
 519                 SD.F[_ix] |= _bx;
 520                 ++r;
 521             }
 522             ++nbV;
 523         END_FORALL_PVertex_in_ParamPolyhedron;
 524         assert(r == (dim-row)+1);
 525
 526         s = volume_simplex(PP, &SD, dim, matrix, point, row, options);
 527         if (!vol)
 528             vol = s;
 529         else {
 530             eadd(s, vol);
 531             free_evalue_refs(s);
 532             free(s);
 533         }
 534     }
 535     Polyhedron_Free(L);
 536     Matrix_Free(FixedRays);
 537     value_clear(tmp);
 538     free(SD.F);
 539
 540     return vol;
 541 }
 542
 543 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 544                                 unsigned dim, evalue ***matrix,
 545                                 struct parameter_point *point,
 546                                 int row, Polyhedron *F,
 547                                 struct barvinok_options *options)
 548 {
 549     int nbV;
 550     Param_Vertices *V;
 551     evalue *vol;
 552
 553     assert(point);
 554
 555     nbV = 0;
 556     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 557         ++nbV;
 558     END_FORALL_PVertex_in_ParamPolyhedron;
 559
 560     if (nbV > (dim-row) + 1) {
 561         if (options->volume_triangulate == BV_VOL_LIFT)
 562             vol = volume_triangulate_lift(PP, D, dim, matrix, point,
 563                                           row, options);
 564         else
 565             vol = volume_triangulate(PP, D, dim, matrix, point,
 566                                      row, F, options);
 567     } else {
 568         assert(nbV == (dim-row) + 1);
 569         vol = volume_simplex(PP, D, dim, matrix, point, row, options);
 570     }
 571
 572     return vol;
 573 }
 574
 575 evalue* Param_Polyhedron_Volume(Polyhedron *P, Polyhedron* C,
 576                                 struct barvinok_options *options)
 577 {
 578     evalue ***matrix;
 579     unsigned nparam = C->Dimension;
 580     unsigned nvar = P->Dimension - C->Dimension;
 581     Param_Polyhedron *PP;
 582     unsigned MaxRays;
 583     int i, j;
 584     Value fact;
 585     evalue *vol;
 586     int nd;
 587     struct section { Polyhedron *D; evalue *E; } *s;
 588     Param_Domain *D;
 589     Polyhedron *TC;
 590
 591     if (options->polynomial_approximation == BV_APPROX_SIGN_NONE)
 592         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 593
 594     if (options->polynomial_approximation != BV_APPROX_SIGN_APPROX) {
 595         int pa = options->polynomial_approximation;
 596         assert(pa == BV_APPROX_SIGN_UPPER || pa == BV_APPROX_SIGN_LOWER);
 597
 598         P = Polyhedron_Flate(P, nparam, pa == BV_APPROX_SIGN_UPPER,
 599                              options->MaxRays);
 600
 601         /* Don't deflate/inflate again (on this polytope) */
 602         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 603         vol = barvinok_enumerate_with_options(P, C, options);
 604         options->polynomial_approximation = pa;
 605
 606         Polyhedron_Free(P);
 607         return vol;
 608     }
 609
 610     TC = true_context(P, C, options->MaxRays);
 611
 612     MaxRays = options->MaxRays;
 613     POL_UNSET(options->MaxRays, POL_INTEGER);
 614
 615     value_init(fact);
 616     Factorial(nvar, &fact);
 617
 618     PP = Polyhedron2Param_Polyhedron(P, C, options);
 619
 620     for (nd = 0, D = PP->D; D; ++nd, D = D->next);
 621     s = ALLOCN(struct section, nd);
 622
 623     matrix = ALLOCN(evalue **, nvar+1);
 624     for (i = 0; i < nvar+1; ++i)
 625         matrix[i] = ALLOCN(evalue *, nvar);
 626
 627     FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
 628         Polyhedron *CA, *F;
 629         struct parameter_point *point;
 630
 631         CA = align_context(D->Domain, P->Dimension, MaxRays);
 632         F = DomainIntersection(P, CA, options->MaxRays);
 633         Domain_Free(CA);
 634
 635         point = non_empty_point(D);
 636         s[i].D = rVD;
 637         s[i].E = volume_in_domain(PP, D, nvar, matrix, point, 0, F, options);
 638         Domain_Free(F);
 639         parameter_point_free(point);
 640         evalue_div(s[i].E, fact);
 641     END_FORALL_REDUCED_DOMAIN
 642     options->MaxRays = MaxRays;
 643     Polyhedron_Free(TC);
 644
 645     vol = ALLOC(evalue);
 646     value_init(vol->d);
 647     value_set_si(vol->d, 0);
 648
 649     if (nd == 0)
 650         evalue_set_si(vol, 0, 1);
 651     else {
 652         vol->x.p = new_enode(partition, 2*nd, C->Dimension);
 653         for (i = 0; i < nd; ++i) {
 654             EVALUE_SET_DOMAIN(vol->x.p->arr[2*i], s[i].D);
 655             value_clear(vol->x.p->arr[2*i+1].d);
 656             vol->x.p->arr[2*i+1] = *s[i].E;
 657             free(s[i].E);
 658         }
 659     }
 660     free(s);
 661
 662     for (i = 0; i < nvar+1; ++i)
 663         free(matrix[i]);
 664     free(matrix);
 665     Param_Polyhedron_Free(PP);
 666     value_clear(fact);
 667
 668     return vol;
 669 }