barvinok_enumerate_e.cc

   1 #include <unistd.h>
   2 #include <stdlib.h>
   3 #include <assert.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/barvinok.h>
   6 #include "argp.h"
   7 #include "progname.h"
   8 #include "error.h"
   9 #include "config.h"
  10 #ifdef HAVE_OMEGA
  11 #include "omega/convert.h"
  12 #endif
  13 #ifdef USE_PARKER
  14 #include "parker/count_solutions.h"
  15 #endif
  16 #include "skewed_genfun.h"
  17 #include "verify.h"
  18 #include "verif_ehrhart.h"
  19 #include "verify_series.h"
  20 #include "evalue_convert.h"
  21 #include "normalization.h"
  22
  23 /* The input of this example program is a polytope in combined
  24  * data and parameter space followed by two lines indicating
  25  * the number of existential variables and parameters respectively.
  26  * The first lines starts with "E ", followed by a number.
  27  * The second lines starts with "P ", followed by a number.
  28  * These two lines are (optionally) followed by the names of the parameters.
  29  * The polytope is in PolyLib notation.
  30  */
  31
  32 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
  33
  34 struct argp_option argp_options[] = {
  35     { "omega",              'o',    0,      0 },
  36     { "parker",             'P',    0,      0 },
  37     { "pip",                'p',    0,      0 },
  38     { "series",             's',    0,      0 },
  39     { "series",             's', 0, 0, "compute rational generating function" },
  40     { "explicit",           'e', 0, 0, "convert rgf to psp" },
  41     { "scarf",              'S',    0,      0 },
  42     { 0 }
  43 };
  44
  45 struct arguments {
  46     struct verify_options    verify;
  47     struct convert_options   convert;
  48     int omega;
  49     int parker;
  50     int pip;
  51     int scarf;
  52     int series;
  53     int function;
  54 };
  55
  56 error_t parse_opt(int key, char *arg, struct argp_state *state)
  57 {
  58     struct arguments *arguments = (struct arguments *)(state->input);
  59
  60     switch (key) {
  61     case ARGP_KEY_INIT:
  62         state->child_inputs[0] = arguments->verify.barvinok;
  63         state->child_inputs[1] = &arguments->verify;
  64         state->child_inputs[2] = &arguments->convert;
  65         break;
  66     case 'e':
  67         arguments->function = 1;
  68         /* fall through */
  69     case 's':
  70         arguments->series = 1;
  71         break;
  72     case 'S':
  73         arguments->scarf = 1;
  74         break;
  75     case 'o':
  76 #ifdef HAVE_OMEGA
  77         arguments->omega = 1;
  78 #else
  79         error(0, 0, "--omega option not supported");
  80 #endif
  81         break;
  82     case 'P':
  83 #ifdef USE_PARKER
  84         arguments->parker = 1;
  85 #else
  86         error(0, 0, "--parker option not supported");
  87 #endif
  88         break;
  89     case 'p':
  90         arguments->pip = 1;
  91         break;
  92     default:
  93         return ARGP_ERR_UNKNOWN;
  94     }
  95     return 0;
  96 }
  97
  98 #ifdef HAVE_OMEGA
  99
 100 Polyhedron *Omega_simplify(Polyhedron *P,
 101                             unsigned exist, unsigned nparam, char **parms)
 102 {
 103     varvector varv;
 104     varvector paramv;
 105     Relation r = Polyhedron2relation(P, exist, nparam, parms);
 106     Polyhedron_Free(P);
 107     return relation2Domain(r, varv, paramv);
 108 }
 109 #else
 110 Polyhedron *Omega_simplify(Polyhedron *P,
 111                             unsigned exist, unsigned nparam, char **parms)
 112 {
 113     return P;
 114 }
 115 #endif
 116
 117 #ifdef USE_PARKER
 118 /*
 119  * Use parker's method to compute the number of integer points in P.
 120  * Since this method assumes all variables are non-negative,
 121  * we have to transform the input polytope first.
 122  */
 123 evalue *barvinok_enumerate_parker(Polyhedron *P,
 124                                         unsigned exist, unsigned nparam,
 125                                         unsigned MaxRays)
 126 {
 127     Polyhedron *R;
 128     evalue *res;
 129
 130     assert(nparam == 0);
 131     R = skew_to_positive_orthant(P, P->Dimension-exist, MaxRays);
 132     Relation r = Polyhedron2relation(R, exist, 0, NULL);
 133     Polyhedron_Free(R);
 134     double d = count_solutions(r);
 135     ALLOC(evalue, res);
 136     value_init(res->d);
 137     value_set_si(res->d, 0);
 138     res->x.p = new_enode(partition, 2, 0);
 139     EVALUE_SET_DOMAIN(res->x.p->arr[0], Universe_Polyhedron(0));
 140     value_set_si(res->x.p->arr[1].d, 1);
 141     value_init(res->x.p->arr[1].x.n);
 142     value_set_double(res->x.p->arr[1].x.n, d);
 143     return res;
 144 }
 145 #else
 146 evalue *barvinok_enumerate_parker(Polyhedron *P,
 147                                         unsigned exist, unsigned nparam,
 148                                         unsigned MaxRays)
 149 {
 150     assert(0);
 151 }
 152 #endif
 153
 154 static void verify_results(Polyhedron *P, evalue *EP, gen_fun *gf,
 155                            int exist, int nparam,
 156                            arguments *options);
 157
 158 int main(int argc, char **argv)
 159 {
 160     Polyhedron *A;
 161     Matrix *MA;
 162     char **param_name;
 163     int exist, nparam, nvar;
 164     char s[128];
 165     evalue *EP = NULL;
 166     gen_fun *gf = NULL;
 167     int print_solution = 1;
 168     struct arguments arguments;
 169     static struct argp_child argp_children[] = {
 170         { &barvinok_argp,       0,      0,                      0 },
 171         { &verify_argp,         0,      "verification",         BV_GRP_LAST+1 },
 172         { &convert_argp,        0,      "output conversion",    BV_GRP_LAST+2 },
 173         { 0 }
 174     };
 175     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 176     struct barvinok_options *options = barvinok_options_new_with_defaults();
 177
 178     arguments.verify.barvinok = options;
 179     arguments.omega = 0;
 180     arguments.parker = 0;
 181     arguments.pip = 0;
 182     arguments.scarf = 0;
 183     arguments.series = 0;
 184     arguments.function = 0;
 185
 186     set_program_name(argv[0]);
 187     argp_parse(&argp, argc, argv, 0, 0, &arguments);
 188
 189     MA = Matrix_Read();
 190     A = Constraints2Polyhedron(MA, options->MaxRays);
 191     Matrix_Free(MA);
 192
 193     fgets(s, 128, stdin);
 194     while ((*s=='#') || (sscanf(s, "E %d", &exist)<1))
 195         fgets(s, 128, stdin);
 196
 197     fgets(s, 128, stdin);
 198     while ((*s=='#') || (sscanf(s, "P %d", &nparam)<1))
 199         fgets(s, 128, stdin);
 200
 201     /******* Read the options: initialize Min and Max ********/
 202
 203     if (arguments.verify.verify) {
 204         verify_options_set_range(&arguments.verify, A->Dimension);
 205         if (!options->verbose)
 206             print_solution = 0;
 207     }
 208
 209     if (print_solution && options->verbose) {
 210         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 211         printf("exist: %d, nparam: %d\n", exist, nparam);
 212     }
 213     param_name = Read_ParamNames(stdin, nparam);
 214     nvar = A->Dimension - exist - nparam;
 215     if (arguments.omega) {
 216         A = Omega_simplify(A, exist, nparam, param_name);
 217         assert(!A->next);
 218         exist = A->Dimension - nvar - nparam;
 219     }
 220     if (arguments.series) {
 221         if (arguments.scarf)
 222             gf = barvinok_enumerate_scarf_series(A, exist, nparam, options);
 223         else
 224             gf = barvinok_enumerate_e_series(A, exist, nparam, options);
 225         if (print_solution) {
 226             gf->print(std::cout, nparam, param_name);
 227             puts("");
 228         }
 229         if (arguments.function) {
 230             EP = *gf;
 231             if (print_solution)
 232                 print_evalue(stdout, EP, param_name);
 233         }
 234     } else {
 235         if (arguments.parker)
 236             EP = barvinok_enumerate_parker(A, exist, nparam, options->MaxRays);
 237         else if (arguments.scarf)
 238             EP = barvinok_enumerate_scarf(A, exist, nparam, options);
 239         else if (arguments.pip && exist > 0)
 240             EP = barvinok_enumerate_pip_with_options(A, exist, nparam, options);
 241         else
 242             EP = barvinok_enumerate_e_with_options(A, exist, nparam, options);
 243         reduce_evalue(EP);
 244         if (evalue_convert(EP, &arguments.convert, options->verbose, nparam,
 245                            param_name))
 246             print_solution = 0;
 247         if (print_solution)
 248             print_evalue(stdout, EP, param_name);
 249     }
 250     if (arguments.verify.verify) {
 251         arguments.verify.params = param_name;
 252         verify_results(A, EP, gf, exist, nparam, &arguments);
 253     }
 254     if (gf)
 255         delete gf;
 256     if (EP)
 257         evalue_free(EP);
 258
 259     if (options->print_stats)
 260         barvinok_stats_print(options->stats, stdout);
 261
 262     Free_ParamNames(param_name, nparam);
 263     Polyhedron_Free(A);
 264     barvinok_options_free(options);
 265     return 0;
 266 }
 267
 268 void verify_results(Polyhedron *P, evalue *EP, gen_fun *gf,
 269                        int exist, int nparam,
 270                        arguments *options)
 271 {
 272     int i;
 273     int res = 0;
 274     Vector *p;
 275     Value tmp;
 276     Polyhedron *S, *CS;
 277     unsigned MaxRays = options->verify.barvinok->MaxRays;
 278     Polyhedron *C = NULL;
 279     value_init(tmp);
 280
 281     p = Vector_Alloc(P->Dimension+2);
 282     value_set_si(p->p[P->Dimension+1], 1);
 283
 284     CS = check_poly_context_scan(P, &C, nparam, &options->verify);
 285     if (!C)
 286         C = Universe_Polyhedron(nparam);
 287
 288     /* S = scanning list of polyhedra */
 289     S = Polyhedron_Scan(P, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 290
 291     check_poly_init(C, &options->verify);
 292
 293     /******* CHECK NOW *********/
 294     if (S) {
 295         if (!options->series || options->function) {
 296             if (!check_poly_EP(S, CS, EP, exist, nparam, 0, p->p,
 297                                 &options->verify))
 298                 res = -1;
 299         } else {
 300             skewed_gen_fun *sgf = new skewed_gen_fun(new gen_fun(gf));
 301             if (!check_poly_gf(S, CS, sgf, exist, nparam, 0, p->p,
 302                                 &options->verify))
 303                 res = -1;
 304             delete sgf;
 305         }
 306     }
 307
 308     if (!options->verify.print_all)
 309         printf( "\n" );
 310
 311     Vector_Free(p);
 312     value_clear(tmp);
 313     Domain_Free(S);
 314     Polyhedron_Free(C);
 315     if (CS)
 316         Domain_Free(CS);
 317 }