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