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