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