update isl to version 0.11
[barvinok.git] / hull.c
blobdece5f0ba6ed4f001bc0d513811b7dff1273faae
1 #include <assert.h>
2 #include <barvinok/options.h>
3 #include <barvinok/sample.h>
4 #include <barvinok/util.h>
5 #include "hilbert.h"
6 #include "hull.h"
7 #include "ilp.h"
8 #include "polysign.h"
9 #include "remove_equalities.h"
10 #include "config.h"
12 #ifndef USE_ZSOLVE
13 Matrix *Cone_Hilbert_Integer_Hull(Polyhedron *C,
14 struct barvinok_options *options)
16 assert(0);
18 #endif
20 struct integer_hull {
21 Polyhedron *P; /* Original polyhedron or cone */
22 Polyhedron *init; /* Initial (under)approximation of integer hull */
24 Matrix *F; /* Set of already computed facets */
25 int n_F; /* The number of computed facets */
27 /* Computes a lower bound for the objective function obj over
28 * the integer hull, possibly exploiting the already computed
29 * facets of the integer hull given in hull->F.
31 void (*set_lower_bound)(struct integer_hull *hull, Value *obj,
32 Value *lower, struct barvinok_options *options);
35 static int matrix_add(Matrix *M, int n, Value *v)
37 if (n >= M->NbRows)
38 Matrix_Extend(M, 3*(M->NbRows+10)/2);
39 Vector_Copy(v, M->p[n], M->NbColumns);
40 return n+1;
43 /* Select "best" constraint from P over which to optimize,
44 * where P is the current approximation of the integer hull
45 * and where a constraint is better if the difference between
46 * the value at the constraint and the lower bound (based on
47 * the original polytope (the linear relaxation) is smaller.
49 static int select_best(struct integer_hull *hull,
50 Polyhedron *P,
51 Value *min, Value *max,
52 struct barvinok_options *options)
54 int i, j, k;
55 Vector *lower = Vector_Alloc(P->NbConstraints);
56 Vector *upper = Vector_Alloc(P->NbConstraints);
57 Vector *distances = Vector_Alloc(P->NbConstraints);
58 int i_min;
59 int non_zero = 0;
61 i_min = -1;
62 for (i = 0; i < P->NbConstraints; ++i) {
63 if (First_Non_Zero(P->Constraint[i]+1, P->Dimension) == -1)
64 continue;
65 for (j = 0; j < hull->n_F; ++j)
66 if (Vector_Equal(hull->F->p[j]+1, P->Constraint[i]+1, P->Dimension))
67 break;
68 if (j < hull->n_F)
69 continue;
70 hull->set_lower_bound(hull, P->Constraint[i]+1, &lower->p[i], options);
71 value_oppose(upper->p[i], P->Constraint[i][1+P->Dimension]);
72 value_subtract(distances->p[i], upper->p[i], lower->p[i]);
73 if (value_zero_p(distances->p[i]))
74 continue;
75 if (i_min == -1 || value_lt(distances->p[i], distances->p[i_min]))
76 i_min = i;
77 non_zero++;
80 if (i_min != -1) {
81 value_decrement(*max, upper->p[i_min]);
82 value_assign(*min, lower->p[i_min]);
85 Vector_Free(lower);
86 Vector_Free(upper);
87 Vector_Free(distances);
89 return i_min;
92 /* Return the (integer) vertices of P */
93 static Matrix *Polyhedron_Vertices(Polyhedron *P)
95 int i, j;
96 Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+1);
98 for (i = 0, j = 0; i < P->NbRays; ++i) {
99 if (value_zero_p(P->Ray[i][1+P->Dimension]))
100 continue;
101 assert(value_one_p(P->Ray[i][1+P->Dimension]));
102 value_set_si(M->p[j][P->Dimension], 1);
103 Vector_Copy(P->Ray[i]+1, M->p[j++], P->Dimension);
105 M->NbRows = j;
107 return M;
110 /* Set INCLUDE_KNOWN_FACETS_IN_LP to 1 to compute a lower bound for
111 * the minimal value wrt a facet based on the facets that have
112 * already been determined to be facets of the integer hull,
113 * rather than simply setting the lower bound to 0 or 1.
114 * This does not appear to be profitable, probably because
115 * the lower bound is recomputed every time we optimize
116 * wrt a facet. If we were to cache the lower bounds for
117 * facets that do not change when a new vertex is added,
118 * this may still turn out to be useful.
120 #define INCLUDE_KNOWN_FACETS_IN_LP 0
121 /* Set INCLUDE_KNOWN_FACETS_IN_ILP to 1 to take into account the
122 * facets that have already been determined to be facets
123 * of the integer hull while optimizing wrt some other
124 * facet. The idea is that the polytope over which to
125 * optimize could be flatter. However, in practice, adding
126 * these constraints seems to slow down the computation
127 * (on at least one example).
129 #define INCLUDE_KNOWN_FACETS_IN_ILP 0
131 Polyhedron *add_facets(Polyhedron *P, Matrix *facets, int n,
132 struct barvinok_options *options)
134 int i;
135 Polyhedron *R;
136 Matrix *M;
138 M = Matrix_Alloc(P->NbConstraints + n, facets->NbColumns);
139 for (i = 0; i < n; ++i)
140 Vector_Copy(facets->p[i], M->p[i], facets->NbColumns);
141 for (i = 0; i < P->NbConstraints; ++i)
142 Vector_Copy(P->Constraint[i], M->p[n + i], M->NbColumns);
143 R = Constraints2Polyhedron(M, options->MaxRays);
144 Matrix_Free(M);
145 return R;
148 /* Add optimal and suboptimal points found by Polyhedron_Integer_Minimum
149 * to Q, replacing (and destroying) the old Q.
150 * rays is a temporary array for use in AddRays that can be reused
151 * over multiple calls and that is extended as needed.
153 Polyhedron *add_points(Polyhedron *Q, Vector *opt, Matrix *subopt,
154 int n_subopt, Matrix *rays,
155 struct barvinok_options *options)
157 int i;
158 Polyhedron *R;
160 if (1 + n_subopt > rays->NbRows) {
161 int n = rays->NbRows;
162 Matrix_Extend(rays, 1+n_subopt);
163 for (i = n; i < rays->NbRows; ++i) {
164 value_set_si(rays->p[i][0], 1);
165 value_set_si(rays->p[i][1+Q->Dimension], 1);
169 Vector_Copy(opt->p, rays->p[0]+1, Q->Dimension);
170 for (i = 0; i < n_subopt; ++i)
171 Vector_Copy(subopt->p[i], rays->p[1+i]+1, Q->Dimension);
173 R = AddRays(rays->p[0], 1+n_subopt, Q, options->MaxRays);
174 Polyhedron_Free(Q);
176 return R;
179 /* Extends an initial approximation hull->init of the integer hull
180 * of hull->P using generalized basis reduction.
181 * In particular, it considers each facet of the current
182 * approximation add optimizes in the direction of the normal,
183 * adding the optimal point to the approximation as well as
184 * all other points that may have been found along the way,
185 * until no facet yields any more new points.
186 * Returns a list of the vertices of this integer hull.
188 static Matrix *gbr_hull_extend(struct integer_hull *hull,
189 struct barvinok_options *options)
191 Value min, max;
192 Polyhedron *Q = hull->init;
193 Matrix *rays = Matrix_Alloc(1, 2+Q->Dimension);
194 Matrix *candidates;
195 int n_candidates = 0;
196 Matrix *vertices;
197 int i_min;
199 hull->F= Matrix_Alloc(Q->NbConstraints, 2+Q->Dimension);
200 hull->n_F = 0;
202 candidates = Matrix_Alloc(0, Q->Dimension);
204 value_init(min);
205 value_init(max);
207 value_set_si(rays->p[0][0], 1);
208 value_set_si(rays->p[0][1+Q->Dimension], 1);
210 while ((i_min = select_best(hull, Q, &min, &max, options)) != -1) {
211 Polyhedron *R;
212 int i;
213 Vector *opt;
215 if (INCLUDE_KNOWN_FACETS_IN_ILP)
216 R = add_facets(hull->P, hull->F, hull->n_F, options);
217 else
218 R = hull->P;
220 n_candidates = 0;
221 opt = Polyhedron_Integer_Minimum(R, Q->Constraint[i_min]+1,
222 min, max, candidates, &n_candidates,
223 options);
224 candidates->NbRows = n_candidates;
225 if (INCLUDE_KNOWN_FACETS_IN_ILP)
226 Polyhedron_Free(R);
228 hull->n_F = matrix_add(hull->F, hull->n_F, Q->Constraint[i_min]);
230 if (!opt)
231 continue;
233 Inner_Product(hull->F->p[hull->n_F-1]+1, opt->p, Q->Dimension,
234 &hull->F->p[hull->n_F-1][1+Q->Dimension]);
235 value_oppose(hull->F->p[hull->n_F-1][1+Q->Dimension],
236 hull->F->p[hull->n_F-1][1+Q->Dimension]);
238 Q = add_points(Q, opt, candidates, n_candidates, rays, options);
240 Vector_Free(opt);
243 vertices = Polyhedron_Vertices(Q);
245 Polyhedron_Free(Q);
246 hull->init = NULL;
248 value_clear(min);
249 value_clear(max);
251 Matrix_Free(rays);
252 Matrix_Free(hull->F);
253 hull->F = NULL;
254 Matrix_Free(candidates);
256 return vertices;
259 /* Returns the Minkowski sum of the cone and the polytope
260 * that is the convex hull of its (integer) extremal rays.
262 static Polyhedron *truncate_cone(Polyhedron *C,
263 struct barvinok_options *options)
265 int i;
266 Matrix *V;
267 Polyhedron *P;
269 V = Matrix_Alloc(2*C->NbRays, 2+C->Dimension);
270 for (i = 0; i < C->NbRays; ++i) {
271 if (value_notzero_p(C->Ray[i][1+C->Dimension]))
272 continue;
273 Vector_Copy(C->Ray[i], V->p[C->NbRays+i], 1+C->Dimension+1);
274 Vector_Copy(C->Ray[i], V->p[i], 1+C->Dimension);
275 value_set_si(V->p[i][1+C->Dimension], 1);
277 P = Rays2Polyhedron(V, options->MaxRays);
278 Matrix_Free(V);
279 return P;
282 /* Frees original list and transformation matrix CV */
283 static Matrix *transform_list_of_vertices(Matrix *list, Matrix *CV)
285 Matrix *T, *M;
286 T = Transpose(CV);
287 M = Matrix_Alloc(list->NbRows, T->NbColumns);
288 Matrix_Product(list, T, M);
289 Matrix_Free(list);
290 Matrix_Free(CV);
291 Matrix_Free(T);
293 return M;
296 /* Set the lower bound for optimizing along the normal of a facet
297 * in case of computing the integer hull of a truncated cone
298 * (i.e., a cone with the origin removed).
299 * In particular, if the constraint is one bounding the original
300 * cone, then the lower bound is zero; if it is a new constraint,
301 * then the lower bound is one.
302 * A more accurate bound can be obtained by looking at the
303 * already computed facets of the integer hull.
305 static void set_to_one(struct integer_hull *hull, Value *obj,
306 Value *lower, struct barvinok_options *options)
308 if (INCLUDE_KNOWN_FACETS_IN_LP) {
309 int i;
310 enum lp_result res;
311 Value one;
313 if (hull->F->NbRows < hull->n_F+hull->P->NbConstraints)
314 Matrix_Extend(hull->F, hull->n_F+hull->P->NbConstraints);
315 for (i = 0; i < hull->P->NbConstraints; ++i)
316 Vector_Copy(hull->P->Constraint[i], hull->F->p[hull->n_F+i],
317 hull->F->NbColumns);
319 value_init(one);
320 value_set_si(one, 1);
322 res = constraints_opt(hull->F, obj, one, lp_min, lower, options);
323 value_subtract(*lower, *lower, obj[hull->P->Dimension]);
324 assert(res == lp_ok);
326 if (value_notzero_p(*lower))
327 return;
330 if (value_zero_p(obj[hull->P->Dimension]))
331 value_set_si(*lower, 0);
332 else
333 value_set_si(*lower, 1);
336 /* Add to Q those points that the calling function already knows about,
337 * replacing (and destroying) the original Q.
339 static Polyhedron *add_known_points(Polyhedron *Q,
340 Matrix *c, int n_c, unsigned MaxRays)
342 Polyhedron *R;
344 R = AddRays(c->p[0], n_c, Q, MaxRays);
345 Polyhedron_Free(Q);
346 return R;
349 /* Computes the integer hull of the (truncated) cone C
350 * using generalized basis reduction, where c contains a list
351 * of vertices that are known to be part of this integer hull.
353 static Matrix *gbr_cone_hull(Polyhedron *C, Matrix *c, int n_c,
354 struct barvinok_options *options)
356 Matrix *vertices;
357 Matrix *CV = NULL;
358 struct integer_hull hull;
360 POL_ENSURE_VERTICES(C);
361 remove_all_equalities(&C, NULL, NULL, &CV, 0, options->MaxRays);
363 POL_ENSURE_VERTICES(C);
365 hull.P = C;
366 hull.init = truncate_cone(C, options);
367 hull.set_lower_bound = set_to_one;
368 if (!CV && c && n_c)
369 hull.init = add_known_points(hull.init, c, n_c, options->MaxRays);
370 vertices = gbr_hull_extend(&hull, options);
372 if (CV) {
373 vertices = transform_list_of_vertices(vertices, CV);
374 Polyhedron_Free(C);
377 return vertices;
380 /* Computes the integer hull of the cone C with the origin removed. */
381 Matrix *Cone_Integer_Hull(Polyhedron *C, Matrix *candidates,
382 int n_candidates, struct barvinok_options *options)
384 switch(options->integer_hull) {
385 case BV_HULL_GBR:
386 return gbr_cone_hull(C, candidates, n_candidates, options);
387 case BV_HULL_HILBERT:
388 return Cone_Hilbert_Integer_Hull(C, options);
389 default:
390 assert(0);
394 /* Computes initial full-dimensional approximation of the integer hull of P,
395 * or discovers an implicit equality of P.
396 * We start off with the integer vertices of P itself, if any.
397 * If there are no such vertices (available), then we simply
398 * take an arbitrary integer point in the polytope.
399 * Then we keep optimizing over normals of equalities in the description
400 * of the current approximation, until we find an equality that holds
401 * for the entire integer hull of P or until we have obtained a
402 * full-dimensional approximation.
404 static Polyhedron *internal_polytope(Polyhedron *P,
405 struct barvinok_options *options)
407 int i, j;
408 Polyhedron *init;
409 Matrix *vertices = Matrix_Alloc(P->NbRays, P->Dimension+2);
411 for (i = 0, j = 0; i < P->NbRays; ++i) {
412 if (value_notone_p(P->Ray[i][1+P->Dimension]))
413 continue;
414 Vector_Copy(P->Ray[i], vertices->p[j++], 2+P->Dimension);
416 vertices->NbRows = j;
417 if (!j) {
418 Vector *sample = Polyhedron_Sample(P, options);
419 if (sample) {
420 value_set_si(vertices->p[0][1], 1);
421 Vector_Copy(sample->p, vertices->p[0]+1, P->Dimension+1);
422 Vector_Free(sample);
423 vertices->NbRows = 1;
426 init = Rays2Polyhedron(vertices, options->MaxRays);
427 Matrix_Free(vertices);
429 while (!emptyQ(init) && init->NbEq) {
430 Vector *sample;
431 Vector *v = Vector_Alloc(init->Dimension+2);
432 Polyhedron *Q;
434 value_set_si(v->p[0], 1);
435 Vector_Copy(init->Constraint[0]+1, v->p+1, init->Dimension+1);
436 value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
438 Q = AddConstraints(v->p, 1, P, options->MaxRays);
439 sample = Polyhedron_Sample(Q, options);
440 Polyhedron_Free(Q);
441 if (!sample) {
442 Vector_Oppose(init->Constraint[0]+1, v->p+1, init->Dimension+1);
443 value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
445 Q = AddConstraints(v->p, 1, P, options->MaxRays);
446 sample = Polyhedron_Sample(Q, options);
447 Polyhedron_Free(Q);
449 if (!sample) {
450 P = AddConstraints(init->Constraint[0], 1, P, options->MaxRays);
451 Polyhedron_Free(init);
452 Vector_Free(v);
453 return P;
455 assert(sample);
457 Vector_Copy(sample->p, v->p+1, init->Dimension+1);
458 Q = AddRays(v->p, 1, init, options->MaxRays);
459 Polyhedron_Free(init);
460 init = Q;
462 Vector_Free(sample);
463 Vector_Free(v);
466 return init;
469 /* Set the lower bound for optimizing along the normal of a facet
470 * in case of computing the integer hull of a polyhedron.
471 * Although obj contains the full affine objective function,
472 * the calling function expect a lower bound on the linear part only.
473 * We therefore need to subtract the constant from the computed
474 * optimum of the linear relaxation.
476 static void set_lower(struct integer_hull *hull, Value *obj,
477 Value *lower, struct barvinok_options *options)
479 enum lp_result res;
480 Value one;
482 value_init(one);
483 value_set_si(one, 1);
485 res = polyhedron_opt(hull->P, obj, one, lp_min, lower, options);
486 value_subtract(*lower, *lower, obj[hull->P->Dimension]);
487 assert(res == lp_ok);
489 value_clear(one);
492 /* Computes the vertices of the integer hull of P
494 Matrix *Polyhedron_Integer_Hull(Polyhedron *P,
495 struct barvinok_options *options)
497 Matrix *vertices;
498 Matrix *CV = NULL;
499 Polyhedron *init;
501 POL_ENSURE_VERTICES(P);
502 remove_all_equalities(&P, NULL, NULL, &CV, 0, options->MaxRays);
504 POL_ENSURE_VERTICES(P);
506 if (P->Dimension == 0) {
507 vertices = Matrix_Alloc(1, 1);
508 value_set_si(vertices->p[0][0], 1);
509 } else {
510 init = internal_polytope(P, options);
511 if (emptyQ(init)) {
512 vertices = Matrix_Alloc(0, P->Dimension+1);
513 Polyhedron_Free(init);
514 } else if (init->NbEq) {
515 vertices = Polyhedron_Integer_Hull(init, options);
516 Polyhedron_Free(init);
517 } else {
518 struct integer_hull hull;
519 hull.P = P;
520 hull.init = init;
521 hull.set_lower_bound = set_lower;
522 vertices = gbr_hull_extend(&hull, options);
526 if (CV) {
527 vertices = transform_list_of_vertices(vertices, CV);
528 Polyhedron_Free(P);
531 return vertices;