?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
add isl_pw_qpolynomial_upper_bound
[barvinok.git] / test_approx.c
blob9106df5f0ca2001cb8175be685e61a7d13061c0b
1 #include <assert.h>
2 #include <limits.h>
3 #include <barvinok/barvinok.h>
4 #include "verify.h"
5 #include "argp.h"
6 #include "progname.h"
7
8 #ifdef HAVE_SYS_TIMES_H
9
10 #include <sys/times.h>
11
12 typedef clock_t my_clock_t;
13
14 static my_clock_t time_diff(struct tms *before, struct tms *after)
15 {
16 return after->tms_utime - before->tms_utime;
17 }
18
19 #else
20
21 typedef int my_clock_t;
22
23 struct tms { int dummy; };
24 static void times(struct tms* time)
25 {
26 }
27 static my_clock_t time_diff(struct tms *before, struct tms *after)
28 {
29 return 0;
30 }
31
32 #endif
33
34 struct {
35 int sign;
36 int method;
37 int flags;
38 } methods[] = {
39 { BV_APPROX_SIGN_NONE, BV_APPROX_NONE, 0 },
40 { BV_APPROX_SIGN_APPROX, BV_APPROX_BERNOULLI, 0 },
41 { BV_APPROX_SIGN_APPROX, BV_APPROX_DROP, 0 },
42 { BV_APPROX_SIGN_APPROX, BV_APPROX_VOLUME, BV_VOL_LIFT },
43 { BV_APPROX_SIGN_APPROX, BV_APPROX_VOLUME, BV_VOL_VERTEX },
44 { BV_APPROX_SIGN_APPROX, BV_APPROX_VOLUME, BV_VOL_BARYCENTER },
45 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, 0 },
46 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, BV_APPROX_SCALE_CHAMBER },
47 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST },
48 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
49 { BV_APPROX_SIGN_LOWER, BV_APPROX_DROP, 0 },
50 { BV_APPROX_SIGN_LOWER, BV_APPROX_VOLUME, BV_VOL_LIFT },
51 { BV_APPROX_SIGN_LOWER, BV_APPROX_VOLUME, BV_VOL_VERTEX },
52 { BV_APPROX_SIGN_LOWER, BV_APPROX_VOLUME, BV_VOL_BARYCENTER },
53 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, 0 },
54 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_CHAMBER },
55 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST },
56 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
57 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW },
58 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER},
59 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
60 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
61 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 },
62 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
63 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
64 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
65 { BV_APPROX_SIGN_UPPER, BV_APPROX_DROP, 0 },
66 { BV_APPROX_SIGN_UPPER, BV_APPROX_VOLUME, BV_VOL_LIFT },
67 { BV_APPROX_SIGN_UPPER, BV_APPROX_VOLUME, BV_VOL_VERTEX },
68 { BV_APPROX_SIGN_UPPER, BV_APPROX_VOLUME, BV_VOL_BARYCENTER },
69 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, 0 },
70 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_CHAMBER },
71 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST },
72 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
73 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW },
74 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
75 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
76 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
77 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 },
78 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
79 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
80 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
81 };
82
83 #define nr_methods (sizeof(methods)/sizeof(*methods))
84
85 struct argp_option argp_options[] = {
86 { "quiet", 'q' },
87 { 0 }
88 };
89
90 struct options {
91 int quiet;
92 struct verify_options verify;
93 };
94
95 static error_t parse_opt(int key, char *arg, struct argp_state *state)
96 {
97 struct options *options = (struct options*) state->input;
98
99 switch (key) {
100 case ARGP_KEY_INIT:
101 state->child_inputs[0] = options->verify.barvinok;
102 state->child_inputs[1] = &options->verify;
103 options->quiet = 0;
104 break;
105 case 'q':
106 options->quiet = 1;
107 break;
108 default:
109 return ARGP_ERR_UNKNOWN;
110 }
111 return 0;
112 }
113
114 struct result_data {
115 Value n;
116 double RE_sum[nr_methods];
117
118 my_clock_t ticks[nr_methods];
119 size_t size[nr_methods];
120 };
121
122 void result_data_init(struct result_data *result)
123 {
124 int i;
125 for (i = 0; i < nr_methods; ++i) {
126 result->RE_sum[i] = 0;
127 result->ticks[i] = 0;
128 result->size[i] = 0;
129 }
130 value_init(result->n);
131 }
132
133 void result_data_clear(struct result_data *result)
134 {
135 int i;
136 value_clear(result->n);
137 }
138
139 void result_data_print(struct result_data *result, int n)
140 {
141 int i;
142
143 fprintf(stderr, "%d", result->ticks[0]/n);
144 for (i = 1; i < nr_methods; ++i)
145 fprintf(stderr, ", %d", result->ticks[i]/n);
146 fprintf(stderr, "\n");
147
148 fprintf(stderr, "%d", result->size[0]/n);
149 for (i = 1; i < nr_methods; ++i)
150 fprintf(stderr, ", %d", result->size[i]/n);
151 fprintf(stderr, "\n");
152
153 fprintf(stderr, "%g\n", VALUE_TO_DOUBLE(result->n));
154 fprintf(stderr, "%g", result->RE_sum[0]/VALUE_TO_DOUBLE(result->n));
155 for (i = 1; i < nr_methods; ++i)
156 fprintf(stderr, ", %g", result->RE_sum[i]/VALUE_TO_DOUBLE(result->n));
157 fprintf(stderr, "\n");
158 }
159
160 struct test_approx_data {
161 struct check_poly_data cp;
162 evalue **EP;
163 struct result_data *result;
164 };
165
166 static void eval(const evalue *EP, Value *z, int sign, Value *v)
167 {
168 evalue *res;
169
170 res = evalue_eval(EP, z);
171 if (sign == BV_APPROX_SIGN_LOWER)
172 mpz_cdiv_q(*v, res->x.n, res->d);
173 else if (sign == BV_APPROX_SIGN_UPPER)
174 mpz_fdiv_q(*v, res->x.n, res->d);
175 else if (sign == BV_APPROX_SIGN_APPROX)
176 mpz_tdiv_q(*v, res->x.n, res->d);
177 else {
178 assert(value_one_p(res->d));
179 value_assign(*v, res->x.n);
180 }
181 evalue_free(res);
182 }
183
184 static int test_approx(const struct check_poly_data *data, int nparam, Value *z,
185 const struct verify_options *options)
186 {
187 struct test_approx_data* ta_data = (struct test_approx_data*) data;
188 Value exact, approx;
189 int i;
190
191 value_init(exact);
192 value_init(approx);
193
194 eval(ta_data->EP[0], z, BV_APPROX_SIGN_NONE, &exact);
195
196 /*
197 value_print(stderr, VALUE_FMT, exact);
198 */
199
200 value_increment(ta_data->result->n, ta_data->result->n);
201 for (i = 1; i < nr_methods; ++i) {
202 double error;
203 eval(ta_data->EP[i], z, methods[i].sign, &approx);
204 /*
205 fprintf(stderr, ", ");
206 value_print(stderr, VALUE_FMT, approx);
207 */
208 if (methods[i].sign == BV_APPROX_SIGN_LOWER)
209 assert(value_le(approx, exact));
210 if (methods[i].sign == BV_APPROX_SIGN_UPPER)
211 assert(value_ge(approx, exact));
212 value_subtract(approx, approx, exact);
213 if (value_zero_p(exact))
214 error = abs(VALUE_TO_DOUBLE(approx));
215 else
216 error = abs(VALUE_TO_DOUBLE(approx)) / VALUE_TO_DOUBLE(exact);
217 ta_data->result->RE_sum[i] += error;
218 }
219
220 /*
221 fprintf(stderr, "\n");
222 */
223
224 value_clear(exact);
225 value_clear(approx);
226 return 1;
227 }
228
229 static void test(Polyhedron *P, Polyhedron *C, evalue **EP,
230 struct result_data *result,
231 struct verify_options *options)
232 {
233 Polyhedron *CS;
234 Vector *p;
235 unsigned nparam = C->Dimension;
236 struct test_approx_data data;
237
238 CS = check_poly_context_scan(P, &C, C->Dimension, options);
239
240 p = Vector_Alloc(P->Dimension+2);
241 value_set_si(p->p[P->Dimension+1], 1);
242
243 check_poly_init(C, options);
244
245 data.cp.z = p->p;
246 data.cp.check = test_approx;
247 data.EP = EP;
248 data.result = result;
249 check_poly(CS, &data.cp, nparam, 0, p->p+P->Dimension-nparam+1,
250 options);
251 if (!options->print_all)
252 printf("\n");
253
254 Vector_Free(p);
255 if (CS) {
256 Domain_Free(CS);
257 Domain_Free(C);
258 }
259 }
260
261 void Matrix_File_Read_Input(FILE *in, Matrix *Mat)
262 {
263 Value *p;
264 int i,j,n;
265 char *c, s[1024],str[1024];
266
267 p = Mat->p_Init;
268 for (i=0;i<Mat->NbRows;i++) {
269 do {
270 c = fgets(s, 1024, in);
271 while(isspace(*c) && *c!='\n')
272 ++c;
273 } while(c && (*c =='#' || *c== '\n'));
274
275 if (!c) {
276 errormsg1( "Matrix_Read", "baddim", "not enough rows" );
277 break;
278 }
279 for (j=0;j<Mat->NbColumns;j++) {
280 if(!c || *c=='\n' || *c=='#') {
281 errormsg1("Matrix_Read", "baddim", "not enough columns");
282 break;
283 }
284 if (sscanf(c,"%s%n",str,&n) == 0) {
285 errormsg1( "Matrix_Read", "baddim", "not enough columns" );
286 break;
287 }
288 value_read(*(p++),str);
289 c += n;
290 }
291 }
292 } /* Matrix_Read_Input */
293
294 /*
295 * Read the contents of the matrix 'Mat' from standard input.
296 * A '#' in the first column is a comment line
297 */
298 Matrix *Matrix_File_Read(FILE *in)
299 {
300 Matrix *Mat;
301 unsigned NbRows, NbColumns;
302 char s[1024];
303
304 fgets(s, 1024, in);
305 while ((*s=='#' || *s=='\n') ||
306 (sscanf(s, "%d %d", &NbRows, &NbColumns)<2))
307 fgets(s, 1024, in);
308 Mat = Matrix_Alloc(NbRows,NbColumns);
309 if(!Mat) {
310 errormsg1("Matrix_Read", "outofmem", "out of memory space");
311 return(NULL);
312 }
313 Matrix_File_Read_Input(in, Mat);
314 return Mat;
315 } /* Matrix_Read */
316
317 void handle(FILE *in, struct result_data *result, struct verify_options *options)
318 {
319 int i;
320 Polyhedron *A, *C;
321 Matrix *M;
322 const char **param_name;
323 evalue *EP[nr_methods];
324
325 M = Matrix_File_Read(in);
326 A = Constraints2Polyhedron(M, options->barvinok->MaxRays);
327 Matrix_Free(M);
328 M = Matrix_File_Read(in);
329 C = Constraints2Polyhedron(M, options->barvinok->MaxRays);
330 Matrix_Free(M);
331 param_name = Read_ParamNames(in, C->Dimension);
332
333 for (i = 0; i < nr_methods; ++i) {
334 struct tms st_cpu;
335 struct tms en_cpu;
336 options->barvinok->polynomial_approximation = methods[i].sign;
337 options->barvinok->approximation_method = methods[i].method;
338 if (methods[i].method == BV_APPROX_SCALE)
339 options->barvinok->scale_flags = methods[i].flags;
340 else if (methods[i].method == BV_APPROX_VOLUME)
341 options->barvinok->volume_triangulate = methods[i].flags;
342
343 times(&st_cpu);
344 EP[i] = barvinok_enumerate_with_options(A, C, options->barvinok);
345 times(&en_cpu);
346 result->ticks[i] = time_diff(&en_cpu, &st_cpu);
347 /*
348 print_evalue(stdout, EP[i], param_name);
349 */
350 }
351 for (i = 0; i < nr_methods; ++i)
352 result->size[i] = evalue_size(EP[i])/4;
353 test(A, C, EP, result, options);
354 for (i = 0; i < nr_methods; ++i)
355 evalue_free(EP[i]);
356
357 Free_ParamNames(param_name, C->Dimension);
358 Polyhedron_Free(A);
359 Polyhedron_Free(C);
360 }
361
362 int main(int argc, char **argv)
363 {
364 struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
365 char path[PATH_MAX+1];
366 struct result_data all_result;
367 int n = 0;
368 static struct argp_child argp_children[] = {
369 { &barvinok_argp, 0, 0, 0 },
370 { &verify_argp, 0, "verification", BV_GRP_LAST+1 },
371 { 0 }
372 };
373 static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
374 struct options options;
375
376 options.verify.barvinok = bv_options;
377 set_program_name(argv[0]);
378 argp_parse(&argp, argc, argv, 0, 0, &options);
379
380 if (options.verify.M == INT_MIN)
381 options.verify.M = 100;
382 if (options.verify.m == INT_MAX)
383 options.verify.m = -100;
384
385 result_data_init(&all_result);
386
387 while (fgets(path, sizeof(path), stdin)) {
388 struct result_data result;
389 FILE *in;
390 int i;
391
392 ++n;
393 result_data_init(&result);
394 fprintf(stderr, "%s", path);
395 *strchr(path, '\n') = '\0';
396 in = fopen(path, "r");
397 assert(in);
398 handle(in, &result, &options.verify);
399 fclose(in);
400
401 if (!options.quiet)
402 result_data_print(&result, 1);
403
404 value_addto(all_result.n, all_result.n, result.n);
405 for (i = 0; i < nr_methods; ++i) {
406 all_result.RE_sum[i] += result.RE_sum[i];
407 all_result.ticks[i] += result.ticks[i];
408 all_result.size[i] += result.size[i];
409 }
410
411 result_data_clear(&result);
412
413 if (!options.quiet) {
414 fprintf(stderr, "average\n");
415 result_data_print(&all_result, n);
416 }
417 }
418
419 result_data_clear(&all_result);
420
421 barvinok_options_free(bv_options);
422
423 return 0;
424 }
lattice point enumeration library