add gen_fun::add(short_rat *r)
[barvinok.git] / barvinok_enumerate_e.cc
blob26df3f24be6da23d0f754839829d2cec8660462b
1 #include <unistd.h>
2 #include <stdlib.h>
3 #include <assert.h>
4 #include <barvinok/util.h>
5 #include <barvinok/barvinok.h>
6 #include "argp.h"
7 #include "error.h"
8 #include "config.h"
9 #ifdef HAVE_OMEGA
10 #include "omega/convert.h"
11 #endif
12 #include "verify.h"
14 /* The input of this example program is a polytope in combined
15 * data and parameter space followed by two lines indicating
16 * the number of existential variables and parameters respectively.
17 * The first lines starts with "E ", followed by a number.
18 * The second lines starts with "P ", followed by a number.
19 * These two lines are (optionally) followed by the names of the parameters.
20 * The polytope is in PolyLib notation.
23 struct argp_option argp_options[] = {
24 { "omega", 'o', 0, 0 },
25 { "pip", 'p', 0, 0 },
26 { "series", 's', 0, 0 },
27 { "scarf", 'S', 0, 0 },
28 { "convert", 'c', 0, 0 },
29 { "floor", 'f', 0, 0 },
30 { "range-reduction", 'R', 0, 0 },
31 { 0 }
34 struct arguments {
35 struct barvinok_options *options;
36 struct verify_options verify;
37 int range;
38 int convert;
39 int omega;
40 int pip;
41 int scarf;
42 int series;
43 int floor;
46 error_t parse_opt(int key, char *arg, struct argp_state *state)
48 struct arguments *arguments = (struct arguments *)(state->input);
50 switch (key) {
51 case ARGP_KEY_INIT:
52 state->child_inputs[0] = arguments->options;
53 state->child_inputs[1] = &arguments->verify;
54 break;
55 case 's':
56 arguments->series = 1;
57 break;
58 case 'S':
59 arguments->scarf = 1;
60 break;
61 case 'o':
62 #ifdef HAVE_OMEGA
63 arguments->omega = 1;
64 #else
65 error(0, 0, "--omega option not supported");
66 #endif
67 break;
68 case 'p':
69 #ifdef HAVE_PIPLIB
70 arguments->pip = 1;
71 #else
72 error(0, 0, "--pip option not supported");
73 #endif
74 break;
75 case 'f':
76 arguments->floor = 1;
77 break;
78 case 'c':
79 arguments->convert = 1;
80 break;
81 case 'R':
82 arguments->range = 1;
83 break;
84 default:
85 return ARGP_ERR_UNKNOWN;
87 return 0;
90 #ifdef HAVE_OMEGA
92 Polyhedron *Omega_simplify(Polyhedron *P,
93 unsigned exist, unsigned nparam, char **parms)
95 varvector varv;
96 varvector paramv;
97 Relation r = Polyhedron2relation(P, exist, nparam, parms);
98 Polyhedron_Free(P);
99 return relation2Domain(r, varv, paramv);
101 #else
102 Polyhedron *Omega_simplify(Polyhedron *P,
103 unsigned exist, unsigned nparam, char **parms)
105 return P;
107 #endif
109 /* define this to continue the test after first error found */
110 /* #define DONT_BREAK_ON_ERROR */
112 static Value min_val, max_val;
114 static char **params;
116 static int st;
118 static int check_poly(Polyhedron *S, Polyhedron *C, evalue *EP,
119 int exist, int nparam, int pos, Value *z, int print_all);
120 static void verify_results(Polyhedron *P, evalue *EP, int exist, int nparam,
121 int m, int M, int print_all, unsigned MaxRays);
123 int main(int argc, char **argv)
125 Polyhedron *A;
126 Matrix *MA;
127 char **param_name;
128 int exist, nparam, nvar;
129 char s[128];
130 evalue *EP;
131 int print_solution = 1;
132 struct arguments arguments;
133 static struct argp_child argp_children[] = {
134 { &barvinok_argp, 0, 0, 0 },
135 { &verify_argp, 0, "verification", 1 },
136 { 0 }
138 static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
139 struct barvinok_options *options = barvinok_options_new_with_defaults();
141 arguments.options = options;
142 arguments.range = 0;
143 arguments.convert = 0;
144 arguments.omega = 0;
145 arguments.pip = 0;
146 arguments.scarf = 0;
147 arguments.series = 0;
148 arguments.floor = 0;
150 argp_parse(&argp, argc, argv, 0, 0, &arguments);
152 if (arguments.series && !arguments.scarf) {
153 fprintf(stderr,
154 "--series currently only available if --scarf is specified\n");
155 exit(1);
158 MA = Matrix_Read();
159 A = Constraints2Polyhedron(MA, options->MaxRays);
160 Matrix_Free(MA);
162 fgets(s, 128, stdin);
163 while ((*s=='#') || (sscanf(s, "E %d", &exist)<1))
164 fgets(s, 128, stdin);
166 fgets(s, 128, stdin);
167 while ((*s=='#') || (sscanf(s, "P %d", &nparam)<1))
168 fgets(s, 128, stdin);
170 /******* Read the options: initialize Min and Max ********/
171 verify_options_set_range(&arguments.verify, A);
173 if (arguments.verify.verify)
174 print_solution = 0;
176 if (print_solution) {
177 Polyhedron_Print(stdout, P_VALUE_FMT, A);
178 printf("exist: %d, nparam: %d\n", exist, nparam);
180 param_name = Read_ParamNames(stdin, nparam);
181 nvar = A->Dimension - exist - nparam;
182 if (arguments.omega) {
183 A = Omega_simplify(A, exist, nparam, param_name);
184 assert(!A->next);
185 exist = A->Dimension - nvar - nparam;
187 if (arguments.series) {
188 gen_fun *gf;
189 barvinok_options *options = barvinok_options_new_with_defaults();
190 assert(arguments.scarf);
191 gf = barvinok_enumerate_scarf_series(A, exist, nparam, options);
192 if (print_solution) {
193 gf->print(std::cout, nparam, param_name);
194 puts("");
196 delete gf;
197 free(options);
198 } else {
199 if (arguments.scarf) {
200 barvinok_options *options = barvinok_options_new_with_defaults();
201 EP = barvinok_enumerate_scarf(A, exist, nparam, options);
202 free(options);
203 } else if (arguments.pip && exist > 0)
204 EP = barvinok_enumerate_pip_with_options(A, exist, nparam, options);
205 else
206 EP = barvinok_enumerate_e_with_options(A, exist, nparam, options);
207 reduce_evalue(EP);
208 evalue_combine(EP);
209 if (arguments.range)
210 evalue_range_reduction(EP);
211 if (print_solution)
212 print_evalue(stdout, EP, param_name);
213 if (arguments.floor) {
214 fprintf(stderr, "WARNING: floor conversion not supported\n");
215 evalue_frac2floor2(EP, 0);
216 if (print_solution)
217 print_evalue(stdout, EP, param_name);
218 } else if (arguments.convert) {
219 evalue_mod2table(EP, nparam);
220 if (print_solution)
221 print_evalue(stdout, EP, param_name);
223 if (arguments.verify.verify)
224 verify_results(A, EP, exist, nparam, arguments.verify.m,
225 arguments.verify.M, arguments.verify.print_all,
226 options->MaxRays);
227 free_evalue_refs(EP);
228 free(EP);
230 Free_ParamNames(param_name, nparam);
231 Polyhedron_Free(A);
232 return 0;
235 void verify_results(Polyhedron *P, evalue *EP, int exist, int nparam, int m, int M,
236 int print_all, unsigned MaxRays)
238 int i;
239 int res;
240 Value *p, tmp;
241 Polyhedron *S;
242 Polyhedron *C = Polyhedron_Project(P, nparam);
243 value_init(min_val);
244 value_init(max_val);
245 value_set_si(min_val,m);
246 value_set_si(max_val,M);
247 value_init(tmp);
249 p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
250 for(i=0;i<=P->Dimension;i++) {
251 value_init(p[i]);
252 value_set_si(p[i],0);
254 value_init(p[i]);
255 value_set_si(p[i],1);
257 /* S = scanning list of polyhedra */
258 S = Polyhedron_Scan(P, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
260 if (!print_all) {
261 if (C->Dimension > 0) {
262 value_subtract(tmp,max_val,min_val);
263 if (VALUE_TO_INT(tmp) > 80)
264 st = 1+(VALUE_TO_INT(tmp))/80;
265 else
266 st=1;
267 for(i=VALUE_TO_INT(min_val);i<=VALUE_TO_INT(max_val);i+=st)
268 printf(".");
269 printf( "\r" );
270 fflush(stdout);
274 /******* CHECK NOW *********/
275 res = 0;
276 if(S && !check_poly(S, C, EP, exist, nparam, 0, p, print_all)) {
277 fprintf(stderr,"Check failed !\n");
278 res = -1;
281 if (!print_all)
282 printf( "\n" );
284 for(i=0;i<=(P->Dimension+1);i++)
285 value_clear(p[i]);
286 free(p);
287 value_clear(tmp);
288 Domain_Free(S);
289 Polyhedron_Free(C);
292 /****************************************************/
293 /* function check_poly : */
294 /* scans the parameter space from min to max (all */
295 /* directions). Computes the number of points in */
296 /* the polytope using both methods, and compare them*/
297 /* returns 1 on success */
298 /****************************************************/
300 int check_poly(Polyhedron *S, Polyhedron *C, evalue *EP,
301 int exist, int nparam, int pos, Value *z, int print_all)
303 int k;
304 Value c,tmp;
306 value_init(c); value_init(tmp);
308 if(pos == nparam) {
310 /* Computes the ehrhart polynomial */
311 value_set_double(c, compute_evalue(EP,&z[S->Dimension-nparam+1])+.25);
312 /* if c=0 we may be out of context. */
313 /* scanning is useless in this case*/
314 if(!in_domain(C,&z[S->Dimension-nparam+1])) {
316 /* ok */ ;
318 else {
320 if (print_all) {
321 printf("EP( ");
322 value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
323 for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
324 printf(", ");
325 value_print(stdout,VALUE_FMT,z[k]);
327 printf(" ) = ");
328 value_print(stdout,VALUE_FMT,c);
329 printf(" ");
332 /* Manually count the number of points */
333 count_points_e(1, S, exist, nparam, z, &tmp);
334 if (print_all) {
335 printf(", count = ");
336 value_print(stdout, P_VALUE_FMT, tmp);
337 printf(". ");
340 if(value_ne(tmp,c)) {
341 printf("\n");
342 fflush(stdout);
343 fprintf(stderr,"Error !\n");
344 fprintf(stderr,"EP( ");
345 value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
346 for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
347 fprintf(stderr,", ");
348 value_print(stderr,VALUE_FMT,z[k]);
350 fprintf(stderr," ) should be ");
351 value_print(stderr,VALUE_FMT,tmp);
352 fprintf(stderr,", while EP eval gives ");
353 value_print(stderr,VALUE_FMT,c);
354 fprintf(stderr,".\n");
355 print_evalue(stderr, EP, params);
356 #ifndef DONT_BREAK_ON_ERROR
357 value_clear(c); value_clear(tmp);
358 return(0);
359 #endif
361 else if (print_all)
362 printf("OK.\n");
365 else
366 for(value_assign(tmp,min_val); value_le(tmp,max_val); value_increment(tmp,tmp)) {
367 if (!print_all) {
368 k = VALUE_TO_INT(tmp);
369 if(!pos && !(k%st)) {
370 printf("o");
371 fflush(stdout);
375 value_assign(z[pos+S->Dimension-nparam+1],tmp);
376 if(!check_poly(S, C, EP, exist, nparam, pos+1, z, print_all)) {
377 value_clear(c); value_clear(tmp);
378 return(0);
381 value_clear(c); value_clear(tmp);
382 return(1);
383 } /* check_poly */