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