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