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 "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 #define PRINT_STATS         (BV_OPT_LAST+1)
  30
  31 struct argp_option argp_options[] = {
  32     { "size",      'S' },
  33     { "series",    's', 0, 0, "compute rational generating function" },
  34     { "explicit",  'e', 0, 0, "convert rgf to psp" },
  35     { "print-stats",        PRINT_STATS,  0,    0 },
  36     { 0 }
  37 };
  38
  39 struct arguments {
  40     int size;
  41     int series;
  42     int function;
  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->print_stats = 0;
  61         break;
  62     case PRINT_STATS:
  63         options->print_stats = 1;
  64         break;
  65     case 'S':
  66         options->size = 1;
  67         break;
  68     case 'e':
  69         options->function = 1;
  70         /* fall through */
  71     case 's':
  72         options->series = 1;
  73         break;
  74     default:
  75         return ARGP_ERR_UNKNOWN;
  76     }
  77     return 0;
  78 }
  79
  80 static int verify(Polyhedron *P, Polyhedron *C, evalue *EP, skewed_gen_fun *gf,
  81                    arguments *options)
  82 {
  83     Polyhedron *CS, *S;
  84     Vector *p;
  85     int result = 0;
  86
  87     CS = check_poly_context_scan(P, &C, C->Dimension, &options->verify);
  88
  89     p = Vector_Alloc(P->Dimension+2);
  90     value_set_si(p->p[P->Dimension+1], 1);
  91
  92     /* S = scanning list of polyhedra */
  93     S = Polyhedron_Scan(P, C, options->verify.barvinok->MaxRays);
  94
  95     check_poly_init(C, &options->verify);
  96
  97     /******* CHECK NOW *********/
  98     if (S) {
  99         if (!options->series || options->function) {
 100             if (!check_poly_EP(S, CS, EP, 0, C->Dimension, 0, p->p,
 101                                 &options->verify))
 102                 result = -1;
 103         } else {
 104             if (!check_poly_gf(S, CS, gf, 0, C->Dimension, 0, p->p,
 105                                 &options->verify))
 106                 result = -1;
 107         }
 108         Domain_Free(S);
 109     }
 110
 111     if (result == -1)
 112         fprintf(stderr,"Check failed !\n");
 113
 114     if (!options->verify.print_all)
 115         printf( "\n" );
 116
 117     Vector_Free(p);
 118     if (CS) {
 119         Domain_Free(CS);
 120         Domain_Free(C);
 121     }
 122
 123     return result;
 124 }
 125
 126 /* frees M and Minv */
 127 static void apply_transformation(Polyhedron **P, Polyhedron **C,
 128                                  bool free_P, bool free_C,
 129                                  Matrix *M, Matrix *Minv, Matrix **inv,
 130                                  barvinok_options *options)
 131 {
 132     Polyhedron *T;
 133     Matrix *M2;
 134
 135     M2 = align_matrix(M, (*P)->Dimension + 1);
 136     T = *P;
 137     *P = Polyhedron_Preimage(*P, M2, options->MaxRays);
 138     if (free_P)
 139         Polyhedron_Free(T);
 140     Matrix_Free(M2);
 141
 142     T = *C;
 143     *C = Polyhedron_Preimage(*C, M, options->MaxRays);
 144     if (free_C)
 145         Polyhedron_Free(T);
 146
 147     Matrix_Free(M);
 148
 149     if (*inv) {
 150         Matrix *T = *inv;
 151         *inv = Matrix_Alloc(Minv->NbRows, T->NbColumns);
 152         Matrix_Product(Minv, T, *inv);
 153         Matrix_Free(T);
 154         Matrix_Free(Minv);
 155     } else
 156         *inv = Minv;
 157 }
 158
 159 /* Since we have "compressed" the parameters (in case there were
 160  * any equalities), the result is independent of the coordinates in the
 161  * coordinate subspace spanned by the lines.  We can therefore assume
 162  * these coordinates are zero and compute the inverse image of the map
 163  * from a lower dimensional space that adds zeros in the appropriate
 164  * places.
 165  */
 166 static void remove_lines(Polyhedron *C, Matrix **M, Matrix **Minv)
 167 {
 168     Matrix *L = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
 169     for (int r = 0; r < C->NbBid; ++r)
 170         Vector_Copy(C->Ray[r]+1, L->p[r], C->Dimension);
 171     unimodular_complete(L, C->NbBid);
 172     assert(value_one_p(L->p[C->Dimension][C->Dimension]));
 173     assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
 174     Matrix_Transposition(L);
 175     assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
 176
 177     *M = Matrix_Alloc(C->Dimension+1, C->Dimension-C->NbBid+1);
 178     for (int i = 0; i < C->Dimension+1; ++i)
 179         Vector_Copy(L->p[i]+C->NbBid, (*M)->p[i], C->Dimension-C->NbBid+1);
 180
 181     Matrix *Linv = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
 182     int ok = Matrix_Inverse(L, Linv);
 183     assert(ok);
 184     Matrix_Free(L);
 185
 186     *Minv = Matrix_Alloc(C->Dimension-C->NbBid+1, C->Dimension+1);
 187     for (int i = C->NbBid; i < C->Dimension+1; ++i)
 188         Vector_AntiScale(Linv->p[i], (*Minv)->p[i-C->NbBid],
 189                          Linv->p[C->Dimension][C->Dimension], C->Dimension+1);
 190     Matrix_Free(Linv);
 191 }
 192
 193 static skewed_gen_fun *series(Polyhedron *P, Polyhedron* C,
 194                                 barvinok_options *options)
 195 {
 196     Polyhedron *C1, *C2;
 197     gen_fun *gf;
 198     Matrix *inv = NULL;
 199     Matrix *eq = NULL;
 200     Matrix *div = NULL;
 201     Polyhedron *PT = P;
 202
 203     /* Compute true context */
 204     C1 = Polyhedron_Project(P, C->Dimension);
 205     C2 = DomainIntersection(C, C1, options->MaxRays);
 206     Polyhedron_Free(C1);
 207
 208     POL_ENSURE_VERTICES(C2);
 209     if (C2->NbBid != 0) {
 210         Polyhedron *T;
 211         Matrix *M, *Minv, *M2;
 212         Matrix *CP;
 213         if (C2->NbEq || P->NbEq) {
 214             /* We remove all equalities to be sure all lines are unit vectors */
 215             Polyhedron *CT = C2;
 216             remove_all_equalities(&PT, &CT, &CP, NULL, C2->Dimension,
 217                                   options->MaxRays);
 218             if (CT != C2) {
 219                 Polyhedron_Free(C2);
 220                 C2 = CT;
 221             }
 222             if (CP) {
 223                 inv = left_inverse(CP, &eq);
 224                 Matrix_Free(CP);
 225
 226                 int d = 0;
 227                 Value tmp;
 228                 value_init(tmp);
 229                 div = Matrix_Alloc(inv->NbRows-1, inv->NbColumns+1);
 230                 for (int i = 0; i < inv->NbRows-1; ++i) {
 231                     Vector_Gcd(inv->p[i], inv->NbColumns, &tmp);
 232                     if (mpz_divisible_p(tmp,
 233                                         inv->p[inv->NbRows-1][inv->NbColumns-1]))
 234                         continue;
 235                     Vector_Copy(inv->p[i], div->p[d], inv->NbColumns);
 236                     value_assign(div->p[d][inv->NbColumns],
 237                                  inv->p[inv->NbRows-1][inv->NbColumns-1]);
 238                     ++d;
 239                 }
 240                 value_clear(tmp);
 241
 242                 if (!d) {
 243                     Matrix_Free(div);
 244                     div = NULL;
 245                 } else
 246                     div->NbRows = d;
 247             }
 248         }
 249         POL_ENSURE_VERTICES(C2);
 250
 251         if (C2->NbBid) {
 252             Matrix *M, *Minv;
 253             remove_lines(C2, &M, &Minv);
 254             apply_transformation(&PT, &C2, PT != P, C2 != C, M, Minv, &inv,
 255                                  options);
 256         }
 257     }
 258     POL_ENSURE_VERTICES(C2);
 259     if (!Polyhedron_has_revlex_positive_rays(C2, C2->Dimension)) {
 260         Polyhedron *T;
 261         Matrix *Constraints;
 262         Matrix *H, *Q, *U;
 263         Constraints = Matrix_Alloc(C2->NbConstraints, C2->Dimension+1);
 264         for (int i = 0; i < C2->NbConstraints; ++i)
 265             Vector_Copy(C2->Constraint[i]+1, Constraints->p[i], C2->Dimension);
 266         left_hermite(Constraints, &H, &Q, &U);
 267         Matrix_Free(Constraints);
 268         /* flip rows of Q */
 269         for (int i = 0; i < C2->Dimension/2; ++i)
 270             Vector_Exchange(Q->p[i], Q->p[C2->Dimension-1-i], C2->Dimension);
 271         Matrix_Free(H);
 272         Matrix_Free(U);
 273         Matrix *M = Matrix_Alloc(C2->Dimension+1, C2->Dimension+1);
 274         U = Matrix_Copy(Q);
 275         int ok = Matrix_Inverse(U, M);
 276         assert(ok);
 277         Matrix_Free(U);
 278
 279         apply_transformation(&PT, &C2, PT != P, C2 != C, M, Q, &inv, options);
 280     }
 281     gf = barvinok_series_with_options(PT, C2, options);
 282     Polyhedron_Free(C2);
 283     if (PT != P)
 284         Polyhedron_Free(PT);
 285     return new skewed_gen_fun(gf, inv, eq, div);
 286 }
 287
 288 int main(int argc, char **argv)
 289 {
 290     Polyhedron *A, *C;
 291     Matrix *M;
 292     evalue *EP = NULL;
 293     skewed_gen_fun *gf = NULL;
 294     char **param_name;
 295     int print_solution = 1;
 296     int result = 0;
 297     struct arguments options;
 298     static struct argp_child argp_children[] = {
 299         { &barvinok_argp,       0,      0,              0 },
 300         { &verify_argp,         0,      "verification",         BV_GRP_LAST+1 },
 301         { &convert_argp,        0,      "output conversion",    BV_GRP_LAST+2 },
 302         { 0 }
 303     };
 304     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 305     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 306
 307     options.verify.barvinok = bv_options;
 308     set_program_name(argv[0]);
 309     argp_parse(&argp, argc, argv, 0, 0, &options);
 310
 311     M = Matrix_Read();
 312     assert(M);
 313     A = Constraints2Polyhedron(M, bv_options->MaxRays);
 314     Matrix_Free(M);
 315     M = Matrix_Read();
 316     assert(M);
 317     C = Constraints2Polyhedron(M, bv_options->MaxRays);
 318     Matrix_Free(M);
 319     param_name = Read_ParamNames(stdin, C->Dimension);
 320
 321     if (options.verify.verify) {
 322         verify_options_set_range(&options.verify, A->Dimension);
 323         if (!bv_options->verbose)
 324             print_solution = 0;
 325     }
 326
 327     if (print_solution && bv_options->verbose) {
 328         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 329         Polyhedron_Print(stdout, P_VALUE_FMT, C);
 330     }
 331
 332     if (options.series) {
 333         gf = series(A, C, bv_options);
 334         if (print_solution) {
 335             gf->print(cout, C->Dimension, param_name);
 336             puts("");
 337         }
 338         if (options.function) {
 339             EP = *gf;
 340             if (print_solution)
 341                 print_evalue(stdout, EP, param_name);
 342         }
 343     } else {
 344         EP = barvinok_enumerate_with_options(A, C, bv_options);
 345         if (evalue_convert(EP, &options.convert, bv_options->verbose,
 346                            C->Dimension, param_name))
 347             print_solution = 0;
 348         if (options.size)
 349             printf("\nSize: %d\n", evalue_size(EP));
 350         if (print_solution)
 351             print_evalue(stdout, EP, param_name);
 352     }
 353
 354     if (options.verify.verify) {
 355         options.verify.params = param_name;
 356         result = verify(A, C, EP, gf, &options);
 357     }
 358
 359     if (gf)
 360         delete gf;
 361     if (EP)
 362         evalue_free(EP);
 363
 364     if (options.print_stats)
 365         barvinok_stats_print(options.verify.barvinok->stats, stdout);
 366
 367     Free_ParamNames(param_name, C->Dimension);
 368     Polyhedron_Free(A);
 369     Polyhedron_Free(C);
 370     barvinok_options_free(bv_options);
 371     return result;
 372 }