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