isl_union_map_apply_union_pw_qpolynomial: rename "dim" variable to "space"
[barvinok.git] / volume.c
blob6153108e1f56b9dc14c20cc4d8bdef2bc1e8b470
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"
12 #define ALLOC(type) (type*)malloc(sizeof(type))
13 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
15 /* Computes an evalue representation of a coordinate
16 * in a ParamVertices.
18 static evalue *vertex2evalue(Value *vertex, int nparam)
20 return affine2evalue(vertex, vertex[nparam+1], nparam);
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)
28 int i, j;
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");
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.
43 static evalue *determinant_cols(evalue ***matrix, int dim, int *cols)
45 evalue *det, *tmp;
46 evalue mone;
47 int j;
48 int *newcols;
50 if (dim == 1) {
51 det = ALLOC(evalue);
52 value_init(det->d);
53 evalue_copy(det, matrix[0][cols[0]]);
54 return det;
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);
77 free(newcols);
78 free_evalue_refs(&mone);
80 return det;
83 static evalue *determinant(evalue ***matrix, int dim)
85 int i;
86 int *cols = ALLOCN(int, dim);
87 evalue *det;
89 for (i = 0; i < dim; ++i)
90 cols[i] = i;
92 det = determinant_cols(matrix, dim, cols);
94 free(cols);
96 return det;
99 /* Compute the facet of P that saturates constraint c.
101 static Polyhedron *facet(Polyhedron *P, int c, unsigned MaxRays)
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;
111 /* Substitute parameters by the corresponding element in subs
113 static evalue *evalue_substitute_new(evalue *e, evalue **subs)
115 evalue *res = NULL;
116 evalue *c;
117 int i;
119 if (value_notzero_p(e->d)) {
120 res = ALLOC(evalue);
121 value_init(res->d);
122 evalue_copy(res, e);
123 return res;
125 assert(e->x.p->type == polynomial);
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);
134 c = evalue_substitute_new(&e->x.p->arr[0], subs);
135 eadd(c, res);
136 evalue_free(c);
138 return res;
141 struct parameter_point {
142 Vector *coord;
143 evalue **e;
146 struct parameter_point *parameter_point_new(unsigned nparam)
148 struct parameter_point *point = ALLOC(struct parameter_point);
149 point->coord = Vector_Alloc(nparam+1);
150 point->e = NULL;
151 return point;
154 evalue **parameter_point_evalue(struct parameter_point *point)
156 int j;
157 unsigned nparam = point->coord->Size-1;
159 if (point->e)
160 return point->e;
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]);
169 return point->e;
172 void parameter_point_free(struct parameter_point *point)
174 int i;
175 unsigned nparam = point->coord->Size-1;
177 Vector_Free(point->coord);
179 if (point->e) {
180 for (i = 0; i < nparam; ++i)
181 evalue_free(point->e[i]);
182 free(point->e);
184 free(point);
187 /* Computes point in pameter space where polyhedron is non-empty.
189 static struct parameter_point *non_empty_point(Param_Domain *D)
191 unsigned nparam = D->Domain->Dimension;
192 struct parameter_point *point;
193 Vector *v;
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);
200 return point;
203 static Matrix *barycenter(Param_Polyhedron *PP, Param_Domain *D)
205 int nbV;
206 Matrix *center = NULL;
207 Value denom;
208 Value fc, fv;
209 unsigned nparam;
210 int i;
211 Param_Vertices *V;
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);
233 END_FORALL_PVertex_in_ParamPolyhedron;
234 value_clear(fc);
235 value_clear(fv);
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);
243 value_clear(denom);
245 return center;
248 static Matrix *triangulation_vertex(Param_Polyhedron *PP, Param_Domain *D,
249 Polyhedron *F)
251 Param_Vertices *V;
253 FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
254 return V->Vertex;
255 END_FORALL_PVertex_in_ParamPolyhedron;
257 assert(0);
258 return NULL;
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.
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.
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);
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)
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;
292 POL_UNSET(cut_MaxRays, POL_INTEGER);
294 value_init(mone.d);
295 evalue_set_si(&mone, -1, 1);
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);
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]);
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;
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);
338 Polyhedron_Free(FF);
339 Param_Domain_Free(FD);
342 if (options->approx->volume_triangulate != BV_VOL_BARYCENTER)
343 Vector_Free(v);
345 for (j = 0; j < dim; ++j)
346 evalue_free(matrix[row][j]);
348 free_evalue_refs(&mone);
349 return vol;
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)
357 evalue mone;
358 Param_Vertices *V;
359 evalue *vol, *val;
360 int i, j;
362 options->stats->volume_simplices++;
364 value_init(mone.d);
365 evalue_set_si(&mone, -1, 1);
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]);
377 ++i;
378 END_FORALL_PVertex_in_ParamPolyhedron;
380 vol = determinant(matrix+1, dim);
382 val = evalue_substitute_new(vol, parameter_point_evalue(point));
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);
389 evalue_free(val);
391 for (i = row; i < dim+1; ++i)
392 for (j = 0; j < dim; ++j)
393 evalue_free(matrix[i][j]);
395 free_evalue_refs(&mone);
397 return vol;
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)
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;
416 SD.Domain = 0;
417 SD.next = 0;
418 nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
419 SD.F = ALLOCN(unsigned, nv);
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]);
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;
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;
452 assert(t <= MAX_TRY);
453 vol = NULL;
455 for (i = 0; i < nbV; ++i)
456 value_set_si(FixedRays->p[i][1+dim], random_int((t+1)*dim*nbV)+1);
458 M = Matrix_Copy(FixedRays);
459 L = Rays2Polyhedron(M, options->MaxRays);
460 Matrix_Free(M);
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;
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;
480 ++nbV;
481 END_FORALL_PVertex_in_ParamPolyhedron;
482 assert(r == (dim-row)+1);
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);
492 Polyhedron_Free(L);
493 Matrix_Free(FixedRays);
494 value_clear(tmp);
495 free(SD.F);
497 return vol;
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)
506 int nbV;
507 Param_Vertices *V;
508 evalue *vol;
510 assert(point);
512 nbV = 0;
513 FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
514 ++nbV;
515 END_FORALL_PVertex_in_ParamPolyhedron;
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);
529 return vol;
532 evalue* Param_Polyhedron_Volume(Polyhedron *P, Polyhedron* C,
533 struct barvinok_options *options)
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;
548 if (options->approx->approximation == BV_APPROX_SIGN_NONE)
549 return NULL;
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);
555 P = Polyhedron_Flate(P, nparam, pa == BV_APPROX_SIGN_UPPER,
556 options->MaxRays);
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;
563 Polyhedron_Free(P);
564 return vol;
567 TC = true_context(P, C, options->MaxRays);
569 MaxRays = options->MaxRays;
570 POL_UNSET(options->MaxRays, POL_INTEGER);
572 value_init(fact);
573 Factorial(nvar, &fact);
575 PP = Polyhedron2Param_Polyhedron(P, C, options);
577 for (nd = 0, D = PP->D; D; ++nd, D = D->next);
578 s = ALLOCN(struct evalue_section, nd);
580 matrix = ALLOCN(evalue **, nvar+1);
581 for (i = 0; i < nvar+1; ++i)
582 matrix[i] = ALLOCN(evalue *, nvar);
584 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
585 Polyhedron *CA, *F;
586 struct parameter_point *point;
588 CA = align_context(D->Domain, P->Dimension, MaxRays);
589 F = DomainIntersection(P, CA, options->MaxRays);
590 Domain_Free(CA);
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);
602 vol = evalue_from_section_array(s, nd);
603 free(s);
605 for (i = 0; i < nvar+1; ++i)
606 free(matrix[i]);
607 free(matrix);
608 Param_Polyhedron_Free(PP);
609 value_clear(fact);
611 return vol;