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