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/bernstein.h>
   8 #include "argp.h"
   9 #include "progname.h"
  10 #include "evalue_convert.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 ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  23
  24 #define OPT_VARS            (BV_OPT_LAST+1)
  25 #define OPT_SPLIT           (BV_OPT_LAST+2)
  26 #define OPT_LOWER           (BV_OPT_LAST+3)
  27 #define OPT_METHOD          (BV_OPT_LAST+4)
  28
  29 struct argp_option argp_options[] = {
  30     { "split",              OPT_SPLIT,  "int" },
  31     { "variables",          OPT_VARS,   "list", 0,
  32         "comma separated list of variables over which to compute a bound" },
  33     { "lower",              OPT_LOWER,  0, 0,   "compute lower bound instead of upper bound"},
  34     { "optimization-method",    OPT_METHOD,     "bernstein|propagation",
  35         0, "optimization method to use" },
  36     { 0 }
  37 };
  38
  39 struct options {
  40     struct convert_options   convert;
  41     struct verify_options    verify;
  42     char* var_list;
  43     int split;
  44     int lower;
  45 #define METHOD_BERNSTEIN        0
  46 #define METHOD_PROPAGATION      1
  47     int method;
  48 };
  49
  50 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  51 {
  52     struct options *options = (struct options*) state->input;
  53
  54     switch (key) {
  55     case ARGP_KEY_INIT:
  56         state->child_inputs[0] = &options->convert;
  57         state->child_inputs[1] = &options->verify;
  58         state->child_inputs[2] = options->verify.barvinok;
  59         options->var_list = NULL;
  60         options->split = 0;
  61         options->lower = 0;
  62         options->method = METHOD_BERNSTEIN;
  63         break;
  64     case OPT_VARS:
  65         options->var_list = strdup(arg);
  66         break;
  67     case OPT_SPLIT:
  68         options->split = atoi(arg);
  69         break;
  70     case OPT_LOWER:
  71         options->lower = 1;
  72         break;
  73     case OPT_METHOD:
  74         if (!strcmp(arg, "bernstein"))
  75             options->method = METHOD_BERNSTEIN;
  76         else if (!strcmp(arg, "propagation"))
  77             options->method = METHOD_PROPAGATION;
  78         else
  79             argp_error(state, "unknown value for --optimization-method option");
  80         break;
  81     default:
  82         return ARGP_ERR_UNKNOWN;
  83     }
  84     return 0;
  85 }
  86
  87 static int check_poly_max(const struct check_poly_data *data,
  88                           int nparam, Value *z,
  89                           const struct verify_options *options);
  90
  91 struct check_poly_max_data : public check_EP_data {
  92     piecewise_lst        *pl;
  93
  94     check_poly_max_data(evalue *EP, piecewise_lst *pl) : pl(pl) {
  95         this->EP = EP;
  96         this->cp.check = check_poly_max;
  97     }
  98 };
  99
 100 static int check_poly_max(const struct check_poly_data *data,
 101                           int nparam, Value *z,
 102                           const struct verify_options *options)
 103 {
 104     int k;
 105     int ok;
 106     const check_poly_max_data *max_data;
 107     max_data = (const check_poly_max_data *)data;
 108     const char *minmax;
 109     Value m, n, d;
 110     value_init(m);
 111     value_init(n);
 112     value_init(d);
 113     int sign;
 114
 115     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
 116         minmax = "max";
 117         sign = 1;
 118     } else {
 119         minmax = "min";
 120         sign = -1;
 121     }
 122
 123     max_data->pl->evaluate(nparam, z, &n, &d);
 124     if (sign > 0)
 125         mpz_fdiv_q(m, n, d);
 126     else
 127         mpz_cdiv_q(m, n, d);
 128
 129     if (options->print_all) {
 130         printf("%s(", minmax);
 131         value_print(stdout, VALUE_FMT, z[0]);
 132         for (k = 1; k < nparam; ++k) {
 133             printf(", ");
 134             value_print(stdout, VALUE_FMT, z[k]);
 135         }
 136         printf(") = ");
 137         value_print(stdout, VALUE_FMT, n);
 138         if (value_notone_p(d)) {
 139             printf("/");
 140             value_print(stdout, VALUE_FMT, d);
 141         }
 142         printf(" (");
 143         value_print(stdout, VALUE_FMT, m);
 144         printf(")");
 145     }
 146
 147     evalue_optimum(max_data, &n, sign);
 148
 149     if (sign > 0)
 150         ok = value_ge(m, n);
 151     else
 152         ok = value_le(m, n);
 153
 154     if (options->print_all) {
 155         printf(", %s(EP) = ", minmax);
 156         value_print(stdout, VALUE_FMT, n);
 157         printf(". ");
 158     }
 159
 160     if (!ok) {
 161         printf("\n");
 162         fflush(stdout);
 163         fprintf(stderr, "Error !\n");
 164         fprintf(stderr, "%s(", minmax);
 165         value_print(stderr, VALUE_FMT, z[0]);
 166         for (k = 1; k < nparam; ++k) {
 167             fprintf(stderr,", ");
 168             value_print(stderr, VALUE_FMT, z[k]);
 169         }
 170         fprintf(stderr, ") should be ");
 171         if (sign > 0)
 172             fprintf(stderr, "greater");
 173         else
 174             fprintf(stderr, "smaller");
 175         fprintf(stderr, " than or equal to ");
 176         value_print(stderr, VALUE_FMT, n);
 177         fprintf(stderr, ", while pl eval gives ");
 178         value_print(stderr, VALUE_FMT, m);
 179         fprintf(stderr, ".\n");
 180         cerr << *max_data->pl << endl;
 181     } else if (options->print_all)
 182         printf("OK.\n");
 183
 184     value_clear(m);
 185     value_clear(n);
 186     value_clear(d);
 187
 188     return ok;
 189 }
 190
 191 static int verify(piecewise_lst *pl, evalue *EP, unsigned nvar, unsigned nparam,
 192                   struct verify_options *options)
 193 {
 194     check_poly_max_data data(EP, pl);
 195     return !check_EP(&data, nvar, nparam, options);
 196 }
 197
 198 static int optimize(evalue *EP, char **all_vars, unsigned nvar, unsigned nparam,
 199                     struct options *options)
 200 {
 201     Polyhedron *U;
 202     piecewise_lst *pl = NULL;
 203     U = Universe_Polyhedron(nparam);
 204     int print_solution = 1;
 205     int result = 0;
 206
 207     exvector params;
 208     params = constructParameterVector(all_vars+nvar, nparam);
 209
 210     if (options->verify.verify) {
 211         verify_options_set_range(&options->verify, nvar+nparam);
 212         if (!options->verify.barvinok->verbose)
 213             print_solution = 0;
 214     }
 215
 216     if (options->lower)
 217         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
 218     else
 219         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
 220     if (options->method == METHOD_BERNSTEIN) {
 221         pl = evalue_bernstein_coefficients(NULL, EP, U, params,
 222                                            options->verify.barvinok);
 223         if (options->lower)
 224             pl->minimize();
 225         else
 226             pl->maximize();
 227     } else
 228         pl = evalue_range_propagation(NULL, EP, params,
 229                                       options->verify.barvinok);
 230     assert(pl);
 231     if (print_solution)
 232         cout << *pl << endl;
 233     if (options->verify.verify) {
 234         result = verify(pl, EP, nvar, nparam, &options->verify);
 235     }
 236     delete pl;
 237
 238     Polyhedron_Free(U);
 239
 240     return result;
 241 }
 242
 243 int main(int argc, char **argv)
 244 {
 245     evalue *EP;
 246     char **all_vars = NULL;
 247     unsigned nvar;
 248     unsigned nparam;
 249     struct options options;
 250     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 251     static struct argp_child argp_children[] = {
 252         { &convert_argp,        0,      "input conversion",     1 },
 253         { &verify_argp,         0,      "verification",         2 },
 254         { &barvinok_argp,       0,      "barvinok options",     3 },
 255         { 0 }
 256     };
 257     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 258     int result = 0;
 259
 260     options.verify.barvinok = bv_options;
 261     set_program_name(argv[0]);
 262     argp_parse(&argp, argc, argv, 0, 0, &options);
 263
 264     EP = evalue_read_from_file(stdin, options.var_list, &all_vars,
 265                                &nvar, &nparam, bv_options->MaxRays);
 266     assert(EP);
 267
 268     if (options.split)
 269         evalue_split_periods(EP, options.split, bv_options->MaxRays);
 270
 271     evalue_convert(EP, &options.convert, bv_options->verbose, nparam, all_vars);
 272
 273     if (EVALUE_IS_ZERO(*EP))
 274         print_evalue(stdout, EP, all_vars);
 275     else
 276         result = optimize(EP, all_vars, nvar, nparam, &options);
 277
 278     evalue_free(EP);
 279
 280     if (options.var_list)
 281         free(options.var_list);
 282     Free_ParamNames(all_vars, nvar+nparam);
 283     barvinok_options_free(bv_options);
 284     return result;
 285 }