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