Polyhedron_Sample: handle empty and 0D polyhedra
[barvinok.git] / sample.c
blob8305a2e9475cb6ae9adfcfef9e26fc6fff625da9
1 #include <polylib/polylibgmp.h>
2 #include <barvinok/util.h>
3 #include "basis_reduction.h"
4 #include "sample.h"
6 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
8 /* If P has no rays, then we return NULL.
9 * Otherwise, look for the coordinate axis with the smallest maximal non-zero
10 * coefficient over all rays and a constraint that bounds the values on
11 * this axis to the maximal value over the vertices plus the above maximal
12 * non-zero coefficient minus 1.
13 * Any integer point outside this region should be the sum of a point inside
14 * the region and an integer multiple of the rays.
16 static Polyhedron *remove_ray(Polyhedron *P, unsigned MaxRays)
18 int r = 0;
19 Vector *min, *max, *c;
20 int i;
21 Value s, v, tmp;
22 int pos;
23 Polyhedron *R;
25 if (P->NbBid == 0)
26 for (; r < P->NbRays; ++r)
27 if (value_zero_p(P->Ray[r][P->Dimension+1]))
28 break;
29 if (P->NbBid == 0 && r == P->NbRays)
30 return NULL;
32 max = Vector_Alloc(P->Dimension);
33 min = Vector_Alloc(P->Dimension);
34 for (r = 0; r < P->NbBid; ++r)
35 for (i = 0 ; i < P->Dimension; ++i)
36 if (value_abs_gt(P->Ray[r][1+i], max->p[i]))
37 value_absolute(max->p[i], P->Ray[r][1+i]);
39 for (i = 0 ; i < P->Dimension; ++i)
40 value_oppose(min->p[i], max->p[i]);
42 for (r = P->NbBid; r < P->NbRays; ++r) {
43 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
44 continue;
45 for (i = 0 ; i < P->Dimension; ++i) {
46 if (value_gt(P->Ray[r][1+i], max->p[i]))
47 value_assign(max->p[i], P->Ray[r][1+i]);
48 if (value_lt(P->Ray[r][1+i], min->p[i]))
49 value_assign(min->p[i], P->Ray[r][1+i]);
53 value_init(s);
54 value_init(v);
55 value_init(tmp);
57 for (i = 0 ; i < P->Dimension; ++i) {
58 if (value_notzero_p(min->p[i]) &&
59 (value_zero_p(s) || value_abs_lt(min->p[i], s))) {
60 value_assign(s, min->p[i]);
61 pos = i;
63 if (value_notzero_p(max->p[i]) &&
64 (value_zero_p(s) || value_abs_lt(max->p[i], s))) {
65 value_assign(s, max->p[i]);
66 pos = i;
70 for (r = P->NbBid; r < P->NbRays; ++r)
71 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
72 break;
74 if (value_pos_p(s))
75 mpz_cdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
76 else
77 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
79 for ( ; r < P->NbRays; ++r) {
80 if (value_zero_p(P->Ray[r][P->Dimension+1]))
81 continue;
83 if (value_pos_p(s)) {
84 mpz_cdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
85 if (value_gt(tmp, v))
86 value_assign(v, tmp);
87 } else {
88 mpz_fdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
89 if (value_lt(tmp, v))
90 value_assign(v, tmp);
94 c = Vector_Alloc(1+P->Dimension+1);
96 value_addto(v, v, s);
97 value_set_si(c->p[0], 1);
98 if (value_pos_p(s)) {
99 value_set_si(c->p[1+pos], -1);
100 value_assign(c->p[1+P->Dimension], v);
101 } else {
102 value_set_si(c->p[1+pos], 1);
103 value_oppose(c->p[1+P->Dimension], v);
105 value_decrement(c->p[1+P->Dimension], c->p[1+P->Dimension]);
107 R = AddConstraints(c->p, 1, P, MaxRays);
109 Vector_Free(c);
111 Vector_Free(min);
112 Vector_Free(max);
114 value_clear(tmp);
115 value_clear(s);
116 value_clear(v);
118 return R;
121 static void print_minmax(Polyhedron *P)
123 int i, j;
124 POL_ENSURE_VERTICES(P);
125 Polyhedron_Print(stderr, P_VALUE_FMT, P);
126 for (i = 0; i < P->Dimension; ++i) {
127 Value min, max, tmp;
128 value_init(min);
129 value_init(max);
130 value_init(tmp);
132 mpz_cdiv_q(min, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
133 mpz_fdiv_q(max, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
135 for (j = 1; j < P->NbRays; ++j) {
136 mpz_cdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
137 if (value_lt(tmp, min))
138 value_assign(min, tmp);
139 mpz_fdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
140 if (value_gt(tmp, max))
141 value_assign(max, tmp);
143 fprintf(stderr, "i: %d, min: ", i);
144 value_print(stderr, VALUE_FMT, min);
145 fprintf(stderr, ", max: ");
146 value_print(stderr, VALUE_FMT, max);
147 fprintf(stderr, "\n");
149 value_clear(min);
150 value_clear(max);
151 value_clear(tmp);
155 /* Remove coordinates that have a fixed value and return the matrix
156 * that adds these fixed coordinates again through T.
158 static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T)
160 int i, j, k, n;
161 int dim = P->Dimension;
162 int *remove = ALLOCN(int, dim);
163 Polyhedron *Q;
164 int NbEq;
166 assert(POL_HAS(P, POL_INEQUALITIES));
167 for (i = 0; i < dim; ++i)
168 remove[i] = 0;
169 NbEq = 0;
170 for (i = 0; i < P->NbEq; ++i) {
171 int pos = First_Non_Zero(P->Constraint[i]+1, dim);
172 if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) != -1)
173 continue;
174 remove[pos] = 1;
175 ++NbEq;
177 assert(NbEq > 0);
178 Q = Polyhedron_Alloc(P->Dimension-NbEq, P->NbConstraints-NbEq, P->NbRays);
179 for (i = 0, k = 0; i < P->NbConstraints; ++i) {
180 if (i < P->NbEq) {
181 int pos = First_Non_Zero(P->Constraint[i]+1, dim);
182 if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) == -1)
183 continue;
185 value_assign(Q->Constraint[k][0], P->Constraint[i][0]);
186 for (j = 0, n = 0; j < P->Dimension; ++j) {
187 if (remove[j])
188 ++n;
189 else
190 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
192 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
193 ++k;
195 for (i = 0; i < Q->NbRays; ++i) {
196 value_assign(Q->Ray[i][0], P->Ray[i][0]);
197 for (j = 0, n = 0; j < P->Dimension; ++j) {
198 if (remove[j])
199 ++n;
200 else
201 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
203 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
205 *T = Matrix_Alloc(P->Dimension+1, Q->Dimension+1);
206 for (i = 0, n = 0; i < P->Dimension; ++i) {
207 if (remove[i]) {
208 value_oppose((*T)->p[i][Q->Dimension], P->Constraint[n][1+P->Dimension]);
209 ++n;
210 } else
211 value_set_si((*T)->p[i][i-n], 1);
213 value_set_si((*T)->p[i][i-n], 1);
214 POL_SET(Q, POL_VALID);
215 if (POL_HAS(P, POL_INEQUALITIES))
216 POL_SET(Q, POL_INEQUALITIES);
217 if (POL_HAS(P, POL_FACETS))
218 POL_SET(Q, POL_FACETS);
219 if (POL_HAS(P, POL_POINTS))
220 POL_SET(Q, POL_POINTS);
221 if (POL_HAS(P, POL_VERTICES))
222 POL_SET(Q, POL_VERTICES);
223 free(remove);
224 return Q;
227 /* This function implements the algorithm described in
228 * "An Implementation of the Generalized Basis Reduction Algorithm
229 * for Integer Programming" of Cook el al. to find an integer point
230 * in a polyhedron.
231 * If the polyhedron is unbounded, we first remove its rays.
233 Vector *Polyhedron_Sample(Polyhedron *P, unsigned MaxRays)
235 int i, j;
236 Vector *sample = NULL;
237 Polyhedron *Q;
238 Matrix *T, *inv, *M;
239 Value min, max, tmp;
240 Vector *v;
241 int ok;
243 POL_ENSURE_VERTICES(P);
244 if (emptyQ(P))
245 return NULL;
247 if (P->Dimension == 0) {
248 sample = Vector_Alloc(1);
249 value_set_si(sample->p[0], 1);
250 return sample;
253 for (i = 0; i < P->NbRays; ++i)
254 if (value_one_p(P->Ray[i][1+P->Dimension])) {
255 sample = Vector_Alloc(P->Dimension+1);
256 Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
257 return sample;
260 /* for now */
261 assert(P->NbEq == 0);
263 Q = remove_ray(P, MaxRays);
264 if (Q) {
265 sample = Polyhedron_Sample(Q, MaxRays);
266 Polyhedron_Free(Q);
267 return sample;
270 Matrix *basis = reduced_basis(P);
272 T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
273 inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
274 for (i = 0; i < P->Dimension; ++i)
275 for (j = 0; j < P->Dimension; ++j)
276 value_assign(T->p[i][j], basis->p[i][j]);
277 value_set_si(T->p[P->Dimension][P->Dimension], 1);
278 Matrix_Free(basis);
280 M = Matrix_Copy(T);
281 ok = Matrix_Inverse(M, inv);
282 assert(ok);
283 Matrix_Free(M);
285 Q = Polyhedron_Image(P, T, MaxRays);
287 value_init(min);
288 value_init(max);
289 value_init(tmp);
291 mpz_cdiv_q(min, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
292 mpz_fdiv_q(max, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
294 for (j = 1; j < Q->NbRays; ++j) {
295 mpz_cdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
296 if (value_lt(tmp, min))
297 value_assign(min, tmp);
298 mpz_fdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
299 if (value_gt(tmp, max))
300 value_assign(max, tmp);
303 v = Vector_Alloc(1+Q->Dimension+1);
304 value_set_si(v->p[1], -1);
306 for (value_assign(tmp, min); value_le(tmp, max); value_increment(tmp, tmp)) {
307 Polyhedron *R, *S;
308 Matrix *T;
309 Vector *S_sample;
310 value_assign(v->p[1+Q->Dimension], tmp);
312 R = AddConstraints(v->p, 1, Q, MaxRays);
313 S = Polyhedron_RemoveFixedColumns(R, &T);
314 Polyhedron_Free(R);
315 S_sample = Polyhedron_Sample(S, MaxRays);
316 Polyhedron_Free(S);
317 if (S_sample) {
318 Vector *Q_sample = Vector_Alloc(Q->Dimension + 1);
319 Matrix_Vector_Product(T, S_sample->p, Q_sample->p);
320 Matrix_Free(T);
321 Vector_Free(S_sample);
322 sample = Vector_Alloc(P->Dimension + 1);
323 Matrix_Vector_Product(inv, Q_sample->p, sample->p);
324 Vector_Free(Q_sample);
325 break;
327 Matrix_Free(T);
330 Matrix_Free(T);
331 Matrix_Free(inv);
332 Polyhedron_Free(Q);
333 Vector_Free(v);
335 value_clear(min);
336 value_clear(max);
337 value_clear(tmp);
339 return sample;