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/%d", result->size[0], n);
  92     for (i = 1; i < nr_methods; ++i)
  93         fprintf(stderr, ", %d/%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 static int number_of_polynomials(piecewise_lst *pl)
 191 {
 192     int n = 0;
 193     for (int i = 0; i < pl->list.size(); ++i)
 194         n += pl->list[i].second.nops();
 195     return n;
 196 }
 197
 198 void handle(FILE *in, struct result_data *result, struct verify_options *options)
 199 {
 200     evalue *EP, *upper, *lower;
 201     char **all_vars = NULL;
 202     unsigned nvar;
 203     unsigned nparam;
 204     Polyhedron *U;
 205     exvector params;
 206     piecewise_lst *pl[2*nr_methods];
 207
 208     EP = evalue_read_from_file(in, NULL, &all_vars,
 209                                &nvar, &nparam, options->barvinok->MaxRays);
 210     assert(EP);
 211     if (EVALUE_IS_ZERO(*EP)) {
 212         evalue_free(EP);
 213         Free_ParamNames(all_vars, nvar+nparam);
 214         return;
 215     }
 216
 217     upper = evalue_dup(EP);
 218     lower = evalue_dup(EP);
 219     evalue_frac2polynomial(upper, 1, options->barvinok->MaxRays);
 220     evalue_frac2polynomial(lower, -1, options->barvinok->MaxRays);
 221
 222     U = Universe_Polyhedron(nparam);
 223     params = constructParameterVector(all_vars+nvar, nparam);
 224
 225     for (int i = 0; i < nr_methods; ++i) {
 226         struct tms st_cpu;
 227         struct tms en_cpu;
 228
 229         times(&st_cpu);
 230         for (int j = 0; j < 2; ++j) {
 231             evalue *poly = j == 0 ? upper : lower;
 232             int sign = j == 0 ? BV_BERNSTEIN_MAX : BV_BERNSTEIN_MIN;
 233             options->barvinok->bernstein_optimize = sign;
 234             if (methods[i].method == METHOD_BERNSTEIN) {
 235                 pl[2*i+j] = evalue_bernstein_coefficients(NULL, poly, U, params,
 236                                                    options->barvinok);
 237                 if (sign == BV_BERNSTEIN_MIN)
 238                     pl[2*i+j]->minimize();
 239                 else
 240                     pl[2*i+j]->maximize();
 241             } else {
 242                 pl[2*i+j] = evalue_range_propagation(NULL, poly, params,
 243                                               options->barvinok);
 244             }
 245         }
 246         times(&en_cpu);
 247         result->ticks[i] = en_cpu.tms_utime - st_cpu.tms_utime;
 248         if (options->barvinok->verbose)
 249             for (int j = 0; j < 2; ++j)
 250                 cerr << *pl[2*i+j] << endl;
 251         result->size[i] = number_of_polynomials(pl[2*i]);
 252         result->size[i] += number_of_polynomials(pl[2*i+1]);
 253     }
 254     test(EP, nvar, nparam, pl, result, options);
 255
 256     for (int i = 0; i < 2*nr_methods; ++i)
 257         delete pl[i];
 258     Polyhedron_Free(U);
 259     evalue_free(EP);
 260     evalue_free(lower);
 261     evalue_free(upper);
 262     Free_ParamNames(all_vars, nvar+nparam);
 263 }
 264
 265 int main(int argc, char **argv)
 266 {
 267     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 268     char path[PATH_MAX+1];
 269     struct result_data all_result;
 270     int n = 0;
 271     static struct argp_child argp_children[] = {
 272         { &barvinok_argp,       0,      0,              0 },
 273         { &verify_argp,         0,      "verification",         BV_GRP_LAST+1 },
 274         { 0 }
 275     };
 276     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 277     struct options options;
 278
 279     options.verify.barvinok = bv_options;
 280     set_program_name(argv[0]);
 281     argp_parse(&argp, argc, argv, 0, 0, &options);
 282
 283     if (options.verify.M == INT_MIN)
 284         options.verify.M = 100;
 285     if (options.verify.m == INT_MAX)
 286         options.verify.m = -100;
 287
 288     result_data_init(&all_result);
 289
 290     while (fgets(path, sizeof(path), stdin)) {
 291         struct result_data result;
 292         FILE *in;
 293         int i;
 294
 295         ++n;
 296         result_data_init(&result);
 297         fprintf(stderr, "%s", path);
 298         *strchr(path, '\n') = '\0';
 299         in = fopen(path, "r");
 300         assert(in);
 301         handle(in, &result, &options.verify);
 302         fclose(in);
 303
 304         result_data_print(&result, 1);
 305
 306         value_addto(all_result.n, all_result.n, result.n);
 307         for (i = 0; i < nr_methods; ++i) {
 308             all_result.RE_sum[i] += result.RE_sum[i];
 309             all_result.ticks[i] += result.ticks[i];
 310             all_result.size[i] += result.size[i];
 311         }
 312
 313         result_data_clear(&result);
 314
 315         fprintf(stderr, "average\n");
 316         result_data_print(&all_result, n);
 317     }
 318
 319     result_data_clear(&all_result);
 320
 321     barvinok_options_free(bv_options);
 322
 323     return 0;
 324 }