barvinok_enumerate_e.cc

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