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_ITERATE         (BV_OPT_LAST+4)
  29
  30 struct argp_option argp_options[] = {
  31     { "split",              OPT_SPLIT,  "int" },
  32     { "variables",          OPT_VARS,   "list", 0,
  33         "comma separated list of variables over which to compute a bound" },
  34     { "lower",              OPT_LOWER,  0, 0,   "compute lower bound instead of upper bound"},
  35     { "iterate",            OPT_ITERATE,        "int", 1,
  36         "exact result by iterating over domain (of specified maximal size)"},
  37     { 0 }
  38 };
  39
  40 struct options {
  41     struct convert_options   convert;
  42     struct verify_options    verify;
  43     char* var_list;
  44     int split;
  45     int lower;
  46     int iterate;
  47 };
  48
  49 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  50 {
  51     struct options *options = (struct options*) state->input;
  52
  53     switch (key) {
  54     case ARGP_KEY_INIT:
  55         state->child_inputs[0] = &options->convert;
  56         state->child_inputs[1] = &options->verify;
  57         state->child_inputs[2] = options->verify.barvinok;
  58         options->var_list = NULL;
  59         options->split = 0;
  60         options->lower = 0;
  61         options->iterate = 0;
  62         break;
  63     case OPT_VARS:
  64         options->var_list = strdup(arg);
  65         break;
  66     case OPT_SPLIT:
  67         options->split = atoi(arg);
  68         break;
  69     case OPT_LOWER:
  70         options->lower = 1;
  71         break;
  72     case OPT_ITERATE:
  73         if (!arg)
  74             options->iterate = -1;
  75         else
  76             options->iterate = atoi(arg);
  77         break;
  78     default:
  79         return ARGP_ERR_UNKNOWN;
  80     }
  81     return 0;
  82 }
  83
  84 struct verify_point_bound {
  85         struct verify_point_data vpd;
  86         isl_pw_qpolynomial *pwqp;
  87         isl_pw_qpolynomial_fold *pwf;
  88 };
  89
  90 static int verify_point(__isl_take isl_point *pnt, void *user)
  91 {
  92         int i;
  93         unsigned nparam;
  94         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
  95         isl_int v, n, d, b, t;
  96         isl_pw_qpolynomial *pwqp;
  97         isl_qpolynomial *bound;
  98         isl_qpolynomial *opt;
  99         const char *minmax;
 100         int sign;
 101         int ok;
 102         int cst;
 103         FILE *out = vpb->vpd.options->print_all ? stdout : stderr;
 104
 105         vpb->vpd.n--;
 106
 107         if (vpb->vpd.options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
 108                 minmax = "max";
 109                 sign = 1;
 110         } else {
 111                 minmax = "min";
 112                 sign = -1;
 113         }
 114
 115         isl_int_init(b);
 116         isl_int_init(t);
 117         isl_int_init(v);
 118         isl_int_init(n);
 119         isl_int_init(d);
 120
 121         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
 122
 123         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
 124         for (i = 0; i < nparam; ++i) {
 125                 isl_point_get_coordinate(pnt, isl_dim_param, i, &v);
 126                 pwqp = isl_pw_qpolynomial_fix_dim(pwqp, isl_dim_param, i, v);
 127         }
 128
 129         bound = isl_pw_qpolynomial_fold_eval(isl_pw_qpolynomial_fold_copy(vpb->pwf),
 130                                                 isl_point_copy(pnt));
 131
 132         if (sign > 0)
 133                 opt = isl_pw_qpolynomial_max(pwqp);
 134         else
 135                 opt = isl_pw_qpolynomial_min(pwqp);
 136
 137         cst = isl_qpolynomial_is_cst(opt, &n, &d);
 138         if (cst != 1)
 139                 goto error;
 140         if (sign > 0)
 141                 isl_int_fdiv_q(v, n, d);
 142         else
 143                 isl_int_cdiv_q(v, n, d);
 144
 145         cst = isl_qpolynomial_is_cst(bound, &n, &d);
 146         if (cst != 1)
 147                 goto error;
 148         if (sign > 0)
 149                 isl_int_fdiv_q(b, n, d);
 150         else
 151                 isl_int_cdiv_q(b, n, d);
 152
 153         if (sign > 0)
 154                 ok = value_ge(b, v);
 155         else
 156                 ok = value_le(b, v);
 157
 158         if (vpb->vpd.options->print_all || !ok) {
 159                 fprintf(out, "%s(", minmax);
 160                 for (i = 0; i < nparam; ++i) {
 161                         if (i)
 162                                 fprintf(out, ", ");
 163                         isl_point_get_coordinate(pnt, isl_dim_param, i, &t);
 164                         isl_int_print(out, t, 0);
 165                 }
 166                 fprintf(out, ") = ");
 167                 isl_int_print(out, n, 0);
 168                 if (!isl_int_is_one(d)) {
 169                         fprintf(out, "/");
 170                         isl_int_print(out, d, 0);
 171                         fprintf(out, " (");
 172                         isl_int_print(out, b, 0);
 173                         fprintf(out, ")");
 174                 }
 175                 fprintf(out, ", %s(EP) = ", minmax);
 176                 isl_int_print(out, v, 0);
 177                 if (ok)
 178                         fprintf(out, ". OK\n");
 179                 else
 180                         fprintf(out, ". NOT OK\n");
 181         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 182                 printf("o");
 183                 fflush(stdout);
 184         }
 185
 186         if (0) {
 187 error:
 188                 ok = 0;
 189         }
 190
 191         isl_qpolynomial_free(bound);
 192         isl_qpolynomial_free(opt);
 193         isl_point_free(pnt);
 194
 195         isl_int_clear(t);
 196         isl_int_clear(d);
 197         isl_int_clear(n);
 198         isl_int_clear(v);
 199         isl_int_clear(b);
 200
 201         if (!ok)
 202                 vpb->vpd.error = 1;
 203
 204         if (vpb->vpd.options->continue_on_error)
 205                 ok = 1;
 206
 207         return (vpb->vpd.n >= 1 && ok) ? 0 : -1;
 208 }
 209
 210 static int verify(piecewise_lst *pl, evalue *EP, unsigned nvar,
 211                   const GiNaC::exvector &params,
 212                   struct verify_options *options)
 213 {
 214         struct verify_point_bound vpb = { { options } };
 215         isl_ctx *ctx = isl_ctx_alloc();
 216         isl_dim *dim;
 217         isl_set *context;
 218         int r;
 219
 220         dim = isl_dim_set_alloc(ctx, params.size(), 0);
 221         vpb.pwf = isl_pw_qpolynomial_fold_from_ginac(dim, pl, params);
 222         dim = isl_dim_set_alloc(ctx, nvar + params.size(), 0);
 223         vpb.pwqp = isl_pw_qpolynomial_from_evalue(dim, EP);
 224         vpb.pwqp = isl_pw_qpolynomial_move(vpb.pwqp, isl_dim_set, 0,
 225                                                 isl_dim_param, 0, nvar);
 226         context = isl_pw_qpolynomial_fold_domain(
 227                                         isl_pw_qpolynomial_fold_copy(vpb.pwf));
 228         context = verify_context_set_bounds(context, options);
 229
 230         r = verify_point_data_init(&vpb.vpd, context);
 231
 232         if (r == 0)
 233                 isl_set_foreach_point(context, verify_point, &vpb);
 234         if (vpb.vpd.error)
 235                 r = -1;
 236
 237         isl_set_free(context);
 238         isl_pw_qpolynomial_free(vpb.pwqp);
 239         isl_pw_qpolynomial_fold_free(vpb.pwf);
 240
 241         isl_ctx_free(ctx);
 242
 243         verify_point_data_fini(&vpb.vpd);
 244
 245         return r;
 246 }
 247
 248 static piecewise_lst *iterate(evalue *EP, unsigned nvar, Polyhedron *C,
 249                                 struct options *options)
 250 {
 251     assert(C->Dimension == 0);
 252     Vector *p = Vector_Alloc(nvar+2);
 253     value_set_si(p->p[nvar+1], 1);
 254     Value opt;
 255
 256     value_init(opt);
 257     struct check_EP_data data;
 258     data.cp.z = p->p;
 259     data.cp.check = NULL;
 260
 261     data.EP = EP;
 262     check_EP_set_scan(&data, C, options->verify.barvinok->MaxRays);
 263     int sign = options->verify.barvinok->bernstein_optimize;
 264     evalue_optimum(&data, &opt, sign);
 265     check_EP_clear_scan(&data);
 266
 267     exvector params;
 268     piecewise_lst *pl = new piecewise_lst(params, sign);
 269     pl->add_guarded_lst(Polyhedron_Copy(C), lst(value2numeric(opt)));
 270
 271     value_clear(opt);
 272     Vector_Free(p);
 273
 274     return pl;
 275 }
 276
 277 /*
 278  * Split (partition) EP into a partition with (sub)domains containing
 279  * size integer points or less and a partition with (sub)domains
 280  * containing more integer points.
 281  */
 282 static void split_on_domain_size(evalue *EP, evalue **EP_less, evalue **EP_more,
 283                                  int size, barvinok_options *options)
 284 {
 285     assert(value_zero_p(EP->d));
 286     assert(EP->x.p->type == partition);
 287     assert(EP->x.p->size >= 2);
 288
 289     struct evalue_section *s_less = new evalue_section[EP->x.p->size/2];
 290     struct evalue_section *s_more = new evalue_section[EP->x.p->size/2];
 291
 292     int n_less = 0;
 293     int n_more = 0;
 294
 295     Value c;
 296     value_init(c);
 297
 298     for (int i = 0; i < EP->x.p->size/2; ++i) {
 299         Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 300         Polyhedron *D_less = NULL;
 301         Polyhedron *D_more = NULL;
 302         Polyhedron **next_less = &D_less;
 303         Polyhedron **next_more = &D_more;
 304
 305         for (Polyhedron *P = D; P; P = P->next) {
 306             Polyhedron *next = P->next;
 307             P->next = NULL;
 308             barvinok_count_with_options(P, &c, options);
 309             P->next = next;
 310
 311             if (value_zero_p(c))
 312                 continue;
 313
 314             if (value_pos_p(c) && value_cmp_si(c, size) <= 0) {
 315                 *next_less = Polyhedron_Copy(P);
 316                 next_less = &(*next_less)->next;
 317             } else {
 318                 *next_more = Polyhedron_Copy(P);
 319                 next_more = &(*next_more)->next;
 320             }
 321         }
 322
 323         if (D_less) {
 324             s_less[n_less].D = D_less;
 325             s_less[n_less].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 326             n_less++;
 327         } else {
 328             s_more[n_more].D = D_more;
 329             s_more[n_more].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 330             n_more++;
 331         }
 332     }
 333
 334     value_clear(c);
 335
 336     *EP_less = evalue_from_section_array(s_less, n_less);
 337     *EP_more = evalue_from_section_array(s_more, n_more);
 338
 339     delete [] s_less;
 340     delete [] s_more;
 341 }
 342
 343 static piecewise_lst *optimize(evalue *EP, unsigned nvar, Polyhedron *C,
 344                                 const exvector& params,
 345                                 struct options *options)
 346 {
 347     piecewise_lst *pl = NULL;
 348     if (options->iterate > 0) {
 349         evalue *EP_less = NULL;
 350         evalue *EP_more = NULL;
 351         piecewise_lst *pl_more = NULL;
 352         split_on_domain_size(EP, &EP_less, &EP_more, options->iterate,
 353                                 options->verify.barvinok);
 354         if (!EVALUE_IS_ZERO(*EP_less)) {
 355             options->iterate = -1;
 356             pl = optimize(EP_less, nvar, C, params, options);
 357         }
 358         if (!EVALUE_IS_ZERO(*EP_more)) {
 359             options->iterate = 0;
 360             pl_more = optimize(EP_more, nvar, C, params, options);
 361         }
 362         evalue_free(EP_less);
 363         evalue_free(EP_more);
 364         if (!pl)
 365             return pl_more;
 366         if (!pl_more)
 367             return pl;
 368         pl->combine(*pl_more);
 369         delete pl_more;
 370         return pl;
 371     }
 372     if (options->iterate)
 373         pl = iterate(EP, nvar, C, options);
 374     else if (options->verify.barvinok->bound == BV_BOUND_BERNSTEIN) {
 375         pl = evalue_bernstein_coefficients(NULL, EP, C, params,
 376                                            options->verify.barvinok);
 377         if (options->lower)
 378             pl->minimize();
 379         else
 380             pl->maximize();
 381     } else
 382         pl = evalue_range_propagation(NULL, EP, params,
 383                                       options->verify.barvinok);
 384     return pl;
 385 }
 386
 387 static int optimize(evalue *EP, const char **all_vars,
 388                     unsigned nvar, unsigned nparam, struct options *options)
 389 {
 390     Polyhedron *U;
 391     piecewise_lst *pl = NULL;
 392     U = Universe_Polyhedron(nparam);
 393     int print_solution = 1;
 394     int result = 0;
 395
 396     exvector params;
 397     params = constructParameterVector(all_vars+nvar, nparam);
 398
 399     if (options->verify.verify) {
 400         verify_options_set_range(&options->verify, nvar+nparam);
 401         if (!options->verify.barvinok->verbose)
 402             print_solution = 0;
 403     }
 404
 405     if (options->lower)
 406         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
 407     else
 408         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
 409     pl = optimize(EP, nvar, U, params, options);
 410     assert(pl);
 411     if (print_solution)
 412         cout << *pl << endl;
 413     if (options->verify.verify) {
 414         result = verify(pl, EP, nvar, params, &options->verify);
 415     }
 416     delete pl;
 417
 418     Polyhedron_Free(U);
 419
 420     return result;
 421 }
 422
 423 int main(int argc, char **argv)
 424 {
 425     evalue *EP;
 426     const char **all_vars = NULL;
 427     unsigned nvar;
 428     unsigned nparam;
 429     struct options options;
 430     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 431     static struct argp_child argp_children[] = {
 432         { &convert_argp,        0,      "input conversion",     1 },
 433         { &verify_argp,         0,      "verification",         2 },
 434         { &barvinok_argp,       0,      "barvinok options",     3 },
 435         { 0 }
 436     };
 437     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 438     int result = 0;
 439
 440     options.verify.barvinok = bv_options;
 441     set_program_name(argv[0]);
 442     argp_parse(&argp, argc, argv, 0, 0, &options);
 443
 444     EP = evalue_read_from_file(stdin, options.var_list, &all_vars,
 445                                &nvar, &nparam, bv_options->MaxRays);
 446     assert(EP);
 447
 448     if (options.split)
 449         evalue_split_periods(EP, options.split, bv_options->MaxRays);
 450
 451     evalue_convert(EP, &options.convert, bv_options->verbose,
 452                         nvar+nparam, all_vars);
 453
 454     if (EVALUE_IS_ZERO(*EP))
 455         print_evalue(stdout, EP, all_vars);
 456     else
 457         result = optimize(EP, all_vars, nvar, nparam, &options);
 458
 459     evalue_free(EP);
 460
 461     if (options.var_list)
 462         free(options.var_list);
 463     Free_ParamNames(all_vars, nvar+nparam);
 464     barvinok_options_free(bv_options);
 465     return result;
 466 }