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 static __isl_give isl_set *set_bounds(__isl_take isl_set *set,
  74         const struct verify_options *options)
  75 {
  76         int i;
  77         unsigned nparam;
  78         isl_point *pt, *pt2;
  79         isl_set *box;
  80
  81         nparam = isl_set_dim(set, isl_dim_param);
  82
  83         if (options->r > 0) {
  84                 pt = isl_set_sample_point(isl_set_copy(set));
  85                 pt2 = isl_point_copy(pt);
  86
  87                 for (i = 0; i < nparam; ++i) {
  88                         pt = isl_point_add_ui(pt, isl_dim_param, i, options->r);
  89                         pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, options->r);
  90                 }
  91         } else {
  92                 isl_int v;
  93
  94                 isl_int_init(v);
  95                 pt = isl_point_zero(isl_set_get_dim(set));
  96                 isl_int_set_si(v, options->m);
  97                 for (i = 0; i < nparam; ++i)
  98                         pt = isl_point_set_coordinate(pt, isl_dim_param, i, v);
  99
 100                 pt2 = isl_point_zero(isl_set_get_dim(set));
 101                 isl_int_set_si(v, options->M);
 102                 for (i = 0; i < nparam; ++i)
 103                         pt2 = isl_point_set_coordinate(pt2, isl_dim_param, i, v);
 104
 105                 isl_int_clear(v);
 106         }
 107
 108         box = isl_set_box_from_points(pt, pt2);
 109
 110         return isl_set_intersect(set, box);
 111 }
 112
 113 struct verify_point_data {
 114         const struct verify_options *options;
 115         isl_set *set;
 116         isl_pw_qpolynomial *pwqp;
 117         int n;
 118         int s;
 119         int error;
 120 };
 121
 122 static int verify_point(__isl_take isl_point *pnt, void *user)
 123 {
 124         struct verify_point_data *vpd = (struct verify_point_data *) user;
 125         isl_set *set;
 126         int i;
 127         unsigned nparam;
 128         isl_int v, n, d;
 129         isl_qpolynomial *cnt = NULL;
 130         int pa = vpd->options->barvinok->polynomial_approximation;
 131         int cst;
 132         int ok;
 133         FILE *out = vpd->options->print_all ? stdout : stderr;
 134
 135         vpd->n--;
 136
 137         isl_int_init(v);
 138         isl_int_init(n);
 139         isl_int_init(d);
 140         set = isl_set_copy(vpd->set);
 141         nparam = isl_set_dim(set, isl_dim_param);
 142         for (i = 0; i < nparam; ++i) {
 143                 isl_point_get_coordinate(pnt, isl_dim_param, i, &v);
 144                 set = isl_set_fix(set, isl_dim_param, i, v);
 145         }
 146
 147         if (isl_set_count(set, &v) < 0)
 148                 goto error;
 149
 150         cnt = isl_pw_qpolynomial_eval(isl_pw_qpolynomial_copy(vpd->pwqp),
 151                                         isl_point_copy(pnt));
 152
 153         cst = isl_qpolynomial_is_cst(cnt, &n, &d);
 154         if (cst != 1)
 155                 goto error;
 156
 157         if (pa == BV_APPROX_SIGN_LOWER)
 158                 isl_int_cdiv_q(n, n, d);
 159         else if (pa == BV_APPROX_SIGN_UPPER)
 160                 isl_int_fdiv_q(n, n, d);
 161         else
 162                 isl_int_tdiv_q(n, n, d);
 163
 164         if (pa == BV_APPROX_SIGN_APPROX)
 165                 /* just accept everything */
 166                 ok = 1;
 167         else if (pa == BV_APPROX_SIGN_LOWER)
 168                 ok = isl_int_le(n, v);
 169         else if (pa == BV_APPROX_SIGN_UPPER)
 170                 ok = isl_int_ge(n, v);
 171         else
 172                 ok = isl_int_eq(n, v);
 173
 174         if (vpd->options->print_all || !ok) {
 175                 fprintf(out, "EP(");
 176                 for (i = 0; i < nparam; ++i) {
 177                         if (i)
 178                                 fprintf(out, ", ");
 179                         isl_point_get_coordinate(pnt, isl_dim_param, i, &d);
 180                         isl_int_print(out, d, 0);
 181                 }
 182                 fprintf(out, ") = ");
 183                 isl_int_print(out, n, 0);
 184                 fprintf(out, ", count = ");
 185                 isl_int_print(out, v, 0);
 186                 if (ok)
 187                         fprintf(out, ". OK\n");
 188                 else
 189                         fprintf(out, ". NOT OK\n");
 190         } else if ((vpd->n % vpd->s) == 0) {
 191                 printf("o");
 192                 fflush(stdout);
 193         }
 194
 195         if (0) {
 196 error:
 197                 ok = 0;
 198         }
 199         isl_set_free(set);
 200         isl_qpolynomial_free(cnt);
 201         isl_int_clear(v);
 202         isl_int_clear(n);
 203         isl_int_clear(d);
 204         isl_point_free(pnt);
 205
 206         if (!ok)
 207                 vpd->error = 1;
 208
 209         if (vpd->options->continue_on_error)
 210                 ok = 1;
 211
 212         return (vpd->n >= 1 && ok) ? 0 : -1;
 213 }
 214
 215 static int verify_isl(Polyhedron *P, Polyhedron *C,
 216                 evalue *EP, const struct verify_options *options)
 217 {
 218         struct verify_point_data vpd = { options };
 219         int i;
 220         isl_ctx *ctx = isl_ctx_alloc();
 221         isl_dim *dim;
 222         isl_set *set;
 223         isl_set *set_C;
 224         isl_int v;
 225         int r;
 226
 227         dim = isl_dim_set_alloc(ctx, C->Dimension, P->Dimension - C->Dimension);
 228         for (i = 0; i < C->Dimension; ++i)
 229                 dim = isl_dim_set_name(dim, isl_dim_param, i, options->params[i]);
 230         set = isl_set_new_from_polylib(P, isl_dim_copy(dim));
 231         dim = isl_dim_drop(dim, isl_dim_set, 0, P->Dimension - C->Dimension);
 232         set_C = isl_set_new_from_polylib(C, dim);
 233         set_C = isl_set_intersect(isl_set_copy(set), set_C);
 234         set_C = isl_set_remove(set_C, isl_dim_set, 0, P->Dimension - C->Dimension);
 235
 236         set_C = set_bounds(set_C, options);
 237
 238         isl_int_init(v);
 239         r = isl_set_count(set_C, &v);
 240         vpd.n = isl_int_cmp_si(v, 200) < 0 ? isl_int_get_si(v) : 200;
 241         isl_int_clear(v);
 242
 243         if (!options->print_all) {
 244                 vpd.s = vpd.n < 80 ? 1 : 1 + vpd.n/80;
 245                 for (i = 0; i < vpd.n; i += vpd.s)
 246                         printf(".");
 247                 printf("\r");
 248                 fflush(stdout);
 249         }
 250
 251         vpd.set = set;
 252         vpd.pwqp = isl_pw_qpolynomial_from_evalue(isl_set_get_dim(set_C), EP);
 253         vpd.error = 0;
 254         if (r == 0)
 255                 isl_set_foreach_point(set_C, verify_point, &vpd);
 256         if (vpd.error)
 257                 r = -1;
 258
 259         isl_pw_qpolynomial_free(vpd.pwqp);
 260         isl_set_free(set);
 261         isl_set_free(set_C);
 262
 263         isl_ctx_free(ctx);
 264
 265         if (!options->print_all)
 266                 printf("\n");
 267
 268         if (r < 0)
 269                 fprintf(stderr, "Check failed !\n");
 270
 271         return r;
 272 }
 273
 274 static int verify(Polyhedron *P, Polyhedron *C, evalue *EP, skewed_gen_fun *gf,
 275                    arguments *options)
 276 {
 277     Polyhedron *CS, *S;
 278     Vector *p;
 279     int result = 0;
 280
 281     if (!options->series || options->function)
 282         return verify_isl(P, C, EP, &options->verify);
 283
 284     CS = check_poly_context_scan(P, &C, C->Dimension, &options->verify);
 285
 286     p = Vector_Alloc(P->Dimension+2);
 287     value_set_si(p->p[P->Dimension+1], 1);
 288
 289     /* S = scanning list of polyhedra */
 290     S = Polyhedron_Scan(P, C, options->verify.barvinok->MaxRays);
 291
 292     check_poly_init(C, &options->verify);
 293
 294     /******* CHECK NOW *********/
 295     if (S) {
 296         if (!options->series || options->function) {
 297             if (!check_poly_EP(S, CS, EP, 0, C->Dimension, 0, p->p,
 298                                 &options->verify))
 299                 result = -1;
 300         } else {
 301             if (!check_poly_gf(S, CS, gf, 0, C->Dimension, 0, p->p,
 302                                 &options->verify))
 303                 result = -1;
 304         }
 305         Domain_Free(S);
 306     }
 307
 308     if (result == -1)
 309         fprintf(stderr,"Check failed !\n");
 310
 311     if (!options->verify.print_all)
 312         printf( "\n" );
 313
 314     Vector_Free(p);
 315     if (CS) {
 316         Domain_Free(CS);
 317         Domain_Free(C);
 318     }
 319
 320     return result;
 321 }
 322
 323 /* frees M and Minv */
 324 static void apply_transformation(Polyhedron **P, Polyhedron **C,
 325                                  bool free_P, bool free_C,
 326                                  Matrix *M, Matrix *Minv, Matrix **inv,
 327                                  barvinok_options *options)
 328 {
 329     Polyhedron *T;
 330     Matrix *M2;
 331
 332     M2 = align_matrix(M, (*P)->Dimension + 1);
 333     T = *P;
 334     *P = Polyhedron_Preimage(*P, M2, options->MaxRays);
 335     if (free_P)
 336         Polyhedron_Free(T);
 337     Matrix_Free(M2);
 338
 339     T = *C;
 340     *C = Polyhedron_Preimage(*C, M, options->MaxRays);
 341     if (free_C)
 342         Polyhedron_Free(T);
 343
 344     Matrix_Free(M);
 345
 346     if (*inv) {
 347         Matrix *T = *inv;
 348         *inv = Matrix_Alloc(Minv->NbRows, T->NbColumns);
 349         Matrix_Product(Minv, T, *inv);
 350         Matrix_Free(T);
 351         Matrix_Free(Minv);
 352     } else
 353         *inv = Minv;
 354 }
 355
 356 /* Since we have "compressed" the parameters (in case there were
 357  * any equalities), the result is independent of the coordinates in the
 358  * coordinate subspace spanned by the lines.  We can therefore assume
 359  * these coordinates are zero and compute the inverse image of the map
 360  * from a lower dimensional space that adds zeros in the appropriate
 361  * places.
 362  */
 363 static void remove_lines(Polyhedron *C, Matrix **M, Matrix **Minv)
 364 {
 365     Matrix *L = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
 366     for (int r = 0; r < C->NbBid; ++r)
 367         Vector_Copy(C->Ray[r]+1, L->p[r], C->Dimension);
 368     unimodular_complete(L, C->NbBid);
 369     assert(value_one_p(L->p[C->Dimension][C->Dimension]));
 370     assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
 371     Matrix_Transposition(L);
 372     assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
 373
 374     *M = Matrix_Alloc(C->Dimension+1, C->Dimension-C->NbBid+1);
 375     for (int i = 0; i < C->Dimension+1; ++i)
 376         Vector_Copy(L->p[i]+C->NbBid, (*M)->p[i], C->Dimension-C->NbBid+1);
 377
 378     Matrix *Linv = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
 379     int ok = Matrix_Inverse(L, Linv);
 380     assert(ok);
 381     Matrix_Free(L);
 382
 383     *Minv = Matrix_Alloc(C->Dimension-C->NbBid+1, C->Dimension+1);
 384     for (int i = C->NbBid; i < C->Dimension+1; ++i)
 385         Vector_AntiScale(Linv->p[i], (*Minv)->p[i-C->NbBid],
 386                          Linv->p[C->Dimension][C->Dimension], C->Dimension+1);
 387     Matrix_Free(Linv);
 388 }
 389
 390 static skewed_gen_fun *series(Polyhedron *P, Polyhedron* C,
 391                                 barvinok_options *options)
 392 {
 393     Polyhedron *C1, *C2;
 394     gen_fun *gf;
 395     Matrix *inv = NULL;
 396     Matrix *eq = NULL;
 397     Matrix *div = NULL;
 398     Polyhedron *PT = P;
 399
 400     /* Compute true context */
 401     C1 = Polyhedron_Project(P, C->Dimension);
 402     C2 = DomainIntersection(C, C1, options->MaxRays);
 403     Polyhedron_Free(C1);
 404
 405     POL_ENSURE_VERTICES(C2);
 406     if (C2->NbBid != 0) {
 407         Polyhedron *T;
 408         Matrix *M, *Minv, *M2;
 409         Matrix *CP;
 410         if (C2->NbEq || P->NbEq) {
 411             /* We remove all equalities to be sure all lines are unit vectors */
 412             Polyhedron *CT = C2;
 413             remove_all_equalities(&PT, &CT, &CP, NULL, C2->Dimension,
 414                                   options->MaxRays);
 415             if (CT != C2) {
 416                 Polyhedron_Free(C2);
 417                 C2 = CT;
 418             }
 419             if (CP) {
 420                 inv = left_inverse(CP, &eq);
 421                 Matrix_Free(CP);
 422
 423                 int d = 0;
 424                 Value tmp;
 425                 value_init(tmp);
 426                 div = Matrix_Alloc(inv->NbRows-1, inv->NbColumns+1);
 427                 for (int i = 0; i < inv->NbRows-1; ++i) {
 428                     Vector_Gcd(inv->p[i], inv->NbColumns, &tmp);
 429                     if (mpz_divisible_p(tmp,
 430                                         inv->p[inv->NbRows-1][inv->NbColumns-1]))
 431                         continue;
 432                     Vector_Copy(inv->p[i], div->p[d], inv->NbColumns);
 433                     value_assign(div->p[d][inv->NbColumns],
 434                                  inv->p[inv->NbRows-1][inv->NbColumns-1]);
 435                     ++d;
 436                 }
 437                 value_clear(tmp);
 438
 439                 if (!d) {
 440                     Matrix_Free(div);
 441                     div = NULL;
 442                 } else
 443                     div->NbRows = d;
 444             }
 445         }
 446         POL_ENSURE_VERTICES(C2);
 447
 448         if (C2->NbBid) {
 449             Matrix *M, *Minv;
 450             remove_lines(C2, &M, &Minv);
 451             apply_transformation(&PT, &C2, PT != P, C2 != C, M, Minv, &inv,
 452                                  options);
 453         }
 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         Matrix_Free(Constraints);
 465         /* flip rows of Q */
 466         for (int i = 0; i < C2->Dimension/2; ++i)
 467             Vector_Exchange(Q->p[i], Q->p[C2->Dimension-1-i], C2->Dimension);
 468         Matrix_Free(H);
 469         Matrix_Free(U);
 470         Matrix *M = Matrix_Alloc(C2->Dimension+1, C2->Dimension+1);
 471         U = Matrix_Copy(Q);
 472         int ok = Matrix_Inverse(U, M);
 473         assert(ok);
 474         Matrix_Free(U);
 475
 476         apply_transformation(&PT, &C2, PT != P, C2 != C, M, Q, &inv, options);
 477     }
 478     gf = barvinok_series_with_options(PT, C2, options);
 479     Polyhedron_Free(C2);
 480     if (PT != P)
 481         Polyhedron_Free(PT);
 482     return new skewed_gen_fun(gf, inv, eq, div);
 483 }
 484
 485 int main(int argc, char **argv)
 486 {
 487     Polyhedron *A, *C;
 488     Matrix *M;
 489     evalue *EP = NULL;
 490     skewed_gen_fun *gf = NULL;
 491     const char **param_name;
 492     int print_solution = 1;
 493     int result = 0;
 494     struct arguments options;
 495     static struct argp_child argp_children[] = {
 496         { &barvinok_argp,       0,      0,              0 },
 497         { &verify_argp,         0,      "verification",         BV_GRP_LAST+1 },
 498         { &convert_argp,        0,      "output conversion",    BV_GRP_LAST+2 },
 499         { 0 }
 500     };
 501     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 502     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 503
 504     options.verify.barvinok = bv_options;
 505     set_program_name(argv[0]);
 506     argp_parse(&argp, argc, argv, 0, 0, &options);
 507
 508     M = Matrix_Read();
 509     assert(M);
 510     A = Constraints2Polyhedron(M, bv_options->MaxRays);
 511     Matrix_Free(M);
 512     M = Matrix_Read();
 513     assert(M);
 514     C = Constraints2Polyhedron(M, bv_options->MaxRays);
 515     Matrix_Free(M);
 516     assert(A->Dimension >= C->Dimension);
 517     param_name = Read_ParamNames(stdin, C->Dimension);
 518
 519     if (options.verify.verify) {
 520         verify_options_set_range(&options.verify, A->Dimension);
 521         if (!bv_options->verbose)
 522             print_solution = 0;
 523     }
 524
 525     if (print_solution && bv_options->verbose) {
 526         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 527         Polyhedron_Print(stdout, P_VALUE_FMT, C);
 528     }
 529
 530     if (options.series) {
 531         gf = series(A, C, bv_options);
 532         if (print_solution) {
 533             gf->print(cout, C->Dimension, param_name);
 534             puts("");
 535         }
 536         if (options.function) {
 537             EP = *gf;
 538             if (print_solution)
 539                 print_evalue(stdout, EP, param_name);
 540         }
 541     } else {
 542         EP = barvinok_enumerate_with_options(A, C, bv_options);
 543         assert(EP);
 544         if (evalue_convert(EP, &options.convert, bv_options->verbose,
 545                            C->Dimension, param_name))
 546             print_solution = 0;
 547         if (options.size)
 548             printf("\nSize: %d\n", evalue_size(EP));
 549         if (print_solution)
 550             print_evalue(stdout, EP, param_name);
 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         evalue_free(EP);
 562
 563     if (options.verify.barvinok->print_stats)
 564         barvinok_stats_print(options.verify.barvinok->stats, stdout);
 565
 566     Free_ParamNames(param_name, C->Dimension);
 567     Polyhedron_Free(A);
 568     Polyhedron_Free(C);
 569     barvinok_options_free(bv_options);
 570     return result;
 571 }