?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
doc: clean up "exponential substitution" section
[barvinok/uuh.git] / volume.c
blob8cf760f593b3781fba1ba24e0e6286efa66927fd
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 * const *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->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->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->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->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->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->polynomial_approximation == BV_APPROX_SIGN_NONE)
546 return NULL;
547
548 if (options->polynomial_approximation != BV_APPROX_SIGN_APPROX) {
549 int pa = options->polynomial_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->polynomial_approximation = BV_APPROX_SIGN_APPROX;
557 vol = barvinok_enumerate_with_options(P, C, options);
558 options->polynomial_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 }
patches by uuh to consider upstream