?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
verify.c: extract some helper functions for isl based verification
[barvinok/uuh.git] / verify.c
blobfdf21105f850ce1a4d2293729a0bd0b0b1a6ed66
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 }
516
517 __isl_give isl_set *verify_context_set_bounds(__isl_take isl_set *set,
518 const struct verify_options *options)
519 {
520 int i;
521 unsigned nparam;
522 isl_point *pt, *pt2;
523 isl_set *box;
524
525 nparam = isl_set_dim(set, isl_dim_param);
526
527 if (options->r > 0) {
528 pt = isl_set_sample_point(isl_set_copy(set));
529 pt2 = isl_point_copy(pt);
530
531 for (i = 0; i < nparam; ++i) {
532 pt = isl_point_add_ui(pt, isl_dim_param, i, options->r);
533 pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, options->r);
534 }
535 } else {
536 isl_int v;
537
538 isl_int_init(v);
539 pt = isl_point_zero(isl_set_get_dim(set));
540 isl_int_set_si(v, options->m);
541 for (i = 0; i < nparam; ++i)
542 pt = isl_point_set_coordinate(pt, isl_dim_param, i, v);
543
544 pt2 = isl_point_zero(isl_set_get_dim(set));
545 isl_int_set_si(v, options->M);
546 for (i = 0; i < nparam; ++i)
547 pt2 = isl_point_set_coordinate(pt2, isl_dim_param, i, v);
548
549 isl_int_clear(v);
550 }
551
552 box = isl_set_box_from_points(pt, pt2);
553
554 return isl_set_intersect(set, box);
555 }
556
557 int verify_point_data_init(struct verify_point_data *vpd,
558 __isl_keep isl_set *context)
559 {
560 isl_int v;
561 int i;
562 int r;
563
564 isl_int_init(v);
565 r = isl_set_count(context, &v);
566 vpd->n = isl_int_cmp_si(v, 200) < 0 ? isl_int_get_si(v) : 200;
567 isl_int_clear(v);
568
569 if (!vpd->options->print_all) {
570 vpd->s = vpd->n < 80 ? 1 : 1 + vpd->n/80;
571 for (i = 0; i < vpd->n; i += vpd->s)
572 printf(".");
573 printf("\r");
574 fflush(stdout);
575 }
576
577 vpd->error = r < 0 ? -1 : 0;
578
579 return r;
580 }
581
582 void verify_point_data_fini(struct verify_point_data *vpd)
583 {
584 if (!vpd->options->print_all)
585 printf("\n");
586
587 if (vpd->error)
588 fprintf(stderr, "Check failed !\n");
589 }
patches by uuh to consider upstream