barvinok_enumerate.cc

   1 #include <unistd.h>
   2 #include <stdlib.h>
   3 #include <barvinok/evalue.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/barvinok.h>
   6 #include "argp.h"
   7 #include "verify.h"
   8 #include "verif_ehrhart.h"
   9
  10 /* The input of this example program is the same as that of testehrhart
  11  * in the PolyLib distribution, i.e., a polytope in combined
  12  * data and parameter space, a context polytope in parameter space
  13  * and (optionally) the names of the parameters.
  14  * Both polytopes are in PolyLib notation.
  15  */
  16
  17 struct argp_option argp_options[] = {
  18     { "convert",   'c', 0, 0, "convert fractionals to periodics" },
  19     { "floor",     'f', 0, 0, "convert fractionals to floorings" },
  20     { "size",      'S' },
  21     { "series",    's', 0, 0, "compute rational generating function" },
  22     { "explicit",  'e', 0, 0, "convert rgf to psp" },
  23     { "verbose",    'v' },
  24     { 0 }
  25 };
  26
  27 struct arguments {
  28     struct barvinok_options *barvinok;
  29     int convert;
  30     int floor;
  31     int size;
  32     int series;
  33     int function;
  34     int verbose;
  35     struct verify_options    verify;
  36 };
  37
  38 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  39 {
  40     struct arguments *options = (struct arguments*) state->input;
  41
  42     switch (key) {
  43     case ARGP_KEY_INIT:
  44         state->child_inputs[0] = options->barvinok;
  45         state->child_inputs[1] = &options->verify;
  46         options->convert = 0;
  47         options->floor = 0;
  48         options->size = 0;
  49         options->series = 0;
  50         options->function = 0;
  51         options->verbose = 0;
  52         break;
  53     case 'c':
  54         options->convert = 1;
  55         break;
  56     case 'f':
  57         options->floor = 1;
  58         break;
  59     case 'S':
  60         options->size = 1;
  61         break;
  62     case 'e':
  63         options->function = 1;
  64         /* fall through */
  65     case 's':
  66         options->series = 1;
  67         break;
  68     case 'v':
  69         options->verbose = 1;
  70         break;
  71     default:
  72         return ARGP_ERR_UNKNOWN;
  73     }
  74     return 0;
  75 }
  76
  77 static int check_series(Polyhedron *S, Polyhedron *CS, gen_fun *gf,
  78                         int nparam, int pos, Value *z, int print_all)
  79 {
  80     int k;
  81     Value c, tmp;
  82     Value LB, UB;
  83
  84     value_init(c);
  85     value_init(tmp);
  86     value_init(LB);
  87     value_init(UB);
  88
  89     if (pos == nparam) {
  90         /* Computes the coefficient */
  91         gf->coefficient(&z[S->Dimension-nparam+1], &c);
  92
  93         /* if c=0 we may be out of context. */
  94         /* scanning is useless in this case*/
  95
  96         if (print_all) {
  97             printf("EP( ");
  98             value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  99             for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 100                 printf(", ");
 101                 value_print(stdout,VALUE_FMT,z[k]);
 102             }
 103             printf(" ) = ");
 104             value_print(stdout,VALUE_FMT,c);
 105             printf(" ");
 106         }
 107
 108         /* Manually count the number of points */
 109         count_points(1,S,z,&tmp);
 110         if (print_all) {
 111             printf(", count = ");
 112             value_print(stdout, P_VALUE_FMT, tmp);
 113             printf(". ");
 114         }
 115
 116         if (value_ne(tmp,c)) {
 117             printf("\n");
 118             fflush(stdout);
 119             fprintf(stderr,"Error !\n");
 120             fprintf(stderr,"EP( ");
 121             value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
 122             for (k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 123                 fprintf(stderr,", ");
 124                 value_print(stderr,VALUE_FMT,z[k]);
 125             }
 126             fprintf(stderr," ) should be ");
 127             value_print(stderr,VALUE_FMT,tmp);
 128             fprintf(stderr,", while EP eval gives ");
 129             value_print(stderr,VALUE_FMT,c);
 130             fprintf(stderr,".\n");
 131 #ifndef DONT_BREAK_ON_ERROR
 132             value_clear(c); value_clear(tmp);
 133             return 0;
 134 #endif
 135         } else if (print_all)
 136             printf("OK.\n");
 137     } else {
 138         int ok =
 139           !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 140         assert(ok);
 141         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 142             if (!print_all) {
 143                 k = VALUE_TO_INT(tmp);
 144                 if(!pos && !(k%st)) {
 145                     printf("o");
 146                     fflush(stdout);
 147                 }
 148             }
 149             value_assign(z[pos+S->Dimension-nparam+1],tmp);
 150             if (!check_series(S, CS->next, gf, nparam, pos+1, z, print_all)) {
 151                 value_clear(c); value_clear(tmp);
 152                 value_clear(LB);
 153                 value_clear(UB);
 154                 return(0);
 155             }
 156         }
 157         value_set_si(z[pos+S->Dimension-nparam+1],0);
 158     }
 159
 160     value_clear(c);
 161     value_clear(tmp);
 162     value_clear(LB);
 163     value_clear(UB);
 164     return 1;
 165 }
 166
 167 static int verify(Polyhedron *P, Polyhedron **C, Enumeration *en, gen_fun *gf,
 168                    arguments *options)
 169 {
 170     Polyhedron *CC, *PP, *CS, *S, *U;
 171     Matrix *C1, *MM;
 172     Vector *p;
 173     int result = 0;
 174
 175     /******* Compute true context *******/
 176     CC = align_context(*C, P->Dimension, options->barvinok->MaxRays);
 177     PP = DomainIntersection(P, CC, options->barvinok->MaxRays);
 178     Domain_Free(CC);
 179     C1 = Matrix_Alloc((*C)->Dimension+1, P->Dimension+1);
 180
 181     for (int i = 0; i < C1->NbRows; i++)
 182         for (int j = 0; j < C1->NbColumns; j++)
 183             if (i == j-P->Dimension+(*C)->Dimension)
 184                 value_set_si(C1->p[i][j], 1);
 185             else
 186                 value_set_si(C1->p[i][j], 0);
 187     CC = Polyhedron_Image(PP, C1, options->barvinok->MaxRays);
 188     Matrix_Free(C1);
 189     Domain_Free(PP);
 190     Domain_Free(*C);
 191     *C = CC;
 192
 193     /* Intersect context with range */
 194     if ((*C)->Dimension > 0) {
 195         MM = Matrix_Alloc(2*(*C)->Dimension, (*C)->Dimension+2);
 196         for (int i = 0; i < (*C)->Dimension; ++i) {
 197             value_set_si(MM->p[2*i][0], 1);
 198             value_set_si(MM->p[2*i][1+i], 1);
 199             value_set_si(MM->p[2*i][1+(*C)->Dimension], -options->verify.m);
 200             value_set_si(MM->p[2*i+1][0], 1);
 201             value_set_si(MM->p[2*i+1][1+i], -1);
 202             value_set_si(MM->p[2*i+1][1+(*C)->Dimension], options->verify.M);
 203         }
 204         CC = AddConstraints(MM->p[0], 2*(*C)->Dimension, *C,
 205                             options->barvinok->MaxRays);
 206         U = Universe_Polyhedron(0);
 207         CS = Polyhedron_Scan(CC, U, options->barvinok->MaxRays);
 208         Polyhedron_Free(U);
 209         Polyhedron_Free(CC);
 210         Matrix_Free(MM);
 211     } else
 212         CS = NULL;
 213
 214     p = Vector_Alloc(P->Dimension+2);
 215     value_set_si(p->p[P->Dimension+1], 1);
 216
 217     /* S = scanning list of polyhedra */
 218     S = Polyhedron_Scan(P, *C, options->barvinok->MaxRays);
 219
 220     if (!options->verify.print_all)
 221         if ((*C)->Dimension > 0) {
 222             int d = options->verify.M - options->verify.m;
 223             if (d > 80)
 224                 st = 1+d/80;
 225             else
 226                 st=1;
 227             for (int i = options->verify.m; i <= options->verify.M; i += st)
 228                 printf(".");
 229             printf( "\r" );
 230             fflush(stdout);
 231         }
 232
 233     /******* CHECK NOW *********/
 234     if (S) {
 235         if (!options->series || options->function) {
 236             if (!check_poly(S, CS, en, (*C)->Dimension, 0, p->p,
 237                             options->verify.print_all))
 238                 result = -1;
 239         } else {
 240             if (!check_series(S, CS, gf, (*C)->Dimension, 0, p->p,
 241                               options->verify.print_all))
 242                 result = -1;
 243         }
 244         Domain_Free(S);
 245     }
 246
 247     if (result == -1)
 248         fprintf(stderr,"Check failed !\n");
 249
 250     if (!options->verify.print_all)
 251         printf( "\n" );
 252
 253     Vector_Free(p);
 254     if (CS)
 255         Domain_Free(CS);
 256
 257     return result;
 258 }
 259
 260 int main(int argc, char **argv)
 261 {
 262     Polyhedron *A, *C;
 263     Matrix *M;
 264     evalue *EP = NULL;
 265     Enumeration *en = NULL;
 266     gen_fun *gf = NULL;
 267     char **param_name;
 268     int print_solution = 1;
 269     int result = 0;
 270     struct arguments options;
 271     static struct argp_child argp_children[] = {
 272         { &barvinok_argp,       0,      0,              0 },
 273         { &verify_argp,         0,      "verification", 1 },
 274         { 0 }
 275     };
 276     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 277     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 278
 279     options.barvinok = bv_options;
 280     argp_parse(&argp, argc, argv, 0, 0, &options);
 281
 282     M = Matrix_Read();
 283     A = Constraints2Polyhedron(M, bv_options->MaxRays);
 284     Matrix_Free(M);
 285     M = Matrix_Read();
 286     C = Constraints2Polyhedron(M, bv_options->MaxRays);
 287     Matrix_Free(M);
 288     param_name = Read_ParamNames(stdin, C->Dimension);
 289
 290     if (options.verify.verify) {
 291         verify_options_set_range(&options.verify, A);
 292         if (!options.verbose)
 293             print_solution = 0;
 294     }
 295
 296     if (print_solution) {
 297         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 298         Polyhedron_Print(stdout, P_VALUE_FMT, C);
 299     }
 300
 301     if (options.series) {
 302         gf = barvinok_series_with_options(A, C, bv_options);
 303         if (print_solution) {
 304             gf->print(std::cout, C->Dimension, param_name);
 305             puts("");
 306         }
 307         if (options.function) {
 308             EP = *gf;
 309             if (print_solution)
 310                 print_evalue(stdout, EP, param_name);
 311         }
 312     } else {
 313         EP = barvinok_enumerate_with_options(A, C, bv_options);
 314         if (print_solution)
 315             print_evalue(stdout, EP, param_name);
 316         if (options.size)
 317             printf("\nSize: %d\n", evalue_size(EP));
 318         if (options.floor) {
 319             fprintf(stderr, "WARNING: floor conversion not supported\n");
 320             evalue_frac2floor2(EP, 0);
 321             print_evalue(stdout, EP, param_name);
 322         } else if (options.convert) {
 323             evalue_mod2table(EP, C->Dimension);
 324             print_evalue(stdout, EP, param_name);
 325             if (options.size)
 326                 printf("\nSize: %d\n", evalue_size(EP));
 327         }
 328     }
 329
 330     if (options.verify.verify) {
 331         if (EP) {
 332             en = partition2enumeration(EP);
 333             EP = NULL;
 334         }
 335         result = verify(A, &C, en, gf, &options);
 336     }
 337
 338     if (en)
 339         Enumeration_Free(en);
 340     if (gf)
 341         delete gf;
 342     if (EP) {
 343         free_evalue_refs(EP);
 344         free(EP);
 345     }
 346
 347     Free_ParamNames(param_name, C->Dimension);
 348     Polyhedron_Free(A);
 349     Polyhedron_Free(C);
 350     free(bv_options);
 351     return result;
 352 }