hull.c

   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
  10 struct integer_hull {
  11     Polyhedron  *P;     /* Original polyhedron or cone */
  12     Polyhedron  *init;  /* Initial (under)approximation of integer hull */
  13
  14     Matrix      *F;     /* Set of already computed facets */
  15     int         n_F;    /* The number of computed facets  */
  16
  17     /* Computes a lower bound for the objective function obj over
  18      * the integer hull, possibly exploiting the already computed
  19      * facets of the integer hull given in hull->F.
  20      */
  21     void (*set_lower_bound)(struct integer_hull *hull, Value *obj,
  22                             Value *lower, struct barvinok_options *options);
  23 };
  24
  25 static int matrix_add(Matrix *M, int n, Value *v)
  26 {
  27     if (n >= M->NbRows)
  28         Matrix_Extend(M, 3*(M->NbRows+10)/2);
  29     Vector_Copy(v, M->p[n], M->NbColumns);
  30     return n+1;
  31 }
  32
  33 static int select_best(struct integer_hull *hull,
  34                        Polyhedron *P, Matrix *candidates,
  35                        int *n_candidates, int *candidate,
  36                        Value *min, Value *max,
  37                        struct barvinok_options *options)
  38 {
  39     int i, j, k;
  40     Matrix *M = candidates;
  41     Vector *lower = Vector_Alloc(P->NbConstraints);
  42     Vector *upper = Vector_Alloc(P->NbConstraints+2);
  43     Vector *distances = Vector_Alloc(P->NbConstraints);
  44     int i_min;
  45     int non_zero = 0;
  46
  47     i_min = -1;
  48     for (i = 0; i < P->NbConstraints; ++i) {
  49         if (First_Non_Zero(P->Constraint[i]+1, P->Dimension) == -1)
  50             continue;
  51         for (j = 0; j < hull->n_F; ++j)
  52             if (Vector_Equal(hull->F->p[j]+1, P->Constraint[i]+1, P->Dimension))
  53                 break;
  54         if (j < hull->n_F)
  55             continue;
  56         hull->set_lower_bound(hull, P->Constraint[i]+1, &lower->p[i], options);
  57         value_oppose(upper->p[i], P->Constraint[i][1+P->Dimension]);
  58         value_subtract(distances->p[i], upper->p[i], lower->p[i]);
  59         if (value_zero_p(distances->p[i]))
  60             continue;
  61         if (i_min == -1 || value_lt(distances->p[i], distances->p[i_min]))
  62             i_min = i;
  63         non_zero++;
  64     }
  65     if (i_min == -1) {
  66         Vector_Free(lower);
  67         Vector_Free(upper);
  68         Vector_Free(distances);
  69         return -1;
  70     }
  71
  72     *candidate = -1;
  73     for (j = 0, k = 0; j < *n_candidates; ++j) {
  74         int keep = 0;
  75         for (i = 0; i < P->NbConstraints; ++i) {
  76             if (value_zero_p(distances->p[i]))
  77                 continue;
  78             Inner_Product(M->p[j], P->Constraint[i]+1, P->Dimension,
  79                           &upper->p[P->NbConstraints]);
  80             value_addto(upper->p[P->NbConstraints+1],
  81                         upper->p[P->NbConstraints],
  82                         P->Constraint[i][1+P->Dimension]);
  83             if (value_neg_p(upper->p[P->NbConstraints+1]))
  84                 keep = 1;
  85             if (value_lt(upper->p[P->NbConstraints], upper->p[i])) {
  86                 value_assign(upper->p[i], upper->p[P->NbConstraints]);
  87                 value_subtract(distances->p[i], upper->p[i], lower->p[i]);
  88                 if (i_min == i)
  89                     *candidate = k;
  90                 else if (value_lt(distances->p[i], distances->p[i_min])) {
  91                     i_min = i;
  92                     *candidate = k;
  93                 }
  94             }
  95         }
  96         if (keep) {
  97             if (k != j)
  98                 Vector_Exchange(M->p[j], M->p[k], M->NbColumns);
  99             ++k;
 100         }
 101     }
 102     *n_candidates = k;
 103
 104     value_decrement(*max, upper->p[i_min]);
 105     value_assign(*min, lower->p[i_min]);
 106
 107     Vector_Free(lower);
 108     Vector_Free(upper);
 109     Vector_Free(distances);
 110
 111     return i_min;
 112 }
 113
 114 /* Return the (integer) vertices of P */
 115 static Matrix *Polyhedron_Vertices(Polyhedron *P)
 116 {
 117     int i, j;
 118     Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+1);
 119
 120     for (i = 0, j = 0; i < P->NbRays; ++i) {
 121         if (value_zero_p(P->Ray[i][1+P->Dimension]))
 122             continue;
 123         assert(value_one_p(P->Ray[i][1+P->Dimension]));
 124         value_set_si(M->p[j][P->Dimension], 1);
 125         Vector_Copy(P->Ray[i]+1, M->p[j++], P->Dimension);
 126     }
 127     M->NbRows = j;
 128
 129     return M;
 130 }
 131
 132 /* Extends an initial approximation hull->init of the integer hull
 133  * of hull->P using generalized basis reduction.
 134  * In particular, it considers each facet of the current
 135  * approximation add optimizes in the direction of the normal,
 136  * adding the optimal point to the approximation as well as
 137  * all other points that may have been found along the way,
 138  * until no facet yields any more new points.
 139  * Returns a list of the vertices of this integer hull.
 140  */
 141 static Matrix *gbr_hull_extend(struct integer_hull *hull,
 142                                 struct barvinok_options *options)
 143 {
 144     Value min, max;
 145     Polyhedron *Q = hull->init;
 146     Vector *ray = Vector_Alloc(2+Q->Dimension);
 147     Matrix *candidates;
 148     int n_candidates = 0;
 149     Matrix *vertices;
 150
 151     hull->F= Matrix_Alloc(Q->NbConstraints, 2+Q->Dimension);
 152     hull->n_F = 0;
 153
 154     candidates = Matrix_Alloc(0, Q->Dimension);
 155
 156     value_init(min);
 157     value_init(max);
 158
 159     value_set_si(ray->p[0], 1);
 160     value_set_si(ray->p[1+Q->Dimension], 1);
 161
 162     for (;;) {
 163         Polyhedron *R;
 164         int i, i_min, candidate;
 165         Vector *opt;
 166         Value *vertex = NULL;
 167
 168         i_min = select_best(hull, Q, candidates, &n_candidates, &candidate,
 169                                 &min, &max, options);
 170         if (i_min == -1)
 171             break;
 172
 173         opt = Polyhedron_Integer_Minimum(hull->P, Q->Constraint[i_min]+1,
 174                                          min, max, candidates, &n_candidates,
 175                                          options);
 176         candidates->NbRows = n_candidates;
 177
 178         hull->n_F = matrix_add(hull->F, hull->n_F, Q->Constraint[i_min]);
 179
 180         if (opt)
 181             vertex = opt->p;
 182         else if (candidate != -1)
 183             vertex = candidates->p[candidate];
 184
 185         if (!vertex)
 186             continue;
 187
 188         Inner_Product(hull->F->p[hull->n_F-1]+1, vertex, Q->Dimension,
 189                       &hull->F->p[hull->n_F-1][1+Q->Dimension]);
 190         value_oppose(hull->F->p[hull->n_F-1][1+Q->Dimension],
 191                      hull->F->p[hull->n_F-1][1+Q->Dimension]);
 192
 193         Vector_Copy(vertex, ray->p+1, Q->Dimension);
 194         R = AddRays(ray->p, 1, Q, options->MaxRays);
 195         Polyhedron_Free(Q);
 196         Q = R;
 197
 198         if (opt)
 199             Vector_Free(opt);
 200     }
 201
 202     vertices = Polyhedron_Vertices(Q);
 203
 204     Polyhedron_Free(Q);
 205     hull->init = NULL;
 206
 207     value_clear(min);
 208     value_clear(max);
 209
 210     Vector_Free(ray);
 211     Matrix_Free(hull->F);
 212     hull->F = NULL;
 213     Matrix_Free(candidates);
 214
 215     return vertices;
 216 }
 217
 218 /* Returns the Minkowski sum of the cone and the polytope
 219  * that is the convex hull of its (integer) extremal rays.
 220  */
 221 static Polyhedron *truncate_cone(Polyhedron *C,
 222                                  struct barvinok_options *options)
 223 {
 224     int i;
 225     Matrix *V;
 226     Polyhedron *P;
 227
 228     V = Matrix_Alloc(2*C->NbRays, 2+C->Dimension);
 229     for (i = 0; i < C->NbRays; ++i) {
 230         if (value_notzero_p(C->Ray[i][1+C->Dimension]))
 231             continue;
 232         Vector_Copy(C->Ray[i], V->p[C->NbRays+i], 1+C->Dimension+1);
 233         Vector_Copy(C->Ray[i], V->p[i], 1+C->Dimension);
 234         value_set_si(V->p[i][1+C->Dimension], 1);
 235     }
 236     P = Rays2Polyhedron(V, options->MaxRays);
 237     Matrix_Free(V);
 238     return P;
 239 }
 240
 241 /* Frees original list and transformation matrix CV */
 242 static Matrix *transform_list_of_vertices(Matrix *list, Matrix *CV)
 243 {
 244     Matrix *T, *M;
 245     T = Transpose(CV);
 246     M = Matrix_Alloc(list->NbRows, T->NbColumns);
 247     Matrix_Product(list, T, M);
 248     Matrix_Free(list);
 249     Matrix_Free(CV);
 250     Matrix_Free(T);
 251
 252     return M;
 253 }
 254
 255 /* Set the lower bound for optimizing along the normal of a facet
 256  * in case of computing the integer hull of a truncated cone
 257  * (i.e., a cone with the origin removed).
 258  * In particular, if the constraint is one bounding the original
 259  * cone, then the lower bound is zero; if it is a new constraint,
 260  * then the lower bound is one.
 261  * A more accurate bound can be obtained by looking at the
 262  * already computed facets of the integer hull.
 263  */
 264 static void set_to_one(struct integer_hull *hull, Value *obj,
 265                        Value *lower, struct barvinok_options *options)
 266 {
 267     if (value_zero_p(obj[hull->P->Dimension]))
 268         value_set_si(*lower, 0);
 269     else
 270         value_set_si(*lower, 1);
 271 }
 272
 273 static Matrix *gbr_hull(Polyhedron *C, struct barvinok_options *options)
 274 {
 275     Matrix *vertices;
 276     Matrix *CV = NULL;
 277     struct integer_hull hull;
 278
 279     POL_ENSURE_VERTICES(C);
 280     remove_all_equalities(&C, NULL, NULL, &CV, 0, options->MaxRays);
 281
 282     POL_ENSURE_VERTICES(C);
 283
 284     hull.P = C;
 285     hull.init = truncate_cone(C, options);
 286     hull.set_lower_bound = set_to_one;
 287     vertices = gbr_hull_extend(&hull, options);
 288
 289     if (CV) {
 290         vertices = transform_list_of_vertices(vertices, CV);
 291         Polyhedron_Free(C);
 292     }
 293
 294     return vertices;
 295 }
 296
 297 Matrix *Cone_Integer_Hull(Polyhedron *C, struct barvinok_options *options)
 298 {
 299     switch(options->integer_hull) {
 300     case BV_HULL_GBR:
 301         return gbr_hull(C, options);
 302     case BV_HULL_HILBERT:
 303         return Cone_Hilbert_Integer_Hull(C, options);
 304     default:
 305         assert(0);
 306     }
 307 }
 308
 309 /* Computes initial full-dimensional approximation of the integer hull of P,
 310  * or discovers an implicit equality of P.
 311  * We start off with the integer vertices of P itself, if any.
 312  * If there are no such vertices (available), then we simply
 313  * take an arbitrary integer point in the polytope.
 314  * Then we keep optimizing over normals of equalities in the description
 315  * of the current approximation, until we find an equality that holds
 316  * for the entire integer hull of P or until we have obtained a
 317  * full-dimensional approximation.
 318  */
 319 static Polyhedron *internal_polytope(Polyhedron *P,
 320                                      struct barvinok_options *options)
 321 {
 322     int i, j;
 323     Polyhedron *init;
 324     Matrix *vertices = Matrix_Alloc(P->NbRays, P->Dimension+2);
 325
 326     for (i = 0, j = 0; i < P->NbRays; ++i) {
 327         if (value_notone_p(P->Ray[i][1+P->Dimension]))
 328             continue;
 329         Vector_Copy(P->Ray[i], vertices->p[j++], 2+P->Dimension);
 330     }
 331     vertices->NbRows = j;
 332     if (!j) {
 333         Vector *sample = Polyhedron_Sample(P, options);
 334         if (sample) {
 335             value_set_si(vertices->p[0][1], 1);
 336             Vector_Copy(sample->p, vertices->p[0]+1, P->Dimension+1);
 337             Vector_Free(sample);
 338             vertices->NbRows = 1;
 339         }
 340     }
 341     init = Rays2Polyhedron(vertices, options->MaxRays);
 342     Matrix_Free(vertices);
 343
 344     while (!emptyQ(init) && init->NbEq) {
 345         Vector *sample;
 346         Vector *v = Vector_Alloc(init->Dimension+2);
 347         Polyhedron *Q;
 348
 349         value_set_si(v->p[0], 1);
 350         Vector_Copy(init->Constraint[0]+1, v->p+1, init->Dimension+1);
 351         value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
 352
 353         Q = AddConstraints(v->p, 1, P, options->MaxRays);
 354         sample = Polyhedron_Sample(Q, options);
 355         Polyhedron_Free(Q);
 356         if (!sample) {
 357             Vector_Oppose(init->Constraint[0]+1, v->p+1, init->Dimension+1);
 358             value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
 359
 360             Q = AddConstraints(v->p, 1, P, options->MaxRays);
 361             sample = Polyhedron_Sample(Q, options);
 362             Polyhedron_Free(Q);
 363         }
 364         if (!sample) {
 365             P = AddConstraints(init->Constraint[0], 1, P, options->MaxRays);
 366             Polyhedron_Free(init);
 367             Vector_Free(v);
 368             return P;
 369         }
 370         assert(sample);
 371
 372         Vector_Copy(sample->p, v->p+1, init->Dimension+1);
 373         Q = AddRays(v->p, 1, init, options->MaxRays);
 374         Polyhedron_Free(init);
 375         init = Q;
 376
 377         Vector_Free(sample);
 378         Vector_Free(v);
 379     }
 380
 381     return init;
 382 }
 383
 384 /* Set the lower bound for optimizing along the normal of a facet
 385  * in case of computing the integer hull of a polyhedron.
 386  * Although obj contains the full affine objective function,
 387  * the calling function expect a lower bound on the linear part only.
 388  * We therefore need to subtract the constant from the computed
 389  * optimum of the linear relaxation.
 390  */
 391 static void set_lower(struct integer_hull *hull, Value *obj,
 392                       Value *lower, struct barvinok_options *options)
 393 {
 394     enum lp_result res;
 395     Value one;
 396
 397     value_init(one);
 398     value_set_si(one, 1);
 399
 400     res = polyhedron_opt(hull->P, obj, one, lp_min, lower, options);
 401     value_subtract(*lower, *lower, obj[hull->P->Dimension]);
 402     assert(res == lp_ok);
 403
 404     value_clear(one);
 405 }
 406
 407 /* Computes the vertices of the integer hull of P
 408  */
 409 Matrix *Polyhedron_Integer_Hull(Polyhedron *P,
 410                                 struct barvinok_options *options)
 411 {
 412     Matrix *vertices;
 413     Matrix *CV = NULL;
 414     Polyhedron *init;
 415
 416     POL_ENSURE_VERTICES(P);
 417     remove_all_equalities(&P, NULL, NULL, &CV, 0, options->MaxRays);
 418
 419     POL_ENSURE_VERTICES(P);
 420
 421     if (P->Dimension == 0) {
 422         vertices = Matrix_Alloc(1, 1);
 423         value_set_si(vertices->p[0][0], 1);
 424     } else {
 425         init = internal_polytope(P, options);
 426         if (emptyQ(init)) {
 427             vertices = Matrix_Alloc(0, P->Dimension+1);
 428             Polyhedron_Free(init);
 429         } else if (init->NbEq) {
 430             vertices = Polyhedron_Integer_Hull(init, options);
 431             Polyhedron_Free(init);
 432         } else {
 433             struct integer_hull hull;
 434             hull.P = P;
 435             hull.init = init;
 436             hull.set_lower_bound = set_lower;
 437             vertices = gbr_hull_extend(&hull, options);
 438         }
 439     }
 440
 441     if (CV) {
 442         vertices = transform_list_of_vertices(vertices, CV);
 443         Polyhedron_Free(P);
 444     }
 445
 446     return vertices;
 447 }