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