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 "error.h"
   8 #include "config.h"
   9 #ifdef HAVE_OMEGA
  10 #include "omega/convert.h"
  11 #endif
  12 #include "verify.h"
  13 #include "verif_ehrhart.h"
  14 #include "evalue_convert.h"
  15
  16 /* The input of this example program is a polytope in combined
  17  * data and parameter space followed by two lines indicating
  18  * the number of existential variables and parameters respectively.
  19  * The first lines starts with "E ", followed by a number.
  20  * The second lines starts with "P ", followed by a number.
  21  * These two lines are (optionally) followed by the names of the parameters.
  22  * The polytope is in PolyLib notation.
  23  */
  24
  25 struct argp_option argp_options[] = {
  26     { "omega",              'o',    0,      0 },
  27     { "pip",                'p',    0,      0 },
  28     { "series",             's',    0,      0 },
  29     { "scarf",              'S',    0,      0 },
  30     { "verbose",            'v' },
  31     { 0 }
  32 };
  33
  34 struct arguments {
  35     struct verify_options    verify;
  36     struct convert_options   convert;
  37     int omega;
  38     int pip;
  39     int scarf;
  40     int series;
  41     int verbose;
  42 };
  43
  44 error_t parse_opt(int key, char *arg, struct argp_state *state)
  45 {
  46     struct arguments *arguments = (struct arguments *)(state->input);
  47
  48     switch (key) {
  49     case ARGP_KEY_INIT:
  50         state->child_inputs[0] = arguments->verify.barvinok;
  51         state->child_inputs[1] = &arguments->verify;
  52         state->child_inputs[2] = &arguments->convert;
  53         break;
  54     case 's':
  55         arguments->series = 1;
  56         break;
  57     case 'S':
  58         arguments->scarf = 1;
  59         break;
  60     case 'o':
  61 #ifdef HAVE_OMEGA
  62         arguments->omega = 1;
  63 #else
  64         error(0, 0, "--omega option not supported");
  65 #endif
  66         break;
  67     case 'p':
  68 #ifdef HAVE_PIPLIB
  69         arguments->pip = 1;
  70 #else
  71         error(0, 0, "--pip option not supported");
  72 #endif
  73         break;
  74     case 'v':
  75         arguments->verbose = 1;
  76         break;
  77     default:
  78         return ARGP_ERR_UNKNOWN;
  79     }
  80     return 0;
  81 }
  82
  83 #ifdef HAVE_OMEGA
  84
  85 Polyhedron *Omega_simplify(Polyhedron *P,
  86                             unsigned exist, unsigned nparam, char **parms)
  87 {
  88     varvector varv;
  89     varvector paramv;
  90     Relation r = Polyhedron2relation(P, exist, nparam, parms);
  91     Polyhedron_Free(P);
  92     return relation2Domain(r, varv, paramv);
  93 }
  94 #else
  95 Polyhedron *Omega_simplify(Polyhedron *P,
  96                             unsigned exist, unsigned nparam, char **parms)
  97 {
  98     return P;
  99 }
 100 #endif
 101
 102 static void verify_results(Polyhedron *P, evalue *EP, int exist, int nparam,
 103                            verify_options *options);
 104
 105 int main(int argc, char **argv)
 106 {
 107     Polyhedron *A;
 108     Matrix *MA;
 109     char **param_name;
 110     int exist, nparam, nvar;
 111     char s[128];
 112     evalue *EP;
 113     int print_solution = 1;
 114     struct arguments arguments;
 115     static struct argp_child argp_children[] = {
 116         { &barvinok_argp,       0,      0,                      0 },
 117         { &verify_argp,         0,      "verification",         1 },
 118         { &convert_argp,        0,      "output conversion",    2 },
 119         { 0 }
 120     };
 121     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 122     struct barvinok_options *options = barvinok_options_new_with_defaults();
 123
 124     arguments.verify.barvinok = options;
 125     arguments.omega = 0;
 126     arguments.pip = 0;
 127     arguments.scarf = 0;
 128     arguments.series = 0;
 129     arguments.verbose = 0;
 130
 131     argp_parse(&argp, argc, argv, 0, 0, &arguments);
 132
 133     if (arguments.series && !arguments.scarf) {
 134         fprintf(stderr,
 135                 "--series currently only available if --scarf is specified\n");
 136         exit(1);
 137     }
 138
 139     MA = Matrix_Read();
 140     A = Constraints2Polyhedron(MA, options->MaxRays);
 141     Matrix_Free(MA);
 142
 143     fgets(s, 128, stdin);
 144     while ((*s=='#') || (sscanf(s, "E %d", &exist)<1))
 145         fgets(s, 128, stdin);
 146
 147     fgets(s, 128, stdin);
 148     while ((*s=='#') || (sscanf(s, "P %d", &nparam)<1))
 149         fgets(s, 128, stdin);
 150
 151     /******* Read the options: initialize Min and Max ********/
 152
 153     if (arguments.verify.verify) {
 154         verify_options_set_range(&arguments.verify, A->Dimension);
 155         if (!arguments.verbose)
 156             print_solution = 0;
 157     }
 158
 159     if (print_solution && arguments.verbose) {
 160         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 161         printf("exist: %d, nparam: %d\n", exist, nparam);
 162     }
 163     param_name = Read_ParamNames(stdin, nparam);
 164     nvar = A->Dimension - exist - nparam;
 165     if (arguments.omega) {
 166         A = Omega_simplify(A, exist, nparam, param_name);
 167         assert(!A->next);
 168         exist = A->Dimension - nvar - nparam;
 169     }
 170     if (arguments.series) {
 171         gen_fun *gf;
 172         assert(arguments.scarf);
 173         gf = barvinok_enumerate_scarf_series(A, exist, nparam, options);
 174         if (print_solution) {
 175             gf->print(std::cout, nparam, param_name);
 176             puts("");
 177         }
 178         delete gf;
 179     } else {
 180         if (arguments.scarf)
 181             EP = barvinok_enumerate_scarf(A, exist, nparam, options);
 182         else if (arguments.pip && exist > 0)
 183             EP = barvinok_enumerate_pip_with_options(A, exist, nparam, options);
 184         else
 185             EP = barvinok_enumerate_e_with_options(A, exist, nparam, options);
 186         reduce_evalue(EP);
 187         evalue_convert(EP, &arguments.convert, nparam,
 188                        arguments.verbose ? param_name : NULL);
 189         if (print_solution && !arguments.verbose)
 190             print_evalue(stdout, EP, param_name);
 191         if (arguments.verify.verify) {
 192             arguments.verify.params = param_name;
 193             verify_results(A, EP, exist, nparam, &arguments.verify);
 194         }
 195         free_evalue_refs(EP);
 196         free(EP);
 197     }
 198     Free_ParamNames(param_name, nparam);
 199     Polyhedron_Free(A);
 200     barvinok_options_free(options);
 201     return 0;
 202 }
 203
 204 void verify_results(Polyhedron *P, evalue *EP, int exist, int nparam,
 205                     verify_options *options)
 206 {
 207     int i;
 208     int res;
 209     Value *p, tmp;
 210     Polyhedron *S, *CS;
 211     unsigned MaxRays = options->barvinok->MaxRays;
 212     Polyhedron *C = NULL;
 213     value_init(tmp);
 214
 215     p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 216     for(i=0;i<=P->Dimension;i++) {
 217       value_init(p[i]);
 218       value_set_si(p[i],0);
 219     }
 220     value_init(p[i]);
 221     value_set_si(p[i],1);
 222
 223     CS = check_poly_context_scan(P, &C, nparam, options);
 224
 225     /* S = scanning list of polyhedra */
 226     S = Polyhedron_Scan(P, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 227
 228     check_poly_init(C, options);
 229
 230     /******* CHECK NOW *********/
 231     res = 0;
 232     if (S && !check_poly_EP(S, CS, EP, exist, nparam, 0, p, options)) {
 233       fprintf(stderr,"Check failed !\n");
 234       res = -1;
 235     }
 236
 237     if (!options->print_all)
 238         printf( "\n" );
 239
 240     for(i=0;i<=(P->Dimension+1);i++)
 241       value_clear(p[i]);
 242     free(p);
 243     value_clear(tmp);
 244     Domain_Free(S);
 245     Polyhedron_Free(C);
 246     if (CS)
 247         Domain_Free(CS);
 248 }