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 "fdstream.h"
   7 #include "argp.h"
   8 #include "verify.h"
   9 #include "verif_ehrhart.h"
  10 #include "remove_equalities.h"
  11
  12 #undef CS   /* for Solaris 10 */
  13
  14 /* The input of this example program is the same as that of testehrhart
  15  * in the PolyLib distribution, i.e., a polytope in combined
  16  * data and parameter space, a context polytope in parameter space
  17  * and (optionally) the names of the parameters.
  18  * Both polytopes are in PolyLib notation.
  19  */
  20
  21 struct argp_option argp_options[] = {
  22     { "convert",   'c', 0, 0, "convert fractionals to periodics" },
  23     { "floor",     'f', 0, 0, "convert fractionals to floorings" },
  24     { "size",      'S' },
  25     { "series",    's', 0, 0, "compute rational generating function" },
  26     { "explicit",  'e', 0, 0, "convert rgf to psp" },
  27     { "verbose",    'v' },
  28     { 0 }
  29 };
  30
  31 struct arguments {
  32     struct barvinok_options *barvinok;
  33     int convert;
  34     int floor;
  35     int size;
  36     int series;
  37     int function;
  38     int verbose;
  39     struct verify_options    verify;
  40 };
  41
  42 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  43 {
  44     struct arguments *options = (struct arguments*) state->input;
  45
  46     switch (key) {
  47     case ARGP_KEY_INIT:
  48         state->child_inputs[0] = options->barvinok;
  49         state->child_inputs[1] = &options->verify;
  50         options->convert = 0;
  51         options->floor = 0;
  52         options->size = 0;
  53         options->series = 0;
  54         options->function = 0;
  55         options->verbose = 0;
  56         break;
  57     case 'c':
  58         options->convert = 1;
  59         break;
  60     case 'f':
  61         options->floor = 1;
  62         break;
  63     case 'S':
  64         options->size = 1;
  65         break;
  66     case 'e':
  67         options->function = 1;
  68         /* fall through */
  69     case 's':
  70         options->series = 1;
  71         break;
  72     case 'v':
  73         options->verbose = 1;
  74         break;
  75     default:
  76         return ARGP_ERR_UNKNOWN;
  77     }
  78     return 0;
  79 }
  80
  81 struct skewed_gen_fun {
  82     gen_fun *gf;
  83     /* maps original space to space in which gf is defined */
  84     Matrix  *T;
  85     /* equalities in the original space that need to be satisfied for
  86      * gf to be valid
  87      */
  88     Matrix  *eq;
  89     /* divisibilities in the original space that need to be satisfied for
  90      * gf to be valid
  91      */
  92     Matrix  *div;
  93
  94     skewed_gen_fun(gen_fun *gf, Matrix *T, Matrix *eq, Matrix *div) :
  95                     gf(gf), T(T), eq(eq), div(div) {}
  96     ~skewed_gen_fun() {
  97         if (T)
  98             Matrix_Free(T);
  99         if (eq)
 100             Matrix_Free(eq);
 101         if (div)
 102             Matrix_Free(div);
 103         delete gf;
 104     }
 105
 106     void print(FILE *out, unsigned int nparam, char **param_name) const;
 107     operator evalue *() const {
 108         assert(T == NULL && eq == NULL); /* other cases not supported for now */
 109         return *gf;
 110     }
 111     void coefficient(Value* params, Value* c, barvinok_options *options) const;
 112 };
 113
 114 void skewed_gen_fun::print(FILE *out, unsigned int nparam,
 115                             char **param_name) const
 116 {
 117     fdostream os(dup(fileno(out)));
 118     if (T) {
 119         fprintf(out, "T:\n");
 120         Matrix_Print(out, P_VALUE_FMT, T);
 121     }
 122     if (eq) {
 123         fprintf(out, "eq:\n");
 124         Matrix_Print(out, P_VALUE_FMT, eq);
 125     }
 126     if (div) {
 127         fprintf(out, "div:\n");
 128         Matrix_Print(out, P_VALUE_FMT, div);
 129     }
 130     gf->print(os, nparam, param_name);
 131 }
 132
 133 void skewed_gen_fun::coefficient(Value* params, Value* c,
 134                                  barvinok_options *options) const
 135 {
 136     if (eq) {
 137         for (int i = 0; i < eq->NbRows; ++i) {
 138             Inner_Product(eq->p[i]+1, params, eq->NbColumns-2, eq->p[i]);
 139             if (value_notzero_p(eq->p[i][0])) {
 140                 value_set_si(*c, 0);
 141                 return;
 142             }
 143         }
 144     }
 145     if (div) {
 146         Value tmp;
 147         value_init(tmp);
 148         for (int i = 0; i < div->NbRows; ++i) {
 149             Inner_Product(div->p[i], params, div->NbColumns-1, &tmp);
 150             if (!mpz_divisible_p(tmp, div->p[i][div->NbColumns-1])) {
 151                 value_set_si(*c, 0);
 152                 return;
 153             }
 154         }
 155         value_clear(tmp);
 156     }
 157
 158     ZZ coeff;
 159     if (!T)
 160         coeff = gf->coefficient(params, options);
 161     else {
 162         Vector *p2 = Vector_Alloc(T->NbRows);
 163         Matrix_Vector_Product(T, params, p2->p);
 164         if (value_notone_p(p2->p[T->NbRows-1]))
 165             Vector_AntiScale(p2->p, p2->p, p2->p[T->NbRows-1], T->NbRows);
 166         coeff = gf->coefficient(p2->p, options);
 167         Vector_Free(p2);
 168     }
 169
 170     zz2value(coeff, *c);
 171 }
 172
 173 static int check_series(Polyhedron *S, Polyhedron *CS, skewed_gen_fun *gf,
 174                         int nparam, int pos, Value *z, int print_all,
 175                         barvinok_options *options)
 176 {
 177     int k;
 178     Value c, tmp;
 179     Value LB, UB;
 180
 181     value_init(c);
 182     value_init(tmp);
 183     value_init(LB);
 184     value_init(UB);
 185
 186     if (pos == nparam) {
 187         /* Computes the coefficient */
 188         gf->coefficient(&z[S->Dimension-nparam+1], &c, options);
 189
 190         /* if c=0 we may be out of context. */
 191         /* scanning is useless in this case*/
 192
 193         if (print_all) {
 194             printf("EP( ");
 195             value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
 196             for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 197                 printf(", ");
 198                 value_print(stdout,VALUE_FMT,z[k]);
 199             }
 200             printf(" ) = ");
 201             value_print(stdout,VALUE_FMT,c);
 202             printf(" ");
 203         }
 204
 205         /* Manually count the number of points */
 206         count_points(1,S,z,&tmp);
 207         if (print_all) {
 208             printf(", count = ");
 209             value_print(stdout, P_VALUE_FMT, tmp);
 210             printf(". ");
 211         }
 212
 213         if (value_ne(tmp,c)) {
 214             printf("\n");
 215             fflush(stdout);
 216             fprintf(stderr,"Error !\n");
 217             fprintf(stderr,"EP( ");
 218             value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
 219             for (k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 220                 fprintf(stderr,", ");
 221                 value_print(stderr,VALUE_FMT,z[k]);
 222             }
 223             fprintf(stderr," ) should be ");
 224             value_print(stderr,VALUE_FMT,tmp);
 225             fprintf(stderr,", while EP eval gives ");
 226             value_print(stderr,VALUE_FMT,c);
 227             fprintf(stderr,".\n");
 228 #ifndef DONT_BREAK_ON_ERROR
 229             value_clear(c); value_clear(tmp);
 230             return 0;
 231 #endif
 232         } else if (print_all)
 233             printf("OK.\n");
 234     } else {
 235         int ok =
 236           !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 237         assert(ok);
 238         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 239             if (!print_all) {
 240                 k = VALUE_TO_INT(tmp);
 241                 if(!pos && !(k%st)) {
 242                     printf("o");
 243                     fflush(stdout);
 244                 }
 245             }
 246             value_assign(z[pos+S->Dimension-nparam+1],tmp);
 247             if (!check_series(S, CS->next, gf, nparam, pos+1, z, print_all,
 248                               options)) {
 249                 value_clear(c); value_clear(tmp);
 250                 value_clear(LB);
 251                 value_clear(UB);
 252                 return(0);
 253             }
 254         }
 255         value_set_si(z[pos+S->Dimension-nparam+1],0);
 256     }
 257
 258     value_clear(c);
 259     value_clear(tmp);
 260     value_clear(LB);
 261     value_clear(UB);
 262     return 1;
 263 }
 264
 265 static int verify(Polyhedron *P, Polyhedron **C, Enumeration *en, skewed_gen_fun *gf,
 266                    arguments *options)
 267 {
 268     Polyhedron *CC, *PP, *CS, *S, *U;
 269     Matrix *C1, *MM;
 270     Vector *p;
 271     int result = 0;
 272
 273     /******* Compute true context *******/
 274     CC = align_context(*C, P->Dimension, options->barvinok->MaxRays);
 275     PP = DomainIntersection(P, CC, options->barvinok->MaxRays);
 276     Domain_Free(CC);
 277     C1 = Matrix_Alloc((*C)->Dimension+1, P->Dimension+1);
 278
 279     for (int i = 0; i < C1->NbRows; i++)
 280         for (int j = 0; j < C1->NbColumns; j++)
 281             if (i == j-P->Dimension+(*C)->Dimension)
 282                 value_set_si(C1->p[i][j], 1);
 283             else
 284                 value_set_si(C1->p[i][j], 0);
 285     CC = Polyhedron_Image(PP, C1, options->barvinok->MaxRays);
 286     Matrix_Free(C1);
 287     Domain_Free(PP);
 288     Domain_Free(*C);
 289     *C = CC;
 290
 291     /* Intersect context with range */
 292     if ((*C)->Dimension > 0) {
 293         MM = Matrix_Alloc(2*(*C)->Dimension, (*C)->Dimension+2);
 294         for (int i = 0; i < (*C)->Dimension; ++i) {
 295             value_set_si(MM->p[2*i][0], 1);
 296             value_set_si(MM->p[2*i][1+i], 1);
 297             value_set_si(MM->p[2*i][1+(*C)->Dimension], -options->verify.m);
 298             value_set_si(MM->p[2*i+1][0], 1);
 299             value_set_si(MM->p[2*i+1][1+i], -1);
 300             value_set_si(MM->p[2*i+1][1+(*C)->Dimension], options->verify.M);
 301         }
 302         CC = AddConstraints(MM->p[0], 2*(*C)->Dimension, *C,
 303                             options->barvinok->MaxRays);
 304         U = Universe_Polyhedron(0);
 305         CS = Polyhedron_Scan(CC, U, options->barvinok->MaxRays);
 306         Polyhedron_Free(U);
 307         Polyhedron_Free(CC);
 308         Matrix_Free(MM);
 309     } else
 310         CS = NULL;
 311
 312     p = Vector_Alloc(P->Dimension+2);
 313     value_set_si(p->p[P->Dimension+1], 1);
 314
 315     /* S = scanning list of polyhedra */
 316     S = Polyhedron_Scan(P, *C, options->barvinok->MaxRays);
 317
 318     if (!options->verify.print_all)
 319         if ((*C)->Dimension > 0) {
 320             int d = options->verify.M - options->verify.m;
 321             if (d > 80)
 322                 st = 1+d/80;
 323             else
 324                 st=1;
 325             for (int i = options->verify.m; i <= options->verify.M; i += st)
 326                 printf(".");
 327             printf( "\r" );
 328             fflush(stdout);
 329         }
 330
 331     /******* CHECK NOW *********/
 332     if (S) {
 333         if (!options->series || options->function) {
 334             if (!check_poly(S, CS, en, (*C)->Dimension, 0, p->p,
 335                             options->verify.print_all))
 336                 result = -1;
 337         } else {
 338             if (!check_series(S, CS, gf, (*C)->Dimension, 0, p->p,
 339                               options->verify.print_all, options->barvinok))
 340                 result = -1;
 341         }
 342         Domain_Free(S);
 343     }
 344
 345     if (result == -1)
 346         fprintf(stderr,"Check failed !\n");
 347
 348     if (!options->verify.print_all)
 349         printf( "\n" );
 350
 351     Vector_Free(p);
 352     if (CS)
 353         Domain_Free(CS);
 354
 355     return result;
 356 }
 357
 358 static void unimodular_complete(Matrix *M, int row)
 359 {
 360     Matrix *H, *Q, *U;
 361     left_hermite(M, &H, &Q, &U);
 362     Matrix_Free(H);
 363     Matrix_Free(U);
 364     for (int r = row; r < M->NbRows; ++r)
 365         Vector_Copy(Q->p[r], M->p[r], M->NbColumns);
 366     Matrix_Free(Q);
 367 }
 368
 369 static skewed_gen_fun *series(Polyhedron *P, Polyhedron* C,
 370                                 barvinok_options *options)
 371 {
 372     Polyhedron *C1, *C2;
 373     gen_fun *gf;
 374     Matrix *inv = NULL;
 375     Matrix *eq = NULL;
 376     Matrix *div = NULL;
 377     Polyhedron *PT = P;
 378
 379     /* Compute true context */
 380     C1 = Polyhedron_Project(P, C->Dimension);
 381     C2 = DomainIntersection(C, C1, options->MaxRays);
 382     Polyhedron_Free(C1);
 383
 384     POL_ENSURE_VERTICES(C2);
 385     if (C2->NbBid != 0) {
 386         Polyhedron *T;
 387         Matrix *M, *Minv, *M2;
 388         Matrix *CP;
 389         if (C2->NbEq || P->NbEq) {
 390             /* We remove all equalities to be sure all lines are unit vectors */
 391             Polyhedron *CT = C2;
 392             remove_all_equalities(&PT, &CT, &CP, NULL, C2->Dimension,
 393                                   options->MaxRays);
 394             if (CT != C2) {
 395                 Polyhedron_Free(C2);
 396                 C2 = CT;
 397             }
 398             if (CP) {
 399                 inv = left_inverse(CP, &eq);
 400                 Matrix_Free(CP);
 401
 402                 int d = 0;
 403                 Value tmp;
 404                 value_init(tmp);
 405                 div = Matrix_Alloc(inv->NbRows-1, inv->NbColumns+1);
 406                 for (int i = 0; i < inv->NbRows-1; ++i) {
 407                     Vector_Gcd(inv->p[i], inv->NbColumns, &tmp);
 408                     if (mpz_divisible_p(tmp,
 409                                         inv->p[inv->NbRows-1][inv->NbColumns-1]))
 410                         continue;
 411                     Vector_Copy(inv->p[i], div->p[d], inv->NbColumns);
 412                     value_assign(div->p[d][inv->NbColumns],
 413                                  inv->p[inv->NbRows-1][inv->NbColumns-1]);
 414                     ++d;
 415                 }
 416                 value_clear(tmp);
 417
 418                 if (!d) {
 419                     Matrix_Free(div);
 420                     div = NULL;
 421                 } else
 422                     div->NbRows = d;
 423             }
 424         }
 425         POL_ENSURE_VERTICES(C2);
 426
 427         /* Since we have "compressed" the parameters (in case there were
 428          * any equalities), the result is independent of the coordinates in the
 429          * coordinate subspace spanned by the lines.  We can therefore assume
 430          * these coordinates are zero and compute the inverse image of the map
 431          * from a lower dimensional space that adds zeros in the appropriate
 432          * places.
 433          */
 434         M = Matrix_Alloc(C2->Dimension+1, C2->Dimension-C2->NbBid+1);
 435         int k = 0;
 436         for (int i = 0; i < C2->NbBid; ++i) {
 437             int j = First_Non_Zero(C2->Ray[i]+1, C2->Dimension);
 438             assert(First_Non_Zero(C2->Ray[i]+1+j+1, C2->Dimension-j-1) == -1);
 439             for ( ; k < j; k++)
 440                 value_set_si(M->p[k+i][k], 1);
 441         }
 442         for ( ; k < C2->Dimension-C2->NbBid+1; k++)
 443             value_set_si(M->p[k+C2->NbBid][k], 1);
 444         Minv = Transpose(M);
 445         M2 = align_matrix(M, PT->Dimension + 1);
 446
 447         T = PT;
 448         PT = Polyhedron_Preimage(PT, M2, options->MaxRays);
 449         if (T != P)
 450             Polyhedron_Free(T);
 451
 452         T = C2;
 453         C2 = Polyhedron_Preimage(C2, M, options->MaxRays);
 454         if (T != C)
 455             Polyhedron_Free(T);
 456
 457         Matrix_Free(M);
 458         Matrix_Free(M2);
 459
 460         if (inv) {
 461             Matrix *T = inv;
 462             inv = Matrix_Alloc(Minv->NbRows, T->NbColumns);
 463             Matrix_Product(Minv, T, inv);
 464             Matrix_Free(T);
 465             Matrix_Free(Minv);
 466         } else
 467             inv = Minv;
 468     }
 469     gf = barvinok_series_with_options(PT, C2, options);
 470     Polyhedron_Free(C2);
 471     if (PT != P)
 472         Polyhedron_Free(PT);
 473     return new skewed_gen_fun(gf, inv, eq, div);
 474 }
 475
 476 int main(int argc, char **argv)
 477 {
 478     Polyhedron *A, *C;
 479     Matrix *M;
 480     evalue *EP = NULL;
 481     Enumeration *en = NULL;
 482     skewed_gen_fun *gf = NULL;
 483     char **param_name;
 484     int print_solution = 1;
 485     int result = 0;
 486     struct arguments options;
 487     static struct argp_child argp_children[] = {
 488         { &barvinok_argp,       0,      0,              0 },
 489         { &verify_argp,         0,      "verification", 1 },
 490         { 0 }
 491     };
 492     static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
 493     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
 494
 495     options.barvinok = bv_options;
 496     argp_parse(&argp, argc, argv, 0, 0, &options);
 497
 498     M = Matrix_Read();
 499     A = Constraints2Polyhedron(M, bv_options->MaxRays);
 500     Matrix_Free(M);
 501     M = Matrix_Read();
 502     C = Constraints2Polyhedron(M, bv_options->MaxRays);
 503     Matrix_Free(M);
 504     param_name = Read_ParamNames(stdin, C->Dimension);
 505
 506     if (options.verify.verify) {
 507         verify_options_set_range(&options.verify, A);
 508         if (!options.verbose)
 509             print_solution = 0;
 510     }
 511
 512     if (print_solution) {
 513         Polyhedron_Print(stdout, P_VALUE_FMT, A);
 514         Polyhedron_Print(stdout, P_VALUE_FMT, C);
 515     }
 516
 517     if (options.series) {
 518         gf = series(A, C, bv_options);
 519         if (print_solution) {
 520             gf->print(stdout, C->Dimension, param_name);
 521             puts("");
 522         }
 523         if (options.function) {
 524             EP = *gf;
 525             if (print_solution)
 526                 print_evalue(stdout, EP, param_name);
 527         }
 528     } else {
 529         EP = barvinok_enumerate_with_options(A, C, bv_options);
 530         if (print_solution)
 531             print_evalue(stdout, EP, param_name);
 532         if (options.size)
 533             printf("\nSize: %d\n", evalue_size(EP));
 534         if (options.floor) {
 535             fprintf(stderr, "WARNING: floor conversion not supported\n");
 536             evalue_frac2floor2(EP, 0);
 537             print_evalue(stdout, EP, param_name);
 538         } else if (options.convert) {
 539             evalue_mod2table(EP, C->Dimension);
 540             print_evalue(stdout, EP, param_name);
 541             if (options.size)
 542                 printf("\nSize: %d\n", evalue_size(EP));
 543         }
 544     }
 545
 546     if (options.verify.verify) {
 547         if (EP) {
 548             en = partition2enumeration(EP);
 549             EP = NULL;
 550         }
 551         result = verify(A, &C, en, gf, &options);
 552     }
 553
 554     if (en)
 555         Enumeration_Free(en);
 556     if (gf)
 557         delete gf;
 558     if (EP) {
 559         free_evalue_refs(EP);
 560         free(EP);
 561     }
 562
 563     Free_ParamNames(param_name, C->Dimension);
 564     Polyhedron_Free(A);
 565     Polyhedron_Free(C);
 566     free(bv_options);
 567     return result;
 568 }