?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
barvinok_sum_over_polytope: special case 0D polytopes
[barvinok.git] / test_approx.c
bloba995bb5b430abf25f16eb6565c5e6936bc627515
1 #include <assert.h>
2 #include <ctype.h>
3 #include <limits.h>
4 #include <barvinok/barvinok.h>
5 #include "verify.h"
6 #include "config.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 options {
86 int quiet;
87 struct verify_options *verify;
88 };
89
90 ISL_ARGS_START(struct options, options_args)
91 ISL_ARG_CHILD(struct options, verify, NULL, &verify_options_args, NULL)
92 ISL_ARG_BOOL(struct options, quiet, 'q', "quiet", 0, NULL)
93 ISL_ARGS_END
94
95 ISL_ARG_DEF(options, struct options, options_args)
96
97 struct result_data {
98 Value n;
99 double RE_sum[nr_methods];
100
101 my_clock_t ticks[nr_methods];
102 size_t size[nr_methods];
103 };
104
105 void result_data_init(struct result_data *result)
106 {
107 int i;
108 for (i = 0; i < nr_methods; ++i) {
109 result->RE_sum[i] = 0;
110 result->ticks[i] = 0;
111 result->size[i] = 0;
112 }
113 value_init(result->n);
114 }
115
116 void result_data_clear(struct result_data *result)
117 {
118 int i;
119 value_clear(result->n);
120 }
121
122 void result_data_print(struct result_data *result, int n)
123 {
124 int i;
125
126 fprintf(stderr, "%g", (double)result->ticks[0]/n);
127 for (i = 1; i < nr_methods; ++i)
128 fprintf(stderr, ", %g", (double)result->ticks[i]/n);
129 fprintf(stderr, "\n");
130
131 fprintf(stderr, "%zd", result->size[0]/n);
132 for (i = 1; i < nr_methods; ++i)
133 fprintf(stderr, ", %zd", result->size[i]/n);
134 fprintf(stderr, "\n");
135
136 fprintf(stderr, "%g\n", VALUE_TO_DOUBLE(result->n));
137 fprintf(stderr, "%g", result->RE_sum[0]/VALUE_TO_DOUBLE(result->n));
138 for (i = 1; i < nr_methods; ++i)
139 fprintf(stderr, ", %g", result->RE_sum[i]/VALUE_TO_DOUBLE(result->n));
140 fprintf(stderr, "\n");
141 }
142
143 struct test_approx_data {
144 struct check_poly_data cp;
145 evalue **EP;
146 struct result_data *result;
147 };
148
149 static void eval(const evalue *EP, Value *z, int sign, Value *v)
150 {
151 evalue *res;
152
153 res = evalue_eval(EP, z);
154 if (sign == BV_APPROX_SIGN_LOWER)
155 mpz_cdiv_q(*v, res->x.n, res->d);
156 else if (sign == BV_APPROX_SIGN_UPPER)
157 mpz_fdiv_q(*v, res->x.n, res->d);
158 else if (sign == BV_APPROX_SIGN_APPROX)
159 mpz_tdiv_q(*v, res->x.n, res->d);
160 else {
161 assert(value_one_p(res->d));
162 value_assign(*v, res->x.n);
163 }
164 evalue_free(res);
165 }
166
167 static int test_approx(const struct check_poly_data *data, int nparam, Value *z,
168 const struct verify_options *options)
169 {
170 struct test_approx_data* ta_data = (struct test_approx_data*) data;
171 Value exact, approx;
172 int i;
173
174 value_init(exact);
175 value_init(approx);
176
177 eval(ta_data->EP[0], z, BV_APPROX_SIGN_NONE, &exact);
178
179 /*
180 value_print(stderr, VALUE_FMT, exact);
181 */
182
183 value_increment(ta_data->result->n, ta_data->result->n);
184 for (i = 1; i < nr_methods; ++i) {
185 double error;
186 eval(ta_data->EP[i], z, methods[i].sign, &approx);
187 /*
188 fprintf(stderr, ", ");
189 value_print(stderr, VALUE_FMT, approx);
190 */
191 if (methods[i].sign == BV_APPROX_SIGN_LOWER)
192 assert(value_le(approx, exact));
193 if (methods[i].sign == BV_APPROX_SIGN_UPPER)
194 assert(value_ge(approx, exact));
195 value_subtract(approx, approx, exact);
196 if (value_zero_p(exact))
197 error = abs(VALUE_TO_DOUBLE(approx));
198 else
199 error = abs(VALUE_TO_DOUBLE(approx)) / VALUE_TO_DOUBLE(exact);
200 ta_data->result->RE_sum[i] += error;
201 }
202
203 /*
204 fprintf(stderr, "\n");
205 */
206
207 value_clear(exact);
208 value_clear(approx);
209 return 1;
210 }
211
212 static void test(Polyhedron *P, Polyhedron *C, evalue **EP,
213 struct result_data *result,
214 struct verify_options *options)
215 {
216 Polyhedron *CS;
217 Vector *p;
218 unsigned nparam = C->Dimension;
219 struct test_approx_data data;
220
221 CS = check_poly_context_scan(P, &C, C->Dimension, options);
222
223 p = Vector_Alloc(P->Dimension+2);
224 value_set_si(p->p[P->Dimension+1], 1);
225
226 check_poly_init(C, options);
227
228 data.cp.z = p->p;
229 data.cp.check = test_approx;
230 data.EP = EP;
231 data.result = result;
232 check_poly(CS, &data.cp, nparam, 0, p->p+P->Dimension-nparam+1,
233 options);
234 if (!options->print_all)
235 printf("\n");
236
237 Vector_Free(p);
238 if (CS) {
239 Domain_Free(CS);
240 Domain_Free(C);
241 }
242 }
243
244 void Matrix_File_Read_Input(FILE *in, Matrix *Mat)
245 {
246 Value *p;
247 int i,j,n;
248 char *c, s[1024],str[1024];
249
250 p = Mat->p_Init;
251 for (i=0;i<Mat->NbRows;i++) {
252 do {
253 c = fgets(s, 1024, in);
254 while(isspace(*c) && *c!='\n')
255 ++c;
256 } while(c && (*c =='#' || *c== '\n'));
257
258 if (!c) {
259 errormsg1( "Matrix_Read", "baddim", "not enough rows" );
260 break;
261 }
262 for (j=0;j<Mat->NbColumns;j++) {
263 if(!c || *c=='\n' || *c=='#') {
264 errormsg1("Matrix_Read", "baddim", "not enough columns");
265 break;
266 }
267 if (sscanf(c,"%s%n",str,&n) == 0) {
268 errormsg1( "Matrix_Read", "baddim", "not enough columns" );
269 break;
270 }
271 value_read(*(p++),str);
272 c += n;
273 }
274 }
275 } /* Matrix_Read_Input */
276
277 /*
278 * Read the contents of the matrix 'Mat' from standard input.
279 * A '#' in the first column is a comment line
280 */
281 Matrix *Matrix_File_Read(FILE *in)
282 {
283 Matrix *Mat;
284 unsigned NbRows, NbColumns;
285 char s[1024];
286
287 fgets(s, 1024, in);
288 while ((*s=='#' || *s=='\n') ||
289 (sscanf(s, "%d %d", &NbRows, &NbColumns)<2))
290 fgets(s, 1024, in);
291 Mat = Matrix_Alloc(NbRows,NbColumns);
292 if(!Mat) {
293 errormsg1("Matrix_Read", "outofmem", "out of memory space");
294 return(NULL);
295 }
296 Matrix_File_Read_Input(in, Mat);
297 return Mat;
298 } /* Matrix_Read */
299
300 void handle(FILE *in, struct result_data *result, struct verify_options *options)
301 {
302 int i;
303 Polyhedron *A, *C;
304 Matrix *M;
305 const char **param_name;
306 evalue *EP[nr_methods];
307
308 M = Matrix_File_Read(in);
309 A = Constraints2Polyhedron(M, options->barvinok->MaxRays);
310 Matrix_Free(M);
311 M = Matrix_File_Read(in);
312 C = Constraints2Polyhedron(M, options->barvinok->MaxRays);
313 Matrix_Free(M);
314 param_name = Read_ParamNames(in, C->Dimension);
315
316 for (i = 0; i < nr_methods; ++i) {
317 struct tms st_cpu;
318 struct tms en_cpu;
319 options->barvinok->approx->approximation = methods[i].sign;
320 options->barvinok->approx->method = methods[i].method;
321 if (methods[i].method == BV_APPROX_SCALE)
322 options->barvinok->approx->scale_flags = methods[i].flags;
323 else if (methods[i].method == BV_APPROX_VOLUME)
324 options->barvinok->approx->volume_triangulate = methods[i].flags;
325
326 times(&st_cpu);
327 EP[i] = barvinok_enumerate_with_options(A, C, options->barvinok);
328 times(&en_cpu);
329 result->ticks[i] = time_diff(&en_cpu, &st_cpu);
330 /*
331 print_evalue(stdout, EP[i], param_name);
332 */
333 }
334 for (i = 0; i < nr_methods; ++i)
335 result->size[i] = evalue_size(EP[i])/4;
336 test(A, C, EP, result, options);
337 for (i = 0; i < nr_methods; ++i)
338 evalue_free(EP[i]);
339
340 Free_ParamNames(param_name, C->Dimension);
341 Polyhedron_Free(A);
342 Polyhedron_Free(C);
343 }
344
345 int main(int argc, char **argv)
346 {
347 char path[PATH_MAX+1];
348 struct result_data all_result;
349 int n = 0;
350 struct options *options = options_new_with_defaults();
351
352 argc = options_parse(options, argc, argv, ISL_ARG_ALL);
353
354 if (options->verify->M == INT_MIN)
355 options->verify->M = 100;
356 if (options->verify->m == INT_MAX)
357 options->verify->m = -100;
358
359 result_data_init(&all_result);
360
361 while (fgets(path, sizeof(path), stdin)) {
362 struct result_data result;
363 FILE *in;
364 int i;
365
366 ++n;
367 result_data_init(&result);
368 fprintf(stderr, "%s", path);
369 *strchr(path, '\n') = '\0';
370 in = fopen(path, "r");
371 assert(in);
372 handle(in, &result, options->verify);
373 fclose(in);
374
375 if (!options->quiet)
376 result_data_print(&result, 1);
377
378 value_addto(all_result.n, all_result.n, result.n);
379 for (i = 0; i < nr_methods; ++i) {
380 all_result.RE_sum[i] += result.RE_sum[i];
381 all_result.ticks[i] += result.ticks[i];
382 all_result.size[i] += result.size[i];
383 }
384
385 result_data_clear(&result);
386
387 if (!options->quiet) {
388 fprintf(stderr, "average\n");
389 result_data_print(&all_result, n);
390 }
391 }
392
393 result_data_clear(&all_result);
394
395 options_free(options);
396
397 return 0;
398 }
lattice point enumeration library