?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
lexmin: optionally use cdd during basis reduction
[barvinok.git] / sample.c
blobe63a9e06594f2a482680fe740536c30d6a79d3f4
1 #include <polylib/polylibgmp.h>
2 #include <barvinok/util.h>
3 #include <barvinok/basis_reduction.h>
4 #include <barvinok/sample.h>
5 #include <barvinok/options.h>
6
7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
8
9 /* If P has no rays, then we return NULL.
10 * Otherwise, look for the coordinate axis with the smallest maximal non-zero
11 * coefficient over all rays and a constraint that bounds the values on
12 * this axis to the maximal value over the vertices plus the above maximal
13 * non-zero coefficient times the number of rays minus 1.
14 * Any integer point x outside this region is the sum of a point inside
15 * the region and an integer multiple of the rays.
16 * Write x = \sum_i a_i v_i + \sum_j b_j r_j
17 * with \sum_i a_i = 1.
18 * Then x = \sum_i a_i v_i + \sum_j {b_j} r_j + \sum_j [b_j] r_j
19 * with y = \sum_i a_i v_i + \sum_j {b_j} r_j a point inside the region.
20 */
21 static Polyhedron *remove_ray(Polyhedron *P, unsigned MaxRays)
22 {
23 int r = 0;
24 Vector *min, *max, *c;
25 int i;
26 Value s, v, tmp;
27 int pos;
28 Polyhedron *R;
29 int rays;
30
31 if (P->NbBid == 0)
32 for (; r < P->NbRays; ++r)
33 if (value_zero_p(P->Ray[r][P->Dimension+1]))
34 break;
35 if (P->NbBid == 0 && r == P->NbRays)
36 return NULL;
37
38 max = Vector_Alloc(P->Dimension);
39 min = Vector_Alloc(P->Dimension);
40 for (r = 0; r < P->NbBid; ++r)
41 for (i = 0 ; i < P->Dimension; ++i)
42 if (value_abs_gt(P->Ray[r][1+i], max->p[i]))
43 value_absolute(max->p[i], P->Ray[r][1+i]);
44
45 for (i = 0 ; i < P->Dimension; ++i)
46 value_oppose(min->p[i], max->p[i]);
47
48 rays = P->NbBid;
49 for (r = P->NbBid; r < P->NbRays; ++r) {
50 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
51 continue;
52 for (i = 0 ; i < P->Dimension; ++i) {
53 if (value_gt(P->Ray[r][1+i], max->p[i]))
54 value_assign(max->p[i], P->Ray[r][1+i]);
55 if (value_lt(P->Ray[r][1+i], min->p[i]))
56 value_assign(min->p[i], P->Ray[r][1+i]);
57 }
58 ++rays;
59 }
60
61 value_init(s);
62 value_init(v);
63 value_init(tmp);
64
65 for (i = 0 ; i < P->Dimension; ++i) {
66 if (value_notzero_p(min->p[i]) &&
67 (value_zero_p(s) || value_abs_lt(min->p[i], s))) {
68 value_assign(s, min->p[i]);
69 pos = i;
70 }
71 if (value_notzero_p(max->p[i]) &&
72 (value_zero_p(s) || value_abs_lt(max->p[i], s))) {
73 value_assign(s, max->p[i]);
74 pos = i;
75 }
76 }
77
78 for (r = P->NbBid; r < P->NbRays; ++r)
79 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
80 break;
81
82 if (value_pos_p(s))
83 mpz_cdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
84 else
85 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
86
87 for ( ; r < P->NbRays; ++r) {
88 if (value_zero_p(P->Ray[r][P->Dimension+1]))
89 continue;
90
91 if (value_pos_p(s)) {
92 mpz_cdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
93 if (value_gt(tmp, v))
94 value_assign(v, tmp);
95 } else {
96 mpz_fdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
97 if (value_lt(tmp, v))
98 value_assign(v, tmp);
99 }
100 }
101
102 c = Vector_Alloc(1+P->Dimension+1);
103
104 value_set_si(tmp, rays);
105 value_addmul(v, tmp, s);
106 value_set_si(c->p[0], 1);
107 if (value_pos_p(s)) {
108 value_set_si(c->p[1+pos], -1);
109 value_assign(c->p[1+P->Dimension], v);
110 } else {
111 value_set_si(c->p[1+pos], 1);
112 value_oppose(c->p[1+P->Dimension], v);
113 }
114 value_decrement(c->p[1+P->Dimension], c->p[1+P->Dimension]);
115
116 R = AddConstraints(c->p, 1, P, MaxRays);
117
118 Vector_Free(c);
119
120 Vector_Free(min);
121 Vector_Free(max);
122
123 value_clear(tmp);
124 value_clear(s);
125 value_clear(v);
126
127 return R;
128 }
129
130 static void print_minmax(Polyhedron *P)
131 {
132 int i, j;
133 POL_ENSURE_VERTICES(P);
134 Polyhedron_Print(stderr, P_VALUE_FMT, P);
135 for (i = 0; i < P->Dimension; ++i) {
136 Value min, max, tmp;
137 value_init(min);
138 value_init(max);
139 value_init(tmp);
140
141 mpz_cdiv_q(min, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
142 mpz_fdiv_q(max, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
143
144 for (j = 1; j < P->NbRays; ++j) {
145 mpz_cdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
146 if (value_lt(tmp, min))
147 value_assign(min, tmp);
148 mpz_fdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
149 if (value_gt(tmp, max))
150 value_assign(max, tmp);
151 }
152 fprintf(stderr, "i: %d, min: ", i);
153 value_print(stderr, VALUE_FMT, min);
154 fprintf(stderr, ", max: ");
155 value_print(stderr, VALUE_FMT, max);
156 fprintf(stderr, "\n");
157
158 value_clear(min);
159 value_clear(max);
160 value_clear(tmp);
161 }
162 }
163
164 /* Remove coordinates that have a fixed value and return the matrix
165 * that adds these fixed coordinates again through T.
166 */
167 static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T)
168 {
169 int i, j, k, n;
170 int dim = P->Dimension;
171 int *remove = ALLOCN(int, dim);
172 Polyhedron *Q;
173 int NbEq;
174
175 assert(POL_HAS(P, POL_INEQUALITIES));
176 for (i = 0; i < dim; ++i)
177 remove[i] = -1;
178 NbEq = 0;
179 for (i = 0; i < P->NbEq; ++i) {
180 int pos = First_Non_Zero(P->Constraint[i]+1, dim);
181 if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) != -1)
182 continue;
183 remove[pos] = i;
184 ++NbEq;
185 }
186 assert(NbEq > 0);
187 Q = Polyhedron_Alloc(P->Dimension-NbEq, P->NbConstraints-NbEq, P->NbRays);
188 for (i = 0, k = 0; i < P->NbConstraints; ++i) {
189 if (i < P->NbEq) {
190 int pos = First_Non_Zero(P->Constraint[i]+1, dim);
191 if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) == -1)
192 continue;
193 }
194 value_assign(Q->Constraint[k][0], P->Constraint[i][0]);
195 for (j = 0, n = 0; j < P->Dimension; ++j) {
196 if (remove[j] != -1)
197 ++n;
198 else
199 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
200 }
201 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
202 ++k;
203 }
204 for (i = 0; i < Q->NbRays; ++i) {
205 value_assign(Q->Ray[i][0], P->Ray[i][0]);
206 for (j = 0, n = 0; j < P->Dimension; ++j) {
207 if (remove[j] != -1)
208 ++n;
209 else
210 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
211 }
212 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
213 }
214 *T = Matrix_Alloc(P->Dimension+1, Q->Dimension+1);
215 for (i = 0, n = 0; i < P->Dimension; ++i) {
216 if (remove[i] != -1) {
217 value_oppose((*T)->p[i][Q->Dimension],
218 P->Constraint[remove[i]][1+P->Dimension]);
219 ++n;
220 } else
221 value_set_si((*T)->p[i][i-n], 1);
222 }
223 value_set_si((*T)->p[i][i-n], 1);
224 POL_SET(Q, POL_VALID);
225 if (POL_HAS(P, POL_INEQUALITIES))
226 POL_SET(Q, POL_INEQUALITIES);
227 if (POL_HAS(P, POL_FACETS))
228 POL_SET(Q, POL_FACETS);
229 if (POL_HAS(P, POL_POINTS))
230 POL_SET(Q, POL_POINTS);
231 if (POL_HAS(P, POL_VERTICES))
232 POL_SET(Q, POL_VERTICES);
233 free(remove);
234 return Q;
235 }
236
237 /* This function implements the algorithm described in
238 * "An Implementation of the Generalized Basis Reduction Algorithm
239 * for Integer Programming" of Cook el al. to find an integer point
240 * in a polyhedron.
241 * If the polyhedron is unbounded, we first remove its rays.
242 */
243 Vector *Polyhedron_Sample(Polyhedron *P, struct barvinok_options *options)
244 {
245 int i, j;
246 Vector *sample = NULL;
247 Polyhedron *Q;
248 Matrix *T, *inv, *M;
249 Value min, max, tmp;
250 Vector *v;
251 int ok;
252
253 POL_ENSURE_VERTICES(P);
254 if (emptyQ(P))
255 return NULL;
256
257 if (P->Dimension == 0) {
258 sample = Vector_Alloc(1);
259 value_set_si(sample->p[0], 1);
260 return sample;
261 }
262
263 for (i = 0; i < P->NbRays; ++i)
264 if (value_one_p(P->Ray[i][1+P->Dimension])) {
265 sample = Vector_Alloc(P->Dimension+1);
266 Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
267 return sample;
268 }
269
270 Q = remove_ray(P, options->MaxRays);
271 if (Q) {
272 sample = Polyhedron_Sample(Q, options);
273 Polyhedron_Free(Q);
274 return sample;
275 }
276
277 Matrix *basis = Polyhedron_Reduced_Basis(P, options);
278
279 T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
280 inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
281 for (i = 0; i < P->Dimension; ++i)
282 for (j = 0; j < P->Dimension; ++j)
283 value_assign(T->p[i][j], basis->p[i][j]);
284 value_set_si(T->p[P->Dimension][P->Dimension], 1);
285 Matrix_Free(basis);
286
287 M = Matrix_Copy(T);
288 ok = Matrix_Inverse(M, inv);
289 assert(ok);
290 Matrix_Free(M);
291
292 Q = Polyhedron_Image(P, T, options->MaxRays);
293
294 POL_ENSURE_VERTICES(Q);
295
296 value_init(min);
297 value_init(max);
298 value_init(tmp);
299
300 mpz_cdiv_q(min, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
301 mpz_fdiv_q(max, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
302
303 for (j = 1; j < Q->NbRays; ++j) {
304 mpz_cdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
305 if (value_lt(tmp, min))
306 value_assign(min, tmp);
307 mpz_fdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
308 if (value_gt(tmp, max))
309 value_assign(max, tmp);
310 }
311
312 v = Vector_Alloc(1+Q->Dimension+1);
313 value_set_si(v->p[1], -1);
314
315 for (value_assign(tmp, min); value_le(tmp, max); value_increment(tmp, tmp)) {
316 Polyhedron *R, *S;
317 Matrix *T;
318 Vector *S_sample;
319 value_assign(v->p[1+Q->Dimension], tmp);
320
321 R = AddConstraints(v->p, 1, Q, options->MaxRays);
322 R = DomainConstraintSimplify(R, options->MaxRays);
323 if (emptyQ(R)) {
324 Polyhedron_Free(R);
325 continue;
326 }
327
328 S = Polyhedron_RemoveFixedColumns(R, &T);
329 Polyhedron_Free(R);
330 S_sample = Polyhedron_Sample(S, options);
331 Polyhedron_Free(S);
332 if (S_sample) {
333 Vector *Q_sample = Vector_Alloc(Q->Dimension + 1);
334 Matrix_Vector_Product(T, S_sample->p, Q_sample->p);
335 Matrix_Free(T);
336 Vector_Free(S_sample);
337 sample = Vector_Alloc(P->Dimension + 1);
338 Matrix_Vector_Product(inv, Q_sample->p, sample->p);
339 Vector_Free(Q_sample);
340 break;
341 }
342 Matrix_Free(T);
343 }
344
345 Matrix_Free(T);
346 Matrix_Free(inv);
347 Polyhedron_Free(Q);
348 Vector_Free(v);
349
350 value_clear(min);
351 value_clear(max);
352 value_clear(tmp);
353
354 return sample;
355 }
lattice point enumeration library