barvinok_enumerate.cc

   1 #include <unistd.h>
   2 #include <stdlib.h>
   3 #include <barvinok/evalue.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/barvinok.h>
   6 #include "fdstream.h"
   7 #include "argp.h"
   8 #include "verify.h"
   9 #include "verif_ehrhart.h"
  10 #include "remove_equalities.h"
  11 #include "evalue_convert.h"
  12
  13 #undef CS   /* for Solaris 10 */
  14
  15 /* The input of this example program is the same as that of testehrhart
  16  * in the PolyLib distribution, i.e., a polytope in combined
  17  * data and parameter space, a context polytope in parameter space
  18  * and (optionally) the names of the parameters.
  19  * Both polytopes are in PolyLib notation.
  20  */
  21
  22 #define PRINT_STATS         (BV_OPT_LAST+1)
  23
  24 struct argp_option argp_options[] = {
  25     { "size",      'S' },
  26     { "series",    's', 0, 0, "compute rational generating function" },
  27     { "explicit",  'e', 0, 0, "convert rgf to psp" },
  28     { "verbose",    'v' },
  29     { "print-stats",        PRINT_STATS,  0,    0 },
  30     { 0 }
  31 };
  32
  33 struct arguments {
  34     int size;
  35     int series;
  36     int function;
  37     int verbose;
  38     int print_stats;
  39     struct verify_options    verify;
  40     struct convert_options   convert;
  41 };
  42
  43 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  44 {
  45     struct arguments *options = (struct arguments*) state->input;
  46
  47     switch (key) {
  48     case ARGP_KEY_INIT:
  49         state->child_inputs[0] = options->verify.barvinok;
  50         state->child_inputs[1] = &options->verify;
  51         state->child_inputs[2] = &options->convert;
  52         options->size = 0;
  53         options->series = 0;
  54         options->function = 0;
  55         options->verbose = 0;
  56         options->print_stats = 0;
  57         break;
  58     case PRINT_STATS:
  59         options->print_stats = 1;
  60         break;
  61     case 'S':
  62         options->size = 1;
  63         break;
  64     case 'e':
  65         options->function = 1;
  66         /* fall through */
  67     case 's':
  68         options->series = 1;
  69         break;
  70     case 'v':
  71         options->verbose = 1;
  72         break;
  73     default:
  74         return ARGP_ERR_UNKNOWN;
  75     }
  76     return 0;
  77 }
  78
  79 struct skewed_gen_fun {
  80     gen_fun *gf;
  81     /* maps original space to space in which gf is defined */
  82     Matrix  *T;
  83     /* equalities in the original space that need to be satisfied for
  84      * gf to be valid
  85      */
  86     Matrix  *eq;
  87     /* divisibilities in the original space that need to be satisfied for
  88      * gf to be valid
  89      */
  90     Matrix  *div;
  91
  92     skewed_gen_fun(gen_fun *gf, Matrix *T, Matrix *eq, Matrix *div) :
  93                     gf(gf), T(T), eq(eq), div(div) {}
  94     ~skewed_gen_fun() {
  95         if (T)
  96             Matrix_Free(T);
  97         if (eq)
  98             Matrix_Free(eq);
  99         if (div)
 100             Matrix_Free(div);
 101         delete gf;
 102     }
 103
 104     void print(FILE *out, unsigned int nparam, char **param_name) const;
 105     operator evalue *() const {
 106         assert(T == NULL && eq == NULL); /* other cases not supported for now */
 107         return *gf;
 108     }
 109     void coefficient(Value* params, Value* c, barvinok_options *options) const;
 110 };
 111
 112 void skewed_gen_fun::print(FILE *out, unsigned int nparam,
 113                             char **param_name) const
 114 {
 115     fdostream os(dup(fileno(out)));
 116     if (T) {
 117         fprintf(out, "T:\n");
 118         Matrix_Print(out, P_VALUE_FMT, T);
 119     }
 120     if (eq) {
 121         fprintf(out, "eq:\n");
 122         Matrix_Print(out, P_VALUE_FMT, eq);
 123     }
 124     if (div) {
 125         fprintf(out, "div:\n");
 126         Matrix_Print(out, P_VALUE_FMT, div);
 127     }
 128     gf->print(os, nparam, param_name);
 129 }
 130
 131 void skewed_gen_fun::coefficient(Value* params, Value* c,
 132                                  barvinok_options *options) const
 133 {
 134     if (eq) {
 135         for (int i = 0; i < eq->NbRows; ++i) {
 136             Inner_Product(eq->p[i]+1, params, eq->NbColumns-2, eq->p[i]);
 137             if (value_notzero_p(eq->p[i][0])) {
 138                 value_set_si(*c, 0);
 139                 return;
 140             }
 141         }
 142     }
 143     if (div) {
 144         Value tmp;
 145         value_init(tmp);
 146         for (int i = 0; i < div->NbRows; ++i) {
 147             Inner_Product(div->p[i], params, div->NbColumns-1, &tmp);
 148             if (!mpz_divisible_p(tmp, div->p[i][div->NbColumns-1])) {
 149                 value_set_si(*c, 0);
 150                 return;
 151             }
 152         }
 153         value_clear(tmp);
 154     }
 155
 156     ZZ coeff;
 157     if (!T)
 158         coeff = gf->coefficient(params, options);
 159     else {
 160         Vector *p2 = Vector_Alloc(T->NbRows);
 161         Matrix_Vector_Product(T, params, p2->p);
 162         if (value_notone_p(p2->p[T->NbRows-1]))
 163             Vector_AntiScale(p2->p, p2->p, p2->p[T->NbRows-1], T->NbRows);
 164         coeff = gf->coefficient(p2->p, options);
 165         Vector_Free(p2);
 166     }
 167
 168     zz2value(coeff, *c);
 169 }
 170
 171 static int check_series(Polyhedron *S, Polyhedron *CS, skewed_gen_fun *gf,
 172                         int nparam, int pos, Value *z, verify_options *options)
 173 {
 174     int k;
 175     Value c, tmp;
 176     Value LB, UB;
 177
 178     value_init(c);
 179     value_init(tmp);
 180     value_init(LB);
 181     value_init(UB);
 182
 183     if (pos == nparam) {
 184         /* Computes the coefficient */
 185         gf->coefficient(&z[S->Dimension-nparam+1], &c, options->barvinok);
 186
 187         /* if c=0 we may be out of context. */
 188         /* scanning is useless in this case*/
 189
 190         if (options->print_all) {
 191             printf("EP( ");
 192             value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
 193             for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 194                 printf(", ");
 195                 value_print(stdout,VALUE_FMT,z[k]);
 196             }
 197             printf(" ) = ");
 198             value_print(stdout,VALUE_FMT,c);
 199             printf(" ");
 200         }
 201
 202         /* Manually count the number of points */
 203         count_points(1,S,z,&tmp);
 204         if (options->print_all) {
 205             printf(", count = ");
 206             value_print(stdout, P_VALUE_FMT, tmp);
 207             printf(". ");
 208         }
 209
 210         if (value_ne(tmp,c)) {
 211             printf("\n");
 212             fflush(stdout);
 213             fprintf(stderr,"Error !\n");
 214             fprintf(stderr,"EP( ");
 215             value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
 216             for (k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 217                 fprintf(stderr,", ");
 218                 value_print(stderr,VALUE_FMT,z[k]);
 219             }
 220             fprintf(stderr," ) should be ");
 221             value_print(stderr,VALUE_FMT,tmp);
 222             fprintf(stderr,", while EP eval gives ");
 223             value_print(stderr,VALUE_FMT,c);
 224             fprintf(stderr,".\n");
 225             if (!options->continue_on_error) {
 226                 value_clear(c); value_clear(tmp);
 227                 return 0;
 228             }
 229         } else if (options->print_all)
 230             printf("OK.\n");
 231     } else {
 232         int ok =
 233           !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 234         assert(ok);
 235         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 236             if (!options->print_all) {
 237                 k = VALUE_TO_INT(tmp);
 238                 if(!pos && !(k % options->st)) {
 239                     printf("o");
 240                     fflush(stdout);
 241                 }
 242             }
 243             value_assign(z[pos+S->Dimension-nparam+1],tmp);
 244             if (!check_series(S, CS->next, gf, nparam, pos+1, z, options)) {
 245                 value_clear(c); value_clear(tmp);
 246                 value_clear(LB);
 247                 value_clear(UB);
 248                 return(0);
 249             }
 250         }
 251         value_set_si(z[pos+S->Dimension-nparam+1],0);
 252     }
 253
 254     value_clear(c);
 255     value_clear(tmp);
 256     value_clear(LB);
 257     value_clear(UB);
 258     return 1;
 259 }
 260
 261 static int verify(Polyhedron *P, Polyhedron **C, evalue *EP, skewed_gen_fun *gf,
 262                    arguments *options)
 263 {
 264     Polyhedron *CC, *PP, *CS, *S;
 265     Matrix *C1;
 266     Vector *p;
 267     int result = 0;
 268
 269     /******* Compute true context *******/
 270     CC = align_context(*C, P->Dimension, options->verify.barvinok->MaxRays);
 271     PP = DomainIntersection(P, CC, options->verify.barvinok->MaxRays);
 272     Domain_Free(CC);
 273     C1 = Matrix_Alloc((*C)->Dimension+1, P->Dimension+1);
 274
 275     for (int i = 0; i < C1->NbRows; i++)
 276         for (int j = 0; j < C1->NbColumns; j++)
 277             if (i == j-P->Dimension+(*C)->Dimension)
 278                 value_set_si(C1->p[i][j], 1);
 279             else
 280                 value_set_si(C1->p[i][j], 0);
 281     CC = Polyhedron_Image(PP, C1, options->verify.barvinok->MaxRays);
 282     Matrix_Free(C1);
 283     Domain_Free(PP);
 284     Domain_Free(*C);
 285     *C = CC;
 286
 287     CS = check_poly_context_scan(*C, &options->verify);
 288
 289     p = Vector_Alloc(P->Dimension+2);
 290     value_set_si(p->p[P->Dimension+1], 1);
 291
 292     /* S = scanning list of polyhedra */
 293     S = Polyhedron_Scan(P, *C, options->verify.barvinok->MaxRays);
 294
 295     check_poly_init(*C, &options->verify);
 296
 297     /******* CHECK NOW *********/
 298     if (S) {
 299         if (!options->series || options->function) {
 300             if (!check_poly_EP(S, CS, EP, 0, (*C)->Dimension, 0, p->p,
 301                                 &options->verify))
 302                 result = -1;
 303         } else {
 304             if (!check_series(S, CS, gf, (*C)->Dimension, 0, p->p, &options->verify))
 305                 result = -1;
 306         }
 307         Domain_Free(S);
 308     }
 309
 310     if (result == -1)
 311         fprintf(stderr,"Check failed !\n");
 312
 313     if (!options->verify.print_all)
 314         printf( "\n" );
 315
 316     Vector_Free(p);
 317     if (CS)
 318         Domain_Free(CS);
 319
 320     return result;
 321 }
 322
 323 static void unimodular_complete(Matrix *M, int row)
 324 {
 325     Matrix *H, *Q, *U;
 326     left_hermite(M, &H, &Q, &U);
 327     Matrix_Free(H);
 328     Matrix_Free(U);
 329     for (int r = row; r < M->NbRows; ++r)
 330         Vector_Copy(Q->p[r], M->p[r], M->NbColumns);
 331     Matrix_Free(Q);
 332 }
 333
 334 /* frees M and Minv */
 335 static void apply_transformation(Polyhedron **P, Polyhedron **C,
 336                                  bool free_P, bool free_C,
 337                                  Matrix *M, Matrix *Minv, Matrix **inv,
 338                                  barvinok_options *options)
 339 {
 340     Polyhedron *T;
 341     Matrix *M2;
 342
 343     M2 = align_matrix(M, (*P)->Dimension + 1);
 344     T = *P;
 345     *P = Polyhedron_Preimage(*P, M2, options->MaxRays);
 346     if (free_P)
 347         Polyhedron_Free(T);
 348     Matrix_Free(M2);
 349
 350     T = *C;
 351     *C = Polyhedron_Preimage(*C, M, options->MaxRays);
 352     if (free_C)
 353         Polyhedron_Free(T);
 354
 355     Matrix_Free(M);
 356
 357     if (*inv) {
 358         Matrix *T = *inv;
 359         *inv = Matrix_Alloc(Minv->NbRows, T->NbColumns);
 360         Matrix_Product(Minv, T, *inv);
 361         Matrix_Free(T);
 362         Matrix_Free(Minv);
 363     } else
 364         *inv = Minv;
 365 }
 366
 367 static skewed_gen_fun *series(Polyhedron *P, Polyhedron* C,
 368                                 barvinok_options *options)
 369 {
 370     Polyhedron *C1, *C2;
 371     gen_fun *gf;
 372     Matrix *inv = NULL;
 373     Matrix *eq = NULL;
 374     Matrix *div = NULL;
 375     Polyhedron *PT = P;
 376
 377     /* Compute true context */
 378     C1 = Polyhedron_Project(P, C->Dimension);
 379     C2 = DomainIntersection(C, C1, options->MaxRays);
 380     Polyhedron_Free(C1);
 381
 382     POL_ENSURE_VERTICES(C2);
 383     if (C2->NbBid != 0) {
 384         Polyhedron *T;
 385         Matrix *M, *Minv, *M2;
 386         Matrix *CP;
 387         if (C2->NbEq || P->NbEq) {
 388             /* We remove all equalities to be sure all lines are unit vectors */
 389             Polyhedron *CT = C2;
 390             remove_all_equalities(&PT, &CT, &CP, NULL, C2->Dimension,
 391                                   options->MaxRays);
 392             if (CT != C2) {
 393                 Polyhedron_Free(C2);
 394                 C2 = CT;
 395             }
 396             if (CP) {
 397                 inv = left_inverse(CP, &eq);
 398                 Matrix_Free(CP);
 399
 400                 int d = 0;
 401                 Value tmp;
 402                 value_init(tmp);
 403                 div = Matrix_Alloc(inv->NbRows-1, inv->NbColumns+1);
 404                 for (int i = 0; i < inv->NbRows-1; ++i) {
 405                     Vector_Gcd(inv->p[i], inv->NbColumns, &tmp);
 406                     if (mpz_divisible_p(tmp,
 407                                         inv->p[inv->NbRows-1][inv->NbColumns-1]))
 408                         continue;
 409                     Vector_Copy(inv->p[i], div->p[d], inv->NbColumns);
 410                     value_assign(div->p[d][inv->NbColumns],
 411                                  inv->p[inv->NbRows-1][inv->NbColumns-1]);
 412                     ++d;
 413                 }
 414                 value_clear(tmp);
 415
 416                 if (!d) {
 417                     Matrix_Free(div);
 418                     div = NULL;
 419                 } else
 420                     div->NbRows = d;
 421             }
 422         }
 423         POL_ENSURE_VERTICES(C2);
 424
 425         /* Since we have "compressed" the parameters (in case there were
 426          * any equalities), the result is independent of the coordinates in the
 427          * coordinate subspace spanned by the lines.  We can therefore assume
 428          * these coordinates are zero and compute the inverse image of the map
 429          * from a lower dimensional space that adds zeros in the appropriate
 430          * places.
 431          */
 432         M = Matrix_Alloc(C2->Dimension+1, C2->Dimension-C2->NbBid+1);
 433         int k = 0;
 434         for (int i = 0; i < C2->NbBid; ++i) {
 435             int j = First_Non_Zero(C2->Ray[i]+1, C2->Dimension);
 436             assert(First_Non_Zero(C2->Ray[i]+1+j+1, C2->Dimension-j-1) == -1);
 437             for ( ; k < j; k++)
 438                 value_set_si(M->p[k+i][k], 1);
 439         }
 440         for ( ; k < C2->Dimension-C2->NbBid+1; k++)
 441             value_set_si(M->p[k+C2->NbBid][k], 1);
 442         Minv = Transpose(M);
 443
 444         apply_transformation(&PT, &C2, PT != P, C2 != C, M, Minv, &inv, options);
 445     }
 446     POL_ENSURE_VERTICES(C2);
 447     if (!Polyhedron_has_revlex_positive_rays(C2, C2->Dimension)) {
 448         Polyhedron *T;
 449         Matrix *Constraints;
 450         Matrix *H, *Q, *U;
 451         Constraints = Matrix_Alloc(C2->NbConstraints, C2->Dimension+1);
 452         for (int i = 0; i < C2->NbConstraints; ++i)
 453             Vector_Copy(C2->Constraint[i]+1, Constraints->p[i], C2->Dimension);
 454         left_hermite(Constraints, &H, &Q, &U);
 455         /* flip rows of Q */
 456         for (int i = 0; i < C2->Dimension/2; ++i)
 457             Vector_Exchange(Q->p[i], Q->p[C2->Dimension-1-i], C2->Dimension);
 458         Matrix_Free(H);
 459         Matrix_Free(U);
 460         Matrix *M = Matrix_Alloc(C2->Dimension+1, C2->Dimension+1);
 461         U = Matrix_Copy(Q);
 462         int ok = Matrix_Inverse(U, M);
 463         assert(ok);
 464         Matrix_Free(U);
 465
 466         apply_transformation(&PT, &C2, PT != P, C2 != C, M, Q, &inv, options);
 467     }
 468     gf = barvinok_series_with_options(PT, C2, options);
 469     Polyhedron_Free(C2);
 470     if (PT != P)
 471         Polyhedron_Free(PT);
 472     return new skewed_gen_fun(gf, inv, eq, div);
 473 }
 474
 475 int main(int argc, char **argv)
 476 {
 477     Polyhedron *A, *C;
 478     Matrix *M;
 479     evalue *EP = NULL;
 480     skewed_gen_fun *gf = NULL;
 481     char **param_name;
 482     int print_solution = 1;
 483     int result = 0;
 484     struct arguments options;
 485     static struct argp_child argp_children[] = {
 486         { &barvinok_argp,       0,      0,              0 },
 487         { &verify_argp,         0,      "verification", 1 },
 488         { &convert_argp,        0,      "output conversion",    2 },
 489         { 0 }
 490     };
 491     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 492     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 493
 494     options.verify.barvinok = bv_options;
 495     argp_parse(&argp, argc, argv, 0, 0, &options);
 496
 497     M = Matrix_Read();
 498     A = Constraints2Polyhedron(M, bv_options->MaxRays);
 499     Matrix_Free(M);
 500     M = Matrix_Read();
 501     C = Constraints2Polyhedron(M, bv_options->MaxRays);
 502     Matrix_Free(M);
 503     param_name = Read_ParamNames(stdin, C->Dimension);
 504
 505     if (options.verify.verify) {
 506         verify_options_set_range(&options.verify, A->Dimension);
 507         if (!options.verbose)
 508             print_solution = 0;
 509     }
 510
 511     if (print_solution && options.verbose) {
 512         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 513         Polyhedron_Print(stdout, P_VALUE_FMT, C);
 514     }
 515
 516     if (options.series) {
 517         gf = series(A, C, bv_options);
 518         if (print_solution) {
 519             gf->print(stdout, C->Dimension, param_name);
 520             puts("");
 521         }
 522         if (options.function) {
 523             EP = *gf;
 524             if (print_solution)
 525                 print_evalue(stdout, EP, param_name);
 526         }
 527     } else {
 528         EP = barvinok_enumerate_with_options(A, C, bv_options);
 529         evalue_convert(EP, &options.convert, C->Dimension,
 530                        options.verbose ? param_name : NULL);
 531         if (options.size)
 532             printf("\nSize: %d\n", evalue_size(EP));
 533     }
 534     if (print_solution && !options.verbose)
 535         print_evalue(stdout, EP, param_name);
 536
 537     if (options.verify.verify) {
 538         options.verify.params = param_name;
 539         result = verify(A, &C, EP, gf, &options);
 540     }
 541
 542     if (gf)
 543         delete gf;
 544     if (EP) {
 545         free_evalue_refs(EP);
 546         free(EP);
 547     }
 548
 549     if (options.print_stats)
 550         barvinok_stats_print(options.verify.barvinok->stats, stdout);
 551
 552     Free_ParamNames(param_name, C->Dimension);
 553     Polyhedron_Free(A);
 554     Polyhedron_Free(C);
 555     barvinok_options_free(bv_options);
 556     return result;
 557 }