test_bound.cc

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