test_bound.cc

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