verify.c: evalue_optimum: allow computation of optimum in 0D domain
[barvinok.git] / verify.c
blob10b59aaf3abf974247ada9e4e60e88fcd06b3895
1 #include <assert.h>
2 #include <stdlib.h>
3 #include <barvinok/options.h>
4 #include <barvinok/util.h>
5 #include "verify.h"
7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
9 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
10 #define RANGE 50
12 /* SRANGE : small range for evalutations */
13 #define SRANGE 15
15 /* if dimension >= BIDDIM, use SRANGE */
16 #define BIGDIM 5
18 /* VSRANGE : very small range for evalutations */
19 #define VSRANGE 5
21 /* if dimension >= VBIDDIM, use VSRANGE */
22 #define VBIGDIM 8
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 }
35 static error_t parse_opt(int key, char *arg, struct argp_state *state)
37 struct verify_options *options = state->input;
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;
80 return 0;
83 void verify_options_set_range(struct verify_options *options, int dim)
85 int r;
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.
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;
103 if (options->verify && options->m > options->M) {
104 fprintf(stderr,"Nothing to do: min > max !\n");
105 exit(0);
109 struct argp verify_argp = {
110 argp_options, parse_opt, 0, 0
113 static Polyhedron *project_on(Polyhedron *P, int i)
115 unsigned dim = P->Dimension;
116 Matrix *T;
117 Polyhedron *I;
119 if (dim == 1)
120 return Polyhedron_Copy(P);
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;
130 static void set_bounds(Polyhedron *C, Value **rows, int i, unsigned nparam,
131 const struct verify_options *options)
133 Value min;
134 Value max;
136 value_init(min);
137 value_init(max);
138 value_set_si(min, options->m);
139 value_set_si(max, options->M);
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);
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);
161 value_clear(min);
162 value_clear(max);
165 Polyhedron *check_poly_context_scan(Polyhedron *P, Polyhedron **C,
166 unsigned nparam,
167 const struct verify_options *options)
169 int i;
170 Matrix *MM;
171 Polyhedron *CC, *CC2, *CS, *U;
172 unsigned MaxRays = options->barvinok->MaxRays;
174 if (nparam <= 0)
175 return NULL;
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;
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;
204 void check_poly_init(Polyhedron *C, struct verify_options *options)
206 int d, i;
208 if (options->print_all)
209 return;
210 if (C->Dimension <= 0)
211 return;
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);
224 static void print_rational(FILE *out, Value n, Value d)
226 value_print(out, VALUE_FMT, n);
227 if (value_notone_p(d)) {
228 fprintf(out, "/");
229 value_print(out, VALUE_FMT, d);
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)
240 int k;
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]);
249 printf(") = ");
250 print_rational(stdout, got_n, got_d);
251 printf(", %s = ", check);
252 print_rational(stdout, want_n, want_d);
253 printf(". ");
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]);
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");
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 /****************************************************/
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)
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);
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;
313 value_set_si(z[pos], 0);
314 value_clear(LB);
315 value_clear(UB);
317 return 1;
318 } /* check_poly */
320 static int check_EP_on_poly(Polyhedron *P,
321 struct check_EP_data *data,
322 unsigned nvar, unsigned nparam,
323 struct verify_options *options)
325 Polyhedron *CS;
326 unsigned MaxRays = options->barvinok->MaxRays;
327 int ok = 1;
328 const evalue *EP = data->EP;
330 CS = check_poly_context_scan(NULL, &P, P->Dimension, options);
332 check_poly_init(P, options);
334 if (!(CS && emptyQ2(CS))) {
335 int i;
336 int n_S = 0;
338 for (i = 0; i < EP->x.p->size/2; ++i) {
339 Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
340 for (; A; A = A->next)
341 ++n_S;
344 data->n_S = n_S;
345 data->S = ALLOCN(Polyhedron *, n_S);
346 n_S = 0;
347 for (i = 0; i < EP->x.p->size/2; ++i) {
348 Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
349 for (; A; A = A->next) {
350 Polyhedron *next = A->next;
351 A->next = NULL;
352 data->S[n_S++] = Polyhedron_Scan(A, P,
353 MaxRays & POL_NO_DUAL ? 0 : MaxRays);
354 A->next = next;
357 ok = check_poly(CS, &data->cp, nparam, 0, data->cp.z+1+nvar, options);
358 for (i = 0; i < data->n_S; ++i)
359 Domain_Free(data->S[i]);
360 free(data->S);
363 if (!options->print_all)
364 printf("\n");
366 if (CS) {
367 Domain_Free(CS);
368 Domain_Free(P);
371 return ok;
375 * Project on final dim dimensions
377 Polyhedron *DomainProject(Polyhedron *D, unsigned dim, unsigned MaxRays)
379 Polyhedron *P;
380 Polyhedron *R;
382 R = Polyhedron_Project(D, dim);
383 for (P = D->next; P; P = P->next) {
384 Polyhedron *R2 = Polyhedron_Project(P, dim);
385 Polyhedron *R3 = DomainUnion(R, R2, MaxRays);
386 Polyhedron_Free(R2);
387 Domain_Free(R);
388 R = R3;
390 return R;
393 static Polyhedron *evalue_parameter_domain(const evalue *e, unsigned nparam,
394 unsigned MaxRays)
396 int i;
397 Polyhedron *U = NULL;
399 if (EVALUE_IS_ZERO(*e))
400 return Universe_Polyhedron(0);
402 assert(value_zero_p(e->d));
403 assert(e->x.p->type == partition);
404 assert(e->x.p->size >= 2);
406 for (i = 0; i < e->x.p->size/2; ++i) {
407 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
408 Polyhedron *P = DomainProject(D, nparam, MaxRays);
409 if (!U) {
410 U = P;
411 } else {
412 Polyhedron *T = U;
413 U = DomainUnion(U, P, MaxRays);
414 Domain_Free(P);
415 Domain_Free(T);
418 return U;
421 int check_EP(struct check_EP_data *data, unsigned nvar, unsigned nparam,
422 struct verify_options *options)
424 Vector *p;
425 int i;
426 int ok = 1;
427 Polyhedron *D, *P;
429 p = Vector_Alloc(nvar+nparam+2);
430 value_set_si(p->p[nvar+nparam+1], 1);
432 data->cp.z = p->p;
434 D = evalue_parameter_domain(data->EP, nparam, options->barvinok->MaxRays);
436 for (P = D; P; P = P->next) {
437 ok = check_EP_on_poly(P, data, nvar, nparam, options);
438 if (!ok && !options->continue_on_error)
439 break;
442 Domain_Free(D);
443 Vector_Free(p);
445 return ok;
448 static void optimum(Polyhedron *S, int pos, const struct check_EP_data *data,
449 Value *opt, int *found, int sign)
451 if (!S) {
452 double d;
453 Value c;
454 value_init(c);
455 /* value_set_double rounds to zero; give it some slop */
456 d = compute_evalue(data->EP, data->cp.z+1);
457 if (d > 0)
458 d += .25;
459 else
460 d -= .25;
461 value_set_double(c, d);
462 if (!*found) {
463 value_assign(*opt, c);
464 *found = 1;
465 } else {
466 if (sign > 0) {
467 if (value_gt(c, *opt))
468 value_assign(*opt, c);
469 } else {
470 if (value_lt(c, *opt))
471 value_assign(*opt, c);
474 value_clear(c);
475 } else {
476 Value LB, UB;
477 int ok;
478 value_init(LB);
479 value_init(UB);
480 ok = !(lower_upper_bounds(1+pos, S, data->cp.z, &LB, &UB));
481 assert(ok);
482 for (; value_le(LB, UB); value_increment(LB, LB)) {
483 value_assign(data->cp.z[1+pos], LB);
484 optimum(S->next, pos+1, data, opt, found, sign);
486 value_set_si(data->cp.z[1+pos], 0);
487 value_clear(LB);
488 value_clear(UB);
492 void evalue_optimum(const struct check_EP_data *data, Value *opt, int sign)
494 int i;
495 int found = 0;
497 for (i = 0; i < data->n_S; ++i)
498 if (!(data->S[i] && emptyQ(data->S[i])))
499 optimum(data->S[i], 0, data, opt, &found, sign);
500 assert(found);