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