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