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