test_bound.c

   1 #include <assert.h>
   2 #include <limits.h>
   3 #include <math.h>
   4 #include <isl/stream.h>
   5 #include <isl/options.h>
   6 #include <barvinok/options.h>
   7 #include <barvinok/util.h>
   8 #include "verify.h"
   9 #include "config.h"
  10
  11 #ifdef HAVE_SYS_TIMES_H
  12
  13 #include <sys/times.h>
  14
  15 typedef clock_t         my_clock_t;
  16
  17 static my_clock_t time_diff(struct tms *before, struct tms *after)
  18 {
  19         return after->tms_utime - before->tms_utime;
  20 }
  21
  22 #else
  23
  24 typedef int             my_clock_t;
  25
  26 struct tms {};
  27 static void times(struct tms* time)
  28 {
  29 }
  30 static my_clock_t time_diff(struct tms *before, struct tms *after)
  31 {
  32         return 0;
  33 }
  34
  35 #endif
  36
  37 static struct {
  38     int     method;
  39 } methods[] = {
  40     { ISL_BOUND_BERNSTEIN },
  41     { ISL_BOUND_RANGE },
  42 };
  43
  44 #define nr_methods (sizeof(methods)/sizeof(*methods))
  45
  46 struct options {
  47     struct verify_options    *verify;
  48     int quiet;
  49 };
  50
  51 ISL_ARGS_START(struct options, options_args)
  52 ISL_ARG_CHILD(struct options, verify, NULL, &verify_options_args, NULL)
  53 ISL_ARG_BOOL(struct options, quiet, 'q', "quiet", 0, NULL)
  54 ISL_ARGS_END
  55
  56 ISL_ARG_DEF(options, struct options, options_args)
  57
  58 struct result_data {
  59     isl_val                 *n;
  60     double                  RE_sum[nr_methods];
  61
  62     my_clock_t              ticks[nr_methods];
  63     size_t                  size[nr_methods];
  64 };
  65
  66 void result_data_init(isl_ctx *ctx, struct result_data *result)
  67 {
  68     int i;
  69     for (i = 0; i < nr_methods; ++i) {
  70         result->RE_sum[i] = 0;
  71         result->ticks[i] = 0;
  72         result->size[i] = 0;
  73     }
  74     result->n = isl_val_zero(ctx);
  75 }
  76
  77 void result_data_clear(struct result_data *result)
  78 {
  79         isl_val_free(result->n);
  80 }
  81
  82 void result_data_print(struct result_data *result, int n)
  83 {
  84     int i;
  85
  86     fprintf(stderr, "%g", (double)result->ticks[0]/n);
  87     for (i = 1; i < nr_methods; ++i)
  88         fprintf(stderr, ", %g", (double)result->ticks[i]/n);
  89     fprintf(stderr, "\n");
  90
  91     fprintf(stderr, "%zd/%d", result->size[0], n);
  92     for (i = 1; i < nr_methods; ++i)
  93         fprintf(stderr, ", %zd/%d", result->size[i], n);
  94     fprintf(stderr, "\n");
  95
  96     fprintf(stderr, "%g\n", isl_val_get_d(result->n));
  97     fprintf(stderr, "%g", result->RE_sum[0]/isl_val_get_d(result->n));
  98     for (i = 1; i < nr_methods; ++i)
  99         fprintf(stderr, ", %g", result->RE_sum[i]/isl_val_get_d(result->n));
 100     fprintf(stderr, "\n");
 101 }
 102
 103 struct verify_point_bound {
 104         struct verify_point_data vpd;
 105         struct result_data *result;
 106         isl_pw_qpolynomial *pwqp;
 107         isl_pw_qpolynomial_fold **pwf;
 108 };
 109
 110 static int verify_point(__isl_take isl_point *pnt, void *user)
 111 {
 112         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
 113         const struct verify_options *options = vpb->vpd.options;
 114         int i;
 115         unsigned nparam;
 116         isl_val *max, *min, *exact, *approx, *t;
 117         isl_pw_qpolynomial *pwqp;
 118         isl_printer *p;
 119
 120         vpb->vpd.n--;
 121
 122         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
 123
 124         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
 125         for (i = 0; i < nparam; ++i) {
 126                 t = isl_point_get_coordinate_val(pnt, isl_dim_param, i);
 127                 pwqp = isl_pw_qpolynomial_fix_val(pwqp, isl_dim_param, i, t);
 128         }
 129
 130         max = isl_pw_qpolynomial_max(isl_pw_qpolynomial_copy(pwqp));
 131         max = isl_val_floor(max);
 132
 133         min = isl_pw_qpolynomial_min(pwqp);
 134         min = isl_val_ceil(min);
 135
 136         exact = isl_val_sub(isl_val_copy(max), isl_val_copy(min));
 137         exact = isl_val_add_ui(exact, 1);
 138
 139         if (options->print_all) {
 140                 p = isl_printer_to_file(isl_point_get_ctx(pnt), stderr);
 141                 p = isl_printer_print_str(p, "max: ");
 142                 p = isl_printer_print_val(p, max);
 143                 p = isl_printer_print_str(p, ", min: ");
 144                 p = isl_printer_print_val(p, min);
 145                 p = isl_printer_print_str(p, ", range: ");
 146                 p = isl_printer_print_val(p, exact);
 147         }
 148         isl_val_free(max);
 149         isl_val_free(min);
 150
 151         for (i = 0; i < nr_methods; ++i) {
 152                 double error;
 153
 154                 max = isl_pw_qpolynomial_fold_eval(
 155                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i]),
 156                                 isl_point_copy(pnt));
 157                 max = isl_val_floor(max);
 158
 159                 min = isl_pw_qpolynomial_fold_eval(
 160                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i + 1]),
 161                                 isl_point_copy(pnt));
 162                 min = isl_val_ceil(min);
 163
 164                 approx = isl_val_sub(max, min);
 165                 approx = isl_val_add_ui(approx, 1);
 166                 if (options->print_all) {
 167                         p = isl_printer_print_str(p, ", ");
 168                         p = isl_printer_print_val(p, approx);
 169                 }
 170
 171                 assert(isl_val_ge(approx, exact));
 172                 approx = isl_val_sub(approx, isl_val_copy(exact));
 173
 174                 error = fabs(isl_val_get_d(approx)) / isl_val_get_d(exact);
 175                 if (options->print_all)
 176                         fprintf(stderr, " (%g)", error);
 177                 vpb->result->RE_sum[i] += error;
 178
 179                 isl_val_free(approx);
 180         }
 181
 182         if (options->print_all) {
 183                 p = isl_printer_end_line(p);
 184                 isl_printer_free(p);
 185         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 186                 printf("o");
 187                 fflush(stdout);
 188         }
 189
 190         isl_val_free(exact);
 191         isl_point_free(pnt);
 192
 193         return vpb->vpd.n >= 1 ? 0 : -1;
 194 }
 195
 196 static void test(__isl_keep isl_pw_qpolynomial *pwqp,
 197         __isl_keep isl_pw_qpolynomial_fold **pwf, struct result_data *result,
 198         struct verify_options *options)
 199 {
 200         struct verify_point_bound vpb = { { options }, result };
 201         isl_set *context;
 202         int r;
 203
 204         vpb.pwf = pwf;
 205         vpb.pwqp = pwqp;
 206         context = isl_pw_qpolynomial_domain(isl_pw_qpolynomial_copy(vpb.pwqp));
 207         context = isl_set_params(context);
 208         context = verify_context_set_bounds(context, options);
 209
 210         r = verify_point_data_init(&vpb.vpd, context);
 211         assert(r == 0);
 212         result->n = isl_val_set_si(result->n, vpb.vpd.n);
 213
 214         isl_set_foreach_point(context, verify_point, &vpb);
 215         assert(!vpb.vpd.error);
 216
 217         isl_set_free(context);
 218
 219         verify_point_data_fini(&vpb.vpd);
 220 }
 221
 222 static void handle(isl_ctx *ctx, FILE *in, struct result_data *result,
 223         struct verify_options *options)
 224 {
 225     int i;
 226     isl_pw_qpolynomial *pwqp, *upper, *lower;
 227     isl_pw_qpolynomial_fold *pwf[2*nr_methods];
 228     struct isl_stream *s;
 229
 230     s = isl_stream_new_file(ctx, in);
 231     pwqp = isl_stream_read_pw_qpolynomial(s);
 232
 233     upper = isl_pw_qpolynomial_to_polynomial(isl_pw_qpolynomial_copy(pwqp), 1);
 234     lower = isl_pw_qpolynomial_to_polynomial(isl_pw_qpolynomial_copy(pwqp), -1);
 235
 236     for (i = 0; i < nr_methods; ++i) {
 237         int j;
 238         struct tms st_cpu;
 239         struct tms en_cpu;
 240
 241         times(&st_cpu);
 242         for (j = 0; j < 2; ++j) {
 243             isl_pw_qpolynomial *poly = j == 0 ? upper : lower;
 244             enum isl_fold type = j == 0 ? isl_fold_max : isl_fold_min;
 245             isl_options_set_bound(ctx, methods[i].method);
 246             poly = isl_pw_qpolynomial_copy(poly);
 247             pwf[2*i+j] = isl_pw_qpolynomial_bound(poly, type, NULL);
 248         }
 249         times(&en_cpu);
 250         result->ticks[i] = time_diff(&en_cpu, &st_cpu);
 251         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i]);
 252         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i+1]);
 253         if (options->barvinok->verbose) {
 254             isl_printer *p = isl_printer_to_file(ctx, stdout);
 255             for (j = 0; j < 2; ++j) {
 256                 p = isl_printer_print_pw_qpolynomial_fold(p, pwf[2*i+j]);
 257                 p = isl_printer_end_line(p);
 258             }
 259             isl_printer_free(p);
 260         }
 261     }
 262     test(pwqp, pwf, result, options);
 263
 264     for (i = 0; i < 2*nr_methods; ++i)
 265         isl_pw_qpolynomial_fold_free(pwf[i]);
 266     isl_pw_qpolynomial_free(pwqp);
 267     isl_pw_qpolynomial_free(lower);
 268     isl_pw_qpolynomial_free(upper);
 269
 270     isl_stream_free(s);
 271 }
 272
 273 int main(int argc, char **argv)
 274 {
 275     char path[PATH_MAX+1];
 276     struct result_data all_result;
 277     int n = 0;
 278     struct options *options = options_new_with_defaults();
 279     isl_ctx *ctx = isl_ctx_alloc();
 280
 281     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 282
 283     if (options->verify->M == INT_MIN)
 284         options->verify->M = 100;
 285     if (options->verify->m == INT_MAX)
 286         options->verify->m = -100;
 287
 288     result_data_init(ctx, &all_result);
 289
 290     while (fgets(path, sizeof(path), stdin)) {
 291         struct result_data result;
 292         FILE *in;
 293         int i;
 294
 295         ++n;
 296         result_data_init(ctx, &result);
 297         fprintf(stderr, "%s", path);
 298         *strchr(path, '\n') = '\0';
 299         in = fopen(path, "r");
 300         assert(in);
 301         handle(ctx, in, &result, options->verify);
 302         fclose(in);
 303
 304         if (!options->quiet)
 305             result_data_print(&result, 1);
 306
 307         all_result.n = isl_val_add(all_result.n, isl_val_copy(result.n));
 308         for (i = 0; i < nr_methods; ++i) {
 309             all_result.RE_sum[i] += result.RE_sum[i];
 310             all_result.ticks[i] += result.ticks[i];
 311             all_result.size[i] += result.size[i];
 312         }
 313
 314         result_data_clear(&result);
 315
 316         if (!options->quiet) {
 317             fprintf(stderr, "average\n");
 318             result_data_print(&all_result, n);
 319         }
 320     }
 321
 322     result_data_clear(&all_result);
 323
 324     options_free(options);
 325     isl_ctx_free(ctx);
 326
 327     return 0;
 328 }