test_bound.cc

   1 #include <assert.h>
   2 #include <limits.h>
   3 #include <math.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/bernstein.h>
   6 #include <barvinok/util.h>
   7 #include <bound_common.h>
   8 #include "argp.h"
   9 #include "progname.h"
  10 #include "evalue_read.h"
  11 #include "verify.h"
  12
  13 #ifdef HAVE_SYS_TIMES_H
  14
  15 #include <sys/times.h>
  16
  17 typedef clock_t         my_clock_t;
  18
  19 static my_clock_t time_diff(struct tms *before, struct tms *after)
  20 {
  21         return after->tms_utime - before->tms_utime;
  22 }
  23
  24 #else
  25
  26 typedef int             my_clock_t;
  27
  28 struct tms {};
  29 static void times(struct tms* time)
  30 {
  31 }
  32 static my_clock_t time_diff(struct tms *before, struct tms *after)
  33 {
  34         return 0;
  35 }
  36
  37 #endif
  38
  39 using std::cout;
  40 using std::cerr;
  41 using std::endl;
  42 using namespace barvinok;
  43
  44 static struct {
  45     int     method;
  46 } methods[] = {
  47     { BV_BOUND_BERNSTEIN },
  48     { BV_BOUND_RANGE },
  49 };
  50
  51 #define nr_methods (sizeof(methods)/sizeof(*methods))
  52
  53 struct argp_option argp_options[] = {
  54     { "quiet",      'q' },
  55     { 0 }
  56 };
  57
  58 struct options {
  59     struct verify_options    verify;
  60     int quiet;
  61 };
  62
  63 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  64 {
  65     struct options *options = (struct options*) state->input;
  66
  67     switch (key) {
  68     case ARGP_KEY_INIT:
  69         state->child_inputs[0] = options->verify.barvinok;
  70         state->child_inputs[1] = &options->verify;
  71         options->quiet = 0;
  72         break;
  73     case 'q':
  74         options->quiet = 1;
  75         break;
  76     default:
  77         return ARGP_ERR_UNKNOWN;
  78     }
  79     return 0;
  80 }
  81
  82 struct result_data {
  83     isl_int                 n;
  84     double                  RE_sum[nr_methods];
  85
  86     my_clock_t              ticks[nr_methods];
  87     size_t                  size[nr_methods];
  88 };
  89
  90 void result_data_init(struct result_data *result)
  91 {
  92     int i;
  93     for (i = 0; i < nr_methods; ++i) {
  94         result->RE_sum[i] = 0;
  95         result->ticks[i] = 0;
  96         result->size[i] = 0;
  97     }
  98     isl_int_init(result->n);
  99 }
 100
 101 void result_data_clear(struct result_data *result)
 102 {
 103     int i;
 104     isl_int_clear(result->n);
 105 }
 106
 107 void result_data_print(struct result_data *result, int n)
 108 {
 109     int i;
 110
 111     fprintf(stderr, "%d", result->ticks[0]/n);
 112     for (i = 1; i < nr_methods; ++i)
 113         fprintf(stderr, ", %d", result->ticks[i]/n);
 114     fprintf(stderr, "\n");
 115
 116     fprintf(stderr, "%zd/%d", result->size[0], n);
 117     for (i = 1; i < nr_methods; ++i)
 118         fprintf(stderr, ", %zd/%d", result->size[i], n);
 119     fprintf(stderr, "\n");
 120
 121     fprintf(stderr, "%g\n", isl_int_get_d(result->n));
 122     fprintf(stderr, "%g", result->RE_sum[0]/isl_int_get_d(result->n));
 123     for (i = 1; i < nr_methods; ++i)
 124         fprintf(stderr, ", %g", result->RE_sum[i]/isl_int_get_d(result->n));
 125     fprintf(stderr, "\n");
 126 }
 127
 128 struct verify_point_bound {
 129         struct verify_point_data vpd;
 130         struct result_data *result;
 131         isl_pw_qpolynomial *pwqp;
 132         isl_pw_qpolynomial_fold **pwf;
 133 };
 134
 135 static int verify_point(__isl_take isl_point *pnt, void *user)
 136 {
 137         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
 138         const struct verify_options *options = vpb->vpd.options;
 139         int i;
 140         unsigned nparam;
 141         isl_int max, min, exact, approx;
 142         isl_int n, d;
 143         isl_qpolynomial *opt;
 144         isl_pw_qpolynomial *pwqp;
 145         int cst;
 146
 147         vpb->vpd.n--;
 148
 149         isl_int_init(max);
 150         isl_int_init(min);
 151         isl_int_init(exact);
 152         isl_int_init(approx);
 153         isl_int_init(n);
 154         isl_int_init(d);
 155
 156         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
 157
 158         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
 159         for (i = 0; i < nparam; ++i) {
 160                 isl_point_get_coordinate(pnt, isl_dim_param, i, &n);
 161                 pwqp = isl_pw_qpolynomial_fix_dim(pwqp, isl_dim_param, i, n);
 162         }
 163
 164         opt = isl_pw_qpolynomial_max(isl_pw_qpolynomial_copy(pwqp));
 165         cst = isl_qpolynomial_is_cst(opt, &n, &d);
 166         isl_qpolynomial_free(opt);
 167         assert(cst == 1);
 168         isl_int_fdiv_q(max, n, d);
 169
 170         opt = isl_pw_qpolynomial_min(pwqp);
 171         cst = isl_qpolynomial_is_cst(opt, &n, &d);
 172         isl_qpolynomial_free(opt);
 173         assert(cst == 1);
 174         isl_int_cdiv_q(min, n, d);
 175
 176         isl_int_sub(exact, max, min);
 177         isl_int_add_ui(exact, exact, 1);
 178
 179         if (options->print_all) {
 180                 fprintf(stderr, "max: "); isl_int_print(stderr, max, 0);
 181                 fprintf(stderr, ", min: "); isl_int_print(stderr, min, 0);
 182                 fprintf(stderr, ", range: "); isl_int_print(stderr, exact, 0);
 183         }
 184
 185         for (int i = 0; i < nr_methods; ++i) {
 186                 double error;
 187
 188                 opt = isl_pw_qpolynomial_fold_eval(
 189                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i]),
 190                                 isl_point_copy(pnt));
 191                 cst = isl_qpolynomial_is_cst(opt, &n, &d);
 192                 isl_qpolynomial_free(opt);
 193                 assert(cst == 1);
 194                 isl_int_fdiv_q(max, n, d);
 195
 196                 opt = isl_pw_qpolynomial_fold_eval(
 197                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i + 1]),
 198                                 isl_point_copy(pnt));
 199                 cst = isl_qpolynomial_is_cst(opt, &n, &d);
 200                 isl_qpolynomial_free(opt);
 201                 assert(cst == 1);
 202                 isl_int_cdiv_q(min, n, d);
 203
 204                 isl_int_sub(approx, max, min);
 205                 isl_int_add_ui(approx, approx, 1);
 206                 if (options->print_all) {
 207                         fprintf(stderr, ", "); isl_int_print(stderr, approx, 0);
 208                 }
 209
 210                 assert(isl_int_ge(approx, exact));
 211                 isl_int_sub(approx, approx, exact);
 212
 213                 error = fabs(isl_int_get_d(approx)) / isl_int_get_d(exact);
 214                 if (options->print_all)
 215                         fprintf(stderr, " (%g)", error);
 216                 vpb->result->RE_sum[i] += error;
 217         }
 218
 219         if (options->print_all) {
 220                 fprintf(stderr, "\n");
 221         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 222                 printf("o");
 223                 fflush(stdout);
 224         }
 225
 226         isl_int_clear(max);
 227         isl_int_clear(min);
 228         isl_int_clear(exact);
 229         isl_int_clear(approx);
 230         isl_int_clear(n);
 231         isl_int_clear(d);
 232
 233         isl_point_free(pnt);
 234
 235         return 0;
 236 }
 237
 238 static void test(evalue *EP, unsigned nvar,
 239         __isl_keep isl_pw_qpolynomial_fold **pwf, struct result_data *result,
 240         struct verify_options *options)
 241 {
 242         struct verify_point_bound vpb = { { options }, result };
 243         isl_ctx *ctx;
 244         isl_dim *dim;
 245         isl_set *context;
 246         int r;
 247         int i;
 248         unsigned nparam;
 249
 250         vpb.pwf = pwf;
 251         ctx = isl_pw_qpolynomial_fold_get_ctx(pwf[0]);
 252         nparam = isl_pw_qpolynomial_fold_dim(pwf[0], isl_dim_param);
 253         dim = isl_dim_set_alloc(ctx, nvar + nparam, 0);
 254         vpb.pwqp = isl_pw_qpolynomial_from_evalue(dim, EP);
 255         vpb.pwqp = isl_pw_qpolynomial_move(vpb.pwqp, isl_dim_set, 0,
 256                                                 isl_dim_param, 0, nvar);
 257         context = isl_pw_qpolynomial_domain(isl_pw_qpolynomial_copy(vpb.pwqp));
 258         context = isl_set_remove(context, isl_dim_set, 0, nvar);
 259         context = verify_context_set_bounds(context, options);
 260
 261         r = verify_point_data_init(&vpb.vpd, context);
 262         assert(r == 0);
 263         isl_int_set_si(result->n, vpb.vpd.n);
 264
 265         isl_set_foreach_point(context, verify_point, &vpb);
 266         assert(!vpb.vpd.error);
 267
 268         isl_set_free(context);
 269         isl_pw_qpolynomial_free(vpb.pwqp);
 270
 271         verify_point_data_fini(&vpb.vpd);
 272 }
 273
 274 void handle(FILE *in, struct result_data *result, struct verify_options *options)
 275 {
 276     evalue *EP, *upper, *lower;
 277     const char **all_vars = NULL;
 278     unsigned nvar;
 279     unsigned nparam;
 280     isl_ctx *ctx = isl_ctx_alloc();
 281     isl_dim *dim;
 282     isl_pw_qpolynomial_fold *pwf[2*nr_methods];
 283
 284     EP = evalue_read_from_file(in, NULL, &all_vars,
 285                                &nvar, &nparam, options->barvinok->MaxRays);
 286     assert(EP);
 287     if (EVALUE_IS_ZERO(*EP)) {
 288         evalue_free(EP);
 289         Free_ParamNames(all_vars, nvar+nparam);
 290         return;
 291     }
 292
 293     upper = evalue_dup(EP);
 294     lower = evalue_dup(EP);
 295     evalue_frac2polynomial(upper, 1, options->barvinok->MaxRays);
 296     evalue_frac2polynomial(lower, -1, options->barvinok->MaxRays);
 297
 298     dim = isl_dim_set_alloc(ctx, nparam, 0);
 299     for (int i = 0; i < nr_methods; ++i) {
 300         struct tms st_cpu;
 301         struct tms en_cpu;
 302
 303         times(&st_cpu);
 304         for (int j = 0; j < 2; ++j) {
 305             evalue *poly = j == 0 ? upper : lower;
 306             int sign = j == 0 ? BV_BERNSTEIN_MAX : BV_BERNSTEIN_MIN;
 307             enum isl_fold type = j == 0 ? isl_fold_max : isl_fold_min;
 308             options->barvinok->bernstein_optimize = sign;
 309             isl_dim *dim_poly;
 310             isl_pw_qpolynomial *pwqp;
 311             dim_poly = isl_dim_insert(isl_dim_copy(dim), isl_dim_param, 0, nvar);
 312             pwqp = isl_pw_qpolynomial_from_evalue(dim_poly, poly);
 313             pwqp = isl_pw_qpolynomial_move(pwqp, isl_dim_set, 0,
 314                                             isl_dim_param, 0, nvar);
 315             pwf[2*i+j] = isl_pw_qpolynomial_bound(pwqp, type, methods[i].method);
 316         }
 317         times(&en_cpu);
 318         result->ticks[i] = time_diff(&en_cpu, &st_cpu);
 319         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i]);
 320         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i+1]);
 321         if (options->barvinok->verbose) {
 322             isl_printer *p = isl_printer_to_file(ctx, stdout);
 323             for (int j = 0; j < 2; ++j) {
 324                 p = isl_printer_print_pw_qpolynomial_fold(p, pwf[2*i+j]);
 325                 p = isl_printer_end_line(p);
 326             }
 327             isl_printer_free(p);
 328         }
 329     }
 330     isl_dim_free(dim);
 331     test(EP, nvar, pwf, result, options);
 332
 333     for (int i = 0; i < 2*nr_methods; ++i)
 334         isl_pw_qpolynomial_fold_free(pwf[i]);
 335     evalue_free(EP);
 336     evalue_free(lower);
 337     evalue_free(upper);
 338     Free_ParamNames(all_vars, nvar+nparam);
 339
 340     isl_ctx_free(ctx);
 341 }
 342
 343 int main(int argc, char **argv)
 344 {
 345     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 346     char path[PATH_MAX+1];
 347     struct result_data all_result;
 348     int n = 0;
 349     static struct argp_child argp_children[] = {
 350         { &barvinok_argp,       0,      0,              0 },
 351         { &verify_argp,         0,      "verification",         BV_GRP_LAST+1 },
 352         { 0 }
 353     };
 354     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 355     struct options options;
 356
 357     options.verify.barvinok = bv_options;
 358     set_program_name(argv[0]);
 359     argp_parse(&argp, argc, argv, 0, 0, &options);
 360
 361     if (options.verify.M == INT_MIN)
 362         options.verify.M = 100;
 363     if (options.verify.m == INT_MAX)
 364         options.verify.m = -100;
 365
 366     result_data_init(&all_result);
 367
 368     while (fgets(path, sizeof(path), stdin)) {
 369         struct result_data result;
 370         FILE *in;
 371         int i;
 372
 373         ++n;
 374         result_data_init(&result);
 375         fprintf(stderr, "%s", path);
 376         *strchr(path, '\n') = '\0';
 377         in = fopen(path, "r");
 378         assert(in);
 379         handle(in, &result, &options.verify);
 380         fclose(in);
 381
 382         if (!options.quiet)
 383             result_data_print(&result, 1);
 384
 385         isl_int_add(all_result.n, all_result.n, result.n);
 386         for (i = 0; i < nr_methods; ++i) {
 387             all_result.RE_sum[i] += result.RE_sum[i];
 388             all_result.ticks[i] += result.ticks[i];
 389             all_result.size[i] += result.size[i];
 390         }
 391
 392         result_data_clear(&result);
 393
 394         if (!options.quiet) {
 395             fprintf(stderr, "average\n");
 396             result_data_print(&all_result, n);
 397         }
 398     }
 399
 400     result_data_clear(&all_result);
 401
 402     barvinok_options_free(bv_options);
 403
 404     return 0;
 405 }