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 "config.h"
   8 #ifdef HAVE_OMEGA
   9 #include "omega/convert.h"
  10 #endif
  11
  12 /* The input of this example program is a polytope in combined
  13  * data and parameter space followed by two lines indicating
  14  * the number of existential variables and parameters respectively.
  15  * The first lines starts with "E ", followed by a number.
  16  * The second lines starts with "P ", followed by a number.
  17  * These two lines are (optionally) followed by the names of the parameters.
  18  * The polytope is in PolyLib notation.
  19  */
  20
  21 struct argp_option argp_options[] = {
  22 #ifdef HAVE_OMEGA
  23     { "omega",              'o',    0,      0 },
  24 #endif
  25 #ifdef HAVE_PIPLIB
  26     { "pip",                'p',    0,      0 },
  27 #endif
  28     { "series",             's',    0,      0 },
  29     { "scarf",              'S',    0,      0 },
  30     { "convert",            'c',    0,      0 },
  31     { "floor",              'f',    0,      0 },
  32     { "range-reduction",    'R',    0,      0 },
  33     { "verify",             'T',    0,      0 },
  34     { "print-all",          'A',    0,      0 },
  35     { "min",                'm',    "int",  0 },
  36     { "max",                'M',    "int",  0 },
  37     { "range",              'r',    "int",  0 },
  38     { 0 }
  39 };
  40
  41 struct arguments {
  42     struct barvinok_options *options;
  43     int range;
  44     int convert;
  45     int omega;
  46     int pip;
  47     int scarf;
  48     int series;
  49     int floor;
  50     int verify;
  51     int print_all;
  52     int m;
  53     int M;
  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 's':
  62         arguments->series = 1;
  63         break;
  64     case 'S':
  65         arguments->scarf = 1;
  66         break;
  67     case 'o':
  68         arguments->omega = 1;
  69         break;
  70     case 'p':
  71         arguments->pip = 1;
  72         break;
  73     case 'f':
  74         arguments->floor = 1;
  75         break;
  76     case 'c':
  77         arguments->convert = 1;
  78         break;
  79     case 'R':
  80         arguments->range = 1;
  81         break;
  82     case 'T':
  83         arguments->verify = 1;
  84         break;
  85     case 'A':
  86         arguments->print_all = 1;
  87         break;
  88     case 'm':
  89         arguments->m = atoi(optarg);
  90         arguments->verify = 1;
  91         break;
  92     case 'M':
  93         arguments->M = atoi(optarg);
  94         arguments->verify = 1;
  95         break;
  96     case 'r':
  97         arguments->M = atoi(optarg);
  98         arguments->m = -arguments->M;
  99         arguments->verify = 1;
 100         break;
 101     default:
 102         return ARGP_ERR_UNKNOWN;
 103     }
 104     return 0;
 105 }
 106
 107 #ifdef HAVE_OMEGA
 108
 109 Polyhedron *Omega_simplify(Polyhedron *P,
 110                             unsigned exist, unsigned nparam, char **parms)
 111 {
 112     varvector varv;
 113     varvector paramv;
 114     Relation r = Polyhedron2relation(P, exist, nparam, parms);
 115     Polyhedron_Free(P);
 116     return relation2Domain(r, varv, paramv);
 117 }
 118 #else
 119 Polyhedron *Omega_simplify(Polyhedron *P,
 120                             unsigned exist, unsigned nparam, char **parms)
 121 {
 122     return P;
 123 }
 124 #endif
 125
 126 /* define this to continue the test after first error found */
 127 /* #define DONT_BREAK_ON_ERROR */
 128
 129 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
 130 #define RANGE 50
 131
 132 /* SRANGE : small range for evalutations */
 133 #define SRANGE 15
 134
 135 /* if dimension >= BIDDIM, use SRANGE */
 136 #define BIGDIM 5
 137
 138 /* VSRANGE : very small range for evalutations */
 139 #define VSRANGE 5
 140
 141 /* if dimension >= VBIDDIM, use VSRANGE */
 142 #define VBIGDIM 8
 143
 144 static Value min_val, max_val;
 145
 146 static char **params;
 147
 148 static int st;
 149
 150 static int check_poly(Polyhedron *S, Polyhedron *C, evalue *EP,
 151                       int exist, int nparam, int pos, Value *z, int print_all);
 152 static void verify_results(Polyhedron *P, evalue *EP, int exist, int nparam,
 153                            int m, int M, int print_all, unsigned MaxRays);
 154
 155 int main(int argc, char **argv)
 156 {
 157     Polyhedron *A;
 158     Matrix *MA;
 159     char **param_name;
 160     int exist, nparam, nvar;
 161     char s[128];
 162     evalue *EP;
 163     int print_solution = 1;
 164     int r;
 165     struct arguments arguments;
 166     static struct argp_child argp_children[] = {
 167         { &barvinok_argp,    0, 0,  0 },
 168         { 0 }
 169     };
 170     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 171     struct barvinok_options *options = barvinok_options_new_with_defaults();
 172
 173     arguments.options = options;
 174     arguments.range = 0;
 175     arguments.convert = 0;
 176     arguments.omega = 0;
 177     arguments.pip = 0;
 178     arguments.scarf = 0;
 179     arguments.series = 0;
 180     arguments.floor = 0;
 181     arguments.verify = 0;
 182     arguments.print_all = 0;
 183     arguments.m = INT_MAX;
 184     arguments.M = INT_MIN;
 185
 186     argp_parse(&argp, argc, argv, ARGP_NO_EXIT, 0, &arguments);
 187
 188     if (arguments.series && !arguments.scarf) {
 189         fprintf(stderr,
 190                 "--series currently only available if --scarf is specified\n");
 191         exit(1);
 192     }
 193
 194     MA = Matrix_Read();
 195     A = Constraints2Polyhedron(MA, options->MaxRays);
 196     Matrix_Free(MA);
 197
 198     fgets(s, 128, stdin);
 199     while ((*s=='#') || (sscanf(s, "E %d", &exist)<1))
 200         fgets(s, 128, stdin);
 201
 202     fgets(s, 128, stdin);
 203     while ((*s=='#') || (sscanf(s, "P %d", &nparam)<1))
 204         fgets(s, 128, stdin);
 205
 206     /******* Read the options: initialize Min and Max ********/
 207     if (A->Dimension >= VBIGDIM)
 208         r = VSRANGE;
 209     else if (A->Dimension >= BIGDIM)
 210         r = SRANGE;
 211     else
 212         r = RANGE;
 213     if (arguments.M == INT_MIN)
 214         arguments.M = r;
 215     if (arguments.m == INT_MAX)
 216         arguments.m = -r;
 217
 218     if (arguments.verify && arguments.m > arguments.M) {
 219         fprintf(stderr,"Nothing to do: min > max !\n");
 220         return(0);
 221     }
 222     if (arguments.verify)
 223         print_solution = 0;
 224
 225     if (print_solution) {
 226         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 227         printf("exist: %d, nparam: %d\n", exist, nparam);
 228     }
 229     param_name = Read_ParamNames(stdin, nparam);
 230     nvar = A->Dimension - exist - nparam;
 231     if (arguments.omega) {
 232         A = Omega_simplify(A, exist, nparam, param_name);
 233         assert(!A->next);
 234         exist = A->Dimension - nvar - nparam;
 235     }
 236     if (arguments.series) {
 237         gen_fun *gf;
 238         barvinok_options *options = barvinok_options_new_with_defaults();
 239         assert(arguments.scarf);
 240         gf = barvinok_enumerate_scarf_series(A, exist, nparam, options);
 241         if (print_solution) {
 242             gf->print(std::cout, nparam, param_name);
 243             puts("");
 244         }
 245         delete gf;
 246         free(options);
 247     } else {
 248         if (arguments.scarf) {
 249             barvinok_options *options = barvinok_options_new_with_defaults();
 250             EP = barvinok_enumerate_scarf(A, exist, nparam, options);
 251             free(options);
 252         } else if (arguments.pip && exist > 0)
 253             EP = barvinok_enumerate_pip(A, exist, nparam, options->MaxRays);
 254         else
 255             EP = barvinok_enumerate_e(A, exist, nparam, options->MaxRays);
 256         reduce_evalue(EP);
 257         evalue_combine(EP);
 258         if (arguments.range)
 259             evalue_range_reduction(EP);
 260         if (print_solution)
 261             print_evalue(stdout, EP, param_name);
 262         if (arguments.floor) {
 263             fprintf(stderr, "WARNING: floor conversion not supported\n");
 264             evalue_frac2floor2(EP, 0);
 265             if (print_solution)
 266                 print_evalue(stdout, EP, param_name);
 267         } else if (arguments.convert) {
 268             evalue_mod2table(EP, nparam);
 269             if (print_solution)
 270                 print_evalue(stdout, EP, param_name);
 271         }
 272         if (arguments.verify)
 273             verify_results(A, EP, exist, nparam, arguments.m, arguments.M,
 274                             arguments.print_all, options->MaxRays);
 275         free_evalue_refs(EP);
 276         free(EP);
 277     }
 278     Free_ParamNames(param_name, nparam);
 279     Polyhedron_Free(A);
 280     return 0;
 281 }
 282
 283 void verify_results(Polyhedron *P, evalue *EP, int exist, int nparam, int m, int M,
 284                     int print_all, unsigned MaxRays)
 285 {
 286     int i;
 287     int res;
 288     Value *p, tmp;
 289     Polyhedron *S;
 290     Polyhedron *C = Polyhedron_Project(P, nparam);
 291     value_init(min_val);
 292     value_init(max_val);
 293     value_set_si(min_val,m);
 294     value_set_si(max_val,M);
 295     value_init(tmp);
 296
 297     p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 298     for(i=0;i<=P->Dimension;i++) {
 299       value_init(p[i]);
 300       value_set_si(p[i],0);
 301     }
 302     value_init(p[i]);
 303     value_set_si(p[i],1);
 304
 305     /* S = scanning list of polyhedra */
 306     S = Polyhedron_Scan(P, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 307
 308     if (!print_all) {
 309         if (C->Dimension > 0) {
 310           value_subtract(tmp,max_val,min_val);
 311           if (VALUE_TO_INT(tmp) > 80)
 312             st = 1+(VALUE_TO_INT(tmp))/80;
 313           else
 314             st=1;
 315           for(i=VALUE_TO_INT(min_val);i<=VALUE_TO_INT(max_val);i+=st)
 316             printf(".");
 317           printf( "\r" );
 318           fflush(stdout);
 319         }
 320     }
 321
 322     /******* CHECK NOW *********/
 323     res = 0;
 324     if(S && !check_poly(S, C, EP, exist, nparam, 0, p, print_all)) {
 325       fprintf(stderr,"Check failed !\n");
 326       res = -1;
 327     }
 328
 329     if (!print_all)
 330         printf( "\n" );
 331
 332     for(i=0;i<=(P->Dimension+1);i++)
 333       value_clear(p[i]);
 334     free(p);
 335     value_clear(tmp);
 336     Domain_Free(S);
 337     Polyhedron_Free(C);
 338 }
 339
 340 /****************************************************/
 341 /* function check_poly :                            */
 342 /* scans the parameter space from min to max (all   */
 343 /* directions). Computes the number of points in    */
 344 /* the polytope using both methods, and compare them*/
 345 /* returns 1 on success                             */
 346 /****************************************************/
 347
 348 int check_poly(Polyhedron *S, Polyhedron *C, evalue *EP,
 349                int exist, int nparam, int pos, Value *z, int print_all)
 350 {
 351   int k;
 352   Value c,tmp;
 353
 354   value_init(c); value_init(tmp);
 355
 356   if(pos == nparam) {
 357
 358     /* Computes the ehrhart polynomial */
 359     value_set_double(c, compute_evalue(EP,&z[S->Dimension-nparam+1])+.25);
 360     /* if c=0 we may be out of context. */
 361     /* scanning is useless in this case*/
 362     if(!in_domain(C,&z[S->Dimension-nparam+1])) {
 363
 364       /* ok */ ;
 365     }
 366     else {
 367
 368     if (print_all) {
 369       printf("EP( ");
 370       value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
 371       for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 372         printf(", ");
 373         value_print(stdout,VALUE_FMT,z[k]);
 374       }
 375       printf(" ) = ");
 376       value_print(stdout,VALUE_FMT,c);
 377       printf(" ");
 378     }
 379
 380       /* Manually count the number of points */
 381       count_points_e(1, S, exist, nparam, z, &tmp);
 382     if (print_all) {
 383         printf(", count = ");
 384         value_print(stdout, P_VALUE_FMT, tmp);
 385         printf(". ");
 386     }
 387
 388       if(value_ne(tmp,c)) {
 389         printf("\n");
 390         fflush(stdout);
 391         fprintf(stderr,"Error !\n");
 392         fprintf(stderr,"EP( ");
 393         value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
 394         for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 395           fprintf(stderr,", ");
 396           value_print(stderr,VALUE_FMT,z[k]);
 397         }
 398         fprintf(stderr," ) should be ");
 399         value_print(stderr,VALUE_FMT,tmp);
 400         fprintf(stderr,", while EP eval gives ");
 401         value_print(stderr,VALUE_FMT,c);
 402         fprintf(stderr,".\n");
 403         print_evalue(stderr, EP, params);
 404 #ifndef DONT_BREAK_ON_ERROR
 405         value_clear(c); value_clear(tmp);
 406         return(0);
 407 #endif
 408       }
 409       else if (print_all)
 410         printf("OK.\n");
 411     }
 412   }
 413   else
 414     for(value_assign(tmp,min_val); value_le(tmp,max_val); value_increment(tmp,tmp)) {
 415       if (!print_all) {
 416           k = VALUE_TO_INT(tmp);
 417           if(!pos && !(k%st)) {
 418             printf("o");
 419             fflush(stdout);
 420           }
 421        }
 422
 423       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 424       if(!check_poly(S, C, EP, exist, nparam, pos+1, z, print_all)) {
 425         value_clear(c); value_clear(tmp);
 426         return(0);
 427       }
 428     }
 429   value_clear(c); value_clear(tmp);
 430   return(1);
 431 } /* check_poly */