bound.cc

   1 #include <assert.h>
   2 #include <iostream>
   3 #include <bernstein/bernstein.h>
   4 #include <barvinok/evalue.h>
   5 #include <barvinok/options.h>
   6 #include <barvinok/util.h>
   7 #include <barvinok/barvinok.h>
   8 #include <barvinok/bernstein.h>
   9 #include "argp.h"
  10 #include "progname.h"
  11 #include "evalue_convert.h"
  12 #include "evalue_read.h"
  13 #include "verify.h"
  14 #include "range.h"
  15
  16 using std::cout;
  17 using std::cerr;
  18 using std::endl;
  19 using namespace GiNaC;
  20 using namespace bernstein;
  21 using namespace barvinok;
  22
  23 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  24
  25 #define OPT_VARS            (BV_OPT_LAST+1)
  26 #define OPT_SPLIT           (BV_OPT_LAST+2)
  27 #define OPT_LOWER           (BV_OPT_LAST+3)
  28 #define OPT_METHOD          (BV_OPT_LAST+4)
  29 #define OPT_ITERATE         (BV_OPT_LAST+5)
  30
  31 struct argp_option argp_options[] = {
  32     { "split",              OPT_SPLIT,  "int" },
  33     { "variables",          OPT_VARS,   "list", 0,
  34         "comma separated list of variables over which to compute a bound" },
  35     { "lower",              OPT_LOWER,  0, 0,   "compute lower bound instead of upper bound"},
  36     { "optimization-method",    OPT_METHOD,     "bernstein|propagation",
  37         0, "optimization method to use" },
  38     { "iterate",            OPT_ITERATE,        "int", 1,
  39         "exact result by iterating over domain (of specified maximal size)"},
  40     { 0 }
  41 };
  42
  43 struct options {
  44     struct convert_options   convert;
  45     struct verify_options    verify;
  46     char* var_list;
  47     int split;
  48     int lower;
  49 #define METHOD_BERNSTEIN        0
  50 #define METHOD_PROPAGATION      1
  51     int method;
  52     int iterate;
  53 };
  54
  55 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  56 {
  57     struct options *options = (struct options*) state->input;
  58
  59     switch (key) {
  60     case ARGP_KEY_INIT:
  61         state->child_inputs[0] = &options->convert;
  62         state->child_inputs[1] = &options->verify;
  63         state->child_inputs[2] = options->verify.barvinok;
  64         options->var_list = NULL;
  65         options->split = 0;
  66         options->lower = 0;
  67         options->method = METHOD_BERNSTEIN;
  68         options->iterate = 0;
  69         break;
  70     case OPT_VARS:
  71         options->var_list = strdup(arg);
  72         break;
  73     case OPT_SPLIT:
  74         options->split = atoi(arg);
  75         break;
  76     case OPT_LOWER:
  77         options->lower = 1;
  78         break;
  79     case OPT_METHOD:
  80         if (!strcmp(arg, "bernstein"))
  81             options->method = METHOD_BERNSTEIN;
  82         else if (!strcmp(arg, "propagation"))
  83             options->method = METHOD_PROPAGATION;
  84         else
  85             argp_error(state, "unknown value for --optimization-method option");
  86         break;
  87     case OPT_ITERATE:
  88         if (!arg)
  89             options->iterate = -1;
  90         else
  91             options->iterate = atoi(arg);
  92         break;
  93     default:
  94         return ARGP_ERR_UNKNOWN;
  95     }
  96     return 0;
  97 }
  98
  99 static int check_poly_max(const struct check_poly_data *data,
 100                           int nparam, Value *z,
 101                           const struct verify_options *options);
 102
 103 struct check_poly_max_data : public check_EP_data {
 104     piecewise_lst        *pl;
 105
 106     check_poly_max_data(evalue *EP, piecewise_lst *pl) : pl(pl) {
 107         this->EP = EP;
 108         this->cp.check = check_poly_max;
 109     }
 110 };
 111
 112 static int check_poly_max(const struct check_poly_data *data,
 113                           int nparam, Value *z,
 114                           const struct verify_options *options)
 115 {
 116     int k;
 117     int ok;
 118     const check_poly_max_data *max_data;
 119     max_data = (const check_poly_max_data *)data;
 120     const char *minmax;
 121     Value m, n, d;
 122     value_init(m);
 123     value_init(n);
 124     value_init(d);
 125     int sign;
 126
 127     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
 128         minmax = "max";
 129         sign = 1;
 130     } else {
 131         minmax = "min";
 132         sign = -1;
 133     }
 134
 135     max_data->pl->evaluate(nparam, z, &n, &d);
 136     if (sign > 0)
 137         mpz_fdiv_q(m, n, d);
 138     else
 139         mpz_cdiv_q(m, n, d);
 140
 141     if (options->print_all) {
 142         printf("%s(", minmax);
 143         value_print(stdout, VALUE_FMT, z[0]);
 144         for (k = 1; k < nparam; ++k) {
 145             printf(", ");
 146             value_print(stdout, VALUE_FMT, z[k]);
 147         }
 148         printf(") = ");
 149         value_print(stdout, VALUE_FMT, n);
 150         if (value_notone_p(d)) {
 151             printf("/");
 152             value_print(stdout, VALUE_FMT, d);
 153         }
 154         printf(" (");
 155         value_print(stdout, VALUE_FMT, m);
 156         printf(")");
 157     }
 158
 159     evalue_optimum(max_data, &n, sign);
 160
 161     if (sign > 0)
 162         ok = value_ge(m, n);
 163     else
 164         ok = value_le(m, n);
 165
 166     if (options->print_all) {
 167         printf(", %s(EP) = ", minmax);
 168         value_print(stdout, VALUE_FMT, n);
 169         printf(". ");
 170     }
 171
 172     if (!ok) {
 173         printf("\n");
 174         fflush(stdout);
 175         fprintf(stderr, "Error !\n");
 176         fprintf(stderr, "%s(", minmax);
 177         value_print(stderr, VALUE_FMT, z[0]);
 178         for (k = 1; k < nparam; ++k) {
 179             fprintf(stderr,", ");
 180             value_print(stderr, VALUE_FMT, z[k]);
 181         }
 182         fprintf(stderr, ") should be ");
 183         if (sign > 0)
 184             fprintf(stderr, "greater");
 185         else
 186             fprintf(stderr, "smaller");
 187         fprintf(stderr, " than or equal to ");
 188         value_print(stderr, VALUE_FMT, n);
 189         fprintf(stderr, ", while pl eval gives ");
 190         value_print(stderr, VALUE_FMT, m);
 191         fprintf(stderr, ".\n");
 192         cerr << *max_data->pl << endl;
 193     } else if (options->print_all)
 194         printf("OK.\n");
 195
 196     value_clear(m);
 197     value_clear(n);
 198     value_clear(d);
 199
 200     return ok;
 201 }
 202
 203 static int verify(piecewise_lst *pl, evalue *EP, unsigned nvar, unsigned nparam,
 204                   struct verify_options *options)
 205 {
 206     check_poly_max_data data(EP, pl);
 207     return !check_EP(&data, nvar, nparam, options);
 208 }
 209
 210 static piecewise_lst *iterate(evalue *EP, unsigned nvar, Polyhedron *C,
 211                                 struct options *options)
 212 {
 213     assert(C->Dimension == 0);
 214     Vector *p = Vector_Alloc(nvar+2);
 215     value_set_si(p->p[nvar+1], 1);
 216     Value opt;
 217
 218     value_init(opt);
 219     struct check_EP_data data;
 220     data.cp.z = p->p;
 221     data.cp.check = NULL;
 222
 223     data.EP = EP;
 224     check_EP_set_scan(&data, C, options->verify.barvinok->MaxRays);
 225     int sign = options->verify.barvinok->bernstein_optimize;
 226     evalue_optimum(&data, &opt, sign);
 227     check_EP_clear_scan(&data);
 228
 229     exvector params;
 230     piecewise_lst *pl = new piecewise_lst(params, sign);
 231     pl->add_guarded_lst(Polyhedron_Copy(C), lst(value2numeric(opt)));
 232
 233     value_clear(opt);
 234     Vector_Free(p);
 235
 236     return pl;
 237 }
 238
 239 /*
 240  * Split (partition) EP into a partition with (sub)domains containing
 241  * size integer points or less and a partition with (sub)domains
 242  * containing more integer points.
 243  */
 244 static void split_on_domain_size(evalue *EP, evalue **EP_less, evalue **EP_more,
 245                                  int size, barvinok_options *options)
 246 {
 247     assert(value_zero_p(EP->d));
 248     assert(EP->x.p->type == partition);
 249     assert(EP->x.p->size >= 2);
 250
 251     struct evalue_section *s_less = new evalue_section[EP->x.p->size/2];
 252     struct evalue_section *s_more = new evalue_section[EP->x.p->size/2];
 253
 254     int n_less = 0;
 255     int n_more = 0;
 256
 257     Value c;
 258     value_init(c);
 259
 260     for (int i = 0; i < EP->x.p->size/2; ++i) {
 261         Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 262         Polyhedron *D_less = NULL;
 263         Polyhedron *D_more = NULL;
 264         Polyhedron **next_less = &D_less;
 265         Polyhedron **next_more = &D_more;
 266
 267         for (Polyhedron *P = D; P; P = P->next) {
 268             Polyhedron *next = P->next;
 269             P->next = NULL;
 270             barvinok_count_with_options(P, &c, options);
 271             P->next = next;
 272
 273             if (value_zero_p(c))
 274                 continue;
 275
 276             if (value_pos_p(c) && value_cmp_si(c, size) <= 0) {
 277                 *next_less = Polyhedron_Copy(P);
 278                 next_less = &(*next_less)->next;
 279             } else {
 280                 *next_more = Polyhedron_Copy(P);
 281                 next_more = &(*next_more)->next;
 282             }
 283         }
 284
 285         if (D_less) {
 286             s_less[n_less].D = D_less;
 287             s_less[n_less].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 288             n_less++;
 289         } else {
 290             s_more[n_more].D = D_more;
 291             s_more[n_more].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 292             n_more++;
 293         }
 294     }
 295
 296     value_clear(c);
 297
 298     *EP_less = evalue_from_section_array(s_less, n_less);
 299     *EP_more = evalue_from_section_array(s_more, n_more);
 300
 301     delete [] s_less;
 302     delete [] s_more;
 303 }
 304
 305 static piecewise_lst *optimize(evalue *EP, unsigned nvar, Polyhedron *C,
 306                                 const exvector& params,
 307                                 struct options *options)
 308 {
 309     piecewise_lst *pl = NULL;
 310     if (options->iterate > 0) {
 311         evalue *EP_less = NULL;
 312         evalue *EP_more = NULL;
 313         piecewise_lst *pl_more = NULL;
 314         split_on_domain_size(EP, &EP_less, &EP_more, options->iterate,
 315                                 options->verify.barvinok);
 316         if (!EVALUE_IS_ZERO(*EP_less)) {
 317             options->iterate = -1;
 318             pl = optimize(EP_less, nvar, C, params, options);
 319         }
 320         if (!EVALUE_IS_ZERO(*EP_more)) {
 321             options->iterate = 0;
 322             pl_more = optimize(EP_more, nvar, C, params, options);
 323         }
 324         evalue_free(EP_less);
 325         evalue_free(EP_more);
 326         if (!pl)
 327             return pl_more;
 328         if (!pl_more)
 329             return pl;
 330         pl->combine(*pl_more);
 331         delete pl_more;
 332         return pl;
 333     }
 334     if (options->iterate)
 335         pl = iterate(EP, nvar, C, options);
 336     else if (options->method == METHOD_BERNSTEIN) {
 337         pl = evalue_bernstein_coefficients(NULL, EP, C, params,
 338                                            options->verify.barvinok);
 339         if (options->lower)
 340             pl->minimize();
 341         else
 342             pl->maximize();
 343     } else
 344         pl = evalue_range_propagation(NULL, EP, params,
 345                                       options->verify.barvinok);
 346     return pl;
 347 }
 348
 349 static int optimize(evalue *EP, const char **all_vars,
 350                     unsigned nvar, unsigned nparam, struct options *options)
 351 {
 352     Polyhedron *U;
 353     piecewise_lst *pl = NULL;
 354     U = Universe_Polyhedron(nparam);
 355     int print_solution = 1;
 356     int result = 0;
 357
 358     exvector params;
 359     params = constructParameterVector(all_vars+nvar, nparam);
 360
 361     if (options->verify.verify) {
 362         verify_options_set_range(&options->verify, nvar+nparam);
 363         if (!options->verify.barvinok->verbose)
 364             print_solution = 0;
 365     }
 366
 367     if (options->lower)
 368         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
 369     else
 370         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
 371     pl = optimize(EP, nvar, U, params, options);
 372     assert(pl);
 373     if (print_solution)
 374         cout << *pl << endl;
 375     if (options->verify.verify) {
 376         result = verify(pl, EP, nvar, nparam, &options->verify);
 377     }
 378     delete pl;
 379
 380     Polyhedron_Free(U);
 381
 382     return result;
 383 }
 384
 385 int main(int argc, char **argv)
 386 {
 387     evalue *EP;
 388     const char **all_vars = NULL;
 389     unsigned nvar;
 390     unsigned nparam;
 391     struct options options;
 392     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 393     static struct argp_child argp_children[] = {
 394         { &convert_argp,        0,      "input conversion",     1 },
 395         { &verify_argp,         0,      "verification",         2 },
 396         { &barvinok_argp,       0,      "barvinok options",     3 },
 397         { 0 }
 398     };
 399     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 400     int result = 0;
 401
 402     options.verify.barvinok = bv_options;
 403     set_program_name(argv[0]);
 404     argp_parse(&argp, argc, argv, 0, 0, &options);
 405
 406     EP = evalue_read_from_file(stdin, options.var_list, &all_vars,
 407                                &nvar, &nparam, bv_options->MaxRays);
 408     assert(EP);
 409
 410     if (options.split)
 411         evalue_split_periods(EP, options.split, bv_options->MaxRays);
 412
 413     evalue_convert(EP, &options.convert, bv_options->verbose,
 414                         nvar+nparam, all_vars);
 415
 416     if (EVALUE_IS_ZERO(*EP))
 417         print_evalue(stdout, EP, all_vars);
 418     else
 419         result = optimize(EP, all_vars, nvar, nparam, &options);
 420
 421     evalue_free(EP);
 422
 423     if (options.var_list)
 424         free(options.var_list);
 425     Free_ParamNames(all_vars, nvar+nparam);
 426     barvinok_options_free(bv_options);
 427     return result;
 428 }