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 /* Substitute parameters by the corresponding element in subs
 110  */
 111 static evalue *evalue_substitute_new(evalue *e, evalue **subs)
 112 {
 113     evalue *res = NULL;
 114     evalue *c;
 115     int i;
 116
 117     if (value_notzero_p(e->d)) {
 118         res = ALLOC(evalue);
 119         value_init(res->d);
 120         evalue_copy(res, e);
 121         return res;
 122     }
 123     assert(e->x.p->type == polynomial);
 124
 125     res = evalue_zero();
 126     for (i = e->x.p->size-1; i > 0; --i) {
 127         c = evalue_substitute_new(&e->x.p->arr[i], subs);
 128         eadd(c, res);
 129         free_evalue_refs(c);
 130         free(c);
 131         emul(subs[e->x.p->pos-1], res);
 132     }
 133     c = evalue_substitute_new(&e->x.p->arr[0], subs);
 134     eadd(c, res);
 135     free_evalue_refs(c);
 136     free(c);
 137
 138     return res;
 139 }
 140
 141 struct parameter_point {
 142     Vector *coord;
 143     evalue **e;
 144 };
 145
 146 struct parameter_point *parameter_point_new(unsigned nparam)
 147 {
 148     struct parameter_point *point = ALLOC(struct parameter_point);
 149     point->coord = Vector_Alloc(nparam+1);
 150     point->e = NULL;
 151     return point;
 152 }
 153
 154 evalue **parameter_point_evalue(struct parameter_point *point)
 155 {
 156     int j;
 157     unsigned nparam = point->coord->Size-1;
 158
 159     if (point->e)
 160         return point->e;
 161
 162     point->e = ALLOCN(evalue *, nparam);
 163     for (j = 0; j < nparam; ++j) {
 164         point->e[j] = ALLOC(evalue);
 165         value_init(point->e[j]->d);
 166         evalue_set(point->e[j], point->coord->p[j], point->coord->p[nparam]);
 167     }
 168
 169     return point->e;
 170 }
 171
 172 void parameter_point_free(struct parameter_point *point)
 173 {
 174     int i;
 175     unsigned nparam = point->coord->Size-1;
 176
 177     Vector_Free(point->coord);
 178
 179     if (point->e) {
 180         for (i = 0; i < nparam; ++i) {
 181             free_evalue_refs(point->e[i]);
 182             free(point->e[i]);
 183         }
 184         free(point->e);
 185     }
 186     free(point);
 187 }
 188
 189 /* Computes point in pameter space where polyhedron is non-empty.
 190  */
 191 static struct parameter_point *non_empty_point(Param_Domain *D)
 192 {
 193     unsigned nparam = D->Domain->Dimension;
 194     struct parameter_point *point;
 195     Vector *v;
 196
 197     v = inner_point(D->Domain);
 198     point = parameter_point_new(nparam);
 199     Vector_Copy(v->p+1, point->coord->p, nparam+1);
 200     Vector_Free(v);
 201
 202     return point;
 203 }
 204
 205 static Matrix *barycenter(Param_Polyhedron *PP, Param_Domain *D)
 206 {
 207     int nbV;
 208     Matrix *center = NULL;
 209     Value denom;
 210     Value fc, fv;
 211     unsigned nparam;
 212     int i;
 213     Param_Vertices *V;
 214
 215     value_init(fc);
 216     value_init(fv);
 217     nbV = 0;
 218     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 219         ++nbV;
 220         if (!center) {
 221             center = Matrix_Copy(V->Vertex);
 222             nparam = center->NbColumns - 2;
 223         } else {
 224             for (i = 0; i < center->NbRows; ++i) {
 225                 value_assign(fc, center->p[i][nparam+1]);
 226                 value_lcm(fc, V->Vertex->p[i][nparam+1],
 227                             &center->p[i][nparam+1]);
 228                 value_division(fc, center->p[i][nparam+1], fc);
 229                 value_division(fv, center->p[i][nparam+1],
 230                                 V->Vertex->p[i][nparam+1]);
 231                 Vector_Combine(center->p[i], V->Vertex->p[i], center->p[i],
 232                                fc, fv, nparam+1);
 233             }
 234         }
 235     END_FORALL_PVertex_in_ParamPolyhedron;
 236     value_clear(fc);
 237     value_clear(fv);
 238
 239     value_init(denom);
 240     value_set_si(denom, nbV);
 241     for (i = 0; i < center->NbRows; ++i) {
 242         value_multiply(center->p[i][nparam+1], center->p[i][nparam+1], denom);
 243         Vector_Normalize(center->p[i], nparam+2);
 244     }
 245     value_clear(denom);
 246
 247     return center;
 248 }
 249
 250 static Matrix *triangulation_vertex(Param_Polyhedron *PP, Param_Domain *D,
 251                                     Polyhedron *F)
 252 {
 253     Param_Vertices *V;
 254
 255     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 256         return V->Vertex;
 257     END_FORALL_PVertex_in_ParamPolyhedron;
 258
 259     assert(0);
 260     return NULL;
 261 }
 262
 263 /* Compute dim! times the volume of polyhedron F in Param_Domain D.
 264  * If F is a simplex, then the volume is computed of a recursive pyramid
 265  * over F with the points already in matrix.
 266  * Otherwise, the barycenter of F is added to matrix and the function
 267  * is called recursively on the facets of F.
 268  *
 269  * The first row of matrix contain the _negative_ of the first point.
 270  * The remaining rows of matrix contain the distance of the corresponding
 271  * point to the first point.
 272  */
 273 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 274                                 unsigned dim, evalue ***matrix,
 275                                 struct parameter_point *point,
 276                                 int row, Polyhedron *F,
 277                                 struct barvinok_options *options);
 278
 279 static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
 280                                   unsigned dim, evalue ***matrix,
 281                                   struct parameter_point *point,
 282                                   int row, Polyhedron *F,
 283                                   struct barvinok_options *options)
 284 {
 285     int j;
 286     evalue *tmp;
 287     evalue *vol;
 288     evalue mone;
 289     Matrix *center;
 290     unsigned cut_MaxRays = options->MaxRays;
 291     unsigned nparam = PP->V->Vertex->NbColumns-2;
 292     Vector *v = NULL;
 293
 294     POL_UNSET(cut_MaxRays, POL_INTEGER);
 295
 296     value_init(mone.d);
 297     evalue_set_si(&mone, -1, 1);
 298
 299     if (options->volume_triangulate == BV_VOL_BARYCENTER)
 300         center = barycenter(PP, D);
 301     else
 302         center = triangulation_vertex(PP, D, F);
 303     for (j = 0; j < dim; ++j)
 304         matrix[row][j] = vertex2evalue(center->p[j], center->NbColumns - 2);
 305     if (options->volume_triangulate == BV_VOL_BARYCENTER)
 306         Matrix_Free(center);
 307     else
 308         v = Vector_Alloc(1+nparam+1);
 309
 310     if (row == 0) {
 311         for (j = 0; j < dim; ++j)
 312             emul(&mone, matrix[row][j]);
 313     } else {
 314         for (j = 0; j < dim; ++j)
 315             eadd(matrix[0][j], matrix[row][j]);
 316     }
 317
 318     vol = NULL;
 319     POL_ENSURE_FACETS(F);
 320     for (j = F->NbEq; j < F->NbConstraints; ++j) {
 321         Polyhedron *FF;
 322         Param_Domain *FD;
 323         if (First_Non_Zero(F->Constraint[j]+1, dim) == -1)
 324             continue;
 325         if (options->volume_triangulate != BV_VOL_BARYCENTER) {
 326             Param_Inner_Product(F->Constraint[j], center, v->p);
 327             if (First_Non_Zero(v->p+1, nparam+1) == -1)
 328                 continue;
 329         }
 330         FF = facet(F, j, options->MaxRays);
 331         FD = Param_Polyhedron_Facet(PP, D, F, j);
 332         tmp = volume_in_domain(PP, FD, dim, matrix, point,
 333                                row+1, FF, options);
 334         if (!vol)
 335             vol = tmp;
 336         else {
 337             eadd(tmp, vol);
 338             free_evalue_refs(tmp);
 339             free(tmp);
 340         }
 341         Polyhedron_Free(FF);
 342         Param_Domain_Free(FD);
 343     }
 344
 345     if (options->volume_triangulate != BV_VOL_BARYCENTER)
 346         Vector_Free(v);
 347
 348     for (j = 0; j < dim; ++j) {
 349         free_evalue_refs(matrix[row][j]);
 350         free(matrix[row][j]);
 351     }
 352
 353     free_evalue_refs(&mone);
 354     return vol;
 355 }
 356
 357 static evalue *volume_simplex(Param_Polyhedron *PP, Param_Domain *D,
 358                                 unsigned dim, evalue ***matrix,
 359                                 struct parameter_point *point,
 360                                 int row, struct barvinok_options *options)
 361 {
 362     evalue mone;
 363     Param_Vertices *V;
 364     evalue *vol, *val;
 365     int i, j;
 366
 367     options->stats->volume_simplices++;
 368
 369     value_init(mone.d);
 370     evalue_set_si(&mone, -1, 1);
 371
 372     i = row;
 373     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal counters */
 374         for (j = 0; j < dim; ++j) {
 375             matrix[i][j] = vertex2evalue(V->Vertex->p[j],
 376                                            V->Vertex->NbColumns - 2);
 377             if (i == 0)
 378                 emul(&mone, matrix[i][j]);
 379             else
 380                 eadd(matrix[0][j], matrix[i][j]);
 381         }
 382         ++i;
 383     END_FORALL_PVertex_in_ParamPolyhedron;
 384
 385     vol = determinant(matrix+1, dim);
 386
 387     val = evalue_substitute_new(vol, parameter_point_evalue(point));
 388
 389     assert(value_notzero_p(val->d));
 390     assert(value_notzero_p(val->x.n));
 391     if (value_neg_p(val->x.n))
 392         emul(&mone, vol);
 393
 394     free_evalue_refs(val);
 395     free(val);
 396
 397     for (i = row; i < dim+1; ++i) {
 398         for (j = 0; j < dim; ++j) {
 399             free_evalue_refs(matrix[i][j]);
 400             free(matrix[i][j]);
 401         }
 402     }
 403
 404     free_evalue_refs(&mone);
 405
 406     return vol;
 407 }
 408
 409 static evalue *volume_triangulate_lift(Param_Polyhedron *PP, Param_Domain *D,
 410                                         unsigned dim, evalue ***matrix,
 411                                         struct parameter_point *point,
 412                                         int row, struct barvinok_options *options)
 413 {
 414     const static int MAX_TRY=10;
 415     Param_Vertices *V;
 416     int nbV, nv;
 417     int i;
 418     int t = 0;
 419     Matrix *FixedRays, *M;
 420     Polyhedron *L;
 421     Param_Domain SD;
 422     Value tmp;
 423     evalue *s, *vol;
 424
 425     SD.Domain = 0;
 426     SD.next = 0;
 427     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 428     SD.F = ALLOCN(unsigned, nv);
 429
 430     FixedRays = Matrix_Alloc(PP->nbV+1, 1+dim+2);
 431     nbV = 0;
 432     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 433         unsigned nparam = V->Vertex->NbColumns - 2;
 434         Param_Vertex_Common_Denominator(V);
 435         for (i = 0; i < V->Vertex->NbRows; ++i) {
 436             value_multiply(FixedRays->p[nbV][1+i], V->Vertex->p[i][nparam],
 437                            point->coord->p[nparam]);
 438             Inner_Product(V->Vertex->p[i], point->coord->p, nparam,
 439                           &FixedRays->p[nbV][1+dim]);
 440             value_addto(FixedRays->p[nbV][1+i], FixedRays->p[nbV][1+i],
 441                         FixedRays->p[nbV][1+dim]);
 442         }
 443         value_multiply(FixedRays->p[nbV][1+dim+1], V->Vertex->p[0][nparam+1],
 444                        point->coord->p[nparam]);
 445         value_set_si(FixedRays->p[nbV][0], 1);
 446         ++nbV;
 447     END_FORALL_PVertex_in_ParamPolyhedron;
 448     value_set_si(FixedRays->p[nbV][0], 1);
 449     value_set_si(FixedRays->p[nbV][1+dim], 1);
 450     FixedRays->NbRows = nbV+1;
 451
 452     value_init(tmp);
 453     if (0) {
 454 try_again:
 455         /* Usually vol should still be NULL */
 456         if (vol) {
 457             free_evalue_refs(vol);
 458             free(vol);
 459         }
 460         Polyhedron_Free(L);
 461         ++t;
 462     }
 463     assert(t <= MAX_TRY);
 464     vol = NULL;
 465
 466     for (i = 0; i < nbV; ++i)
 467         value_set_si(FixedRays->p[i][1+dim], random_int((t+1)*dim*nbV)+1);
 468
 469     M = Matrix_Copy(FixedRays);
 470     L = Rays2Polyhedron(M, options->MaxRays);
 471     Matrix_Free(M);
 472
 473     POL_ENSURE_FACETS(L);
 474     for (i = 0; i < L->NbConstraints; ++i) {
 475         int r;
 476         /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
 477         if (value_negz_p(L->Constraint[i][1+dim]))
 478             continue;
 479
 480         memset(SD.F, 0, nv * sizeof(unsigned));
 481         nbV = 0;
 482         r = 0;
 483         FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal */
 484             Inner_Product(FixedRays->p[nbV]+1, L->Constraint[i]+1, dim+2, &tmp);
 485             if (value_zero_p(tmp)) {
 486                 if (r > dim-row)
 487                     goto try_again;
 488                 SD.F[_ix] |= _bx;
 489                 ++r;
 490             }
 491             ++nbV;
 492         END_FORALL_PVertex_in_ParamPolyhedron;
 493         assert(r == (dim-row)+1);
 494
 495         s = volume_simplex(PP, &SD, dim, matrix, point, row, options);
 496         if (!vol)
 497             vol = s;
 498         else {
 499             eadd(s, vol);
 500             free_evalue_refs(s);
 501             free(s);
 502         }
 503     }
 504     Polyhedron_Free(L);
 505     Matrix_Free(FixedRays);
 506     value_clear(tmp);
 507     free(SD.F);
 508
 509     return vol;
 510 }
 511
 512 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 513                                 unsigned dim, evalue ***matrix,
 514                                 struct parameter_point *point,
 515                                 int row, Polyhedron *F,
 516                                 struct barvinok_options *options)
 517 {
 518     int nbV;
 519     Param_Vertices *V;
 520     evalue *vol;
 521
 522     assert(point);
 523
 524     nbV = 0;
 525     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 526         ++nbV;
 527     END_FORALL_PVertex_in_ParamPolyhedron;
 528
 529     if (nbV > (dim-row) + 1) {
 530         if (options->volume_triangulate == BV_VOL_LIFT)
 531             vol = volume_triangulate_lift(PP, D, dim, matrix, point,
 532                                           row, options);
 533         else
 534             vol = volume_triangulate(PP, D, dim, matrix, point,
 535                                      row, F, options);
 536     } else {
 537         assert(nbV == (dim-row) + 1);
 538         vol = volume_simplex(PP, D, dim, matrix, point, row, options);
 539     }
 540
 541     return vol;
 542 }
 543
 544 evalue* Param_Polyhedron_Volume(Polyhedron *P, Polyhedron* C,
 545                                 struct barvinok_options *options)
 546 {
 547     evalue ***matrix;
 548     unsigned nparam = C->Dimension;
 549     unsigned nvar = P->Dimension - C->Dimension;
 550     Param_Polyhedron *PP;
 551     unsigned MaxRays;
 552     int i, j;
 553     Value fact;
 554     evalue *vol;
 555     int nd;
 556     struct evalue_section *s;
 557     Param_Domain *D;
 558     Polyhedron *TC;
 559
 560     if (options->polynomial_approximation == BV_APPROX_SIGN_NONE)
 561         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 562
 563     if (options->polynomial_approximation != BV_APPROX_SIGN_APPROX) {
 564         int pa = options->polynomial_approximation;
 565         assert(pa == BV_APPROX_SIGN_UPPER || pa == BV_APPROX_SIGN_LOWER);
 566
 567         P = Polyhedron_Flate(P, nparam, pa == BV_APPROX_SIGN_UPPER,
 568                              options->MaxRays);
 569
 570         /* Don't deflate/inflate again (on this polytope) */
 571         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 572         vol = barvinok_enumerate_with_options(P, C, options);
 573         options->polynomial_approximation = pa;
 574
 575         Polyhedron_Free(P);
 576         return vol;
 577     }
 578
 579     TC = true_context(P, C, options->MaxRays);
 580
 581     MaxRays = options->MaxRays;
 582     POL_UNSET(options->MaxRays, POL_INTEGER);
 583
 584     value_init(fact);
 585     Factorial(nvar, &fact);
 586
 587     PP = Polyhedron2Param_Polyhedron(P, C, options);
 588
 589     for (nd = 0, D = PP->D; D; ++nd, D = D->next);
 590     s = ALLOCN(struct evalue_section, nd);
 591
 592     matrix = ALLOCN(evalue **, nvar+1);
 593     for (i = 0; i < nvar+1; ++i)
 594         matrix[i] = ALLOCN(evalue *, nvar);
 595
 596     FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
 597         Polyhedron *CA, *F;
 598         struct parameter_point *point;
 599
 600         CA = align_context(D->Domain, P->Dimension, MaxRays);
 601         F = DomainIntersection(P, CA, options->MaxRays);
 602         Domain_Free(CA);
 603
 604         point = non_empty_point(D);
 605         s[i].D = rVD;
 606         s[i].E = volume_in_domain(PP, D, nvar, matrix, point, 0, F, options);
 607         Domain_Free(F);
 608         parameter_point_free(point);
 609         evalue_div(s[i].E, fact);
 610     END_FORALL_REDUCED_DOMAIN
 611     options->MaxRays = MaxRays;
 612     Polyhedron_Free(TC);
 613
 614     vol = evalue_from_section_array(s, nd);
 615     free(s);
 616
 617     for (i = 0; i < nvar+1; ++i)
 618         free(matrix[i]);
 619     free(matrix);
 620     Param_Polyhedron_Free(PP);
 621     value_clear(fact);
 622
 623     return vol;
 624 }