volume.c

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