?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
volume.c: volume_triangulate: use vertex instead of barycenter by default
[barvinok.git] / barvinok_enumerate.cc
blobc2f9c818379844f36651c6a15f4b75b42194554f
1 #include <unistd.h>
2 #include <stdlib.h>
3 #include <barvinok/evalue.h>
4 #include <barvinok/util.h>
5 #include <barvinok/barvinok.h>
6 #include "fdstream.h"
7 #include "argp.h"
8 #include "progname.h"
9 #include "verify.h"
10 #include "verif_ehrhart.h"
11 #include "remove_equalities.h"
12 #include "evalue_convert.h"
13
14 #undef CS /* for Solaris 10 */
15
16 /* The input of this example program is the same as that of testehrhart
17 * in the PolyLib distribution, i.e., a polytope in combined
18 * data and parameter space, a context polytope in parameter space
19 * and (optionally) the names of the parameters.
20 * Both polytopes are in PolyLib notation.
21 */
22
23 #define PRINT_STATS (BV_OPT_LAST+1)
24
25 struct argp_option argp_options[] = {
26 { "size", 'S' },
27 { "series", 's', 0, 0, "compute rational generating function" },
28 { "explicit", 'e', 0, 0, "convert rgf to psp" },
29 { "verbose", 'v' },
30 { "print-stats", PRINT_STATS, 0, 0 },
31 { 0 }
32 };
33
34 struct arguments {
35 int size;
36 int series;
37 int function;
38 int verbose;
39 int print_stats;
40 struct verify_options verify;
41 struct convert_options convert;
42 };
43
44 static error_t parse_opt(int key, char *arg, struct argp_state *state)
45 {
46 struct arguments *options = (struct arguments*) state->input;
47
48 switch (key) {
49 case ARGP_KEY_INIT:
50 state->child_inputs[0] = options->verify.barvinok;
51 state->child_inputs[1] = &options->verify;
52 state->child_inputs[2] = &options->convert;
53 options->size = 0;
54 options->series = 0;
55 options->function = 0;
56 options->verbose = 0;
57 options->print_stats = 0;
58 break;
59 case PRINT_STATS:
60 options->print_stats = 1;
61 break;
62 case 'S':
63 options->size = 1;
64 break;
65 case 'e':
66 options->function = 1;
67 /* fall through */
68 case 's':
69 options->series = 1;
70 break;
71 case 'v':
72 options->verbose = 1;
73 break;
74 default:
75 return ARGP_ERR_UNKNOWN;
76 }
77 return 0;
78 }
79
80 struct skewed_gen_fun {
81 gen_fun *gf;
82 /* maps original space to space in which gf is defined */
83 Matrix *T;
84 /* equalities in the original space that need to be satisfied for
85 * gf to be valid
86 */
87 Matrix *eq;
88 /* divisibilities in the original space that need to be satisfied for
89 * gf to be valid
90 */
91 Matrix *div;
92
93 skewed_gen_fun(gen_fun *gf, Matrix *T, Matrix *eq, Matrix *div) :
94 gf(gf), T(T), eq(eq), div(div) {}
95 ~skewed_gen_fun() {
96 if (T)
97 Matrix_Free(T);
98 if (eq)
99 Matrix_Free(eq);
100 if (div)
101 Matrix_Free(div);
102 delete gf;
103 }
104
105 void print(FILE *out, unsigned int nparam, char **param_name) const;
106 operator evalue *() const {
107 assert(T == NULL && eq == NULL); /* other cases not supported for now */
108 return *gf;
109 }
110 void coefficient(Value* params, Value* c, barvinok_options *options) const;
111 };
112
113 void skewed_gen_fun::print(FILE *out, unsigned int nparam,
114 char **param_name) const
115 {
116 fdostream os(dup(fileno(out)));
117 if (T) {
118 fprintf(out, "T:\n");
119 Matrix_Print(out, P_VALUE_FMT, T);
120 }
121 if (eq) {
122 fprintf(out, "eq:\n");
123 Matrix_Print(out, P_VALUE_FMT, eq);
124 }
125 if (div) {
126 fprintf(out, "div:\n");
127 Matrix_Print(out, P_VALUE_FMT, div);
128 }
129 gf->print(os, nparam, param_name);
130 }
131
132 void skewed_gen_fun::coefficient(Value* params, Value* c,
133 barvinok_options *options) const
134 {
135 if (eq) {
136 for (int i = 0; i < eq->NbRows; ++i) {
137 Inner_Product(eq->p[i]+1, params, eq->NbColumns-2, eq->p[i]);
138 if (value_notzero_p(eq->p[i][0])) {
139 value_set_si(*c, 0);
140 return;
141 }
142 }
143 }
144 if (div) {
145 Value tmp;
146 value_init(tmp);
147 for (int i = 0; i < div->NbRows; ++i) {
148 Inner_Product(div->p[i], params, div->NbColumns-1, &tmp);
149 if (!mpz_divisible_p(tmp, div->p[i][div->NbColumns-1])) {
150 value_set_si(*c, 0);
151 return;
152 }
153 }
154 value_clear(tmp);
155 }
156
157 ZZ coeff;
158 if (!T)
159 coeff = gf->coefficient(params, options);
160 else {
161 Vector *p2 = Vector_Alloc(T->NbRows);
162 Matrix_Vector_Product(T, params, p2->p);
163 if (value_notone_p(p2->p[T->NbRows-1]))
164 Vector_AntiScale(p2->p, p2->p, p2->p[T->NbRows-1], T->NbRows);
165 coeff = gf->coefficient(p2->p, options);
166 Vector_Free(p2);
167 }
168
169 zz2value(coeff, *c);
170 }
171
172 static int check_series(Polyhedron *S, Polyhedron *CS, skewed_gen_fun *gf,
173 int nparam, int pos, Value *z, verify_options *options)
174 {
175 int k;
176 Value c, tmp;
177 Value LB, UB;
178
179 value_init(c);
180 value_init(tmp);
181 value_init(LB);
182 value_init(UB);
183
184 if (pos == nparam) {
185 /* Computes the coefficient */
186 gf->coefficient(&z[S->Dimension-nparam+1], &c, options->barvinok);
187
188 /* if c=0 we may be out of context. */
189 /* scanning is useless in this case*/
190
191 if (options->print_all) {
192 printf("EP( ");
193 value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
194 for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
195 printf(", ");
196 value_print(stdout,VALUE_FMT,z[k]);
197 }
198 printf(" ) = ");
199 value_print(stdout,VALUE_FMT,c);
200 printf(" ");
201 }
202
203 /* Manually count the number of points */
204 count_points(1,S,z,&tmp);
205 if (options->print_all) {
206 printf(", count = ");
207 value_print(stdout, P_VALUE_FMT, tmp);
208 printf(". ");
209 }
210
211 if (value_ne(tmp,c)) {
212 printf("\n");
213 fflush(stdout);
214 fprintf(stderr,"Error !\n");
215 fprintf(stderr,"EP( ");
216 value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
217 for (k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
218 fprintf(stderr,", ");
219 value_print(stderr,VALUE_FMT,z[k]);
220 }
221 fprintf(stderr," ) should be ");
222 value_print(stderr,VALUE_FMT,tmp);
223 fprintf(stderr,", while EP eval gives ");
224 value_print(stderr,VALUE_FMT,c);
225 fprintf(stderr,".\n");
226 if (!options->continue_on_error) {
227 value_clear(c); value_clear(tmp);
228 return 0;
229 }
230 } else if (options->print_all)
231 printf("OK.\n");
232 } else {
233 int ok =
234 !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
235 assert(ok);
236 for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
237 if (!options->print_all) {
238 k = VALUE_TO_INT(tmp);
239 if(!pos && !(k % options->st)) {
240 printf("o");
241 fflush(stdout);
242 }
243 }
244 value_assign(z[pos+S->Dimension-nparam+1],tmp);
245 if (!check_series(S, CS->next, gf, nparam, pos+1, z, options)) {
246 value_clear(c); value_clear(tmp);
247 value_clear(LB);
248 value_clear(UB);
249 return(0);
250 }
251 }
252 value_set_si(z[pos+S->Dimension-nparam+1],0);
253 }
254
255 value_clear(c);
256 value_clear(tmp);
257 value_clear(LB);
258 value_clear(UB);
259 return 1;
260 }
261
262 static int verify(Polyhedron *P, Polyhedron *C, evalue *EP, skewed_gen_fun *gf,
263 arguments *options)
264 {
265 Polyhedron *CS, *S;
266 Vector *p;
267 int result = 0;
268
269 CS = check_poly_context_scan(P, &C, C->Dimension, &options->verify);
270
271 p = Vector_Alloc(P->Dimension+2);
272 value_set_si(p->p[P->Dimension+1], 1);
273
274 /* S = scanning list of polyhedra */
275 S = Polyhedron_Scan(P, C, options->verify.barvinok->MaxRays);
276
277 check_poly_init(C, &options->verify);
278
279 /******* CHECK NOW *********/
280 if (S) {
281 if (!options->series || options->function) {
282 if (!check_poly_EP(S, CS, EP, 0, C->Dimension, 0, p->p,
283 &options->verify))
284 result = -1;
285 } else {
286 if (!check_series(S, CS, gf, C->Dimension, 0, p->p, &options->verify))
287 result = -1;
288 }
289 Domain_Free(S);
290 }
291
292 if (result == -1)
293 fprintf(stderr,"Check failed !\n");
294
295 if (!options->verify.print_all)
296 printf( "\n" );
297
298 Vector_Free(p);
299 if (CS) {
300 Domain_Free(CS);
301 Domain_Free(C);
302 }
303
304 return result;
305 }
306
307 static void unimodular_complete(Matrix *M, int row)
308 {
309 Matrix *H, *Q, *U;
310 left_hermite(M, &H, &Q, &U);
311 Matrix_Free(H);
312 Matrix_Free(U);
313 for (int r = row; r < M->NbRows; ++r)
314 Vector_Copy(Q->p[r], M->p[r], M->NbColumns);
315 Matrix_Free(Q);
316 }
317
318 /* frees M and Minv */
319 static void apply_transformation(Polyhedron **P, Polyhedron **C,
320 bool free_P, bool free_C,
321 Matrix *M, Matrix *Minv, Matrix **inv,
322 barvinok_options *options)
323 {
324 Polyhedron *T;
325 Matrix *M2;
326
327 M2 = align_matrix(M, (*P)->Dimension + 1);
328 T = *P;
329 *P = Polyhedron_Preimage(*P, M2, options->MaxRays);
330 if (free_P)
331 Polyhedron_Free(T);
332 Matrix_Free(M2);
333
334 T = *C;
335 *C = Polyhedron_Preimage(*C, M, options->MaxRays);
336 if (free_C)
337 Polyhedron_Free(T);
338
339 Matrix_Free(M);
340
341 if (*inv) {
342 Matrix *T = *inv;
343 *inv = Matrix_Alloc(Minv->NbRows, T->NbColumns);
344 Matrix_Product(Minv, T, *inv);
345 Matrix_Free(T);
346 Matrix_Free(Minv);
347 } else
348 *inv = Minv;
349 }
350
351 static skewed_gen_fun *series(Polyhedron *P, Polyhedron* C,
352 barvinok_options *options)
353 {
354 Polyhedron *C1, *C2;
355 gen_fun *gf;
356 Matrix *inv = NULL;
357 Matrix *eq = NULL;
358 Matrix *div = NULL;
359 Polyhedron *PT = P;
360
361 /* Compute true context */
362 C1 = Polyhedron_Project(P, C->Dimension);
363 C2 = DomainIntersection(C, C1, options->MaxRays);
364 Polyhedron_Free(C1);
365
366 POL_ENSURE_VERTICES(C2);
367 if (C2->NbBid != 0) {
368 Polyhedron *T;
369 Matrix *M, *Minv, *M2;
370 Matrix *CP;
371 if (C2->NbEq || P->NbEq) {
372 /* We remove all equalities to be sure all lines are unit vectors */
373 Polyhedron *CT = C2;
374 remove_all_equalities(&PT, &CT, &CP, NULL, C2->Dimension,
375 options->MaxRays);
376 if (CT != C2) {
377 Polyhedron_Free(C2);
378 C2 = CT;
379 }
380 if (CP) {
381 inv = left_inverse(CP, &eq);
382 Matrix_Free(CP);
383
384 int d = 0;
385 Value tmp;
386 value_init(tmp);
387 div = Matrix_Alloc(inv->NbRows-1, inv->NbColumns+1);
388 for (int i = 0; i < inv->NbRows-1; ++i) {
389 Vector_Gcd(inv->p[i], inv->NbColumns, &tmp);
390 if (mpz_divisible_p(tmp,
391 inv->p[inv->NbRows-1][inv->NbColumns-1]))
392 continue;
393 Vector_Copy(inv->p[i], div->p[d], inv->NbColumns);
394 value_assign(div->p[d][inv->NbColumns],
395 inv->p[inv->NbRows-1][inv->NbColumns-1]);
396 ++d;
397 }
398 value_clear(tmp);
399
400 if (!d) {
401 Matrix_Free(div);
402 div = NULL;
403 } else
404 div->NbRows = d;
405 }
406 }
407 POL_ENSURE_VERTICES(C2);
408
409 /* Since we have "compressed" the parameters (in case there were
410 * any equalities), the result is independent of the coordinates in the
411 * coordinate subspace spanned by the lines. We can therefore assume
412 * these coordinates are zero and compute the inverse image of the map
413 * from a lower dimensional space that adds zeros in the appropriate
414 * places.
415 */
416 M = Matrix_Alloc(C2->Dimension+1, C2->Dimension-C2->NbBid+1);
417 int k = 0;
418 for (int i = 0; i < C2->NbBid; ++i) {
419 int j = First_Non_Zero(C2->Ray[i]+1, C2->Dimension);
420 assert(First_Non_Zero(C2->Ray[i]+1+j+1, C2->Dimension-j-1) == -1);
421 for ( ; k < j; k++)
422 value_set_si(M->p[k+i][k], 1);
423 }
424 for ( ; k < C2->Dimension-C2->NbBid+1; k++)
425 value_set_si(M->p[k+C2->NbBid][k], 1);
426 Minv = Transpose(M);
427
428 apply_transformation(&PT, &C2, PT != P, C2 != C, M, Minv, &inv, options);
429 }
430 POL_ENSURE_VERTICES(C2);
431 if (!Polyhedron_has_revlex_positive_rays(C2, C2->Dimension)) {
432 Polyhedron *T;
433 Matrix *Constraints;
434 Matrix *H, *Q, *U;
435 Constraints = Matrix_Alloc(C2->NbConstraints, C2->Dimension+1);
436 for (int i = 0; i < C2->NbConstraints; ++i)
437 Vector_Copy(C2->Constraint[i]+1, Constraints->p[i], C2->Dimension);
438 left_hermite(Constraints, &H, &Q, &U);
439 /* flip rows of Q */
440 for (int i = 0; i < C2->Dimension/2; ++i)
441 Vector_Exchange(Q->p[i], Q->p[C2->Dimension-1-i], C2->Dimension);
442 Matrix_Free(H);
443 Matrix_Free(U);
444 Matrix *M = Matrix_Alloc(C2->Dimension+1, C2->Dimension+1);
445 U = Matrix_Copy(Q);
446 int ok = Matrix_Inverse(U, M);
447 assert(ok);
448 Matrix_Free(U);
449
450 apply_transformation(&PT, &C2, PT != P, C2 != C, M, Q, &inv, options);
451 }
452 gf = barvinok_series_with_options(PT, C2, options);
453 Polyhedron_Free(C2);
454 if (PT != P)
455 Polyhedron_Free(PT);
456 return new skewed_gen_fun(gf, inv, eq, div);
457 }
458
459 int main(int argc, char **argv)
460 {
461 Polyhedron *A, *C;
462 Matrix *M;
463 evalue *EP = NULL;
464 skewed_gen_fun *gf = NULL;
465 char **param_name;
466 int print_solution = 1;
467 int result = 0;
468 struct arguments options;
469 static struct argp_child argp_children[] = {
470 { &barvinok_argp, 0, 0, 0 },
471 { &verify_argp, 0, "verification", BV_GRP_LAST+1 },
472 { &convert_argp, 0, "output conversion", BV_GRP_LAST+2 },
473 { 0 }
474 };
475 static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
476 struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
477
478 options.verify.barvinok = bv_options;
479 set_program_name(argv[0]);
480 argp_parse(&argp, argc, argv, 0, 0, &options);
481
482 M = Matrix_Read();
483 A = Constraints2Polyhedron(M, bv_options->MaxRays);
484 Matrix_Free(M);
485 M = Matrix_Read();
486 C = Constraints2Polyhedron(M, bv_options->MaxRays);
487 Matrix_Free(M);
488 param_name = Read_ParamNames(stdin, C->Dimension);
489
490 if (options.verify.verify) {
491 verify_options_set_range(&options.verify, A->Dimension);
492 if (!options.verbose)
493 print_solution = 0;
494 }
495
496 if (print_solution && options.verbose) {
497 Polyhedron_Print(stdout, P_VALUE_FMT, A);
498 Polyhedron_Print(stdout, P_VALUE_FMT, C);
499 }
500
501 if (options.series) {
502 gf = series(A, C, bv_options);
503 if (print_solution) {
504 gf->print(stdout, C->Dimension, param_name);
505 puts("");
506 }
507 if (options.function) {
508 EP = *gf;
509 if (print_solution)
510 print_evalue(stdout, EP, param_name);
511 }
512 } else {
513 EP = barvinok_enumerate_with_options(A, C, bv_options);
514 if (evalue_convert(EP, &options.convert, options.verbose, C->Dimension,
515 param_name))
516 print_solution = 0;
517 if (options.size)
518 printf("\nSize: %d\n", evalue_size(EP));
519 }
520 if (print_solution)
521 print_evalue(stdout, EP, param_name);
522
523 if (options.verify.verify) {
524 options.verify.params = param_name;
525 result = verify(A, C, EP, gf, &options);
526 }
527
528 if (gf)
529 delete gf;
530 if (EP) {
531 free_evalue_refs(EP);
532 free(EP);
533 }
534
535 if (options.print_stats)
536 barvinok_stats_print(options.verify.barvinok->stats, stdout);
537
538 Free_ParamNames(param_name, C->Dimension);
539 Polyhedron_Free(A);
540 Polyhedron_Free(C);
541 barvinok_options_free(bv_options);
542 return result;
543 }
lattice point enumeration library