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