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