test_bound.c

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