barvinok_enumerate.cc

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