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 /* Set INCLUDE_KNOWN_FACETS_IN_LP to 1 to compute a lower bound for
 133  * the minimal value wrt a facet based on the facets that have
 134  * already been determined to be facets of the integer hull,
 135  * rather than simply setting the lower bound to 0 or 1.
 136  * This does not appear to be profitable, probably because
 137  * the lower bound is recomputed every time we optimize
 138  * wrt a facet.  If we were to cache the lower bounds for
 139  * facets that do not change when a new vertex is added,
 140  * this may still turn out to be useful.
 141  */
 142 #define INCLUDE_KNOWN_FACETS_IN_LP      0
 143 /* Set INCLUDE_KNOWN_FACETS_IN_ILP to 1 to take into account the
 144  * facets that have already been determined to be facets
 145  * of the integer hull while optimizing wrt some other
 146  * facet.  The idea is that the polytope over which to
 147  * optimize could be flatter.  However, in practice, adding
 148  * these constraints seems to slow down the computation
 149  * (on at least one example).
 150  */
 151 #define INCLUDE_KNOWN_FACETS_IN_ILP     0
 152
 153 Polyhedron *add_facets(Polyhedron *P, Matrix *facets, int n,
 154                         struct barvinok_options *options)
 155 {
 156     int i;
 157     Polyhedron *R;
 158     Matrix *M;
 159
 160     M = Matrix_Alloc(P->NbConstraints + n, facets->NbColumns);
 161     for (i = 0; i < n; ++i)
 162         Vector_Copy(facets->p[i], M->p[i], facets->NbColumns);
 163     for (i = 0; i < P->NbConstraints; ++i)
 164         Vector_Copy(P->Constraint[i], M->p[n + i], M->NbColumns);
 165     R = Constraints2Polyhedron(M, options->MaxRays);
 166     Matrix_Free(M);
 167     return R;
 168 }
 169
 170 /* Extends an initial approximation hull->init of the integer hull
 171  * of hull->P using generalized basis reduction.
 172  * In particular, it considers each facet of the current
 173  * approximation add optimizes in the direction of the normal,
 174  * adding the optimal point to the approximation as well as
 175  * all other points that may have been found along the way,
 176  * until no facet yields any more new points.
 177  * Returns a list of the vertices of this integer hull.
 178  */
 179 static Matrix *gbr_hull_extend(struct integer_hull *hull,
 180                                 struct barvinok_options *options)
 181 {
 182     Value min, max;
 183     Polyhedron *Q = hull->init;
 184     Vector *ray = Vector_Alloc(2+Q->Dimension);
 185     Matrix *candidates;
 186     int n_candidates = 0;
 187     Matrix *vertices;
 188
 189     hull->F= Matrix_Alloc(Q->NbConstraints, 2+Q->Dimension);
 190     hull->n_F = 0;
 191
 192     candidates = Matrix_Alloc(0, Q->Dimension);
 193
 194     value_init(min);
 195     value_init(max);
 196
 197     value_set_si(ray->p[0], 1);
 198     value_set_si(ray->p[1+Q->Dimension], 1);
 199
 200     for (;;) {
 201         Polyhedron *R;
 202         int i, i_min, candidate;
 203         Vector *opt;
 204         Value *vertex = NULL;
 205
 206         i_min = select_best(hull, Q, candidates, &n_candidates, &candidate,
 207                                 &min, &max, options);
 208         if (i_min == -1)
 209             break;
 210
 211         if (INCLUDE_KNOWN_FACETS_IN_ILP)
 212             R = add_facets(hull->P, hull->F, hull->n_F, options);
 213         else
 214             R = hull->P;
 215
 216         opt = Polyhedron_Integer_Minimum(R, Q->Constraint[i_min]+1,
 217                                          min, max, candidates, &n_candidates,
 218                                          options);
 219         candidates->NbRows = n_candidates;
 220         if (INCLUDE_KNOWN_FACETS_IN_ILP)
 221             Polyhedron_Free(R);
 222
 223         hull->n_F = matrix_add(hull->F, hull->n_F, Q->Constraint[i_min]);
 224
 225         if (opt)
 226             vertex = opt->p;
 227         else if (candidate != -1)
 228             vertex = candidates->p[candidate];
 229
 230         if (!vertex)
 231             continue;
 232
 233         Inner_Product(hull->F->p[hull->n_F-1]+1, vertex, 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]);
 237
 238         Vector_Copy(vertex, ray->p+1, Q->Dimension);
 239         R = AddRays(ray->p, 1, Q, options->MaxRays);
 240         Polyhedron_Free(Q);
 241         Q = R;
 242
 243         if (opt)
 244             Vector_Free(opt);
 245     }
 246
 247     vertices = Polyhedron_Vertices(Q);
 248
 249     Polyhedron_Free(Q);
 250     hull->init = NULL;
 251
 252     value_clear(min);
 253     value_clear(max);
 254
 255     Vector_Free(ray);
 256     Matrix_Free(hull->F);
 257     hull->F = NULL;
 258     Matrix_Free(candidates);
 259
 260     return vertices;
 261 }
 262
 263 /* Returns the Minkowski sum of the cone and the polytope
 264  * that is the convex hull of its (integer) extremal rays.
 265  */
 266 static Polyhedron *truncate_cone(Polyhedron *C,
 267                                  struct barvinok_options *options)
 268 {
 269     int i;
 270     Matrix *V;
 271     Polyhedron *P;
 272
 273     V = Matrix_Alloc(2*C->NbRays, 2+C->Dimension);
 274     for (i = 0; i < C->NbRays; ++i) {
 275         if (value_notzero_p(C->Ray[i][1+C->Dimension]))
 276             continue;
 277         Vector_Copy(C->Ray[i], V->p[C->NbRays+i], 1+C->Dimension+1);
 278         Vector_Copy(C->Ray[i], V->p[i], 1+C->Dimension);
 279         value_set_si(V->p[i][1+C->Dimension], 1);
 280     }
 281     P = Rays2Polyhedron(V, options->MaxRays);
 282     Matrix_Free(V);
 283     return P;
 284 }
 285
 286 /* Frees original list and transformation matrix CV */
 287 static Matrix *transform_list_of_vertices(Matrix *list, Matrix *CV)
 288 {
 289     Matrix *T, *M;
 290     T = Transpose(CV);
 291     M = Matrix_Alloc(list->NbRows, T->NbColumns);
 292     Matrix_Product(list, T, M);
 293     Matrix_Free(list);
 294     Matrix_Free(CV);
 295     Matrix_Free(T);
 296
 297     return M;
 298 }
 299
 300 /* Set the lower bound for optimizing along the normal of a facet
 301  * in case of computing the integer hull of a truncated cone
 302  * (i.e., a cone with the origin removed).
 303  * In particular, if the constraint is one bounding the original
 304  * cone, then the lower bound is zero; if it is a new constraint,
 305  * then the lower bound is one.
 306  * A more accurate bound can be obtained by looking at the
 307  * already computed facets of the integer hull.
 308  */
 309 static void set_to_one(struct integer_hull *hull, Value *obj,
 310                        Value *lower, struct barvinok_options *options)
 311 {
 312     if (INCLUDE_KNOWN_FACETS_IN_LP) {
 313         int i;
 314         enum lp_result res;
 315         Value one;
 316
 317         if (hull->F->NbRows < hull->n_F+hull->P->NbConstraints)
 318             Matrix_Extend(hull->F, hull->n_F+hull->P->NbConstraints);
 319         for (i = 0; i < hull->P->NbConstraints; ++i)
 320             Vector_Copy(hull->P->Constraint[i], hull->F->p[hull->n_F+i],
 321                         hull->F->NbColumns);
 322
 323         value_init(one);
 324         value_set_si(one, 1);
 325
 326         res = constraints_opt(hull->F, obj, one, lp_min, lower, options);
 327         value_subtract(*lower, *lower, obj[hull->P->Dimension]);
 328         assert(res == lp_ok);
 329
 330         if (value_notzero_p(*lower))
 331             return;
 332     }
 333
 334     if (value_zero_p(obj[hull->P->Dimension]))
 335         value_set_si(*lower, 0);
 336     else
 337         value_set_si(*lower, 1);
 338 }
 339
 340 static Matrix *gbr_hull(Polyhedron *C, struct barvinok_options *options)
 341 {
 342     Matrix *vertices;
 343     Matrix *CV = NULL;
 344     struct integer_hull hull;
 345
 346     POL_ENSURE_VERTICES(C);
 347     remove_all_equalities(&C, NULL, NULL, &CV, 0, options->MaxRays);
 348
 349     POL_ENSURE_VERTICES(C);
 350
 351     hull.P = C;
 352     hull.init = truncate_cone(C, options);
 353     hull.set_lower_bound = set_to_one;
 354     vertices = gbr_hull_extend(&hull, options);
 355
 356     if (CV) {
 357         vertices = transform_list_of_vertices(vertices, CV);
 358         Polyhedron_Free(C);
 359     }
 360
 361     return vertices;
 362 }
 363
 364 Matrix *Cone_Integer_Hull(Polyhedron *C, struct barvinok_options *options)
 365 {
 366     switch(options->integer_hull) {
 367     case BV_HULL_GBR:
 368         return gbr_hull(C, options);
 369     case BV_HULL_HILBERT:
 370         return Cone_Hilbert_Integer_Hull(C, options);
 371     default:
 372         assert(0);
 373     }
 374 }
 375
 376 /* Computes initial full-dimensional approximation of the integer hull of P,
 377  * or discovers an implicit equality of P.
 378  * We start off with the integer vertices of P itself, if any.
 379  * If there are no such vertices (available), then we simply
 380  * take an arbitrary integer point in the polytope.
 381  * Then we keep optimizing over normals of equalities in the description
 382  * of the current approximation, until we find an equality that holds
 383  * for the entire integer hull of P or until we have obtained a
 384  * full-dimensional approximation.
 385  */
 386 static Polyhedron *internal_polytope(Polyhedron *P,
 387                                      struct barvinok_options *options)
 388 {
 389     int i, j;
 390     Polyhedron *init;
 391     Matrix *vertices = Matrix_Alloc(P->NbRays, P->Dimension+2);
 392
 393     for (i = 0, j = 0; i < P->NbRays; ++i) {
 394         if (value_notone_p(P->Ray[i][1+P->Dimension]))
 395             continue;
 396         Vector_Copy(P->Ray[i], vertices->p[j++], 2+P->Dimension);
 397     }
 398     vertices->NbRows = j;
 399     if (!j) {
 400         Vector *sample = Polyhedron_Sample(P, options);
 401         if (sample) {
 402             value_set_si(vertices->p[0][1], 1);
 403             Vector_Copy(sample->p, vertices->p[0]+1, P->Dimension+1);
 404             Vector_Free(sample);
 405             vertices->NbRows = 1;
 406         }
 407     }
 408     init = Rays2Polyhedron(vertices, options->MaxRays);
 409     Matrix_Free(vertices);
 410
 411     while (!emptyQ(init) && init->NbEq) {
 412         Vector *sample;
 413         Vector *v = Vector_Alloc(init->Dimension+2);
 414         Polyhedron *Q;
 415
 416         value_set_si(v->p[0], 1);
 417         Vector_Copy(init->Constraint[0]+1, v->p+1, init->Dimension+1);
 418         value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
 419
 420         Q = AddConstraints(v->p, 1, P, options->MaxRays);
 421         sample = Polyhedron_Sample(Q, options);
 422         Polyhedron_Free(Q);
 423         if (!sample) {
 424             Vector_Oppose(init->Constraint[0]+1, v->p+1, init->Dimension+1);
 425             value_decrement(v->p[1+init->Dimension],v->p[1+init->Dimension]);
 426
 427             Q = AddConstraints(v->p, 1, P, options->MaxRays);
 428             sample = Polyhedron_Sample(Q, options);
 429             Polyhedron_Free(Q);
 430         }
 431         if (!sample) {
 432             P = AddConstraints(init->Constraint[0], 1, P, options->MaxRays);
 433             Polyhedron_Free(init);
 434             Vector_Free(v);
 435             return P;
 436         }
 437         assert(sample);
 438
 439         Vector_Copy(sample->p, v->p+1, init->Dimension+1);
 440         Q = AddRays(v->p, 1, init, options->MaxRays);
 441         Polyhedron_Free(init);
 442         init = Q;
 443
 444         Vector_Free(sample);
 445         Vector_Free(v);
 446     }
 447
 448     return init;
 449 }
 450
 451 /* Set the lower bound for optimizing along the normal of a facet
 452  * in case of computing the integer hull of a polyhedron.
 453  * Although obj contains the full affine objective function,
 454  * the calling function expect a lower bound on the linear part only.
 455  * We therefore need to subtract the constant from the computed
 456  * optimum of the linear relaxation.
 457  */
 458 static void set_lower(struct integer_hull *hull, Value *obj,
 459                       Value *lower, struct barvinok_options *options)
 460 {
 461     enum lp_result res;
 462     Value one;
 463
 464     value_init(one);
 465     value_set_si(one, 1);
 466
 467     res = polyhedron_opt(hull->P, obj, one, lp_min, lower, options);
 468     value_subtract(*lower, *lower, obj[hull->P->Dimension]);
 469     assert(res == lp_ok);
 470
 471     value_clear(one);
 472 }
 473
 474 /* Computes the vertices of the integer hull of P
 475  */
 476 Matrix *Polyhedron_Integer_Hull(Polyhedron *P,
 477                                 struct barvinok_options *options)
 478 {
 479     Matrix *vertices;
 480     Matrix *CV = NULL;
 481     Polyhedron *init;
 482
 483     POL_ENSURE_VERTICES(P);
 484     remove_all_equalities(&P, NULL, NULL, &CV, 0, options->MaxRays);
 485
 486     POL_ENSURE_VERTICES(P);
 487
 488     if (P->Dimension == 0) {
 489         vertices = Matrix_Alloc(1, 1);
 490         value_set_si(vertices->p[0][0], 1);
 491     } else {
 492         init = internal_polytope(P, options);
 493         if (emptyQ(init)) {
 494             vertices = Matrix_Alloc(0, P->Dimension+1);
 495             Polyhedron_Free(init);
 496         } else if (init->NbEq) {
 497             vertices = Polyhedron_Integer_Hull(init, options);
 498             Polyhedron_Free(init);
 499         } else {
 500             struct integer_hull hull;
 501             hull.P = P;
 502             hull.init = init;
 503             hull.set_lower_bound = set_lower;
 504             vertices = gbr_hull_extend(&hull, options);
 505         }
 506     }
 507
 508     if (CV) {
 509         vertices = transform_list_of_vertices(vertices, CV);
 510         Polyhedron_Free(P);
 511     }
 512
 513     return vertices;
 514 }