barvinok_enumerate.cc

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