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