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