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)
 367         return NULL;
 368
 369     if (tok->type == '(') {
 370         token_free(tok);
 371         e = evalue_read_term(s, p);
 372         tok = stream_next_token(s);
 373         if (!tok || tok->type != ')') {
 374             stream_error(s, tok, "expecting \")\"");
 375             if (tok)
 376                 stream_push_token(s, tok);
 377         } else
 378             token_free(tok);
 379     } else if (tok->type == TOKEN_VALUE) {
 380         e = ALLOC(evalue);
 381         value_init(e->d);
 382         value_set_si(e->d, 1);
 383         value_init(e->x.n);
 384         value_assign(e->x.n, tok->u.v);
 385         token_free(tok);
 386         tok = stream_next_token(s);
 387         if (tok && tok->type == '/') {
 388             token_free(tok);
 389             tok = stream_next_token(s);
 390             if (!tok || tok->type != TOKEN_VALUE) {
 391                 stream_error(s, tok, "expecting denominator");
 392                 if (tok)
 393                     stream_push_token(s, tok);
 394                 return NULL;
 395             }
 396             value_assign(e->d, tok->u.v);
 397             token_free(tok);
 398         } else if (tok)
 399             stream_push_token(s, tok);
 400     } else if (tok->type == TOKEN_IDENT) {
 401         int pos = parameter_pos(p, tok->u.s);
 402         int pow = optional_power(s);
 403         token_free(tok);
 404         e = ALLOC(evalue);
 405         value_init(e->d);
 406         e->x.p = new_enode(polynomial, pow+1, pos+1);
 407         evalue_set_si(&e->x.p->arr[pow], 1, 1);
 408         while (--pow >= 0)
 409             evalue_set_si(&e->x.p->arr[pow], 0, 1);
 410     } else if (tok->type == '[' || tok->type == '{')
 411         e = read_fract_like(s, tok, p);
 412
 413     tok = stream_next_token(s);
 414     if (tok && tok->type == '*') {
 415         evalue *e2;
 416         token_free(tok);
 417         e2 = evalue_read_factor(s, p);
 418         if (!e2) {
 419             stream_error(s, NULL, "unexpected EOF");
 420             return NULL;
 421         }
 422         emul(e2, e);
 423         free_evalue_refs(e2);
 424         free(e2);
 425     } else if (tok)
 426         stream_push_token(s, tok);
 427
 428     return e;
 429 }
 430
 431 static evalue *evalue_read_term(struct stream *s, struct parameter **p)
 432 {
 433     struct token *tok;
 434     evalue *e = NULL;
 435
 436     e = evalue_read_factor(s, p);
 437     if (!e)
 438         return NULL;
 439
 440     tok = stream_next_token(s);
 441     if (!tok)
 442         return e;
 443
 444     if (tok->type == '+') {
 445         evalue *e2;
 446         token_free(tok);
 447         e2 = evalue_read_term(s, p);
 448         if (!e2) {
 449             stream_error(s, NULL, "unexpected EOF");
 450             return NULL;
 451         }
 452         eadd(e2, e);
 453         free_evalue_refs(e2);
 454         free(e2);
 455     } else
 456         stream_push_token(s, tok);
 457
 458     return e;
 459 }
 460
 461 struct constraint {
 462     int                  type;
 463     Vector              *v;
 464     struct constraint   *next;
 465 };
 466
 467 struct constraint *constraint_new()
 468 {
 469     struct constraint *c = ALLOC(struct constraint);
 470     c->type = -1;
 471     c->v = Vector_Alloc(16);
 472     c->next = NULL;
 473     return c;
 474 }
 475
 476 void constraint_free(struct constraint *c)
 477 {
 478     while (c) {
 479         struct constraint *next = c->next;
 480         Vector_Free(c->v);
 481         free(c);
 482         c = next;
 483     }
 484 }
 485
 486 void constraint_extend(struct constraint *c, int pos)
 487 {
 488     Vector *v;
 489     if (pos < c->v->Size)
 490         return;
 491
 492     v = Vector_Alloc((3*c->v->Size)/2);
 493     Vector_Copy(c->v->p, v->p, c->v->Size);
 494     Vector_Free(c->v);
 495     c->v = v;
 496 }
 497
 498 static int evalue_read_constraint(struct stream *s, struct parameter **p,
 499                                   struct constraint **constraints)
 500 {
 501     struct token *tok;
 502     struct constraint *c = NULL;
 503
 504     while ((tok = stream_next_token(s))) {
 505         struct token *tok2;
 506         int pos;
 507         if (tok->type == '+')
 508             token_free(tok);
 509         else if (tok->type == TOKEN_IDENT) {
 510             if (!c)
 511                 c = constraint_new();
 512             pos = parameter_pos(p, tok->u.s);
 513             constraint_extend(c, 1+pos);
 514             value_set_si(c->v->p[1+pos], 1);
 515             token_free(tok);
 516         } else if (tok->type == TOKEN_VALUE) {
 517             if (!c)
 518                 c = constraint_new();
 519             tok2 = stream_next_token(s);
 520             if (tok2 && tok2->type == TOKEN_IDENT) {
 521                 pos = parameter_pos(p, tok2->u.s);
 522                 constraint_extend(c, 1+pos);
 523                 value_assign(c->v->p[1+pos], tok->u.v);
 524                 token_free(tok);
 525                 token_free(tok2);
 526             } else {
 527                 if (tok2)
 528                     stream_push_token(s, tok2);
 529                 value_assign(c->v->p[0], tok->u.v);
 530                 token_free(tok);
 531             }
 532         } else if (tok->type == TOKEN_GE || tok->type == '=') {
 533             int type = tok->type == TOKEN_GE;
 534             token_free(tok);
 535             tok = stream_next_token(s);
 536             if (!tok || tok->type != TOKEN_VALUE || value_notzero_p(tok->u.v)) {
 537                 stream_error(s, tok, "expecting \"0\"");
 538                 if (tok)
 539                     stream_push_token(s, tok);
 540                 *constraints = NULL;
 541             } else {
 542                 c->type = type;
 543                 c->next = *constraints;
 544                 *constraints = c;
 545                 token_free(tok);
 546             }
 547             break;
 548         } else {
 549             if (!c)
 550                 stream_push_token(s, tok);
 551             else {
 552                 stream_error(s, tok, "unexpected token");
 553                 *constraints = NULL;
 554             }
 555             return 0;
 556         }
 557     }
 558     return tok != NULL;
 559 }
 560
 561 static struct constraint *evalue_read_domain(struct stream *s, struct parameter **p,
 562                                              unsigned MaxRays)
 563 {
 564     struct constraint *constraints = NULL;
 565     struct token *tok;
 566     int line;
 567
 568     tok = stream_next_token(s);
 569     if (!tok)
 570         return NULL;
 571     stream_push_token(s, tok);
 572
 573     line = tok->line;
 574     while (evalue_read_constraint(s, p, &constraints)) {
 575         tok = stream_next_token(s);
 576         if (tok) {
 577             stream_push_token(s, tok);
 578             /* empty line separates domain from evalue */
 579             if (tok->line > line+1)
 580                 break;
 581             line = tok->line;
 582         }
 583     }
 584     return constraints;
 585 }
 586
 587 struct section {
 588     struct constraint   *constraints;
 589     evalue              *e;
 590     struct section      *next;
 591 };
 592
 593 static char **extract_parameters(struct parameter *p, unsigned *nparam)
 594 {
 595     int i;
 596     char **params;
 597
 598     *nparam = p ? p->pos+1 : 0;
 599     params = ALLOCN(char *, *nparam);
 600     for (i = 0; i < *nparam; ++i) {
 601         struct parameter *next = p->next;
 602         params[p->pos] = p->name;
 603         free(p);
 604         p = next;
 605     }
 606     return params;
 607 }
 608
 609 static evalue *evalue_read_partition(struct stream *s, char ***ppp,
 610                                      unsigned *nparam, unsigned MaxRays)
 611 {
 612     struct section *part = NULL;
 613     struct constraint *constraints;
 614     evalue *e = NULL;
 615     int m = 0;
 616     struct parameter *p = NULL;
 617
 618     while ((constraints = evalue_read_domain(s, &p, MaxRays))) {
 619         evalue *e = evalue_read_term(s, &p);
 620         if (!e) {
 621             stream_error(s, NULL, "missing evalue");
 622             break;
 623         }
 624         struct section *sect = ALLOC(struct section);
 625         sect->constraints = constraints;
 626         sect->e = e;
 627         sect->next = part;
 628         part = sect;
 629         ++m;
 630     }
 631
 632     if (part) {
 633         Polyhedron *D;
 634         Matrix *M;
 635         int j;
 636
 637         *ppp = extract_parameters(p, nparam);
 638         e = ALLOC(evalue);
 639         value_init(e->d);
 640         e->x.p = new_enode(partition, 2*m, *nparam);
 641
 642         for (j = 0; j < m; ++j) {
 643             int n;
 644             struct section *next = part->next;
 645             struct constraint *c;
 646             constraints = part->constraints;
 647             for (n = 0, c = constraints; c; ++n, c = c->next)
 648                 ;
 649             M = Matrix_Alloc(n, 1+*nparam+1);
 650             while (--n >= 0) {
 651                 struct constraint *next = constraints->next;
 652                 Vector_Copy(constraints->v->p+1, M->p[n]+1, *nparam);
 653                 if (constraints->type)
 654                     value_set_si(M->p[n][0], 1);
 655                 value_assign(M->p[n][1+*nparam], constraints->v->p[0]);
 656                 constraints->next = NULL;
 657                 constraint_free(constraints);
 658                 constraints = next;
 659             }
 660             D = Constraints2Polyhedron(M, MaxRays);
 661             Matrix_Free(M);
 662             EVALUE_SET_DOMAIN(e->x.p->arr[2*j], D);
 663             value_clear(e->x.p->arr[2*j+1].d);
 664             e->x.p->arr[2*j+1] = *part->e;
 665             free(part->e);
 666             free(part);
 667             part = next;
 668         }
 669     }
 670     return e;
 671 }
 672
 673 static evalue *evalue_read(FILE *in, char ***ppp, unsigned *nparam,
 674                            unsigned MaxRays)
 675 {
 676     struct stream *s = stream_new(in);
 677     struct token *tok;
 678     evalue *e;
 679     struct parameter *p = NULL;
 680
 681     if (!(tok = stream_next_token(s)))
 682         return NULL;
 683
 684     if (tok->type == TOKEN_VALUE) {
 685         struct token *tok2 = stream_next_token(s);
 686         if (tok2)
 687             stream_push_token(s, tok2);
 688         stream_push_token(s, tok);
 689         if (tok2 && (tok2->type == TOKEN_IDENT || tok2->type == TOKEN_GE))
 690             e = evalue_read_partition(s, ppp, nparam, MaxRays);
 691         else {
 692             e = evalue_read_term(s, &p);
 693             *ppp = extract_parameters(p, nparam);
 694         }
 695     } else if (tok->type == TOKEN_IDENT) {
 696         stream_push_token(s, tok);
 697         e = evalue_read_partition(s, ppp, nparam, MaxRays);
 698     }
 699     stream_free(s);
 700     return e;
 701 }
 702
 703 int main(int argc, char **argv)
 704 {
 705     evalue *EP;
 706     char **all_vars = NULL;
 707     unsigned ntotal;
 708     struct options options;
 709     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 710     static struct argp argp = { argp_options, parse_opt, 0, 0, 0 };
 711     Polyhedron *U;
 712     piecewise_lst *pl = NULL;
 713
 714     argp_parse(&argp, argc, argv, 0, 0, &options);
 715
 716     EP = evalue_read(stdin, &all_vars, &ntotal, bv_options->MaxRays);
 717     assert(EP);
 718
 719     if (options.verbose)
 720         print_evalue(stderr, EP, all_vars);
 721
 722     if (options.nvar == -1)
 723         options.nvar = ntotal;
 724
 725     U = Universe_Polyhedron(ntotal - options.nvar);
 726
 727     exvector params;
 728     params = constructParameterVector(all_vars+options.nvar, ntotal-options.nvar);
 729
 730     pl = evalue_bernstein_coefficients(NULL, EP, U, params, bv_options);
 731     assert(pl);
 732     pl->maximize();
 733     cout << *pl << endl;
 734     delete pl;
 735
 736     free_evalue_refs(EP);
 737     free(EP);
 738
 739     Polyhedron_Free(U);
 740     Free_ParamNames(all_vars, ntotal);
 741     barvinok_options_free(bv_options);
 742     return 0;
 743 }