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