volume.c

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