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 #define METHOD_BERNSTEIN        0
  48 #define METHOD_PROPAGATION      1
  49
  50 static struct {
  51     int     method;
  52 } methods[] = {
  53     { METHOD_BERNSTEIN },
  54     { METHOD_PROPAGATION },
  55 };
  56
  57 #define nr_methods (sizeof(methods)/sizeof(*methods))
  58
  59 struct argp_option argp_options[] = {
  60     { "quiet",      'q' },
  61     { 0 }
  62 };
  63
  64 struct options {
  65     struct verify_options    verify;
  66     int quiet;
  67 };
  68
  69 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  70 {
  71     struct options *options = (struct options*) state->input;
  72
  73     switch (key) {
  74     case ARGP_KEY_INIT:
  75         state->child_inputs[0] = options->verify.barvinok;
  76         state->child_inputs[1] = &options->verify;
  77         options->quiet = 0;
  78         break;
  79     case 'q':
  80         options->quiet = 1;
  81         break;
  82     default:
  83         return ARGP_ERR_UNKNOWN;
  84     }
  85     return 0;
  86 }
  87
  88 struct result_data {
  89     Value                   n;
  90     double                  RE_sum[nr_methods];
  91
  92     my_clock_t              ticks[nr_methods];
  93     size_t                  size[nr_methods];
  94 };
  95
  96 void result_data_init(struct result_data *result)
  97 {
  98     int i;
  99     for (i = 0; i < nr_methods; ++i) {
 100         result->RE_sum[i] = 0;
 101         result->ticks[i] = 0;
 102         result->size[i] = 0;
 103     }
 104     value_init(result->n);
 105 }
 106
 107 void result_data_clear(struct result_data *result)
 108 {
 109     int i;
 110     value_clear(result->n);
 111 }
 112
 113 void result_data_print(struct result_data *result, int n)
 114 {
 115     int i;
 116
 117     fprintf(stderr, "%d", result->ticks[0]/n);
 118     for (i = 1; i < nr_methods; ++i)
 119         fprintf(stderr, ", %d", result->ticks[i]/n);
 120     fprintf(stderr, "\n");
 121
 122     fprintf(stderr, "%zd/%d", result->size[0], n);
 123     for (i = 1; i < nr_methods; ++i)
 124         fprintf(stderr, ", %zd/%d", result->size[i], n);
 125     fprintf(stderr, "\n");
 126
 127     fprintf(stderr, "%g\n", VALUE_TO_DOUBLE(result->n));
 128     fprintf(stderr, "%g", result->RE_sum[0]/VALUE_TO_DOUBLE(result->n));
 129     for (i = 1; i < nr_methods; ++i)
 130         fprintf(stderr, ", %g", result->RE_sum[i]/VALUE_TO_DOUBLE(result->n));
 131     fprintf(stderr, "\n");
 132 }
 133
 134 static int test_bound(const struct check_poly_data *data,
 135                           int nparam, Value *z,
 136                           const struct verify_options *options);
 137
 138 struct test_bound_data : public check_EP_data {
 139     piecewise_lst       **pl;
 140     struct result_data   *result;
 141
 142     test_bound_data(evalue *EP, piecewise_lst **pl, result_data *result) :
 143                         pl(pl), result(result) {
 144         this->EP = EP;
 145         this->cp.check = test_bound;
 146     }
 147 };
 148
 149 static int test_bound(const struct check_poly_data *data,
 150                           int nparam, Value *z,
 151                           const struct verify_options *options)
 152 {
 153     const test_bound_data *tb_data = (const test_bound_data *)data;
 154     Value max, min, exact, approx;
 155     Value n, d;
 156
 157     value_init(exact);
 158     value_init(approx);
 159     value_init(max);
 160     value_init(min);
 161     value_init(n);
 162     value_init(d);
 163
 164     evalue_optimum(tb_data, &max, 1);
 165     evalue_optimum(tb_data, &min, -1);
 166     value_assign(exact, max);
 167     value_subtract(exact, exact, min);
 168     value_add_int(exact, exact, 1);
 169
 170     if (options->print_all) {
 171         value_print(stderr, "max: "VALUE_FMT, max);
 172         value_print(stderr, ", min: "VALUE_FMT, min);
 173         value_print(stderr, ", range: "VALUE_FMT, exact);
 174     }
 175
 176     value_increment(tb_data->result->n, tb_data->result->n);
 177     for (int i = 0; i < nr_methods; ++i) {
 178         double error;
 179
 180         tb_data->pl[2*i]->evaluate(nparam, z, &n, &d);
 181         mpz_fdiv_q(max, n, d);
 182         tb_data->pl[2*i+1]->evaluate(nparam, z, &n, &d);
 183         mpz_cdiv_q(min, n, d);
 184
 185         value_assign(approx, max);
 186         value_subtract(approx, approx, min);
 187         value_add_int(approx, approx, 1);
 188         if (options->print_all)
 189             value_print(stderr, ", "VALUE_FMT, approx);
 190
 191         assert(value_ge(approx, exact));
 192         value_subtract(approx, approx, exact);
 193
 194         error = fabs(VALUE_TO_DOUBLE(approx)) / VALUE_TO_DOUBLE(exact);
 195         if (options->print_all)
 196             fprintf(stderr, " (%g)", error);
 197         tb_data->result->RE_sum[i] += error;
 198     }
 199
 200     if (options->print_all)
 201         fprintf(stderr, "\n");
 202
 203     value_clear(n);
 204     value_clear(d);
 205     value_clear(max);
 206     value_clear(min);
 207     value_clear(exact);
 208     value_clear(approx);
 209
 210     return 1;
 211 }
 212
 213 static void test(evalue *EP, unsigned nvar, unsigned nparam,
 214                  piecewise_lst **pl, struct result_data *result,
 215                  struct verify_options *options)
 216 {
 217     test_bound_data data(EP, pl, result);
 218     check_EP(&data, nvar, nparam, options);
 219 }
 220
 221 static int number_of_polynomials(piecewise_lst *pl)
 222 {
 223     int n = 0;
 224     for (int i = 0; i < pl->list.size(); ++i)
 225         n += pl->list[i].second.nops();
 226     return n;
 227 }
 228
 229 void handle(FILE *in, struct result_data *result, struct verify_options *options)
 230 {
 231     evalue *EP, *upper, *lower;
 232     const char **all_vars = NULL;
 233     unsigned nvar;
 234     unsigned nparam;
 235     Polyhedron *U;
 236     exvector params;
 237     piecewise_lst *pl[2*nr_methods];
 238
 239     EP = evalue_read_from_file(in, NULL, &all_vars,
 240                                &nvar, &nparam, options->barvinok->MaxRays);
 241     assert(EP);
 242     if (EVALUE_IS_ZERO(*EP)) {
 243         evalue_free(EP);
 244         Free_ParamNames(all_vars, nvar+nparam);
 245         return;
 246     }
 247
 248     upper = evalue_dup(EP);
 249     lower = evalue_dup(EP);
 250     evalue_frac2polynomial(upper, 1, options->barvinok->MaxRays);
 251     evalue_frac2polynomial(lower, -1, options->barvinok->MaxRays);
 252
 253     U = Universe_Polyhedron(nparam);
 254     params = constructParameterVector(all_vars+nvar, nparam);
 255
 256     for (int i = 0; i < nr_methods; ++i) {
 257         struct tms st_cpu;
 258         struct tms en_cpu;
 259
 260         times(&st_cpu);
 261         for (int j = 0; j < 2; ++j) {
 262             evalue *poly = j == 0 ? upper : lower;
 263             int sign = j == 0 ? BV_BERNSTEIN_MAX : BV_BERNSTEIN_MIN;
 264             options->barvinok->bernstein_optimize = sign;
 265             if (methods[i].method == METHOD_BERNSTEIN) {
 266                 pl[2*i+j] = evalue_bernstein_coefficients(NULL, poly, U, params,
 267                                                    options->barvinok);
 268                 if (sign == BV_BERNSTEIN_MIN)
 269                     pl[2*i+j]->minimize();
 270                 else
 271                     pl[2*i+j]->maximize();
 272             } else {
 273                 pl[2*i+j] = evalue_range_propagation(NULL, poly, params,
 274                                               options->barvinok);
 275             }
 276         }
 277         times(&en_cpu);
 278         result->ticks[i] = time_diff(&en_cpu, &st_cpu);
 279         if (options->barvinok->verbose)
 280             for (int j = 0; j < 2; ++j)
 281                 cerr << *pl[2*i+j] << endl;
 282         result->size[i] = number_of_polynomials(pl[2*i]);
 283         result->size[i] += number_of_polynomials(pl[2*i+1]);
 284     }
 285     test(EP, nvar, nparam, pl, result, options);
 286
 287     for (int i = 0; i < 2*nr_methods; ++i)
 288         delete pl[i];
 289     Polyhedron_Free(U);
 290     evalue_free(EP);
 291     evalue_free(lower);
 292     evalue_free(upper);
 293     Free_ParamNames(all_vars, nvar+nparam);
 294 }
 295
 296 int main(int argc, char **argv)
 297 {
 298     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 299     char path[PATH_MAX+1];
 300     struct result_data all_result;
 301     int n = 0;
 302     static struct argp_child argp_children[] = {
 303         { &barvinok_argp,       0,      0,              0 },
 304         { &verify_argp,         0,      "verification",         BV_GRP_LAST+1 },
 305         { 0 }
 306     };
 307     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 308     struct options options;
 309
 310     options.verify.barvinok = bv_options;
 311     set_program_name(argv[0]);
 312     argp_parse(&argp, argc, argv, 0, 0, &options);
 313
 314     if (options.verify.M == INT_MIN)
 315         options.verify.M = 100;
 316     if (options.verify.m == INT_MAX)
 317         options.verify.m = -100;
 318
 319     result_data_init(&all_result);
 320
 321     while (fgets(path, sizeof(path), stdin)) {
 322         struct result_data result;
 323         FILE *in;
 324         int i;
 325
 326         ++n;
 327         result_data_init(&result);
 328         fprintf(stderr, "%s", path);
 329         *strchr(path, '\n') = '\0';
 330         in = fopen(path, "r");
 331         assert(in);
 332         handle(in, &result, &options.verify);
 333         fclose(in);
 334
 335         if (!options.quiet)
 336             result_data_print(&result, 1);
 337
 338         value_addto(all_result.n, all_result.n, result.n);
 339         for (i = 0; i < nr_methods; ++i) {
 340             all_result.RE_sum[i] += result.RE_sum[i];
 341             all_result.ticks[i] += result.ticks[i];
 342             all_result.size[i] += result.size[i];
 343         }
 344
 345         result_data_clear(&result);
 346
 347         if (!options.quiet) {
 348             fprintf(stderr, "average\n");
 349             result_data_print(&all_result, n);
 350         }
 351     }
 352
 353     result_data_clear(&all_result);
 354
 355     barvinok_options_free(bv_options);
 356
 357     return 0;
 358 }