volume.c: face_vertices: use Param_Inner_Product
[barvinok.git] / test_approx.c
blob7b58ceb27f5f051c74e3be4f46dc9c4bc7e164cb
1 #include <limits.h>
2 #include <sys/times.h>
3 #include <barvinok/barvinok.h>
4 #include "verify.h"
5 #include "argp.h"
6 #include "progname.h"
8 struct {
9 int sign;
10 int method;
11 int flags;
12 } methods[] = {
13 { BV_APPROX_SIGN_NONE, BV_APPROX_NONE, 0 },
14 { BV_APPROX_SIGN_APPROX, BV_APPROX_DROP, 0 },
15 { BV_APPROX_SIGN_APPROX, BV_APPROX_VOLUME, 0 },
16 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, 0 },
17 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, BV_APPROX_SCALE_CHAMBER },
18 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST },
19 { BV_APPROX_SIGN_APPROX, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
20 { BV_APPROX_SIGN_LOWER, BV_APPROX_DROP, 0 },
21 { BV_APPROX_SIGN_LOWER, BV_APPROX_VOLUME, 0 },
22 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, 0 },
23 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_CHAMBER },
24 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST },
25 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
26 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW },
27 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER},
28 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
29 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
30 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 },
31 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
32 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
33 { BV_APPROX_SIGN_LOWER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
34 { BV_APPROX_SIGN_UPPER, BV_APPROX_DROP, 0 },
35 { BV_APPROX_SIGN_UPPER, BV_APPROX_VOLUME, 0 },
36 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, 0 },
37 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_CHAMBER },
38 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST },
39 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_CHAMBER },
40 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW },
41 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
42 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW },
43 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW | BV_APPROX_SCALE_CHAMBER },
44 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 },
45 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
46 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 },
47 { BV_APPROX_SIGN_UPPER, BV_APPROX_SCALE, BV_APPROX_SCALE_FAST | BV_APPROX_SCALE_NARROW2 | BV_APPROX_SCALE_CHAMBER },
50 #define nr_methods (sizeof(methods)/sizeof(*methods))
52 struct argp_option argp_options[] = {
53 { "quiet", 'q' },
54 { 0 }
57 struct options {
58 int quiet;
59 struct verify_options verify;
62 static error_t parse_opt(int key, char *arg, struct argp_state *state)
64 struct options *options = (struct options*) state->input;
66 switch (key) {
67 case ARGP_KEY_INIT:
68 state->child_inputs[0] = options->verify.barvinok;
69 state->child_inputs[1] = &options->verify;
70 options->quiet = 0;
71 break;
72 case 'q':
73 options->quiet = 1;
74 break;
75 default:
76 return ARGP_ERR_UNKNOWN;
78 return 0;
81 struct result_data {
82 Value n;
83 double RE_sum[nr_methods];
85 clock_t ticks[nr_methods];
86 size_t size[nr_methods];
89 void result_data_init(struct result_data *result)
91 int i;
92 for (i = 0; i < nr_methods; ++i) {
93 result->RE_sum[i] = 0;
94 result->ticks[i] = 0;
95 result->size[i] = 0;
97 value_init(result->n);
100 void result_data_clear(struct result_data *result)
102 int i;
103 value_clear(result->n);
106 void result_data_print(struct result_data *result, int n)
108 int i;
110 fprintf(stderr, "%d", result->ticks[0]/n);
111 for (i = 1; i < nr_methods; ++i)
112 fprintf(stderr, ", %d", result->ticks[i]/n);
113 fprintf(stderr, "\n");
115 fprintf(stderr, "%d", result->size[0]/n);
116 for (i = 1; i < nr_methods; ++i)
117 fprintf(stderr, ", %d", result->size[i]/n);
118 fprintf(stderr, "\n");
120 fprintf(stderr, "%g\n", VALUE_TO_DOUBLE(result->n));
121 fprintf(stderr, "%g", result->RE_sum[0]/VALUE_TO_DOUBLE(result->n));
122 for (i = 1; i < nr_methods; ++i)
123 fprintf(stderr, ", %g", result->RE_sum[i]/VALUE_TO_DOUBLE(result->n));
124 fprintf(stderr, "\n");
127 struct test_approx_data {
128 struct check_poly_data cp;
129 evalue **EP;
130 struct result_data *result;
133 static void eval(const evalue *EP, Value *z, int sign, Value *v)
135 evalue *res;
137 res = evalue_eval(EP, z);
138 if (sign == BV_APPROX_SIGN_LOWER)
139 mpz_cdiv_q(*v, res->x.n, res->d);
140 else if (sign == BV_APPROX_SIGN_UPPER)
141 mpz_fdiv_q(*v, res->x.n, res->d);
142 else if (sign == BV_APPROX_SIGN_APPROX)
143 mpz_tdiv_q(*v, res->x.n, res->d);
144 else {
145 assert(value_one_p(res->d));
146 value_assign(*v, res->x.n);
148 free_evalue_refs(res);
149 free(res);
152 static int test_approx(const struct check_poly_data *data, int nparam, Value *z,
153 const struct verify_options *options)
155 struct test_approx_data* ta_data = (struct test_approx_data*) data;
156 Value exact, approx;
157 int i;
159 value_init(exact);
160 value_init(approx);
162 eval(ta_data->EP[0], z, BV_APPROX_SIGN_NONE, &exact);
165 value_print(stderr, VALUE_FMT, exact);
168 value_increment(ta_data->result->n, ta_data->result->n);
169 for (i = 1; i < nr_methods; ++i) {
170 double error;
171 eval(ta_data->EP[i], z, methods[i].sign, &approx);
173 fprintf(stderr, ", ");
174 value_print(stderr, VALUE_FMT, approx);
176 if (methods[i].sign == BV_APPROX_SIGN_LOWER)
177 assert(value_le(approx, exact));
178 if (methods[i].sign == BV_APPROX_SIGN_UPPER)
179 assert(value_ge(approx, exact));
180 value_subtract(approx, approx, exact);
181 if (value_zero_p(exact))
182 error = abs(VALUE_TO_DOUBLE(approx));
183 else
184 error = abs(VALUE_TO_DOUBLE(approx)) / VALUE_TO_DOUBLE(exact);
185 ta_data->result->RE_sum[i] += error;
189 fprintf(stderr, "\n");
192 value_clear(exact);
193 value_clear(approx);
194 return 1;
197 static void test(Polyhedron *P, Polyhedron *C, evalue **EP,
198 struct result_data *result,
199 struct verify_options *options)
201 Polyhedron *CS;
202 Vector *p;
203 unsigned nparam = C->Dimension;
204 struct test_approx_data data;
206 CS = check_poly_context_scan(P, &C, C->Dimension, options);
208 p = Vector_Alloc(P->Dimension+2);
209 value_set_si(p->p[P->Dimension+1], 1);
211 check_poly_init(C, options);
213 data.cp.z = p->p;
214 data.cp.check = test_approx;
215 data.EP = EP;
216 data.result = result;
217 check_poly(CS, &data.cp, nparam, 0, p->p+P->Dimension-nparam+1,
218 options);
219 if (!options->print_all)
220 printf("\n");
222 Vector_Free(p);
223 if (CS) {
224 Domain_Free(CS);
225 Domain_Free(C);
229 void Matrix_File_Read_Input(FILE *in, Matrix *Mat)
231 Value *p;
232 int i,j,n;
233 char *c, s[1024],str[1024];
235 p = Mat->p_Init;
236 for (i=0;i<Mat->NbRows;i++) {
237 do {
238 c = fgets(s, 1024, in);
239 while(isspace(*c) && *c!='\n')
240 ++c;
241 } while(c && (*c =='#' || *c== '\n'));
243 if (!c) {
244 errormsg1( "Matrix_Read", "baddim", "not enough rows" );
245 break;
247 for (j=0;j<Mat->NbColumns;j++) {
248 if(!c || *c=='\n' || *c=='#') {
249 errormsg1("Matrix_Read", "baddim", "not enough columns");
250 break;
252 if (sscanf(c,"%s%n",str,&n) == 0) {
253 errormsg1( "Matrix_Read", "baddim", "not enough columns" );
254 break;
256 value_read(*(p++),str);
257 c += n;
260 } /* Matrix_Read_Input */
263 * Read the contents of the matrix 'Mat' from standard input.
264 * A '#' in the first column is a comment line
266 Matrix *Matrix_File_Read(FILE *in)
268 Matrix *Mat;
269 unsigned NbRows, NbColumns;
270 char s[1024];
272 fgets(s, 1024, in);
273 while ((*s=='#' || *s=='\n') ||
274 (sscanf(s, "%d %d", &NbRows, &NbColumns)<2))
275 fgets(s, 1024, in);
276 Mat = Matrix_Alloc(NbRows,NbColumns);
277 if(!Mat) {
278 errormsg1("Matrix_Read", "outofmem", "out of memory space");
279 return(NULL);
281 Matrix_File_Read_Input(in, Mat);
282 return Mat;
283 } /* Matrix_Read */
285 void handle(FILE *in, struct result_data *result, struct verify_options *options)
287 int i;
288 Polyhedron *A, *C;
289 Matrix *M;
290 char **param_name;
291 evalue *EP[nr_methods];
293 M = Matrix_File_Read(in);
294 A = Constraints2Polyhedron(M, options->barvinok->MaxRays);
295 Matrix_Free(M);
296 M = Matrix_File_Read(in);
297 C = Constraints2Polyhedron(M, options->barvinok->MaxRays);
298 Matrix_Free(M);
299 param_name = Read_ParamNames(in, C->Dimension);
301 for (i = 0; i < nr_methods; ++i) {
302 struct tms st_cpu;
303 struct tms en_cpu;
304 options->barvinok->polynomial_approximation = methods[i].sign;
305 options->barvinok->approximation_method = methods[i].method;
306 options->barvinok->scale_flags = methods[i].flags;
308 times(&st_cpu);
309 EP[i] = barvinok_enumerate_with_options(A, C, options->barvinok);
310 times(&en_cpu);
311 result->ticks[i] = en_cpu.tms_utime - st_cpu.tms_utime;
313 print_evalue(stdout, EP[i], param_name);
316 for (i = 0; i < nr_methods; ++i)
317 result->size[i] = evalue_size(EP[i])/4;
318 test(A, C, EP, result, options);
319 for (i = 0; i < nr_methods; ++i) {
320 free_evalue_refs(EP[i]);
321 free(EP[i]);
324 Free_ParamNames(param_name, C->Dimension);
325 Polyhedron_Free(A);
326 Polyhedron_Free(C);
329 int main(int argc, char **argv)
331 struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
332 char path[PATH_MAX+1];
333 struct result_data all_result;
334 int n = 0;
335 static struct argp_child argp_children[] = {
336 { &barvinok_argp, 0, 0, 0 },
337 { &verify_argp, 0, "verification", BV_GRP_LAST+1 },
338 { 0 }
340 static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
341 struct options options;
343 options.verify.barvinok = bv_options;
344 set_program_name(argv[0]);
345 argp_parse(&argp, argc, argv, 0, 0, &options);
347 if (options.verify.M == INT_MIN)
348 options.verify.M = 100;
349 if (options.verify.m == INT_MAX)
350 options.verify.m = -100;
352 result_data_init(&all_result);
354 while (fgets(path, sizeof(path), stdin)) {
355 struct result_data result;
356 FILE *in;
357 int i;
359 ++n;
360 result_data_init(&result);
361 fprintf(stderr, "%s", path);
362 *strchr(path, '\n') = '\0';
363 in = fopen(path, "r");
364 assert(in);
365 handle(in, &result, &options.verify);
366 fclose(in);
368 if (!options.quiet)
369 result_data_print(&result, 1);
371 value_addto(all_result.n, all_result.n, result.n);
372 for (i = 0; i < nr_methods; ++i) {
373 all_result.RE_sum[i] += result.RE_sum[i];
374 all_result.ticks[i] += result.ticks[i];
375 all_result.size[i] += result.size[i];
378 result_data_clear(&result);
380 if (!options.quiet) {
381 fprintf(stderr, "average\n");
382 result_data_print(&all_result, n);
386 result_data_clear(&all_result);
388 barvinok_options_free(bv_options);
390 return 0;