?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
configure.in: perform check for GMP support in NTL in C++ context
[barvinok.git] / sample.c
blobc98dcd8065c3bca7947efd6bec152698b2a9b8de
1 #include <assert.h>
2 #include <isl_set_polylib.h>
3 #include <barvinok/util.h>
4 #include <barvinok/basis_reduction.h>
5 #include <barvinok/sample.h>
6 #include <barvinok/options.h>
7 #include "polysign.h"
8
9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
10
11 /* If P has no rays, then we return NULL.
12 * Otherwise, look for the coordinate axis with the smallest maximal non-zero
13 * coefficient over all rays and a constraint that bounds the values on
14 * this axis to the maximal value over the vertices plus the above maximal
15 * non-zero coefficient times the number of rays minus 1.
16 * Any integer point x outside this region is the sum of a point inside
17 * the region and an integer multiple of the rays.
18 * Write x = \sum_i a_i v_i + \sum_j b_j r_j
19 * with \sum_i a_i = 1.
20 * Then x = \sum_i a_i v_i + \sum_j {b_j} r_j + \sum_j [b_j] r_j
21 * with y = \sum_i a_i v_i + \sum_j {b_j} r_j a point inside the region.
22 */
23 static Polyhedron *remove_ray(Polyhedron *P, unsigned MaxRays)
24 {
25 int r = 0;
26 Vector *min, *max, *c;
27 int i;
28 Value s, v, tmp;
29 int pos;
30 Polyhedron *R;
31 int rays;
32
33 POL_ENSURE_VERTICES(P);
34 if (P->NbBid == 0)
35 for (; r < P->NbRays; ++r)
36 if (value_zero_p(P->Ray[r][P->Dimension+1]))
37 break;
38 if (P->NbBid == 0 && r == P->NbRays)
39 return NULL;
40
41 max = Vector_Alloc(P->Dimension);
42 min = Vector_Alloc(P->Dimension);
43 for (r = 0; r < P->NbBid; ++r)
44 for (i = 0 ; i < P->Dimension; ++i)
45 if (value_abs_gt(P->Ray[r][1+i], max->p[i]))
46 value_absolute(max->p[i], P->Ray[r][1+i]);
47
48 for (i = 0 ; i < P->Dimension; ++i)
49 value_oppose(min->p[i], max->p[i]);
50
51 rays = P->NbBid;
52 for (r = P->NbBid; r < P->NbRays; ++r) {
53 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
54 continue;
55 for (i = 0 ; i < P->Dimension; ++i) {
56 if (value_gt(P->Ray[r][1+i], max->p[i]))
57 value_assign(max->p[i], P->Ray[r][1+i]);
58 if (value_lt(P->Ray[r][1+i], min->p[i]))
59 value_assign(min->p[i], P->Ray[r][1+i]);
60 }
61 ++rays;
62 }
63
64 value_init(s);
65 value_init(v);
66 value_init(tmp);
67
68 for (i = 0 ; i < P->Dimension; ++i) {
69 if (value_notzero_p(min->p[i]) &&
70 (value_zero_p(s) || value_abs_lt(min->p[i], s))) {
71 value_assign(s, min->p[i]);
72 pos = i;
73 }
74 if (value_notzero_p(max->p[i]) &&
75 (value_zero_p(s) || value_abs_lt(max->p[i], s))) {
76 value_assign(s, max->p[i]);
77 pos = i;
78 }
79 }
80
81 for (r = P->NbBid; r < P->NbRays; ++r)
82 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
83 break;
84
85 if (value_pos_p(s))
86 mpz_cdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
87 else
88 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
89
90 for ( ; r < P->NbRays; ++r) {
91 if (value_zero_p(P->Ray[r][P->Dimension+1]))
92 continue;
93
94 if (value_pos_p(s)) {
95 mpz_cdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
96 if (value_gt(tmp, v))
97 value_assign(v, tmp);
98 } else {
99 mpz_fdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
100 if (value_lt(tmp, v))
101 value_assign(v, tmp);
102 }
103 }
104
105 c = Vector_Alloc(1+P->Dimension+1);
106
107 value_set_si(tmp, rays);
108 value_addmul(v, tmp, s);
109 value_set_si(c->p[0], 1);
110 if (value_pos_p(s)) {
111 value_set_si(c->p[1+pos], -1);
112 value_assign(c->p[1+P->Dimension], v);
113 } else {
114 value_set_si(c->p[1+pos], 1);
115 value_oppose(c->p[1+P->Dimension], v);
116 }
117 value_decrement(c->p[1+P->Dimension], c->p[1+P->Dimension]);
118
119 R = AddConstraints(c->p, 1, P, MaxRays);
120
121 Vector_Free(c);
122
123 Vector_Free(min);
124 Vector_Free(max);
125
126 value_clear(tmp);
127 value_clear(s);
128 value_clear(v);
129
130 return R;
131 }
132
133 static void print_minmax(Polyhedron *P)
134 {
135 int i, j;
136 POL_ENSURE_VERTICES(P);
137 Polyhedron_Print(stderr, P_VALUE_FMT, P);
138 for (i = 0; i < P->Dimension; ++i) {
139 Value min, max, tmp;
140 value_init(min);
141 value_init(max);
142 value_init(tmp);
143
144 mpz_cdiv_q(min, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
145 mpz_fdiv_q(max, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
146
147 for (j = 1; j < P->NbRays; ++j) {
148 mpz_cdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
149 if (value_lt(tmp, min))
150 value_assign(min, tmp);
151 mpz_fdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
152 if (value_gt(tmp, max))
153 value_assign(max, tmp);
154 }
155 fprintf(stderr, "i: %d, min: ", i);
156 value_print(stderr, VALUE_FMT, min);
157 fprintf(stderr, ", max: ");
158 value_print(stderr, VALUE_FMT, max);
159 fprintf(stderr, "\n");
160
161 value_clear(min);
162 value_clear(max);
163 value_clear(tmp);
164 }
165 }
166
167 /* Remove coordinates that have a fixed value and return the matrix
168 * that adds these fixed coordinates again through T.
169 */
170 static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T)
171 {
172 int i, j, k, n;
173 int dim = P->Dimension;
174 int *remove = ALLOCN(int, dim);
175 Polyhedron *Q;
176 int NbEq;
177
178 assert(POL_HAS(P, POL_INEQUALITIES));
179 for (i = 0; i < dim; ++i)
180 remove[i] = -1;
181 NbEq = 0;
182 for (i = 0; i < P->NbEq; ++i) {
183 int pos = First_Non_Zero(P->Constraint[i]+1, dim);
184 if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) != -1)
185 continue;
186 remove[pos] = i;
187 ++NbEq;
188 }
189 assert(NbEq > 0);
190 Q = Polyhedron_Alloc(P->Dimension-NbEq, P->NbConstraints-NbEq, P->NbRays);
191 Q->NbEq = P->NbEq - NbEq;
192 for (i = 0, k = 0; i < P->NbConstraints; ++i) {
193 if (i < P->NbEq) {
194 int pos = First_Non_Zero(P->Constraint[i]+1, dim);
195 if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) == -1)
196 continue;
197 }
198 value_assign(Q->Constraint[k][0], P->Constraint[i][0]);
199 for (j = 0, n = 0; j < P->Dimension; ++j) {
200 if (remove[j] != -1)
201 ++n;
202 else
203 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
204 }
205 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
206 ++k;
207 }
208 for (i = 0; i < Q->NbRays; ++i) {
209 value_assign(Q->Ray[i][0], P->Ray[i][0]);
210 for (j = 0, n = 0; j < P->Dimension; ++j) {
211 if (remove[j] != -1)
212 ++n;
213 else
214 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
215 }
216 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
217 }
218 *T = Matrix_Alloc(P->Dimension+1, Q->Dimension+1);
219 for (i = 0, n = 0; i < P->Dimension; ++i) {
220 if (remove[i] != -1) {
221 value_oppose((*T)->p[i][Q->Dimension],
222 P->Constraint[remove[i]][1+P->Dimension]);
223 ++n;
224 } else
225 value_set_si((*T)->p[i][i-n], 1);
226 }
227 value_set_si((*T)->p[i][i-n], 1);
228 POL_SET(Q, POL_VALID);
229 if (POL_HAS(P, POL_INEQUALITIES))
230 POL_SET(Q, POL_INEQUALITIES);
231 if (POL_HAS(P, POL_FACETS))
232 POL_SET(Q, POL_FACETS);
233 if (POL_HAS(P, POL_POINTS))
234 POL_SET(Q, POL_POINTS);
235 if (POL_HAS(P, POL_VERTICES))
236 POL_SET(Q, POL_VERTICES);
237 free(remove);
238 return Q;
239 }
240
241 static Polyhedron *remove_all_equalities(Polyhedron *P, Matrix **T,
242 unsigned MaxRays)
243 {
244 Matrix M;
245 Polyhedron_Matrix_View(P, &M, P->NbEq);
246
247 *T = compress_variables(&M, 0);
248
249 if (!*T)
250 return NULL;
251 P = Polyhedron_Preimage(P, *T, MaxRays);
252
253 return P;
254 }
255
256 static Vector *product_sample(Polyhedron *P, Matrix *T,
257 struct barvinok_options *options)
258 {
259 int i;
260 Vector *sample = NULL;
261 Vector *tmp = Vector_Alloc(T->NbRows);
262 i = 0;
263 for (; P; P = P->next) {
264 Vector *P_sample;
265 Polyhedron *next = P->next;
266 P->next = NULL;
267 P_sample = Polyhedron_Sample(P, options);
268 P->next = next;
269 if (!P_sample) {
270 Vector_Free(tmp);
271 tmp = NULL;
272 break;
273 }
274 Vector_Copy(P_sample->p, tmp->p+i, P->Dimension);
275 Vector_Free(P_sample);
276 i += P->Dimension;
277 }
278 if (tmp) {
279 sample = Vector_Alloc(T->NbRows + 1);
280 Matrix_Vector_Product(T, tmp->p, sample->p);
281 value_set_si(sample->p[T->NbRows], 1);
282 Vector_Free(tmp);
283 }
284 return sample;
285 }
286
287 static Vector *isl_Polyhedron_Sample(Polyhedron *P,
288 struct barvinok_options *options)
289 {
290 int i;
291 isl_ctx *ctx = isl_ctx_alloc();
292 isl_dim *dim;
293 int nvar = P->Dimension;
294 isl_basic_set *bset;
295 isl_point *pnt;
296 Vector *sample = NULL;
297
298 dim = isl_dim_set_alloc(ctx, 0, nvar);
299 bset = isl_basic_set_new_from_polylib(P, dim);
300 pnt = isl_basic_set_sample_point(bset);
301
302 if (!isl_point_is_void(pnt)) {
303 isl_int v;
304
305 isl_int_init(v);
306 sample = Vector_Alloc(1 + nvar);
307 assert(sample);
308 for (i = 0; i < nvar; ++i) {
309 isl_point_get_coordinate(pnt, isl_dim_set, i, &v);
310 isl_int_get_gmp(v, sample->p[i]);
311 }
312 value_set_si(sample->p[nvar], 1);
313 isl_int_clear(v);
314 }
315
316 isl_point_free(pnt);
317
318 isl_ctx_free(ctx);
319
320 return sample;
321 }
322
323 /* This function implements the algorithm described in
324 * "An Implementation of the Generalized Basis Reduction Algorithm
325 * for Integer Programming" of Cook el al. to find an integer point
326 * in a polyhedron.
327 * If the polyhedron is unbounded, we first remove its rays.
328 */
329 Vector *Polyhedron_Sample(Polyhedron *P, struct barvinok_options *options)
330 {
331 int i, j;
332 Vector *sample = NULL, *obj = NULL;
333 Polyhedron *Q;
334 Matrix *inv = NULL;
335 Value min, max, tmp;
336 Vector *v;
337 int ok;
338 enum lp_result res;
339 Matrix *T;
340
341 if (options->gbr_lp_solver == BV_GBR_ISL)
342 return isl_Polyhedron_Sample(P, options);
343
344 if (emptyQ2(P))
345 return NULL;
346
347 if (P->Dimension == 0) {
348 sample = Vector_Alloc(1);
349 value_set_si(sample->p[0], 1);
350 return sample;
351 }
352
353 if (P->Dimension == 1)
354 POL_ENSURE_VERTICES(P);
355
356 for (i = 0; i < P->NbRays; ++i)
357 if (value_one_p(P->Ray[i][1+P->Dimension])) {
358 sample = Vector_Alloc(P->Dimension+1);
359 Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
360 return sample;
361 }
362
363 if (P->NbEq > 0) {
364 Vector *Q_sample;
365 Polyhedron *Q = remove_all_equalities(P, &T, options->MaxRays);
366 if (!Q)
367 return NULL;
368 Q_sample = Polyhedron_Sample(Q, options);
369 Polyhedron_Free(Q);
370 if (Q_sample) {
371 sample = Vector_Alloc(P->Dimension + 1);
372 Matrix_Vector_Product(T, Q_sample->p, sample->p);
373 Vector_Free(Q_sample);
374 }
375 Matrix_Free(T);
376 return sample;
377 }
378
379 Q = Polyhedron_Factor(P, 0, &T, options->MaxRays);
380 if (Q) {
381 sample = product_sample(Q, T, options);
382 Domain_Free(Q);
383 Matrix_Free(T);
384 return sample;
385 }
386
387 value_init(min);
388 value_init(max);
389
390 obj = Vector_Alloc(P->Dimension+1);
391 value_set_si(obj->p[0], 1);
392 res = polyhedron_range(P, obj->p, obj->p[0], &min, &max, options);
393 if (res == lp_unbounded) {
394 unbounded:
395 value_clear(min);
396 value_clear(max);
397 Vector_Free(obj);
398
399 Q = remove_ray(P, options->MaxRays);
400 assert(Q);
401 sample = Polyhedron_Sample(Q, options);
402 Polyhedron_Free(Q);
403 return sample;
404 }
405 if (res == lp_empty)
406 goto out;
407 assert(res == lp_ok);
408
409 if (value_eq(min, max)) {
410 Q = P;
411 } else {
412 Matrix *M, *T;
413 Matrix *basis;
414
415 options->gbr_only_first = 1;
416 basis = Polyhedron_Reduced_Basis(P, options);
417 options->gbr_only_first = 0;
418
419 if (!basis)
420 goto unbounded;
421 T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
422 inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
423 for (i = 0; i < P->Dimension; ++i)
424 for (j = 0; j < P->Dimension; ++j)
425 value_assign(T->p[i][j], basis->p[i][j]);
426 value_set_si(T->p[P->Dimension][P->Dimension], 1);
427 Matrix_Free(basis);
428
429 M = Matrix_Copy(T);
430 ok = Matrix_Inverse(M, inv);
431 assert(ok);
432 Matrix_Free(M);
433
434 Q = Polyhedron_Image(P, T, options->MaxRays);
435 res = polyhedron_range(Q, obj->p, obj->p[0], &min, &max, options);
436
437 Matrix_Free(T);
438 if (res == lp_empty)
439 goto out;
440 assert(res == lp_ok);
441 }
442
443 value_init(tmp);
444
445 v = Vector_Alloc(1+Q->Dimension+1);
446 value_set_si(v->p[1], -1);
447
448 for (value_assign(tmp, min); value_le(tmp, max); value_increment(tmp, tmp)) {
449 Polyhedron *R, *S;
450 Matrix *T;
451 Vector *S_sample;
452 value_assign(v->p[1+Q->Dimension], tmp);
453
454 R = AddConstraints(v->p, 1, Q, options->MaxRays);
455 R = DomainConstraintSimplify(R, options->MaxRays);
456 if (emptyQ(R)) {
457 Polyhedron_Free(R);
458 continue;
459 }
460
461 S = Polyhedron_RemoveFixedColumns(R, &T);
462 Polyhedron_Free(R);
463 S_sample = Polyhedron_Sample(S, options);
464 Polyhedron_Free(S);
465 if (S_sample) {
466 Vector *Q_sample = obj;
467 obj = NULL;
468 Matrix_Vector_Product(T, S_sample->p, Q_sample->p);
469 Matrix_Free(T);
470 Vector_Free(S_sample);
471 if (!inv)
472 sample = Q_sample;
473 else {
474 sample = Vector_Alloc(P->Dimension + 1);
475 Matrix_Vector_Product(inv, Q_sample->p, sample->p);
476 Vector_Free(Q_sample);
477 }
478 break;
479 }
480 Matrix_Free(T);
481 }
482 value_clear(tmp);
483
484 Vector_Free(v);
485
486 out:
487 if (inv)
488 Matrix_Free(inv);
489 if (P != Q)
490 Polyhedron_Free(Q);
491 if (obj)
492 Vector_Free(obj);
493 value_clear(min);
494 value_clear(max);
495
496 return sample;
497 }
lattice point enumeration library