maximize.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
  14 using std::cout;
  15 using std::cerr;
  16 using std::endl;
  17 using namespace GiNaC;
  18 using namespace bernstein;
  19 using namespace barvinok;
  20
  21 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  22
  23 #define OPT_VARS            (BV_OPT_LAST+1)
  24 #define OPT_SPLIT           (BV_OPT_LAST+2)
  25 #define OPT_MIN             (BV_OPT_LAST+3)
  26
  27 struct argp_option argp_options[] = {
  28     { "split",              OPT_SPLIT,  "int" },
  29     { "variables",          OPT_VARS,   "list", 0,
  30         "comma separated list of variables over which to maximize" },
  31     { "verbose",            'v',        0,      0, },
  32     { "minimize",           OPT_MIN,    0, 0,   "minimize instead of maximize"},
  33     { 0 }
  34 };
  35
  36 struct options {
  37     struct convert_options   convert;
  38     struct verify_options    verify;
  39     char* var_list;
  40     int verbose;
  41     int split;
  42     int minimize;
  43 };
  44
  45 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  46 {
  47     struct options *options = (struct options*) state->input;
  48
  49     switch (key) {
  50     case ARGP_KEY_INIT:
  51         state->child_inputs[0] = &options->convert;
  52         state->child_inputs[1] = &options->verify;
  53         state->child_inputs[2] = options->verify.barvinok;
  54         options->var_list = NULL;
  55         options->verbose = 0;
  56         options->split = 0;
  57         options->minimize = 0;
  58         break;
  59     case 'v':
  60         options->verbose = 1;
  61         break;
  62     case OPT_VARS:
  63         options->var_list = strdup(arg);
  64         break;
  65     case OPT_SPLIT:
  66         options->split = atoi(arg);
  67         break;
  68     case OPT_MIN:
  69         options->minimize = 1;
  70         break;
  71     default:
  72         return ARGP_ERR_UNKNOWN;
  73     }
  74     return 0;
  75 }
  76
  77 static int check_poly_max(const struct check_poly_data *data,
  78                           int nparam, Value *z,
  79                           const struct verify_options *options);
  80
  81 struct check_poly_max_data : public check_poly_data {
  82     Polyhedron          **S;
  83     evalue               *EP;
  84     piecewise_lst        *pl;
  85
  86     check_poly_max_data(Value *z, evalue *EP, piecewise_lst *pl) :
  87                         EP(EP), pl(pl) {
  88         this->z = z;
  89         this->check = check_poly_max;
  90     }
  91 };
  92
  93 static void optimum(Polyhedron *S, int pos, const check_poly_max_data *data,
  94                     Value *opt, bool& found,
  95                     const struct verify_options *options)
  96 {
  97     if (!S) {
  98         Value c;
  99         value_init(c);
 100         value_set_double(c, compute_evalue(data->EP, data->z+1)+.25);
 101         if (!found) {
 102             value_assign(*opt, c);
 103             found = true;
 104         } else {
 105             if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
 106                 if (value_gt(c, *opt))
 107                     value_assign(*opt, c);
 108             } else {
 109                 if (value_lt(c, *opt))
 110                     value_assign(*opt, c);
 111             }
 112         }
 113         value_clear(c);
 114     } else {
 115         Value LB, UB;
 116         int ok;
 117         value_init(LB);
 118         value_init(UB);
 119         ok = !(lower_upper_bounds(1+pos, S, data->z, &LB, &UB));
 120         assert(ok);
 121         for (; value_le(LB, UB); value_increment(LB, LB)) {
 122             value_assign(data->z[1+pos], LB);
 123             optimum(S->next, pos+1, data, opt, found, options);
 124         }
 125         value_set_si(data->z[1+pos], 0);
 126         value_clear(LB);
 127         value_clear(UB);
 128     }
 129 }
 130
 131 static void optimum(const check_poly_max_data *data, Value *opt,
 132                     const struct verify_options *options)
 133 {
 134     bool found = false;
 135     for (int i = 0; i < data->EP->x.p->size/2; ++i)
 136         if (!emptyQ2(data->S[i]))
 137             optimum(data->S[i], 0, data, opt, found, options);
 138     assert(found);
 139 }
 140
 141 static int check_poly_max(const struct check_poly_data *data,
 142                           int nparam, Value *z,
 143                           const struct verify_options *options)
 144 {
 145     int k;
 146     int ok;
 147     const check_poly_max_data *max_data;
 148     max_data = static_cast<const check_poly_max_data *>(data);
 149     const char *minmax;
 150     Value m, n, d;
 151     value_init(m);
 152     value_init(n);
 153     value_init(d);
 154
 155     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
 156         minmax = "max";
 157     else
 158         minmax = "min";
 159
 160     max_data->pl->evaluate(nparam, z, &n, &d);
 161     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
 162         mpz_fdiv_q(m, n, d);
 163     else
 164         mpz_cdiv_q(m, n, d);
 165
 166     if (options->print_all) {
 167         printf("%s(", minmax);
 168         value_print(stdout, VALUE_FMT, z[0]);
 169         for (k = 1; k < nparam; ++k) {
 170             printf(", ");
 171             value_print(stdout, VALUE_FMT, z[k]);
 172         }
 173         printf(") = ");
 174         value_print(stdout, VALUE_FMT, n);
 175         if (value_notone_p(d)) {
 176             printf("/");
 177             value_print(stdout, VALUE_FMT, d);
 178         }
 179         printf(" (");
 180         value_print(stdout, VALUE_FMT, m);
 181         printf(")");
 182     }
 183
 184     optimum(max_data, &n, options);
 185
 186     if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
 187         ok = value_ge(m, n);
 188     else
 189         ok = value_le(m, n);
 190
 191     if (options->print_all) {
 192         printf(", %s(EP) = ", minmax);
 193         value_print(stdout, VALUE_FMT, n);
 194         printf(". ");
 195     }
 196
 197     if (!ok) {
 198         printf("\n");
 199         fflush(stdout);
 200         fprintf(stderr, "Error !\n");
 201         fprintf(stderr, "%s(", minmax);
 202         value_print(stderr, VALUE_FMT, z[0]);
 203         for (k = 1; k < nparam; ++k) {
 204             fprintf(stderr,", ");
 205             value_print(stderr, VALUE_FMT, z[k]);
 206         }
 207         fprintf(stderr, ") should be ");
 208         if (options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX)
 209             fprintf(stderr, "greater");
 210         else
 211             fprintf(stderr, "smaller");
 212         fprintf(stderr, " than or equal to ");
 213         value_print(stderr, VALUE_FMT, n);
 214         fprintf(stderr, ", while pl eval gives ");
 215         value_print(stderr, VALUE_FMT, m);
 216         fprintf(stderr, ".\n");
 217         cerr << *max_data->pl << endl;
 218     } else if (options->print_all)
 219         printf("OK.\n");
 220
 221     value_clear(m);
 222     value_clear(n);
 223     value_clear(d);
 224
 225     return ok;
 226 }
 227
 228 static int verify(Polyhedron *D, piecewise_lst *pl, evalue *EP,
 229                   unsigned nvar, unsigned nparam, Vector *p,
 230                   struct verify_options *options)
 231 {
 232     Polyhedron *CS, *S;
 233     unsigned MaxRays = options->barvinok->MaxRays;
 234     assert(value_zero_p(EP->d));
 235     assert(EP->x.p->type == partition);
 236     int ok = 1;
 237
 238     CS = check_poly_context_scan(NULL, &D, D->Dimension, options);
 239
 240     check_poly_init(D, options);
 241
 242     if (!(CS && emptyQ2(CS))) {
 243         check_poly_max_data data(p->p, EP, pl);
 244         data.S = ALLOCN(Polyhedron *, EP->x.p->size/2);
 245         for (int i = 0; i < EP->x.p->size/2; ++i) {
 246             Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 247             data.S[i] = Polyhedron_Scan(A, D, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 248         }
 249         ok = check_poly(CS, &data, nparam, 0, p->p+1+nvar, options);
 250         for (int i = 0; i < EP->x.p->size/2; ++i)
 251             Domain_Free(data.S[i]);
 252         free(data.S);
 253     }
 254
 255     if (!options->print_all)
 256         printf("\n");
 257
 258     if (CS) {
 259         Domain_Free(CS);
 260         Domain_Free(D);
 261     }
 262
 263     return ok;
 264 }
 265
 266 static int verify(piecewise_lst *pl, evalue *EP, unsigned nvar, unsigned nparam,
 267                   struct verify_options *options)
 268 {
 269     Vector *p;
 270
 271     p = Vector_Alloc(nvar+nparam+2);
 272     value_set_si(p->p[nvar+nparam+1], 1);
 273
 274     for (int i = 0; i < pl->list.size(); ++i) {
 275         int ok = verify(pl->list[i].first, pl, EP, nvar, nparam, p, options);
 276         if (!ok && !options->continue_on_error)
 277             break;
 278     }
 279
 280     Vector_Free(p);
 281
 282     return 0;
 283 }
 284
 285 static int optimize(evalue *EP, char **all_vars, unsigned nvar, unsigned nparam,
 286                     struct options *options)
 287 {
 288     Polyhedron *U;
 289     piecewise_lst *pl = NULL;
 290     U = Universe_Polyhedron(nparam);
 291     int print_solution = 1;
 292     int result = 0;
 293
 294     exvector params;
 295     params = constructParameterVector(all_vars+nvar, nparam);
 296
 297     if (options->verify.verify) {
 298         verify_options_set_range(&options->verify, nvar+nparam);
 299         if (!options->verbose)
 300             print_solution = 0;
 301     }
 302
 303     if (options->minimize)
 304         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
 305     else
 306         options->verify.barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
 307     pl = evalue_bernstein_coefficients(NULL, EP, U, params,
 308                                        options->verify.barvinok);
 309     assert(pl);
 310     if (options->minimize)
 311         pl->minimize();
 312     else
 313         pl->maximize();
 314     if (print_solution)
 315         cout << *pl << endl;
 316     if (options->verify.verify)
 317         result = verify(pl, EP, nvar, nparam, &options->verify);
 318     delete pl;
 319
 320     Polyhedron_Free(U);
 321
 322     return result;
 323 }
 324
 325 int main(int argc, char **argv)
 326 {
 327     evalue *EP;
 328     char **all_vars = NULL;
 329     unsigned nvar;
 330     unsigned nparam;
 331     struct options options;
 332     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 333     static struct argp_child argp_children[] = {
 334         { &convert_argp,        0,      "input conversion",     1 },
 335         { &verify_argp,         0,      "verification",         2 },
 336         { &barvinok_argp,       0,      "barvinok options",     3 },
 337         { 0 }
 338     };
 339     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 340     int result = 0;
 341
 342     options.verify.barvinok = bv_options;
 343     set_program_name(argv[0]);
 344     argp_parse(&argp, argc, argv, 0, 0, &options);
 345
 346     EP = evalue_read_from_file(stdin, options.var_list, &all_vars,
 347                                &nvar, &nparam, bv_options->MaxRays);
 348     assert(EP);
 349
 350     if (options.split)
 351         evalue_split_periods(EP, options.split, bv_options->MaxRays);
 352
 353     evalue_convert(EP, &options.convert, options.verbose, nparam, all_vars);
 354
 355     if (EVALUE_IS_ZERO(*EP))
 356         print_evalue(stdout, EP, all_vars);
 357     else
 358         result = optimize(EP, all_vars, nvar, nparam, &options);
 359
 360     free_evalue_refs(EP);
 361     free(EP);
 362
 363     if (options.var_list)
 364         free(options.var_list);
 365     Free_ParamNames(all_vars, nvar+nparam);
 366     barvinok_options_free(bv_options);
 367     return result;
 368 }