drop verify_lexsmaller test application
[barvinok.git] / hull.c
blob3896d96afcc2efef83ec9eb0059e946df4fa7e78
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;
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 Vector *opt;
214 if (INCLUDE_KNOWN_FACETS_IN_ILP)
215 R = add_facets(hull->P, hull->F, hull->n_F, options);
216 else
217 R = hull->P;
219 n_candidates = 0;
220 opt = Polyhedron_Integer_Minimum(R, Q->Constraint[i_min]+1,
221 min, max, candidates, &n_candidates,
222 options);
223 candidates->NbRows = n_candidates;
224 if (INCLUDE_KNOWN_FACETS_IN_ILP)
225 Polyhedron_Free(R);
227 hull->n_F = matrix_add(hull->F, hull->n_F, Q->Constraint[i_min]);
229 if (!opt)
230 continue;
232 Inner_Product(hull->F->p[hull->n_F-1]+1, opt->p, Q->Dimension,
233 &hull->F->p[hull->n_F-1][1+Q->Dimension]);
234 value_oppose(hull->F->p[hull->n_F-1][1+Q->Dimension],
235 hull->F->p[hull->n_F-1][1+Q->Dimension]);
237 Q = add_points(Q, opt, candidates, n_candidates, rays, options);
239 Vector_Free(opt);
242 vertices = Polyhedron_Vertices(Q);
244 Polyhedron_Free(Q);
245 hull->init = NULL;
247 value_clear(min);
248 value_clear(max);
250 Matrix_Free(rays);
251 Matrix_Free(hull->F);
252 hull->F = NULL;
253 Matrix_Free(candidates);
255 return vertices;
258 /* Returns the Minkowski sum of the cone and the polytope
259 * that is the convex hull of its (integer) extremal rays.
261 static Polyhedron *truncate_cone(Polyhedron *C,
262 struct barvinok_options *options)
264 int i;
265 Matrix *V;
266 Polyhedron *P;
268 V = Matrix_Alloc(2*C->NbRays, 2+C->Dimension);
269 for (i = 0; i < C->NbRays; ++i) {
270 if (value_notzero_p(C->Ray[i][1+C->Dimension]))
271 continue;
272 Vector_Copy(C->Ray[i], V->p[C->NbRays+i], 1+C->Dimension+1);
273 Vector_Copy(C->Ray[i], V->p[i], 1+C->Dimension);
274 value_set_si(V->p[i][1+C->Dimension], 1);
276 P = Rays2Polyhedron(V, options->MaxRays);
277 Matrix_Free(V);
278 return P;
281 /* Frees original list and transformation matrix CV */
282 static Matrix *transform_list_of_vertices(Matrix *list, Matrix *CV)
284 Matrix *T, *M;
285 T = Transpose(CV);
286 M = Matrix_Alloc(list->NbRows, T->NbColumns);
287 Matrix_Product(list, T, M);
288 Matrix_Free(list);
289 Matrix_Free(CV);
290 Matrix_Free(T);
292 return M;
295 /* Set the lower bound for optimizing along the normal of a facet
296 * in case of computing the integer hull of a truncated cone
297 * (i.e., a cone with the origin removed).
298 * In particular, if the constraint is one bounding the original
299 * cone, then the lower bound is zero; if it is a new constraint,
300 * then the lower bound is one.
301 * A more accurate bound can be obtained by looking at the
302 * already computed facets of the integer hull.
304 static void set_to_one(struct integer_hull *hull, Value *obj,
305 Value *lower, struct barvinok_options *options)
307 if (INCLUDE_KNOWN_FACETS_IN_LP) {
308 int i;
309 enum lp_result res;
310 Value one;
312 if (hull->F->NbRows < hull->n_F+hull->P->NbConstraints)
313 Matrix_Extend(hull->F, hull->n_F+hull->P->NbConstraints);
314 for (i = 0; i < hull->P->NbConstraints; ++i)
315 Vector_Copy(hull->P->Constraint[i], hull->F->p[hull->n_F+i],
316 hull->F->NbColumns);
318 value_init(one);
319 value_set_si(one, 1);
321 res = constraints_opt(hull->F, obj, one, lp_min, lower, options);
322 value_subtract(*lower, *lower, obj[hull->P->Dimension]);
323 assert(res == lp_ok);
325 if (value_notzero_p(*lower))
326 return;
329 if (value_zero_p(obj[hull->P->Dimension]))
330 value_set_si(*lower, 0);
331 else
332 value_set_si(*lower, 1);
335 /* Add to Q those points that the calling function already knows about,
336 * replacing (and destroying) the original Q.
338 static Polyhedron *add_known_points(Polyhedron *Q,
339 Matrix *c, int n_c, unsigned MaxRays)
341 Polyhedron *R;
343 R = AddRays(c->p[0], n_c, Q, MaxRays);
344 Polyhedron_Free(Q);
345 return R;
348 /* Computes the integer hull of the (truncated) cone C
349 * using generalized basis reduction, where c contains a list
350 * of vertices that are known to be part of this integer hull.
352 static Matrix *gbr_cone_hull(Polyhedron *C, Matrix *c, int n_c,
353 struct barvinok_options *options)
355 Matrix *vertices;
356 Matrix *CV = NULL;
357 struct integer_hull hull;
359 POL_ENSURE_VERTICES(C);
360 remove_all_equalities(&C, NULL, NULL, &CV, 0, options->MaxRays);
362 POL_ENSURE_VERTICES(C);
364 hull.P = C;
365 hull.init = truncate_cone(C, options);
366 hull.set_lower_bound = set_to_one;
367 if (!CV && c && n_c)
368 hull.init = add_known_points(hull.init, c, n_c, options->MaxRays);
369 vertices = gbr_hull_extend(&hull, options);
371 if (CV) {
372 vertices = transform_list_of_vertices(vertices, CV);
373 Polyhedron_Free(C);
376 return vertices;
379 /* Computes the integer hull of the cone C with the origin removed. */
380 Matrix *Cone_Integer_Hull(Polyhedron *C, Matrix *candidates,
381 int n_candidates, struct barvinok_options *options)
383 switch(options->integer_hull) {
384 case BV_HULL_GBR:
385 return gbr_cone_hull(C, candidates, n_candidates, options);
386 case BV_HULL_HILBERT:
387 return Cone_Hilbert_Integer_Hull(C, options);
388 default:
389 assert(0);
393 /* Computes initial full-dimensional approximation of the integer hull of P,
394 * or discovers an implicit equality of P.
395 * We start off with the integer vertices of P itself, if any.
396 * If there are no such vertices (available), then we simply
397 * take an arbitrary integer point in the polytope.
398 * Then we keep optimizing over normals of equalities in the description
399 * of the current approximation, until we find an equality that holds
400 * for the entire integer hull of P or until we have obtained a
401 * full-dimensional approximation.
403 static Polyhedron *internal_polytope(Polyhedron *P,
404 struct barvinok_options *options)
406 int i, j;
407 Polyhedron *init;
408 Matrix *vertices = Matrix_Alloc(P->NbRays, P->Dimension+2);
410 for (i = 0, j = 0; i < P->NbRays; ++i) {
411 if (value_notone_p(P->Ray[i][1+P->Dimension]))
412 continue;
413 Vector_Copy(P->Ray[i], vertices->p[j++], 2+P->Dimension);
415 vertices->NbRows = j;
416 if (!j) {
417 Vector *sample = Polyhedron_Sample(P, options);
418 if (sample) {
419 value_set_si(vertices->p[0][1], 1);
420 Vector_Copy(sample->p, vertices->p[0]+1, P->Dimension+1);
421 Vector_Free(sample);
422 vertices->NbRows = 1;
425 init = Rays2Polyhedron(vertices, options->MaxRays);
426 Matrix_Free(vertices);
428 while (!emptyQ(init) && init->NbEq) {
429 Vector *sample;
430 Vector *v = Vector_Alloc(init->Dimension+2);
431 Polyhedron *Q;
433 value_set_si(v->p[0], 1);
434 Vector_Copy(init->Constraint[0]+1, v->p+1, init->Dimension+1);
435 value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
437 Q = AddConstraints(v->p, 1, P, options->MaxRays);
438 sample = Polyhedron_Sample(Q, options);
439 Polyhedron_Free(Q);
440 if (!sample) {
441 Vector_Oppose(init->Constraint[0]+1, v->p+1, init->Dimension+1);
442 value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
444 Q = AddConstraints(v->p, 1, P, options->MaxRays);
445 sample = Polyhedron_Sample(Q, options);
446 Polyhedron_Free(Q);
448 if (!sample) {
449 P = AddConstraints(init->Constraint[0], 1, P, options->MaxRays);
450 Polyhedron_Free(init);
451 Vector_Free(v);
452 return P;
454 assert(sample);
456 Vector_Copy(sample->p, v->p+1, init->Dimension+1);
457 Q = AddRays(v->p, 1, init, options->MaxRays);
458 Polyhedron_Free(init);
459 init = Q;
461 Vector_Free(sample);
462 Vector_Free(v);
465 return init;
468 /* Set the lower bound for optimizing along the normal of a facet
469 * in case of computing the integer hull of a polyhedron.
470 * Although obj contains the full affine objective function,
471 * the calling function expect a lower bound on the linear part only.
472 * We therefore need to subtract the constant from the computed
473 * optimum of the linear relaxation.
475 static void set_lower(struct integer_hull *hull, Value *obj,
476 Value *lower, struct barvinok_options *options)
478 enum lp_result res;
479 Value one;
481 value_init(one);
482 value_set_si(one, 1);
484 res = polyhedron_opt(hull->P, obj, one, lp_min, lower, options);
485 value_subtract(*lower, *lower, obj[hull->P->Dimension]);
486 assert(res == lp_ok);
488 value_clear(one);
491 /* Computes the vertices of the integer hull of P
493 Matrix *Polyhedron_Integer_Hull(Polyhedron *P,
494 struct barvinok_options *options)
496 Matrix *vertices;
497 Matrix *CV = NULL;
498 Polyhedron *init;
500 POL_ENSURE_VERTICES(P);
501 remove_all_equalities(&P, NULL, NULL, &CV, 0, options->MaxRays);
503 POL_ENSURE_VERTICES(P);
505 if (P->Dimension == 0) {
506 vertices = Matrix_Alloc(1, 1);
507 value_set_si(vertices->p[0][0], 1);
508 } else {
509 init = internal_polytope(P, options);
510 if (emptyQ(init)) {
511 vertices = Matrix_Alloc(0, P->Dimension+1);
512 Polyhedron_Free(init);
513 } else if (init->NbEq) {
514 vertices = Polyhedron_Integer_Hull(init, options);
515 Polyhedron_Free(init);
516 } else {
517 struct integer_hull hull;
518 hull.P = P;
519 hull.init = init;
520 hull.set_lower_bound = set_lower;
521 vertices = gbr_hull_extend(&hull, options);
525 if (CV) {
526 vertices = transform_list_of_vertices(vertices, CV);
527 Polyhedron_Free(P);
530 return vertices;