maximize.cc

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