?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
test_bound.c: verify_point: use isl_val
[barvinok.git] / barvinok_enumerate.cc
blobbe98124749fe138355577ef25dbebfe3d6d92c12
1 #include <assert.h>
2 #include <unistd.h>
3 #include <stdlib.h>
4 #include <gmp.h>
5 #include <isl_set_polylib.h>
6 #include <barvinok/evalue.h>
7 #include <barvinok/util.h>
8 #include <barvinok/barvinok.h>
9 #include "barvinok_enumerate_options.h"
10 #include "verify.h"
11 #include "verif_ehrhart.h"
12 #include "verify_series.h"
13 #include "remove_equalities.h"
14 #include "evalue_convert.h"
15 #include "conversion.h"
16 #include "skewed_genfun.h"
17
18 #undef CS /* for Solaris 10 */
19
20 using std::cout;
21 using std::endl;
22
23 /* The input of this example program is the same as that of testehrhart
24 * in the PolyLib distribution, i.e., a polytope in combined
25 * data and parameter space, a context polytope in parameter space
26 * and (optionally) the names of the parameters.
27 * Both polytopes are in PolyLib notation.
28 */
29
30 struct verify_point_enum {
31 struct verify_point_data vpd;
32 isl_set *set;
33 isl_pw_qpolynomial *pwqp;
34 };
35
36 static int verify_point(__isl_take isl_point *pnt, void *user)
37 {
38 struct verify_point_enum *vpe = (struct verify_point_enum *) user;
39 isl_set *set;
40 int i;
41 unsigned nparam;
42 isl_val *v, *n, *t;
43 isl_qpolynomial *cnt;
44 int pa = vpe->vpd.options->barvinok->approx->approximation;
45 int cst;
46 int ok;
47 FILE *out = vpe->vpd.options->print_all ? stdout : stderr;
48
49 vpe->vpd.n--;
50
51 set = isl_set_copy(vpe->set);
52 nparam = isl_set_dim(set, isl_dim_param);
53 for (i = 0; i < nparam; ++i) {
54 v = isl_point_get_coordinate_val(pnt, isl_dim_param, i);
55 set = isl_set_fix_val(set, isl_dim_param, i, v);
56 }
57
58 v = isl_set_count_val(set);
59
60 cnt = isl_pw_qpolynomial_eval(isl_pw_qpolynomial_copy(vpe->pwqp),
61 isl_point_copy(pnt));
62 n = isl_qpolynomial_get_constant_val(cnt);
63 isl_qpolynomial_free(cnt);
64
65 if (pa == BV_APPROX_SIGN_LOWER)
66 n = isl_val_ceil(n);
67 else if (pa == BV_APPROX_SIGN_UPPER)
68 n = isl_val_floor(n);
69 else
70 n = isl_val_trunc(n);
71
72 if (pa == BV_APPROX_SIGN_APPROX)
73 /* just accept everything */
74 ok = 1;
75 else if (pa == BV_APPROX_SIGN_LOWER)
76 ok = isl_val_le(n, v);
77 else if (pa == BV_APPROX_SIGN_UPPER)
78 ok = isl_val_ge(n, v);
79 else
80 ok = isl_val_eq(n, v);
81
82 if (vpe->vpd.options->print_all || !ok) {
83 isl_ctx *ctx = isl_point_get_ctx(pnt);
84 isl_printer *p;
85 p = isl_printer_to_file(ctx, out);
86 p = isl_printer_print_str(p, "EP(");
87 for (i = 0; i < nparam; ++i) {
88 if (i)
89 p = isl_printer_print_str(p, ", ");
90 t = isl_point_get_coordinate_val(pnt, isl_dim_param, i);
91 p = isl_printer_print_val(p, t);
92 isl_val_free(t);
93 }
94 p = isl_printer_print_str(p, ") = ");
95 p = isl_printer_print_val(p, n);
96 p = isl_printer_print_str(p, ", count = ");
97 p = isl_printer_print_val(p, v);
98 if (ok)
99 p = isl_printer_print_str(p, ". OK");
100 else
101 p = isl_printer_print_str(p, ". NOT OK");
102 p = isl_printer_end_line(p);
103 isl_printer_free(p);
104 } else if ((vpe->vpd.n % vpe->vpd.s) == 0) {
105 printf("o");
106 fflush(stdout);
107 }
108
109 if (0) {
110 error:
111 ok = 0;
112 }
113 isl_set_free(set);
114 isl_val_free(v);
115 isl_val_free(n);
116 isl_point_free(pnt);
117
118 if (!ok)
119 vpe->vpd.error = 1;
120
121 if (vpe->vpd.options->continue_on_error)
122 ok = 1;
123
124 return (vpe->vpd.n >= 1 && ok) ? 0 : -1;
125 }
126
127 static int verify_isl(Polyhedron *P, Polyhedron *C,
128 evalue *EP, const struct verify_options *options)
129 {
130 struct verify_point_enum vpe = { { options } };
131 int i;
132 isl_ctx *ctx = isl_ctx_alloc();
133 isl_space *dim;
134 isl_set *set;
135 isl_set *set_C;
136 int r;
137
138 dim = isl_space_set_alloc(ctx, C->Dimension, P->Dimension - C->Dimension);
139 for (i = 0; i < C->Dimension; ++i)
140 dim = isl_space_set_dim_name(dim, isl_dim_param, i, options->params[i]);
141 set = isl_set_new_from_polylib(P, isl_space_copy(dim));
142 dim = isl_space_params(dim);
143 set_C = isl_set_new_from_polylib(C, dim);
144 set_C = isl_set_intersect_params(isl_set_copy(set), set_C);
145 set_C = isl_set_params(set_C);
146
147 set_C = verify_context_set_bounds(set_C, options);
148
149 r = verify_point_data_init(&vpe.vpd, set_C);
150
151 vpe.set = set;
152 vpe.pwqp = isl_pw_qpolynomial_from_evalue(isl_set_get_space(set_C), EP);
153 if (r == 0)
154 isl_set_foreach_point(set_C, verify_point, &vpe);
155 if (vpe.vpd.error)
156 r = -1;
157
158 isl_pw_qpolynomial_free(vpe.pwqp);
159 isl_set_free(set);
160 isl_set_free(set_C);
161
162 isl_ctx_free(ctx);
163
164 verify_point_data_fini(&vpe.vpd);
165
166 return r;
167 }
168
169 static int verify(Polyhedron *P, Polyhedron *C, evalue *EP, skewed_gen_fun *gf,
170 struct enumerate_options *options)
171 {
172 Polyhedron *CS, *S;
173 Vector *p;
174 int result = 0;
175
176 if (!options->series || options->function)
177 return verify_isl(P, C, EP, options->verify);
178
179 CS = check_poly_context_scan(P, &C, C->Dimension, options->verify);
180
181 p = Vector_Alloc(P->Dimension+2);
182 value_set_si(p->p[P->Dimension+1], 1);
183
184 /* S = scanning list of polyhedra */
185 S = Polyhedron_Scan(P, C, options->verify->barvinok->MaxRays);
186
187 check_poly_init(C, options->verify);
188
189 /******* CHECK NOW *********/
190 if (S) {
191 if (!options->series || options->function) {
192 if (!check_poly_EP(S, CS, EP, 0, C->Dimension, 0, p->p,
193 options->verify))
194 result = -1;
195 } else {
196 if (!check_poly_gf(S, CS, gf, 0, C->Dimension, 0, p->p,
197 options->verify))
198 result = -1;
199 }
200 Domain_Free(S);
201 }
202
203 if (result == -1)
204 fprintf(stderr,"Check failed !\n");
205
206 if (!options->verify->print_all)
207 printf( "\n" );
208
209 Vector_Free(p);
210 if (CS) {
211 Domain_Free(CS);
212 Domain_Free(C);
213 }
214
215 return result;
216 }
217
218 /* frees M and Minv */
219 static void apply_transformation(Polyhedron **P, Polyhedron **C,
220 bool free_P, bool free_C,
221 Matrix *M, Matrix *Minv, Matrix **inv,
222 barvinok_options *options)
223 {
224 Polyhedron *T;
225 Matrix *M2;
226
227 M2 = align_matrix(M, (*P)->Dimension + 1);
228 T = *P;
229 *P = Polyhedron_Preimage(*P, M2, options->MaxRays);
230 if (free_P)
231 Polyhedron_Free(T);
232 Matrix_Free(M2);
233
234 T = *C;
235 *C = Polyhedron_Preimage(*C, M, options->MaxRays);
236 if (free_C)
237 Polyhedron_Free(T);
238
239 Matrix_Free(M);
240
241 if (*inv) {
242 Matrix *T = *inv;
243 *inv = Matrix_Alloc(Minv->NbRows, T->NbColumns);
244 Matrix_Product(Minv, T, *inv);
245 Matrix_Free(T);
246 Matrix_Free(Minv);
247 } else
248 *inv = Minv;
249 }
250
251 /* Since we have "compressed" the parameters (in case there were
252 * any equalities), the result is independent of the coordinates in the
253 * coordinate subspace spanned by the lines. We can therefore assume
254 * these coordinates are zero and compute the inverse image of the map
255 * from a lower dimensional space that adds zeros in the appropriate
256 * places.
257 */
258 static void remove_lines(Polyhedron *C, Matrix **M, Matrix **Minv)
259 {
260 Matrix *L = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
261 for (int r = 0; r < C->NbBid; ++r)
262 Vector_Copy(C->Ray[r]+1, L->p[r], C->Dimension);
263 unimodular_complete(L, C->NbBid);
264 assert(value_one_p(L->p[C->Dimension][C->Dimension]));
265 assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
266 Matrix_Transposition(L);
267 assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
268
269 *M = Matrix_Alloc(C->Dimension+1, C->Dimension-C->NbBid+1);
270 for (int i = 0; i < C->Dimension+1; ++i)
271 Vector_Copy(L->p[i]+C->NbBid, (*M)->p[i], C->Dimension-C->NbBid+1);
272
273 Matrix *Linv = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
274 int ok = Matrix_Inverse(L, Linv);
275 assert(ok);
276 Matrix_Free(L);
277
278 *Minv = Matrix_Alloc(C->Dimension-C->NbBid+1, C->Dimension+1);
279 for (int i = C->NbBid; i < C->Dimension+1; ++i)
280 Vector_AntiScale(Linv->p[i], (*Minv)->p[i-C->NbBid],
281 Linv->p[C->Dimension][C->Dimension], C->Dimension+1);
282 Matrix_Free(Linv);
283 }
284
285 static skewed_gen_fun *series(Polyhedron *P, Polyhedron* C,
286 barvinok_options *options)
287 {
288 Polyhedron *C1, *C2;
289 gen_fun *gf;
290 Matrix *inv = NULL;
291 Matrix *eq = NULL;
292 Matrix *div = NULL;
293 Polyhedron *PT = P;
294
295 /* Compute true context */
296 C1 = Polyhedron_Project(P, C->Dimension);
297 C2 = DomainIntersection(C, C1, options->MaxRays);
298 Polyhedron_Free(C1);
299
300 POL_ENSURE_VERTICES(C2);
301 if (C2->NbBid != 0) {
302 Polyhedron *T;
303 Matrix *M, *Minv, *M2;
304 Matrix *CP;
305 if (C2->NbEq || P->NbEq) {
306 /* We remove all equalities to be sure all lines are unit vectors */
307 Polyhedron *CT = C2;
308 remove_all_equalities(&PT, &CT, &CP, NULL, C2->Dimension,
309 options->MaxRays);
310 if (CT != C2) {
311 Polyhedron_Free(C2);
312 C2 = CT;
313 }
314 if (CP) {
315 inv = left_inverse(CP, &eq);
316 Matrix_Free(CP);
317
318 int d = 0;
319 Value tmp;
320 value_init(tmp);
321 div = Matrix_Alloc(inv->NbRows-1, inv->NbColumns+1);
322 for (int i = 0; i < inv->NbRows-1; ++i) {
323 Vector_Gcd(inv->p[i], inv->NbColumns, &tmp);
324 if (mpz_divisible_p(tmp,
325 inv->p[inv->NbRows-1][inv->NbColumns-1]))
326 continue;
327 Vector_Copy(inv->p[i], div->p[d], inv->NbColumns);
328 value_assign(div->p[d][inv->NbColumns],
329 inv->p[inv->NbRows-1][inv->NbColumns-1]);
330 ++d;
331 }
332 value_clear(tmp);
333
334 if (!d) {
335 Matrix_Free(div);
336 div = NULL;
337 } else
338 div->NbRows = d;
339 }
340 }
341 POL_ENSURE_VERTICES(C2);
342
343 if (C2->NbBid) {
344 Matrix *M, *Minv;
345 remove_lines(C2, &M, &Minv);
346 apply_transformation(&PT, &C2, PT != P, C2 != C, M, Minv, &inv,
347 options);
348 }
349 }
350 POL_ENSURE_VERTICES(C2);
351 if (!Polyhedron_has_revlex_positive_rays(C2, C2->Dimension)) {
352 Polyhedron *T;
353 Matrix *Constraints;
354 Matrix *H, *Q, *U;
355 Constraints = Matrix_Alloc(C2->NbConstraints, C2->Dimension+1);
356 for (int i = 0; i < C2->NbConstraints; ++i)
357 Vector_Copy(C2->Constraint[i]+1, Constraints->p[i], C2->Dimension);
358 left_hermite(Constraints, &H, &Q, &U);
359 Matrix_Free(Constraints);
360 /* flip rows of Q */
361 for (int i = 0; i < C2->Dimension/2; ++i)
362 Vector_Exchange(Q->p[i], Q->p[C2->Dimension-1-i], C2->Dimension);
363 Matrix_Free(H);
364 Matrix_Free(U);
365 Matrix *M = Matrix_Alloc(C2->Dimension+1, C2->Dimension+1);
366 U = Matrix_Copy(Q);
367 int ok = Matrix_Inverse(U, M);
368 assert(ok);
369 Matrix_Free(U);
370
371 apply_transformation(&PT, &C2, PT != P, C2 != C, M, Q, &inv, options);
372 }
373 gf = barvinok_series_with_options(PT, C2, options);
374 Polyhedron_Free(C2);
375 if (PT != P)
376 Polyhedron_Free(PT);
377 return new skewed_gen_fun(gf, inv, eq, div);
378 }
379
380 int main(int argc, char **argv)
381 {
382 Polyhedron *A, *C;
383 Matrix *M;
384 evalue *EP = NULL;
385 skewed_gen_fun *gf = NULL;
386 const char **param_name;
387 int print_solution = 1;
388 int result = 0;
389 struct enumerate_options *options = enumerate_options_new_with_defaults();
390
391 argc = enumerate_options_parse(options, argc, argv, ISL_ARG_ALL);
392
393 M = Matrix_Read();
394 assert(M);
395 A = Constraints2Polyhedron(M, options->verify->barvinok->MaxRays);
396 Matrix_Free(M);
397 M = Matrix_Read();
398 assert(M);
399 C = Constraints2Polyhedron(M, options->verify->barvinok->MaxRays);
400 Matrix_Free(M);
401 assert(A->Dimension >= C->Dimension);
402 param_name = Read_ParamNames(stdin, C->Dimension);
403
404 if (options->verify->verify) {
405 verify_options_set_range(options->verify, A->Dimension);
406 if (!options->verify->barvinok->verbose)
407 print_solution = 0;
408 }
409
410 if (print_solution && options->verify->barvinok->verbose) {
411 Polyhedron_Print(stdout, P_VALUE_FMT, A);
412 Polyhedron_Print(stdout, P_VALUE_FMT, C);
413 }
414
415 if (options->series) {
416 gf = series(A, C, options->verify->barvinok);
417 if (print_solution) {
418 gf->print(cout, C->Dimension, param_name);
419 puts("");
420 }
421 if (options->function) {
422 EP = *gf;
423 if (print_solution)
424 print_evalue(stdout, EP, param_name);
425 }
426 } else {
427 EP = barvinok_enumerate_with_options(A, C, options->verify->barvinok);
428 assert(EP);
429 if (evalue_convert(EP, options->convert, options->verify->barvinok->verbose,
430 C->Dimension, param_name))
431 print_solution = 0;
432 if (options->size)
433 printf("\nSize: %zd\n", evalue_size(EP));
434 if (print_solution)
435 print_evalue(stdout, EP, param_name);
436 }
437
438 if (options->verify->verify) {
439 options->verify->params = param_name;
440 result = verify(A, C, EP, gf, options);
441 }
442
443 if (gf)
444 delete gf;
445 if (EP)
446 evalue_free(EP);
447
448 if (options->verify->barvinok->print_stats)
449 barvinok_stats_print(options->verify->barvinok->stats, stdout);
450
451 Free_ParamNames(param_name, C->Dimension);
452 Polyhedron_Free(A);
453 Polyhedron_Free(C);
454 enumerate_options_free(options);
455 return result;
456 }
lattice point enumeration library