volume.c

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