?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
change type of some "sign"s to int as they can only be 1 or -1
[barvinok.git] / sample.c
blob191948677b23f5819368f1fa92cfa5f445f7c022
1 #include <barvinok/util.h>
2 #include <barvinok/basis_reduction.h>
3 #include <barvinok/sample.h>
4 #include <barvinok/options.h>
5 #include "polysign.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 POL_ENSURE_VERTICES(P);
32 if (P->NbBid == 0)
33 for (; r < P->NbRays; ++r)
34 if (value_zero_p(P->Ray[r][P->Dimension+1]))
35 break;
36 if (P->NbBid == 0 && r == P->NbRays)
37 return NULL;
38
39 max = Vector_Alloc(P->Dimension);
40 min = Vector_Alloc(P->Dimension);
41 for (r = 0; r < P->NbBid; ++r)
42 for (i = 0 ; i < P->Dimension; ++i)
43 if (value_abs_gt(P->Ray[r][1+i], max->p[i]))
44 value_absolute(max->p[i], P->Ray[r][1+i]);
45
46 for (i = 0 ; i < P->Dimension; ++i)
47 value_oppose(min->p[i], max->p[i]);
48
49 rays = P->NbBid;
50 for (r = P->NbBid; r < P->NbRays; ++r) {
51 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
52 continue;
53 for (i = 0 ; i < P->Dimension; ++i) {
54 if (value_gt(P->Ray[r][1+i], max->p[i]))
55 value_assign(max->p[i], P->Ray[r][1+i]);
56 if (value_lt(P->Ray[r][1+i], min->p[i]))
57 value_assign(min->p[i], P->Ray[r][1+i]);
58 }
59 ++rays;
60 }
61
62 value_init(s);
63 value_init(v);
64 value_init(tmp);
65
66 for (i = 0 ; i < P->Dimension; ++i) {
67 if (value_notzero_p(min->p[i]) &&
68 (value_zero_p(s) || value_abs_lt(min->p[i], s))) {
69 value_assign(s, min->p[i]);
70 pos = i;
71 }
72 if (value_notzero_p(max->p[i]) &&
73 (value_zero_p(s) || value_abs_lt(max->p[i], s))) {
74 value_assign(s, max->p[i]);
75 pos = i;
76 }
77 }
78
79 for (r = P->NbBid; r < P->NbRays; ++r)
80 if (value_notzero_p(P->Ray[r][P->Dimension+1]))
81 break;
82
83 if (value_pos_p(s))
84 mpz_cdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
85 else
86 mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
87
88 for ( ; r < P->NbRays; ++r) {
89 if (value_zero_p(P->Ray[r][P->Dimension+1]))
90 continue;
91
92 if (value_pos_p(s)) {
93 mpz_cdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
94 if (value_gt(tmp, v))
95 value_assign(v, tmp);
96 } else {
97 mpz_fdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
98 if (value_lt(tmp, v))
99 value_assign(v, tmp);
100 }
101 }
102
103 c = Vector_Alloc(1+P->Dimension+1);
104
105 value_set_si(tmp, rays);
106 value_addmul(v, tmp, s);
107 value_set_si(c->p[0], 1);
108 if (value_pos_p(s)) {
109 value_set_si(c->p[1+pos], -1);
110 value_assign(c->p[1+P->Dimension], v);
111 } else {
112 value_set_si(c->p[1+pos], 1);
113 value_oppose(c->p[1+P->Dimension], v);
114 }
115 value_decrement(c->p[1+P->Dimension], c->p[1+P->Dimension]);
116
117 R = AddConstraints(c->p, 1, P, MaxRays);
118
119 Vector_Free(c);
120
121 Vector_Free(min);
122 Vector_Free(max);
123
124 value_clear(tmp);
125 value_clear(s);
126 value_clear(v);
127
128 return R;
129 }
130
131 static void print_minmax(Polyhedron *P)
132 {
133 int i, j;
134 POL_ENSURE_VERTICES(P);
135 Polyhedron_Print(stderr, P_VALUE_FMT, P);
136 for (i = 0; i < P->Dimension; ++i) {
137 Value min, max, tmp;
138 value_init(min);
139 value_init(max);
140 value_init(tmp);
141
142 mpz_cdiv_q(min, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
143 mpz_fdiv_q(max, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
144
145 for (j = 1; j < P->NbRays; ++j) {
146 mpz_cdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
147 if (value_lt(tmp, min))
148 value_assign(min, tmp);
149 mpz_fdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
150 if (value_gt(tmp, max))
151 value_assign(max, tmp);
152 }
153 fprintf(stderr, "i: %d, min: ", i);
154 value_print(stderr, VALUE_FMT, min);
155 fprintf(stderr, ", max: ");
156 value_print(stderr, VALUE_FMT, max);
157 fprintf(stderr, "\n");
158
159 value_clear(min);
160 value_clear(max);
161 value_clear(tmp);
162 }
163 }
164
165 /* Remove coordinates that have a fixed value and return the matrix
166 * that adds these fixed coordinates again through T.
167 */
168 static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T)
169 {
170 int i, j, k, n;
171 int dim = P->Dimension;
172 int *remove = ALLOCN(int, dim);
173 Polyhedron *Q;
174 int NbEq;
175
176 assert(POL_HAS(P, POL_INEQUALITIES));
177 for (i = 0; i < dim; ++i)
178 remove[i] = -1;
179 NbEq = 0;
180 for (i = 0; i < P->NbEq; ++i) {
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;
184 remove[pos] = i;
185 ++NbEq;
186 }
187 assert(NbEq > 0);
188 Q = Polyhedron_Alloc(P->Dimension-NbEq, P->NbConstraints-NbEq, P->NbRays);
189 Q->NbEq = P->NbEq - NbEq;
190 for (i = 0, k = 0; i < P->NbConstraints; ++i) {
191 if (i < P->NbEq) {
192 int pos = First_Non_Zero(P->Constraint[i]+1, dim);
193 if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) == -1)
194 continue;
195 }
196 value_assign(Q->Constraint[k][0], P->Constraint[i][0]);
197 for (j = 0, n = 0; j < P->Dimension; ++j) {
198 if (remove[j] != -1)
199 ++n;
200 else
201 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
202 }
203 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
204 ++k;
205 }
206 for (i = 0; i < Q->NbRays; ++i) {
207 value_assign(Q->Ray[i][0], P->Ray[i][0]);
208 for (j = 0, n = 0; j < P->Dimension; ++j) {
209 if (remove[j] != -1)
210 ++n;
211 else
212 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
213 }
214 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
215 }
216 *T = Matrix_Alloc(P->Dimension+1, Q->Dimension+1);
217 for (i = 0, n = 0; i < P->Dimension; ++i) {
218 if (remove[i] != -1) {
219 value_oppose((*T)->p[i][Q->Dimension],
220 P->Constraint[remove[i]][1+P->Dimension]);
221 ++n;
222 } else
223 value_set_si((*T)->p[i][i-n], 1);
224 }
225 value_set_si((*T)->p[i][i-n], 1);
226 POL_SET(Q, POL_VALID);
227 if (POL_HAS(P, POL_INEQUALITIES))
228 POL_SET(Q, POL_INEQUALITIES);
229 if (POL_HAS(P, POL_FACETS))
230 POL_SET(Q, POL_FACETS);
231 if (POL_HAS(P, POL_POINTS))
232 POL_SET(Q, POL_POINTS);
233 if (POL_HAS(P, POL_VERTICES))
234 POL_SET(Q, POL_VERTICES);
235 free(remove);
236 return Q;
237 }
238
239 static Polyhedron *remove_all_equalities(Polyhedron *P, Matrix **T,
240 unsigned MaxRays)
241 {
242 /* Matrix "view" of equalities */
243 Matrix M;
244 M.NbRows = P->NbEq;
245 M.NbColumns = P->Dimension+2;
246 M.p_Init = P->p_Init;
247 M.p = P->Constraint;
248
249 *T = compress_variables(&M, 0);
250
251 if (!*T)
252 return NULL;
253 P = Polyhedron_Preimage(P, *T, MaxRays);
254
255 return P;
256 }
257
258 static Vector *product_sample(Polyhedron *P, Matrix *T,
259 struct barvinok_options *options)
260 {
261 int i;
262 Vector *sample = NULL;
263 Vector *tmp = Vector_Alloc(T->NbRows);
264 i = 0;
265 for (; P; P = P->next) {
266 Vector *P_sample;
267 Polyhedron *next = P->next;
268 P->next = NULL;
269 P_sample = Polyhedron_Sample(P, options);
270 if (!P_sample) {
271 Vector_Free(tmp);
272 tmp = NULL;
273 break;
274 }
275 Vector_Copy(P_sample->p, tmp->p+i, P->Dimension);
276 Vector_Free(P_sample);
277 i += P->Dimension;
278 P->next = next;
279 }
280 if (tmp) {
281 sample = Vector_Alloc(T->NbRows + 1);
282 Matrix_Vector_Product(T, tmp->p, sample->p);
283 value_set_si(sample->p[T->NbRows], 1);
284 Vector_Free(tmp);
285 }
286 return sample;
287 }
288
289 /* This function implements the algorithm described in
290 * "An Implementation of the Generalized Basis Reduction Algorithm
291 * for Integer Programming" of Cook el al. to find an integer point
292 * in a polyhedron.
293 * If the polyhedron is unbounded, we first remove its rays.
294 */
295 Vector *Polyhedron_Sample(Polyhedron *P, struct barvinok_options *options)
296 {
297 int i, j;
298 Vector *sample = NULL, *obj = NULL;
299 Polyhedron *Q;
300 Matrix *inv = NULL;
301 Value min, max, tmp;
302 Vector *v;
303 int ok;
304 enum lp_result res;
305 Matrix *T;
306
307 if (emptyQ2(P))
308 return NULL;
309
310 if (P->Dimension == 0) {
311 sample = Vector_Alloc(1);
312 value_set_si(sample->p[0], 1);
313 return sample;
314 }
315
316 if (P->Dimension == 1)
317 POL_ENSURE_VERTICES(P);
318
319 for (i = 0; i < P->NbRays; ++i)
320 if (value_one_p(P->Ray[i][1+P->Dimension])) {
321 sample = Vector_Alloc(P->Dimension+1);
322 Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
323 return sample;
324 }
325
326 if (P->NbEq > 0) {
327 Vector *Q_sample;
328 Polyhedron *Q = remove_all_equalities(P, &T, options->MaxRays);
329 if (!Q)
330 return NULL;
331 Q_sample = Polyhedron_Sample(Q, options);
332 Polyhedron_Free(Q);
333 if (Q_sample) {
334 sample = Vector_Alloc(P->Dimension + 1);
335 Matrix_Vector_Product(T, Q_sample->p, sample->p);
336 Vector_Free(Q_sample);
337 }
338 Matrix_Free(T);
339 return sample;
340 }
341
342 Q = Polyhedron_Factor(P, 0, &T, options->MaxRays);
343 if (Q) {
344 sample = product_sample(Q, T, options);
345 Domain_Free(Q);
346 Matrix_Free(T);
347 return sample;
348 }
349
350 value_init(min);
351 value_init(max);
352
353 obj = Vector_Alloc(P->Dimension+1);
354 value_set_si(obj->p[0], 1);
355 res = polyhedron_range(P, obj->p, obj->p[0], &min, &max, options);
356 if (res == lp_unbounded) {
357 unbounded:
358 value_clear(min);
359 value_clear(max);
360 Vector_Free(obj);
361
362 Q = remove_ray(P, options->MaxRays);
363 assert(Q);
364 sample = Polyhedron_Sample(Q, options);
365 Polyhedron_Free(Q);
366 return sample;
367 }
368 if (res == lp_empty) {
369 value_clear(min);
370 value_clear(max);
371 Vector_Free(obj);
372 return NULL;
373 }
374 assert(res == lp_ok);
375
376 if (value_eq(min, max)) {
377 Q = P;
378 } else {
379 Matrix *M, *T;
380 Matrix *basis = Polyhedron_Reduced_Basis(P, options);
381
382 if (!basis)
383 goto unbounded;
384 T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
385 inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
386 for (i = 0; i < P->Dimension; ++i)
387 for (j = 0; j < P->Dimension; ++j)
388 value_assign(T->p[i][j], basis->p[i][j]);
389 value_set_si(T->p[P->Dimension][P->Dimension], 1);
390 Matrix_Free(basis);
391
392 M = Matrix_Copy(T);
393 ok = Matrix_Inverse(M, inv);
394 assert(ok);
395 Matrix_Free(M);
396
397 Q = Polyhedron_Image(P, T, options->MaxRays);
398 res = polyhedron_range(Q, obj->p, obj->p[0], &min, &max, options);
399 assert(res == lp_ok);
400
401 Matrix_Free(T);
402 }
403
404 value_init(tmp);
405
406 v = Vector_Alloc(1+Q->Dimension+1);
407 value_set_si(v->p[1], -1);
408
409 for (value_assign(tmp, min); value_le(tmp, max); value_increment(tmp, tmp)) {
410 Polyhedron *R, *S;
411 Matrix *T;
412 Vector *S_sample;
413 value_assign(v->p[1+Q->Dimension], tmp);
414
415 R = AddConstraints(v->p, 1, Q, options->MaxRays);
416 R = DomainConstraintSimplify(R, options->MaxRays);
417 if (emptyQ(R)) {
418 Polyhedron_Free(R);
419 continue;
420 }
421
422 S = Polyhedron_RemoveFixedColumns(R, &T);
423 Polyhedron_Free(R);
424 S_sample = Polyhedron_Sample(S, options);
425 Polyhedron_Free(S);
426 if (S_sample) {
427 Vector *Q_sample = obj;
428 obj = NULL;
429 Matrix_Vector_Product(T, S_sample->p, Q_sample->p);
430 Matrix_Free(T);
431 Vector_Free(S_sample);
432 if (!inv)
433 sample = Q_sample;
434 else {
435 sample = Vector_Alloc(P->Dimension + 1);
436 Matrix_Vector_Product(inv, Q_sample->p, sample->p);
437 Vector_Free(Q_sample);
438 }
439 break;
440 }
441 Matrix_Free(T);
442 }
443
444 if (obj)
445 Vector_Free(obj);
446 if (inv)
447 Matrix_Free(inv);
448 if (P != Q)
449 Polyhedron_Free(Q);
450 Vector_Free(v);
451
452 value_clear(min);
453 value_clear(max);
454 value_clear(tmp);
455
456 return sample;
457 }
lattice point enumeration library