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