barvinok_enumerate_e.cc

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