verify.c

   1 #include <assert.h>
   2 #include <stdlib.h>
   3 #include <barvinok/options.h>
   4 #include <barvinok/util.h>
   5 #include "verify.h"
   6
   7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   8
   9 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
  10 #define RANGE 50
  11
  12 /* SRANGE : small range for evalutations */
  13 #define SRANGE 15
  14
  15 /* if dimension >= BIDDIM, use SRANGE */
  16 #define BIGDIM 5
  17
  18 /* VSRANGE : very small range for evalutations */
  19 #define VSRANGE 5
  20
  21 /* if dimension >= VBIDDIM, use VSRANGE */
  22 #define VBIGDIM 8
  23
  24 static struct argp_option argp_options[] = {
  25     { "verify",             'T',    0,      0 },
  26     { "exact",              'E',    0,      0 },
  27     { "print-all",          'A',    0,      0 },
  28     { "continue-on-error",  'C',    0,      0 },
  29     { "min",                'm',    "int",  0 },
  30     { "max",                'M',    "int",  0 },
  31     { "range",              'r',    "int",  0 },
  32     { 0 }
  33 };
  34
  35 static error_t parse_opt(int key, char *arg, struct argp_state *state)
  36 {
  37     struct verify_options *options = state->input;
  38
  39     switch (key) {
  40     case ARGP_KEY_INIT:
  41         options->verify = 0;
  42         options->exact = 0;
  43         options->print_all = 0;
  44         options->continue_on_error = 0;
  45         options->m = INT_MAX;
  46         options->M = INT_MIN;
  47         options->r = -1;
  48         break;
  49     case ARGP_KEY_FINI:
  50         break;
  51     case 'T':
  52         options->verify = 1;
  53         break;
  54     case 'E':
  55         options->exact = 1;
  56         break;
  57     case 'A':
  58         options->print_all = 1;
  59         break;
  60     case 'C':
  61         options->continue_on_error = 1;
  62         break;
  63     case 'm':
  64         options->m = atoi(arg);
  65         options->verify = 1;
  66         break;
  67     case 'M':
  68         options->M = atoi(arg);
  69         options->verify = 1;
  70         break;
  71     case 'r':
  72         options->r = atoi(arg);
  73         options->M = atoi(arg);
  74         options->m = -options->M;
  75         options->verify = 1;
  76         break;
  77     default:
  78         return ARGP_ERR_UNKNOWN;
  79     }
  80     return 0;
  81 }
  82
  83 void verify_options_set_range(struct verify_options *options, int dim)
  84 {
  85     int r;
  86
  87     if (dim >= VBIGDIM)
  88         r = VSRANGE;
  89     else if (dim >= BIGDIM)
  90         r = SRANGE;
  91     else
  92         r = RANGE;
  93     /* If the user didn't set m or M, then we try to adjust the window
  94      * to the context in check_poly_context_scan.
  95      */
  96     if (options->m == INT_MAX && options->M == INT_MIN)
  97         options->r = r;
  98     if (options->M == INT_MIN)
  99         options->M = r;
 100     if (options->m == INT_MAX)
 101         options->m = -r;
 102
 103     if (options->verify && options->m > options->M) {
 104         fprintf(stderr,"Nothing to do: min > max !\n");
 105         exit(0);
 106     }
 107 }
 108
 109 struct argp verify_argp = {
 110     argp_options, parse_opt, 0, 0
 111 };
 112
 113 static Polyhedron *project_on(Polyhedron *P, int i)
 114 {
 115     unsigned dim = P->Dimension;
 116     Matrix *T;
 117     Polyhedron *I;
 118
 119     if (dim == 1)
 120         return Polyhedron_Copy(P);
 121
 122     T = Matrix_Alloc(2, dim+1);
 123     value_set_si(T->p[0][i], 1);
 124     value_set_si(T->p[1][dim], 1);
 125     I = Polyhedron_Image(P, T, P->NbConstraints);
 126     Matrix_Free(T);
 127     return I;
 128 }
 129
 130 static void set_bounds(Polyhedron *C, Value **rows, int i, unsigned nparam,
 131                         const struct verify_options *options)
 132 {
 133     Value min;
 134     Value max;
 135
 136     value_init(min);
 137     value_init(max);
 138     value_set_si(min, options->m);
 139     value_set_si(max, options->M);
 140
 141     if (options->r > 0) {
 142         Polyhedron *I = project_on(C, i);
 143         line_minmax(I, &min, &max);
 144         if (value_cmp_si(min, options->M) >= 0)
 145             value_add_int(max, min, options->r);
 146         else if (value_cmp_si(max, options->m) <= 0)
 147             value_sub_int(min, max, options->r);
 148         else {
 149             value_set_si(min, options->m);
 150             value_set_si(max, options->M);
 151         }
 152     }
 153
 154     value_set_si(rows[0][0], 1);
 155     value_set_si(rows[0][1+i], 1);
 156     value_oppose(rows[0][1+nparam], min);
 157     value_set_si(rows[1][0], 1);
 158     value_set_si(rows[1][1+i], -1);
 159     value_assign(rows[1][1+nparam], max);
 160
 161     value_clear(min);
 162     value_clear(max);
 163 }
 164
 165 Polyhedron *check_poly_context_scan(Polyhedron *P, Polyhedron **C,
 166                                     unsigned nparam,
 167                                     const struct verify_options *options)
 168 {
 169     int i;
 170     Matrix *MM;
 171     Polyhedron *CC, *CC2, *CS, *U;
 172     unsigned MaxRays = options->barvinok->MaxRays;
 173
 174     if (nparam <= 0)
 175         return NULL;
 176
 177     if (!P)
 178         CC = *C;
 179     else {
 180         CC = Polyhedron_Project(P, nparam);
 181         if (*C) {
 182             CC2 = DomainIntersection(*C, CC, MaxRays);
 183             Domain_Free(CC);
 184             CC = CC2;
 185         }
 186     }
 187
 188     /* Intersect context with range */
 189     MM = Matrix_Alloc(2*nparam, nparam+2);
 190     for (i = 0; i < nparam; ++i)
 191         set_bounds(CC, &MM->p[2*i], i, nparam, options);
 192     CC2 = AddConstraints(MM->p[0], 2*nparam, CC, options->barvinok->MaxRays);
 193     if (CC != *C)
 194         Domain_Free(CC);
 195     CC = CC2;
 196     U = Universe_Polyhedron(0);
 197     CS = Polyhedron_Scan(CC, U, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 198     Polyhedron_Free(U);
 199     *C = CC;
 200     Matrix_Free(MM);
 201     return CS;
 202 }
 203
 204 void check_poly_init(Polyhedron *C, struct verify_options *options)
 205 {
 206     int d, i;
 207
 208     if (options->print_all)
 209         return;
 210     if (C->Dimension <= 0)
 211         return;
 212
 213     d = options->M - options->m;
 214     if (d > 80)
 215         options->st = 1+d/80;
 216     else
 217         options->st = 1;
 218     for (i = options->m; i <= options->M; i += options->st)
 219         printf(".");
 220     printf( "\r" );
 221     fflush(stdout);
 222 }
 223
 224 static void print_rational(FILE *out, Value n, Value d)
 225 {
 226     value_print(out, VALUE_FMT, n);
 227     if (value_notone_p(d)) {
 228         fprintf(out, "/");
 229         value_print(out, VALUE_FMT, d);
 230     }
 231 }
 232
 233 void check_poly_print(int ok, int nparam, Value *z,
 234                       Value want_n, Value want_d,
 235                       Value got_n, Value got_d,
 236                       const char *op, const char *check,
 237                       const char *long_op,
 238                       const struct verify_options *options)
 239 {
 240     int k;
 241
 242     if (options->print_all) {
 243         printf("%s(", op);
 244         value_print(stdout, VALUE_FMT, z[0]);
 245         for (k = 1; k < nparam; ++k) {
 246             printf(", ");
 247             value_print(stdout, VALUE_FMT, z[k]);
 248         }
 249         printf(") = ");
 250         print_rational(stdout, got_n, got_d);
 251         printf(", %s = ", check);
 252         print_rational(stdout, want_n, want_d);
 253         printf(". ");
 254     }
 255
 256     if (!ok) {
 257         printf("\n");
 258         fflush(stdout);
 259         fprintf(stderr, "Error !\n");
 260         fprintf(stderr, "%s(", op);
 261         value_print(stderr, VALUE_FMT, z[0]);
 262         for (k = 1; k < nparam; ++k) {
 263             fprintf(stderr,", ");
 264             value_print(stderr, VALUE_FMT, z[k]);
 265         }
 266         fprintf(stderr, ") should be ");
 267         print_rational(stderr, want_n, want_d);
 268         fprintf(stderr, ", while %s gives ", long_op);
 269         print_rational(stderr, got_n, got_d);
 270         fprintf(stderr, ".\n");
 271     } else if (options->print_all)
 272         printf("OK.\n");
 273 }
 274
 275 /****************************************************/
 276 /* function check_poly :                            */
 277 /* scans the parameter space from Min to Max (all   */
 278 /* directions). Computes the number of points in    */
 279 /* the polytope using both methods, and compare them*/
 280 /* returns 1 on success                             */
 281 /****************************************************/
 282
 283 int check_poly(Polyhedron *CS, const struct check_poly_data *data,
 284                int nparam, int pos, Value *z,
 285                const struct verify_options *options)
 286 {
 287     if (pos == nparam) {
 288         if (!data->check(data, nparam, z, options) && !options->continue_on_error)
 289             return 0;
 290     } else {
 291         Value LB, UB;
 292         int ok;
 293         value_init(LB);
 294         value_init(UB);
 295         ok = !(lower_upper_bounds(1+pos, CS, z-1, &LB, &UB));
 296         assert(ok);
 297         for (; value_le(LB, UB); value_increment(LB, LB)) {
 298             if (!options->print_all) {
 299                 int k = VALUE_TO_INT(LB);
 300                 if (!pos && !(k % options->st)) {
 301                     printf("o");
 302                     fflush(stdout);
 303                 }
 304             }
 305
 306             value_assign(z[pos], LB);
 307             if (!check_poly(CS->next, data, nparam, pos+1, z, options)) {
 308                 value_clear(LB);
 309                 value_clear(UB);
 310                 return 0;
 311             }
 312         }
 313         value_set_si(z[pos], 0);
 314         value_clear(LB);
 315         value_clear(UB);
 316     }
 317     return 1;
 318 } /* check_poly */
 319
 320 void check_EP_set_scan(struct check_EP_data *data, Polyhedron *C,
 321                         unsigned MaxRays)
 322 {
 323     const evalue *EP = data->EP;
 324     int i;
 325     int n_S = 0;
 326
 327     for (i = 0; i < EP->x.p->size/2; ++i) {
 328         Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 329         for (; A; A = A->next)
 330             ++n_S;
 331     }
 332
 333     data->n_S = n_S;
 334     data->S = ALLOCN(Polyhedron *, n_S);
 335     n_S = 0;
 336     for (i = 0; i < EP->x.p->size/2; ++i) {
 337         Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 338         for (; A; A = A->next) {
 339             Polyhedron *next = A->next;
 340             A->next = NULL;
 341             data->S[n_S++] = Polyhedron_Scan(A, C,
 342                                     MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 343             A->next = next;
 344         }
 345     }
 346 }
 347
 348 void check_EP_clear_scan(struct check_EP_data *data)
 349 {
 350     int i;
 351
 352     for (i = 0; i < data->n_S; ++i)
 353         Domain_Free(data->S[i]);
 354     free(data->S);
 355 }
 356
 357 static int check_EP_on_poly(Polyhedron *P,
 358                             struct check_EP_data *data,
 359                             unsigned nvar, unsigned nparam,
 360                             struct verify_options *options)
 361 {
 362     Polyhedron *CS;
 363     unsigned MaxRays = options->barvinok->MaxRays;
 364     int ok = 1;
 365     const evalue *EP = data->EP;
 366
 367     CS = check_poly_context_scan(NULL, &P, P->Dimension, options);
 368
 369     check_poly_init(P, options);
 370
 371     if (!(CS && emptyQ2(CS))) {
 372         check_EP_set_scan(data, P, MaxRays);
 373         ok = check_poly(CS, &data->cp, nparam, 0, data->cp.z+1+nvar, options);
 374         check_EP_clear_scan(data);
 375     }
 376
 377     if (!options->print_all)
 378         printf("\n");
 379
 380     if (CS) {
 381         Domain_Free(CS);
 382         Domain_Free(P);
 383     }
 384
 385     return ok;
 386 }
 387
 388 /*
 389  * Project on final dim dimensions
 390  */
 391 Polyhedron *DomainProject(Polyhedron *D, unsigned dim, unsigned MaxRays)
 392 {
 393     Polyhedron *P;
 394     Polyhedron *R;
 395
 396     R = Polyhedron_Project(D, dim);
 397     for (P = D->next; P; P = P->next) {
 398         Polyhedron *R2 = Polyhedron_Project(P, dim);
 399         Polyhedron *R3 = DomainUnion(R, R2, MaxRays);
 400         Polyhedron_Free(R2);
 401         Domain_Free(R);
 402         R = R3;
 403     }
 404     return R;
 405 }
 406
 407 static Polyhedron *evalue_parameter_domain(const evalue *e, unsigned nparam,
 408                                              unsigned MaxRays)
 409 {
 410     int i;
 411     Polyhedron *U = NULL;
 412
 413     if (EVALUE_IS_ZERO(*e))
 414         return Universe_Polyhedron(0);
 415
 416     assert(value_zero_p(e->d));
 417     assert(e->x.p->type == partition);
 418     assert(e->x.p->size >= 2);
 419
 420     for (i = 0; i < e->x.p->size/2; ++i) {
 421         Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 422         Polyhedron *P = DomainProject(D, nparam, MaxRays);
 423         if (!U) {
 424             U = P;
 425         } else {
 426             Polyhedron *T = U;
 427             U = DomainUnion(U, P, MaxRays);
 428             Domain_Free(P);
 429             Domain_Free(T);
 430         }
 431     }
 432     return U;
 433 }
 434
 435 int check_EP(struct check_EP_data *data, unsigned nvar, unsigned nparam,
 436              struct verify_options *options)
 437 {
 438     Vector *p;
 439     int i;
 440     int ok = 1;
 441     Polyhedron *D, *P;
 442
 443     p = Vector_Alloc(nvar+nparam+2);
 444     value_set_si(p->p[nvar+nparam+1], 1);
 445
 446     data->cp.z = p->p;
 447
 448     D = evalue_parameter_domain(data->EP, nparam, options->barvinok->MaxRays);
 449
 450     for (P = D; P; P = P->next) {
 451         ok = check_EP_on_poly(P, data, nvar, nparam, options);
 452         if (!ok && !options->continue_on_error)
 453             break;
 454     }
 455
 456     Domain_Free(D);
 457     Vector_Free(p);
 458
 459     return ok;
 460 }
 461
 462 static void optimum(Polyhedron *S, int pos, const struct check_EP_data *data,
 463                     Value *opt, int *found, int sign)
 464 {
 465     if (!S) {
 466         double d;
 467         Value c;
 468         value_init(c);
 469         /* value_set_double rounds to zero; give it some slop */
 470         d = compute_evalue(data->EP, data->cp.z+1);
 471         if (d > 0)
 472             d += .25;
 473         else
 474             d -= .25;
 475         value_set_double(c, d);
 476         if (!*found) {
 477             value_assign(*opt, c);
 478             *found = 1;
 479         } else {
 480             if (sign > 0) {
 481                 if (value_gt(c, *opt))
 482                     value_assign(*opt, c);
 483             } else {
 484                 if (value_lt(c, *opt))
 485                     value_assign(*opt, c);
 486             }
 487         }
 488         value_clear(c);
 489     } else {
 490         Value LB, UB;
 491         int ok;
 492         value_init(LB);
 493         value_init(UB);
 494         ok = !(lower_upper_bounds(1+pos, S, data->cp.z, &LB, &UB));
 495         assert(ok);
 496         for (; value_le(LB, UB); value_increment(LB, LB)) {
 497             value_assign(data->cp.z[1+pos], LB);
 498             optimum(S->next, pos+1, data, opt, found, sign);
 499         }
 500         value_set_si(data->cp.z[1+pos], 0);
 501         value_clear(LB);
 502         value_clear(UB);
 503     }
 504 }
 505
 506 void evalue_optimum(const struct check_EP_data *data, Value *opt, int sign)
 507 {
 508     int i;
 509     int found = 0;
 510
 511     for (i = 0; i < data->n_S; ++i)
 512         if (!(data->S[i] && emptyQ(data->S[i])))
 513             optimum(data->S[i], 0, data, opt, &found, sign);
 514     assert(found);
 515 }