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