barvinok_enumerate.cc

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