test_bound.c

   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 "evalue_read.h"
   9 #include "verify.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     { BV_BOUND_BERNSTEIN },
  41     { BV_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 struct isl_arg options_arg[] = {
  52 ISL_ARG_CHILD(struct options, verify, NULL, verify_options_arg, NULL)
  53 ISL_ARG_BOOL(struct options, quiet, 'q', "quiet", 0, NULL)
  54 ISL_ARG_END
  55 };
  56
  57 ISL_ARG_DEF(options, struct options, options_arg)
  58
  59 struct result_data {
  60     isl_int                 n;
  61     double                  RE_sum[nr_methods];
  62
  63     my_clock_t              ticks[nr_methods];
  64     size_t                  size[nr_methods];
  65 };
  66
  67 void result_data_init(struct result_data *result)
  68 {
  69     int i;
  70     for (i = 0; i < nr_methods; ++i) {
  71         result->RE_sum[i] = 0;
  72         result->ticks[i] = 0;
  73         result->size[i] = 0;
  74     }
  75     isl_int_init(result->n);
  76 }
  77
  78 void result_data_clear(struct result_data *result)
  79 {
  80     int i;
  81     isl_int_clear(result->n);
  82 }
  83
  84 void result_data_print(struct result_data *result, int n)
  85 {
  86     int i;
  87
  88     fprintf(stderr, "%d", result->ticks[0]/n);
  89     for (i = 1; i < nr_methods; ++i)
  90         fprintf(stderr, ", %d", result->ticks[i]/n);
  91     fprintf(stderr, "\n");
  92
  93     fprintf(stderr, "%zd/%d", result->size[0], n);
  94     for (i = 1; i < nr_methods; ++i)
  95         fprintf(stderr, ", %zd/%d", result->size[i], n);
  96     fprintf(stderr, "\n");
  97
  98     fprintf(stderr, "%g\n", isl_int_get_d(result->n));
  99     fprintf(stderr, "%g", result->RE_sum[0]/isl_int_get_d(result->n));
 100     for (i = 1; i < nr_methods; ++i)
 101         fprintf(stderr, ", %g", result->RE_sum[i]/isl_int_get_d(result->n));
 102     fprintf(stderr, "\n");
 103 }
 104
 105 struct verify_point_bound {
 106         struct verify_point_data vpd;
 107         struct result_data *result;
 108         isl_pw_qpolynomial *pwqp;
 109         isl_pw_qpolynomial_fold **pwf;
 110 };
 111
 112 static int verify_point(__isl_take isl_point *pnt, void *user)
 113 {
 114         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
 115         const struct verify_options *options = vpb->vpd.options;
 116         int i;
 117         unsigned nparam;
 118         isl_int max, min, exact, approx;
 119         isl_int n, d;
 120         isl_qpolynomial *opt;
 121         isl_pw_qpolynomial *pwqp;
 122         int cst;
 123
 124         vpb->vpd.n--;
 125
 126         isl_int_init(max);
 127         isl_int_init(min);
 128         isl_int_init(exact);
 129         isl_int_init(approx);
 130         isl_int_init(n);
 131         isl_int_init(d);
 132
 133         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
 134
 135         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
 136         for (i = 0; i < nparam; ++i) {
 137                 isl_point_get_coordinate(pnt, isl_dim_param, i, &n);
 138                 pwqp = isl_pw_qpolynomial_fix_dim(pwqp, isl_dim_param, i, n);
 139         }
 140
 141         opt = isl_pw_qpolynomial_max(isl_pw_qpolynomial_copy(pwqp));
 142         cst = isl_qpolynomial_is_cst(opt, &n, &d);
 143         isl_qpolynomial_free(opt);
 144         assert(cst == 1);
 145         isl_int_fdiv_q(max, n, d);
 146
 147         opt = isl_pw_qpolynomial_min(pwqp);
 148         cst = isl_qpolynomial_is_cst(opt, &n, &d);
 149         isl_qpolynomial_free(opt);
 150         assert(cst == 1);
 151         isl_int_cdiv_q(min, n, d);
 152
 153         isl_int_sub(exact, max, min);
 154         isl_int_add_ui(exact, exact, 1);
 155
 156         if (options->print_all) {
 157                 fprintf(stderr, "max: "); isl_int_print(stderr, max, 0);
 158                 fprintf(stderr, ", min: "); isl_int_print(stderr, min, 0);
 159                 fprintf(stderr, ", range: "); isl_int_print(stderr, exact, 0);
 160         }
 161
 162         for (i = 0; i < nr_methods; ++i) {
 163                 double error;
 164
 165                 opt = isl_pw_qpolynomial_fold_eval(
 166                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i]),
 167                                 isl_point_copy(pnt));
 168                 cst = isl_qpolynomial_is_cst(opt, &n, &d);
 169                 isl_qpolynomial_free(opt);
 170                 assert(cst == 1);
 171                 isl_int_fdiv_q(max, n, d);
 172
 173                 opt = isl_pw_qpolynomial_fold_eval(
 174                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i + 1]),
 175                                 isl_point_copy(pnt));
 176                 cst = isl_qpolynomial_is_cst(opt, &n, &d);
 177                 isl_qpolynomial_free(opt);
 178                 assert(cst == 1);
 179                 isl_int_cdiv_q(min, n, d);
 180
 181                 isl_int_sub(approx, max, min);
 182                 isl_int_add_ui(approx, approx, 1);
 183                 if (options->print_all) {
 184                         fprintf(stderr, ", "); isl_int_print(stderr, approx, 0);
 185                 }
 186
 187                 assert(isl_int_ge(approx, exact));
 188                 isl_int_sub(approx, approx, exact);
 189
 190                 error = fabs(isl_int_get_d(approx)) / isl_int_get_d(exact);
 191                 if (options->print_all)
 192                         fprintf(stderr, " (%g)", error);
 193                 vpb->result->RE_sum[i] += error;
 194         }
 195
 196         if (options->print_all) {
 197                 fprintf(stderr, "\n");
 198         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 199                 printf("o");
 200                 fflush(stdout);
 201         }
 202
 203         isl_int_clear(max);
 204         isl_int_clear(min);
 205         isl_int_clear(exact);
 206         isl_int_clear(approx);
 207         isl_int_clear(n);
 208         isl_int_clear(d);
 209
 210         isl_point_free(pnt);
 211
 212         return 0;
 213 }
 214
 215 static void test(evalue *EP, unsigned nvar,
 216         __isl_keep isl_pw_qpolynomial_fold **pwf, struct result_data *result,
 217         struct verify_options *options)
 218 {
 219         struct verify_point_bound vpb = { { options }, result };
 220         isl_ctx *ctx;
 221         isl_dim *dim;
 222         isl_set *context;
 223         int r;
 224         int i;
 225         unsigned nparam;
 226
 227         vpb.pwf = pwf;
 228         ctx = isl_pw_qpolynomial_fold_get_ctx(pwf[0]);
 229         nparam = isl_pw_qpolynomial_fold_dim(pwf[0], isl_dim_param);
 230         dim = isl_dim_set_alloc(ctx, nvar + nparam, 0);
 231         vpb.pwqp = isl_pw_qpolynomial_from_evalue(dim, EP);
 232         vpb.pwqp = isl_pw_qpolynomial_move_dims(vpb.pwqp, isl_dim_set, 0,
 233                                                 isl_dim_param, 0, nvar);
 234         context = isl_pw_qpolynomial_domain(isl_pw_qpolynomial_copy(vpb.pwqp));
 235         context = isl_set_remove(context, isl_dim_set, 0, nvar);
 236         context = verify_context_set_bounds(context, options);
 237
 238         r = verify_point_data_init(&vpb.vpd, context);
 239         assert(r == 0);
 240         isl_int_set_si(result->n, vpb.vpd.n);
 241
 242         isl_set_foreach_point(context, verify_point, &vpb);
 243         assert(!vpb.vpd.error);
 244
 245         isl_set_free(context);
 246         isl_pw_qpolynomial_free(vpb.pwqp);
 247
 248         verify_point_data_fini(&vpb.vpd);
 249 }
 250
 251 void handle(FILE *in, struct result_data *result, struct verify_options *options)
 252 {
 253     int i;
 254     evalue *EP, *upper, *lower;
 255     const char **all_vars = NULL;
 256     unsigned nvar;
 257     unsigned nparam;
 258     isl_ctx *ctx = isl_ctx_alloc();
 259     isl_dim *dim;
 260     isl_pw_qpolynomial_fold *pwf[2*nr_methods];
 261
 262     EP = evalue_read_from_file(in, NULL, &all_vars,
 263                                &nvar, &nparam, options->barvinok->MaxRays);
 264     assert(EP);
 265     if (EVALUE_IS_ZERO(*EP)) {
 266         evalue_free(EP);
 267         Free_ParamNames(all_vars, nvar+nparam);
 268         return;
 269     }
 270
 271     upper = evalue_dup(EP);
 272     lower = evalue_dup(EP);
 273     evalue_frac2polynomial(upper, 1, options->barvinok->MaxRays);
 274     evalue_frac2polynomial(lower, -1, options->barvinok->MaxRays);
 275
 276     dim = isl_dim_set_alloc(ctx, nparam, 0);
 277     for (i = 0; i < nr_methods; ++i) {
 278         int j;
 279         struct tms st_cpu;
 280         struct tms en_cpu;
 281
 282         times(&st_cpu);
 283         for (j = 0; j < 2; ++j) {
 284             evalue *poly = j == 0 ? upper : lower;
 285             int sign = j == 0 ? BV_BERNSTEIN_MAX : BV_BERNSTEIN_MIN;
 286             enum isl_fold type = j == 0 ? isl_fold_max : isl_fold_min;
 287             isl_dim *dim_poly;
 288             isl_pw_qpolynomial *pwqp;
 289             options->barvinok->bernstein_optimize = sign;
 290             dim_poly = isl_dim_insert(isl_dim_copy(dim), isl_dim_param, 0, nvar);
 291             pwqp = isl_pw_qpolynomial_from_evalue(dim_poly, poly);
 292             pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_set, 0,
 293                                             isl_dim_param, 0, nvar);
 294             pwf[2*i+j] = isl_pw_qpolynomial_bound(pwqp, type, methods[i].method);
 295         }
 296         times(&en_cpu);
 297         result->ticks[i] = time_diff(&en_cpu, &st_cpu);
 298         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i]);
 299         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i+1]);
 300         if (options->barvinok->verbose) {
 301             isl_printer *p = isl_printer_to_file(ctx, stdout);
 302             for (j = 0; j < 2; ++j) {
 303                 p = isl_printer_print_pw_qpolynomial_fold(p, pwf[2*i+j]);
 304                 p = isl_printer_end_line(p);
 305             }
 306             isl_printer_free(p);
 307         }
 308     }
 309     isl_dim_free(dim);
 310     test(EP, nvar, pwf, result, options);
 311
 312     for (i = 0; i < 2*nr_methods; ++i)
 313         isl_pw_qpolynomial_fold_free(pwf[i]);
 314     evalue_free(EP);
 315     evalue_free(lower);
 316     evalue_free(upper);
 317     Free_ParamNames(all_vars, nvar+nparam);
 318
 319     isl_ctx_free(ctx);
 320 }
 321
 322 int main(int argc, char **argv)
 323 {
 324     char path[PATH_MAX+1];
 325     struct result_data all_result;
 326     int n = 0;
 327     struct options *options = options_new_with_defaults();
 328
 329     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 330
 331     if (options->verify->M == INT_MIN)
 332         options->verify->M = 100;
 333     if (options->verify->m == INT_MAX)
 334         options->verify->m = -100;
 335
 336     result_data_init(&all_result);
 337
 338     while (fgets(path, sizeof(path), stdin)) {
 339         struct result_data result;
 340         FILE *in;
 341         int i;
 342
 343         ++n;
 344         result_data_init(&result);
 345         fprintf(stderr, "%s", path);
 346         *strchr(path, '\n') = '\0';
 347         in = fopen(path, "r");
 348         assert(in);
 349         handle(in, &result, options->verify);
 350         fclose(in);
 351
 352         if (!options->quiet)
 353             result_data_print(&result, 1);
 354
 355         isl_int_add(all_result.n, all_result.n, result.n);
 356         for (i = 0; i < nr_methods; ++i) {
 357             all_result.RE_sum[i] += result.RE_sum[i];
 358             all_result.ticks[i] += result.ticks[i];
 359             all_result.size[i] += result.size[i];
 360         }
 361
 362         result_data_clear(&result);
 363
 364         if (!options->quiet) {
 365             fprintf(stderr, "average\n");
 366             result_data_print(&all_result, n);
 367         }
 368     }
 369
 370     result_data_clear(&all_result);
 371
 372     options_free(options);
 373
 374     return 0;
 375 }