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