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