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 static int check_poly_max(const struct check_poly_data *data,
  85                           int nparam, Value *z,
  86                           const struct verify_options *options);
  87
  88 struct check_poly_max_data : public check_EP_data {
  89     piecewise_lst        *pl;
  90
  91     check_poly_max_data(evalue *EP, piecewise_lst *pl) : pl(pl) {
  92         this->EP = EP;
  93         this->cp.check = check_poly_max;
  94     }
  95 };
  96
  97 static int check_poly_max(const struct check_poly_data *data,
  98                           int nparam, Value *z,
  99                           const struct verify_options *options)
 100 {
 101     int k;
 102     int ok;
 103     const check_poly_max_data *max_data;
 104     max_data = (const check_poly_max_data *)data;
 105     const char *minmax;
 106     Value m, n, d;
 107     value_init(m);
 108     value_init(n);
 109     value_init(d);
 110     int sign;
 111
 112     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
 113         minmax = "max";
 114         sign = 1;
 115     } else {
 116         minmax = "min";
 117         sign = -1;
 118     }
 119
 120     max_data->pl->evaluate(nparam, z, &n, &d);
 121     if (sign > 0)
 122         mpz_fdiv_q(m, n, d);
 123     else
 124         mpz_cdiv_q(m, n, d);
 125
 126     if (options->print_all) {
 127         printf("%s(", minmax);
 128         value_print(stdout, VALUE_FMT, z[0]);
 129         for (k = 1; k < nparam; ++k) {
 130             printf(", ");
 131             value_print(stdout, VALUE_FMT, z[k]);
 132         }
 133         printf(") = ");
 134         value_print(stdout, VALUE_FMT, n);
 135         if (value_notone_p(d)) {
 136             printf("/");
 137             value_print(stdout, VALUE_FMT, d);
 138         }
 139         printf(" (");
 140         value_print(stdout, VALUE_FMT, m);
 141         printf(")");
 142     }
 143
 144     evalue_optimum(max_data, &n, sign);
 145
 146     if (sign > 0)
 147         ok = value_ge(m, n);
 148     else
 149         ok = value_le(m, n);
 150
 151     if (options->print_all) {
 152         printf(", %s(EP) = ", minmax);
 153         value_print(stdout, VALUE_FMT, n);
 154         printf(". ");
 155     }
 156
 157     if (!ok) {
 158         printf("\n");
 159         fflush(stdout);
 160         fprintf(stderr, "Error !\n");
 161         fprintf(stderr, "%s(", minmax);
 162         value_print(stderr, VALUE_FMT, z[0]);
 163         for (k = 1; k < nparam; ++k) {
 164             fprintf(stderr,", ");
 165             value_print(stderr, VALUE_FMT, z[k]);
 166         }
 167         fprintf(stderr, ") should be ");
 168         if (sign > 0)
 169             fprintf(stderr, "greater");
 170         else
 171             fprintf(stderr, "smaller");
 172         fprintf(stderr, " than or equal to ");
 173         value_print(stderr, VALUE_FMT, n);
 174         fprintf(stderr, ", while pl eval gives ");
 175         value_print(stderr, VALUE_FMT, m);
 176         fprintf(stderr, ".\n");
 177         cerr << *max_data->pl << endl;
 178     } else if (options->print_all)
 179         printf("OK.\n");
 180
 181     value_clear(m);
 182     value_clear(n);
 183     value_clear(d);
 184
 185     return ok;
 186 }
 187
 188 static int verify(piecewise_lst *pl, evalue *EP, unsigned nvar, unsigned nparam,
 189                   struct verify_options *options)
 190 {
 191     check_poly_max_data data(EP, pl);
 192     return !check_EP(&data, nvar, nparam, options);
 193 }
 194
 195 static piecewise_lst *iterate(evalue *EP, unsigned nvar, Polyhedron *C,
 196                                 struct options *options)
 197 {
 198     assert(C->Dimension == 0);
 199     Vector *p = Vector_Alloc(nvar+2);
 200     value_set_si(p->p[nvar+1], 1);
 201     Value opt;
 202
 203     value_init(opt);
 204     struct check_EP_data data;
 205     data.cp.z = p->p;
 206     data.cp.check = NULL;
 207
 208     data.EP = EP;
 209     check_EP_set_scan(&data, C, options->verify.barvinok->MaxRays);
 210     int sign = options->verify.barvinok->bernstein_optimize;
 211     evalue_optimum(&data, &opt, sign);
 212     check_EP_clear_scan(&data);
 213
 214     exvector params;
 215     piecewise_lst *pl = new piecewise_lst(params, sign);
 216     pl->add_guarded_lst(Polyhedron_Copy(C), lst(value2numeric(opt)));
 217
 218     value_clear(opt);
 219     Vector_Free(p);
 220
 221     return pl;
 222 }
 223
 224 /*
 225  * Split (partition) EP into a partition with (sub)domains containing
 226  * size integer points or less and a partition with (sub)domains
 227  * containing more integer points.
 228  */
 229 static void split_on_domain_size(evalue *EP, evalue **EP_less, evalue **EP_more,
 230                                  int size, barvinok_options *options)
 231 {
 232     assert(value_zero_p(EP->d));
 233     assert(EP->x.p->type == partition);
 234     assert(EP->x.p->size >= 2);
 235
 236     struct evalue_section *s_less = new evalue_section[EP->x.p->size/2];
 237     struct evalue_section *s_more = new evalue_section[EP->x.p->size/2];
 238
 239     int n_less = 0;
 240     int n_more = 0;
 241
 242     Value c;
 243     value_init(c);
 244
 245     for (int i = 0; i < EP->x.p->size/2; ++i) {
 246         Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 247         Polyhedron *D_less = NULL;
 248         Polyhedron *D_more = NULL;
 249         Polyhedron **next_less = &D_less;
 250         Polyhedron **next_more = &D_more;
 251
 252         for (Polyhedron *P = D; P; P = P->next) {
 253             Polyhedron *next = P->next;
 254             P->next = NULL;
 255             barvinok_count_with_options(P, &c, options);
 256             P->next = next;
 257
 258             if (value_zero_p(c))
 259                 continue;
 260
 261             if (value_pos_p(c) && value_cmp_si(c, size) <= 0) {
 262                 *next_less = Polyhedron_Copy(P);
 263                 next_less = &(*next_less)->next;
 264             } else {
 265                 *next_more = Polyhedron_Copy(P);
 266                 next_more = &(*next_more)->next;
 267             }
 268         }
 269
 270         if (D_less) {
 271             s_less[n_less].D = D_less;
 272             s_less[n_less].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 273             n_less++;
 274         } else {
 275             s_more[n_more].D = D_more;
 276             s_more[n_more].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 277             n_more++;
 278         }
 279     }
 280
 281     value_clear(c);
 282
 283     *EP_less = evalue_from_section_array(s_less, n_less);
 284     *EP_more = evalue_from_section_array(s_more, n_more);
 285
 286     delete [] s_less;
 287     delete [] s_more;
 288 }
 289
 290 static piecewise_lst *optimize(evalue *EP, unsigned nvar, Polyhedron *C,
 291                                 const exvector& params,
 292                                 struct options *options)
 293 {
 294     piecewise_lst *pl = NULL;
 295     if (options->iterate > 0) {
 296         evalue *EP_less = NULL;
 297         evalue *EP_more = NULL;
 298         piecewise_lst *pl_more = NULL;
 299         split_on_domain_size(EP, &EP_less, &EP_more, options->iterate,
 300                                 options->verify.barvinok);
 301         if (!EVALUE_IS_ZERO(*EP_less)) {
 302             options->iterate = -1;
 303             pl = optimize(EP_less, nvar, C, params, options);
 304         }
 305         if (!EVALUE_IS_ZERO(*EP_more)) {
 306             options->iterate = 0;
 307             pl_more = optimize(EP_more, nvar, C, params, options);
 308         }
 309         evalue_free(EP_less);
 310         evalue_free(EP_more);
 311         if (!pl)
 312             return pl_more;
 313         if (!pl_more)
 314             return pl;
 315         pl->combine(*pl_more);
 316         delete pl_more;
 317         return pl;
 318     }
 319     if (options->iterate)
 320         pl = iterate(EP, nvar, C, options);
 321     else if (options->verify.barvinok->bound == BV_BOUND_BERNSTEIN) {
 322         pl = evalue_bernstein_coefficients(NULL, EP, C, params,
 323                                            options->verify.barvinok);
 324         if (options->lower)
 325             pl->minimize();
 326         else
 327             pl->maximize();
 328     } else
 329         pl = evalue_range_propagation(NULL, EP, params,
 330                                       options->verify.barvinok);
 331     return pl;
 332 }
 333
 334 static int optimize(evalue *EP, const char **all_vars,
 335                     unsigned nvar, unsigned nparam, struct options *options)
 336 {
 337     Polyhedron *U;
 338     piecewise_lst *pl = NULL;
 339     U = Universe_Polyhedron(nparam);
 340     int print_solution = 1;
 341     int result = 0;
 342
 343     exvector params;
 344     params = constructParameterVector(all_vars+nvar, nparam);
 345
 346     if (options->verify.verify) {
 347         verify_options_set_range(&options->verify, nvar+nparam);
 348         if (!options->verify.barvinok->verbose)
 349             print_solution = 0;
 350     }
 351
 352     if (options->lower)
 353         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
 354     else
 355         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
 356     pl = optimize(EP, nvar, U, params, options);
 357     assert(pl);
 358     if (print_solution)
 359         cout << *pl << endl;
 360     if (options->verify.verify) {
 361         result = verify(pl, EP, nvar, nparam, &options->verify);
 362     }
 363     delete pl;
 364
 365     Polyhedron_Free(U);
 366
 367     return result;
 368 }
 369
 370 int main(int argc, char **argv)
 371 {
 372     evalue *EP;
 373     const char **all_vars = NULL;
 374     unsigned nvar;
 375     unsigned nparam;
 376     struct options options;
 377     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 378     static struct argp_child argp_children[] = {
 379         { &convert_argp,        0,      "input conversion",     1 },
 380         { &verify_argp,         0,      "verification",         2 },
 381         { &barvinok_argp,       0,      "barvinok options",     3 },
 382         { 0 }
 383     };
 384     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 385     int result = 0;
 386
 387     options.verify.barvinok = bv_options;
 388     set_program_name(argv[0]);
 389     argp_parse(&argp, argc, argv, 0, 0, &options);
 390
 391     EP = evalue_read_from_file(stdin, options.var_list, &all_vars,
 392                                &nvar, &nparam, bv_options->MaxRays);
 393     assert(EP);
 394
 395     if (options.split)
 396         evalue_split_periods(EP, options.split, bv_options->MaxRays);
 397
 398     evalue_convert(EP, &options.convert, bv_options->verbose,
 399                         nvar+nparam, all_vars);
 400
 401     if (EVALUE_IS_ZERO(*EP))
 402         print_evalue(stdout, EP, all_vars);
 403     else
 404         result = optimize(EP, all_vars, nvar, nparam, &options);
 405
 406     evalue_free(EP);
 407
 408     if (options.var_list)
 409         free(options.var_list);
 410     Free_ParamNames(all_vars, nvar+nparam);
 411     barvinok_options_free(bv_options);
 412     return result;
 413 }