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
   9 using std::cout;
  10 using std::cerr;
  11 using std::endl;
  12 using namespace GiNaC;
  13 using namespace bernstein;
  14 using namespace barvinok;
  15
  16 #define OPT_VARS            (BV_OPT_LAST+1)
  17
  18 struct argp_option argp_options[] = {
  19     { "variables",          OPT_VARS,   "int",  0,
  20         "number of variables over which to maximize" },
  21     { "verbose",            'V',        0,      0, },
  22     { 0 }
  23 };
  24
  25 struct options {
  26     int nvar;
  27     int verbose;
  28 };
  29
  30 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  31 {
  32     struct options *options = (struct options*) state->input;
  33
  34     switch (key) {
  35     case ARGP_KEY_INIT:
  36         options->nvar = -1;
  37         options->verbose = 0;
  38         break;
  39     case 'V':
  40         options->verbose = 1;
  41         break;
  42     case OPT_VARS:
  43         options->nvar = atoi(arg);
  44         break;
  45     default:
  46         return ARGP_ERR_UNKNOWN;
  47     }
  48     return 0;
  49 }
  50
  51
  52 #define ALLOC(type) (type*)malloc(sizeof(type))
  53 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  54
  55 enum token_type { TOKEN_UNKNOWN = 256, TOKEN_VALUE, TOKEN_IDENT, TOKEN_GE };
  56
  57 struct token {
  58     enum token_type  type;
  59
  60     int line;
  61     int col;
  62
  63     union {
  64         Value   v;
  65         char    *s;
  66     } u;
  67 };
  68
  69 static struct token *token_new(int line, int col)
  70 {
  71     struct token *tok = ALLOC(struct token);
  72     tok->line = line;
  73     tok->col = col;
  74     return tok;
  75 }
  76
  77 void token_free(struct token *tok)
  78 {
  79     if (tok->type == TOKEN_VALUE)
  80         value_clear(tok->u.v);
  81     else if (tok->type == TOKEN_IDENT)
  82         free(tok->u.s);
  83     free(tok);
  84 }
  85
  86 struct stream {
  87     FILE            *file;
  88     int             line;
  89     int             col;
  90     int             eof;
  91
  92     char            *buffer;
  93     size_t          size;
  94     size_t          len;
  95     int             c;
  96
  97     struct token    *tokens[5];
  98     int             n_token;
  99 };
 100
 101 static struct stream* stream_new(FILE *file)
 102 {
 103     int i;
 104     struct stream *s = ALLOC(struct stream);
 105     s->file = file;
 106     s->size = 256;
 107     s->buffer = (char*)malloc(s->size);
 108     s->len = 0;
 109     s->line = 1;
 110     s->col = 0;
 111     s->eof = 0;
 112     s->c = -1;
 113     for (i = 0; i < 5; ++i)
 114         s->tokens[i] = NULL;
 115     s->n_token = 0;
 116     return s;
 117 }
 118
 119 static void stream_free(struct stream *s)
 120 {
 121     free(s->buffer);
 122     assert(s->n_token == 0);
 123     free(s);
 124 }
 125
 126 static int stream_getc(struct stream *s)
 127 {
 128     int c;
 129     if (s->eof)
 130         return -1;
 131     c = fgetc(s->file);
 132     if (c == -1)
 133         s->eof = 1;
 134     if (s->c != -1) {
 135         if (s->c == '\n') {
 136             s->line++;
 137             s->col = 0;
 138         } else
 139             s->col++;
 140     }
 141     s->c = c;
 142     return c;
 143 }
 144
 145 static void stream_ungetc(struct stream *s, int c)
 146 {
 147     ungetc(c, s->file);
 148     s->c = -1;
 149 }
 150
 151 static void stream_push_char(struct stream *s, int c)
 152 {
 153     if (s->len >= s->size) {
 154         s->size = (3*s->size)/2;
 155         s->buffer = (char*)realloc(s->buffer, s->size);
 156     }
 157     s->buffer[s->len++] = c;
 158 }
 159
 160 static struct token *stream_push_token(struct stream *s, struct token *tok)
 161 {
 162     assert(s->n_token < 5);
 163     s->tokens[s->n_token++] = tok;
 164 }
 165
 166 static struct token *stream_next_token(struct stream *s)
 167 {
 168     int c;
 169     struct token *tok;
 170     int line, col;
 171
 172     if (s->n_token)
 173         return s->tokens[--s->n_token];
 174
 175     s->len = 0;
 176
 177     /* skip spaces */
 178     while ((c = stream_getc(s)) != -1 && isspace(c))
 179         /* nothing */
 180         ;
 181
 182     line = s->line;
 183     col = s->col;
 184
 185     if (c == -1)
 186         return NULL;
 187     if (c == '(' ||
 188         c == ')' ||
 189         c == '+' ||
 190         c == '/' ||
 191         c == '*' ||
 192         c == '^' ||
 193         c == '=' ||
 194         c == '[' ||
 195         c == ']' ||
 196         c == '{' ||
 197         c == '}') {
 198         tok = token_new(line, col);
 199         tok->type = (token_type)c;
 200         return tok;
 201     }
 202     if (c == '-' || isdigit(c)) {
 203         tok = token_new(line, col);
 204         tok->type = TOKEN_VALUE;
 205         value_init(tok->u.v);
 206         stream_push_char(s, c);
 207         while ((c = stream_getc(s)) != -1 && isdigit(c))
 208             stream_push_char(s, c);
 209         if (c != -1)
 210             stream_ungetc(s, c);
 211         if (s->len == 1 && s->buffer[0] == '-')
 212             value_set_si(tok->u.v, -1);
 213         else {
 214             stream_push_char(s, '\0');
 215             mpz_set_str(tok->u.v, s->buffer, 0);
 216         }
 217         return tok;
 218     }
 219     if (isalpha(c)) {
 220         tok = token_new(line, col);
 221         tok->type = TOKEN_IDENT;
 222         stream_push_char(s, c);
 223         while ((c = stream_getc(s)) != -1 && isalpha(c))
 224             stream_push_char(s, c);
 225         if (c != -1)
 226             stream_ungetc(s, c);
 227         stream_push_char(s, '\0');
 228         tok->u.s = strdup(s->buffer);
 229         return tok;
 230     }
 231     if (c == '>') {
 232         if ((c = stream_getc(s)) == '=') {
 233             tok = token_new(line, col);
 234             tok->type = TOKEN_GE;
 235             return tok;
 236         }
 237         if (c != -1)
 238             stream_ungetc(s, c);
 239     }
 240
 241     tok = token_new(line, col);
 242     tok->type = TOKEN_UNKNOWN;
 243     return tok;
 244 }
 245
 246 void stream_error(struct stream *s, struct token *tok, char *msg)
 247 {
 248     int line = tok ? tok->line : s->line;
 249     int col = tok ? tok->col : s->col;
 250     fprintf(stderr, "syntax error (%d, %d): %s\n", line, col, msg);
 251     if (tok) {
 252         if (tok->type < 256)
 253             fprintf(stderr, "got '%c'\n", tok->type);
 254     }
 255 }
 256
 257 struct parameter {
 258     char                *name;
 259     int                  pos;
 260     struct parameter    *next;
 261 };
 262
 263 struct parameter *parameter_new(char *name, int pos, struct parameter *next)
 264 {
 265     struct parameter *p = ALLOC(struct parameter);
 266     p->name = strdup(name);
 267     p->pos = pos;
 268     p->next = next;
 269     return p;
 270 }
 271
 272 static int parameter_pos(struct parameter **p, char *s)
 273 {
 274     int pos = *p ? (*p)->pos+1 : 0;
 275     struct parameter *q;
 276
 277     for (q = *p; q; q = q->next) {
 278         if (strcmp(q->name, s) == 0)
 279             break;
 280     }
 281     if (q)
 282         pos = q->pos;
 283     else
 284         *p = parameter_new(s, pos, *p);
 285     return pos;
 286 }
 287
 288 static int optional_power(struct stream *s)
 289 {
 290     int pow;
 291     struct token *tok;
 292     tok = stream_next_token(s);
 293     if (!tok)
 294         return 1;
 295     if (tok->type != '^') {
 296         stream_push_token(s, tok);
 297         return 1;
 298     }
 299     token_free(tok);
 300     tok = stream_next_token(s);
 301     if (!tok || tok->type != TOKEN_VALUE) {
 302         stream_error(s, tok, "expecting exponent");
 303         if (tok)
 304             stream_push_token(s, tok);
 305         return 1;
 306     }
 307     pow = VALUE_TO_INT(tok->u.v);
 308     token_free(tok);
 309     return pow;
 310 }
 311
 312 static evalue *evalue_read_factor(struct stream *s, struct parameter **p);
 313 static evalue *evalue_read_term(struct stream *s, struct parameter **p);
 314
 315 static evalue *read_fract_like(struct stream *s, struct token *tok,
 316                                struct parameter **p)
 317 {
 318     int pow;
 319     evalue *arg;
 320     evalue *e;
 321     char next;
 322     char *error;
 323     enode_type type;
 324
 325     if (tok->type == '[') {
 326         next = ']';
 327         error = "expecting \"]\"";
 328         type = flooring;
 329     } else if (tok->type == '{') {
 330         next = '}';
 331         error = "expecting \"}\"";
 332         type = fractional;
 333     } else
 334         assert(0);
 335
 336     token_free(tok);
 337     arg = evalue_read_term(s, p);
 338     tok = stream_next_token(s);
 339     if (!tok || tok->type != next) {
 340         stream_error(s, tok, error);
 341         if (tok)
 342             stream_push_token(s, tok);
 343     } else
 344         token_free(tok);
 345     pow = optional_power(s);
 346
 347     e = ALLOC(evalue);
 348     value_init(e->d);
 349     e->x.p = new_enode(type, pow+2, -1);
 350     value_clear(e->x.p->arr[0].d);
 351     e->x.p->arr[0] = *arg;
 352     free(arg);
 353     evalue_set_si(&e->x.p->arr[1+pow], 1, 1);
 354     while (--pow >= 0)
 355         evalue_set_si(&e->x.p->arr[1+pow], 0, 1);
 356
 357     return e;
 358 }
 359
 360 static evalue *evalue_read_factor(struct stream *s, struct parameter **p)
 361 {
 362     struct token *tok;
 363     evalue *e = NULL;
 364
 365     tok = stream_next_token(s);
 366     if (tok->type == '(') {
 367         token_free(tok);
 368         e = evalue_read_term(s, p);
 369         tok = stream_next_token(s);
 370         if (!tok || tok->type != ')') {
 371             stream_error(s, tok, "expecting \")\"");
 372             if (tok)
 373                 stream_push_token(s, tok);
 374         } else
 375             token_free(tok);
 376     } else if (tok->type == TOKEN_VALUE) {
 377         e = ALLOC(evalue);
 378         value_init(e->d);
 379         value_set_si(e->d, 1);
 380         value_init(e->x.n);
 381         value_assign(e->x.n, tok->u.v);
 382         token_free(tok);
 383         tok = stream_next_token(s);
 384         if (tok && tok->type == '/') {
 385             token_free(tok);
 386             tok = stream_next_token(s);
 387             if (!tok || tok->type != TOKEN_VALUE) {
 388                 stream_error(s, tok, "expecting denominator");
 389                 if (tok)
 390                     stream_push_token(s, tok);
 391                 return NULL;
 392             }
 393             value_assign(e->d, tok->u.v);
 394             token_free(tok);
 395         } else if (tok)
 396             stream_push_token(s, tok);
 397     } else if (tok->type == TOKEN_IDENT) {
 398         int pos = parameter_pos(p, tok->u.s);
 399         int pow = optional_power(s);
 400         token_free(tok);
 401         e = ALLOC(evalue);
 402         value_init(e->d);
 403         e->x.p = new_enode(polynomial, pow+1, pos+1);
 404         evalue_set_si(&e->x.p->arr[pow], 1, 1);
 405         while (--pow >= 0)
 406             evalue_set_si(&e->x.p->arr[pow], 0, 1);
 407     } else if (tok->type == '[' || tok->type == '{')
 408         e = read_fract_like(s, tok, p);
 409
 410     tok = stream_next_token(s);
 411     if (tok && tok->type == '*') {
 412         evalue *e2;
 413         token_free(tok);
 414         e2 = evalue_read_factor(s, p);
 415         if (!e2) {
 416             stream_error(s, NULL, "unexpected EOF");
 417             return NULL;
 418         }
 419         emul(e2, e);
 420         free_evalue_refs(e2);
 421         free(e2);
 422     } else if (tok)
 423         stream_push_token(s, tok);
 424
 425     return e;
 426 }
 427
 428 static evalue *evalue_read_term(struct stream *s, struct parameter **p)
 429 {
 430     struct token *tok;
 431     evalue *e = NULL;
 432
 433     e = evalue_read_factor(s, p);
 434     if (!e)
 435         return NULL;
 436
 437     tok = stream_next_token(s);
 438     if (!tok)
 439         return e;
 440
 441     if (tok->type == '+') {
 442         evalue *e2;
 443         token_free(tok);
 444         e2 = evalue_read_term(s, p);
 445         if (!e2) {
 446             stream_error(s, NULL, "unexpected EOF");
 447             return NULL;
 448         }
 449         eadd(e2, e);
 450         free_evalue_refs(e2);
 451         free(e2);
 452     } else
 453         stream_push_token(s, tok);
 454
 455     return e;
 456 }
 457
 458 struct constraint {
 459     int                  type;
 460     Vector              *v;
 461     struct constraint   *next;
 462 };
 463
 464 struct constraint *constraint_new()
 465 {
 466     struct constraint *c = ALLOC(struct constraint);
 467     c->type = -1;
 468     c->v = Vector_Alloc(16);
 469     c->next = NULL;
 470     return c;
 471 }
 472
 473 void constraint_free(struct constraint *c)
 474 {
 475     while (c) {
 476         struct constraint *next = c->next;
 477         Vector_Free(c->v);
 478         free(c);
 479         c = next;
 480     }
 481 }
 482
 483 void constraint_extend(struct constraint *c, int pos)
 484 {
 485     Vector *v;
 486     if (pos < c->v->Size)
 487         return;
 488
 489     v = Vector_Alloc((3*c->v->Size)/2);
 490     Vector_Copy(c->v->p, v->p, c->v->Size);
 491     Vector_Free(c->v);
 492     c->v = v;
 493 }
 494
 495 static int evalue_read_constraint(struct stream *s, struct parameter **p,
 496                                   struct constraint **constraints)
 497 {
 498     struct token *tok;
 499     struct constraint *c = NULL;
 500
 501     while ((tok = stream_next_token(s))) {
 502         struct token *tok2;
 503         int pos;
 504         if (tok->type == '+')
 505             token_free(tok);
 506         else if (tok->type == TOKEN_IDENT) {
 507             if (!c)
 508                 c = constraint_new();
 509             pos = parameter_pos(p, tok->u.s);
 510             constraint_extend(c, 1+pos);
 511             value_set_si(c->v->p[1+pos], 1);
 512             token_free(tok);
 513         } else if (tok->type == TOKEN_VALUE) {
 514             if (!c)
 515                 c = constraint_new();
 516             tok2 = stream_next_token(s);
 517             if (tok2 && tok2->type == TOKEN_IDENT) {
 518                 pos = parameter_pos(p, tok2->u.s);
 519                 constraint_extend(c, 1+pos);
 520                 value_assign(c->v->p[1+pos], tok->u.v);
 521                 token_free(tok);
 522                 token_free(tok2);
 523             } else {
 524                 if (tok2)
 525                     stream_push_token(s, tok2);
 526                 value_assign(c->v->p[0], tok->u.v);
 527                 token_free(tok);
 528             }
 529         } else if (tok->type == TOKEN_GE || tok->type == '=') {
 530             int type = tok->type == TOKEN_GE;
 531             token_free(tok);
 532             tok = stream_next_token(s);
 533             if (!tok || tok->type != TOKEN_VALUE || value_notzero_p(tok->u.v)) {
 534                 stream_error(s, tok, "expecting \"0\"");
 535                 if (tok)
 536                     stream_push_token(s, tok);
 537                 *constraints = NULL;
 538             } else {
 539                 c->type = type;
 540                 c->next = *constraints;
 541                 *constraints = c;
 542                 token_free(tok);
 543             }
 544             break;
 545         } else {
 546             if (!c)
 547                 stream_push_token(s, tok);
 548             else {
 549                 stream_error(s, tok, "unexpected token");
 550                 *constraints = NULL;
 551             }
 552             return 0;
 553         }
 554     }
 555     return tok != NULL;
 556 }
 557
 558 static struct constraint *evalue_read_domain(struct stream *s, struct parameter **p,
 559                                              unsigned MaxRays)
 560 {
 561     struct constraint *constraints = NULL;
 562     struct token *tok;
 563     int line;
 564
 565     tok = stream_next_token(s);
 566     if (!tok)
 567         return NULL;
 568     stream_push_token(s, tok);
 569
 570     line = tok->line;
 571     while (evalue_read_constraint(s, p, &constraints)) {
 572         tok = stream_next_token(s);
 573         if (tok) {
 574             stream_push_token(s, tok);
 575             /* empty line separates domain from evalue */
 576             if (tok->line > line+1)
 577                 break;
 578             line = tok->line;
 579         }
 580     }
 581     return constraints;
 582 }
 583
 584 struct section {
 585     struct constraint   *constraints;
 586     evalue              *e;
 587     struct section      *next;
 588 };
 589
 590 static char **extract_parameters(struct parameter *p, unsigned *nparam)
 591 {
 592     int i;
 593     char **params;
 594
 595     *nparam = p ? p->pos+1 : 0;
 596     params = ALLOCN(char *, *nparam);
 597     for (i = 0; i < *nparam; ++i) {
 598         struct parameter *next = p->next;
 599         params[p->pos] = p->name;
 600         free(p);
 601         p = next;
 602     }
 603     return params;
 604 }
 605
 606 static evalue *evalue_read_partition(struct stream *s, char ***ppp,
 607                                      unsigned *nparam, unsigned MaxRays)
 608 {
 609     struct section *part = NULL;
 610     struct constraint *constraints;
 611     evalue *e = NULL;
 612     int m = 0;
 613     struct parameter *p = NULL;
 614
 615     while ((constraints = evalue_read_domain(s, &p, MaxRays))) {
 616         evalue *e = evalue_read_term(s, &p);
 617         if (!e) {
 618             stream_error(s, NULL, "missing evalue");
 619             break;
 620         }
 621         struct section *sect = ALLOC(struct section);
 622         sect->constraints = constraints;
 623         sect->e = e;
 624         sect->next = part;
 625         part = sect;
 626         ++m;
 627     }
 628
 629     if (part) {
 630         Polyhedron *D;
 631         Matrix *M;
 632         int j;
 633
 634         *ppp = extract_parameters(p, nparam);
 635         e = ALLOC(evalue);
 636         value_init(e->d);
 637         e->x.p = new_enode(partition, 2*m, *nparam);
 638
 639         for (j = 0; j < m; ++j) {
 640             int n;
 641             struct section *next = part->next;
 642             struct constraint *c;
 643             constraints = part->constraints;
 644             for (n = 0, c = constraints; c; ++n, c = c->next)
 645                 ;
 646             M = Matrix_Alloc(n, 1+*nparam+1);
 647             while (--n >= 0) {
 648                 struct constraint *next = constraints->next;
 649                 Vector_Copy(constraints->v->p+1, M->p[n]+1, *nparam);
 650                 if (constraints->type)
 651                     value_set_si(M->p[n][0], 1);
 652                 value_assign(M->p[n][1+*nparam], constraints->v->p[0]);
 653                 constraints->next = NULL;
 654                 constraint_free(constraints);
 655                 constraints = next;
 656             }
 657             D = Constraints2Polyhedron(M, MaxRays);
 658             Matrix_Free(M);
 659             EVALUE_SET_DOMAIN(e->x.p->arr[2*j], D);
 660             value_clear(e->x.p->arr[2*j+1].d);
 661             e->x.p->arr[2*j+1] = *part->e;
 662             free(part->e);
 663             free(part);
 664             part = next;
 665         }
 666     }
 667     return e;
 668 }
 669
 670 static evalue *evalue_read(FILE *in, char ***ppp, unsigned *nparam,
 671                            unsigned MaxRays)
 672 {
 673     struct stream *s = stream_new(in);
 674     struct token *tok;
 675     evalue *e;
 676     struct parameter *p = NULL;
 677
 678     if (!(tok = stream_next_token(s)))
 679         return NULL;
 680
 681     if (tok->type == TOKEN_VALUE) {
 682         struct token *tok2 = stream_next_token(s);
 683         if (tok2)
 684             stream_push_token(s, tok2);
 685         stream_push_token(s, tok);
 686         if (tok2 && (tok2->type == TOKEN_IDENT || tok2->type == TOKEN_GE))
 687             e = evalue_read_partition(s, ppp, nparam, MaxRays);
 688         else {
 689             e = evalue_read_term(s, &p);
 690             *ppp = extract_parameters(p, nparam);
 691         }
 692     } else if (tok->type == TOKEN_IDENT) {
 693         stream_push_token(s, tok);
 694         e = evalue_read_partition(s, ppp, nparam, MaxRays);
 695     }
 696     stream_free(s);
 697     return e;
 698 }
 699
 700 int main(int argc, char **argv)
 701 {
 702     evalue *EP;
 703     char **all_vars = NULL;
 704     unsigned ntotal;
 705     struct options options;
 706     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 707     static struct argp argp = { argp_options, parse_opt, 0, 0, 0 };
 708     Polyhedron *U;
 709     piecewise_lst *pl = NULL;
 710
 711     argp_parse(&argp, argc, argv, 0, 0, &options);
 712
 713     EP = evalue_read(stdin, &all_vars, &ntotal, bv_options->MaxRays);
 714     assert(EP);
 715
 716     if (options.verbose)
 717         print_evalue(stderr, EP, all_vars);
 718
 719     if (options.nvar == -1)
 720         options.nvar = ntotal;
 721
 722     U = Universe_Polyhedron(ntotal - options.nvar);
 723
 724     exvector params;
 725     params = constructParameterVector(all_vars+options.nvar, ntotal-options.nvar);
 726
 727     pl = evalue_bernstein_coefficients(NULL, EP, U, params, bv_options);
 728     assert(pl);
 729     pl->maximize();
 730     cout << *pl << endl;
 731     delete pl;
 732
 733     free_evalue_refs(EP);
 734     free(EP);
 735
 736     Polyhedron_Free(U);
 737     Free_ParamNames(all_vars, ntotal);
 738     barvinok_options_free(bv_options);
 739     return 0;
 740 }