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 #include "config.h"
  10
  11 #ifndef USE_ZSOLVE
  12 Matrix *Cone_Hilbert_Integer_Hull(Polyhedron *C,
  13                                   struct barvinok_options *options)
  14 {
  15     assert(0);
  16 }
  17 #endif
  18
  19 struct integer_hull {
  20     Polyhedron  *P;     /* Original polyhedron or cone */
  21     Polyhedron  *init;  /* Initial (under)approximation of integer hull */
  22
  23     Matrix      *F;     /* Set of already computed facets */
  24     int         n_F;    /* The number of computed facets  */
  25
  26     /* Computes a lower bound for the objective function obj over
  27      * the integer hull, possibly exploiting the already computed
  28      * facets of the integer hull given in hull->F.
  29      */
  30     void (*set_lower_bound)(struct integer_hull *hull, Value *obj,
  31                             Value *lower, struct barvinok_options *options);
  32 };
  33
  34 static int matrix_add(Matrix *M, int n, Value *v)
  35 {
  36     if (n >= M->NbRows)
  37         Matrix_Extend(M, 3*(M->NbRows+10)/2);
  38     Vector_Copy(v, M->p[n], M->NbColumns);
  39     return n+1;
  40 }
  41
  42 /* Select "best" constraint from P over which to optimize,
  43  * where P is the current approximation of the integer hull
  44  * and where a constraint is better if the difference between
  45  * the value at the constraint and the lower bound (based on
  46  * the original polytope (the linear relaxation) is smaller.
  47  */
  48 static int select_best(struct integer_hull *hull,
  49                        Polyhedron *P,
  50                        Value *min, Value *max,
  51                        struct barvinok_options *options)
  52 {
  53     int i, j, k;
  54     Vector *lower = Vector_Alloc(P->NbConstraints);
  55     Vector *upper = Vector_Alloc(P->NbConstraints);
  56     Vector *distances = Vector_Alloc(P->NbConstraints);
  57     int i_min;
  58     int non_zero = 0;
  59
  60     i_min = -1;
  61     for (i = 0; i < P->NbConstraints; ++i) {
  62         if (First_Non_Zero(P->Constraint[i]+1, P->Dimension) == -1)
  63             continue;
  64         for (j = 0; j < hull->n_F; ++j)
  65             if (Vector_Equal(hull->F->p[j]+1, P->Constraint[i]+1, P->Dimension))
  66                 break;
  67         if (j < hull->n_F)
  68             continue;
  69         hull->set_lower_bound(hull, P->Constraint[i]+1, &lower->p[i], options);
  70         value_oppose(upper->p[i], P->Constraint[i][1+P->Dimension]);
  71         value_subtract(distances->p[i], upper->p[i], lower->p[i]);
  72         if (value_zero_p(distances->p[i]))
  73             continue;
  74         if (i_min == -1 || value_lt(distances->p[i], distances->p[i_min]))
  75             i_min = i;
  76         non_zero++;
  77     }
  78
  79     if (i_min != -1) {
  80         value_decrement(*max, upper->p[i_min]);
  81         value_assign(*min, lower->p[i_min]);
  82     }
  83
  84     Vector_Free(lower);
  85     Vector_Free(upper);
  86     Vector_Free(distances);
  87
  88     return i_min;
  89 }
  90
  91 /* Return the (integer) vertices of P */
  92 static Matrix *Polyhedron_Vertices(Polyhedron *P)
  93 {
  94     int i, j;
  95     Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+1);
  96
  97     for (i = 0, j = 0; i < P->NbRays; ++i) {
  98         if (value_zero_p(P->Ray[i][1+P->Dimension]))
  99             continue;
 100         assert(value_one_p(P->Ray[i][1+P->Dimension]));
 101         value_set_si(M->p[j][P->Dimension], 1);
 102         Vector_Copy(P->Ray[i]+1, M->p[j++], P->Dimension);
 103     }
 104     M->NbRows = j;
 105
 106     return M;
 107 }
 108
 109 /* Set INCLUDE_KNOWN_FACETS_IN_LP to 1 to compute a lower bound for
 110  * the minimal value wrt a facet based on the facets that have
 111  * already been determined to be facets of the integer hull,
 112  * rather than simply setting the lower bound to 0 or 1.
 113  * This does not appear to be profitable, probably because
 114  * the lower bound is recomputed every time we optimize
 115  * wrt a facet.  If we were to cache the lower bounds for
 116  * facets that do not change when a new vertex is added,
 117  * this may still turn out to be useful.
 118  */
 119 #define INCLUDE_KNOWN_FACETS_IN_LP      0
 120 /* Set INCLUDE_KNOWN_FACETS_IN_ILP to 1 to take into account the
 121  * facets that have already been determined to be facets
 122  * of the integer hull while optimizing wrt some other
 123  * facet.  The idea is that the polytope over which to
 124  * optimize could be flatter.  However, in practice, adding
 125  * these constraints seems to slow down the computation
 126  * (on at least one example).
 127  */
 128 #define INCLUDE_KNOWN_FACETS_IN_ILP     0
 129
 130 Polyhedron *add_facets(Polyhedron *P, Matrix *facets, int n,
 131                         struct barvinok_options *options)
 132 {
 133     int i;
 134     Polyhedron *R;
 135     Matrix *M;
 136
 137     M = Matrix_Alloc(P->NbConstraints + n, facets->NbColumns);
 138     for (i = 0; i < n; ++i)
 139         Vector_Copy(facets->p[i], M->p[i], facets->NbColumns);
 140     for (i = 0; i < P->NbConstraints; ++i)
 141         Vector_Copy(P->Constraint[i], M->p[n + i], M->NbColumns);
 142     R = Constraints2Polyhedron(M, options->MaxRays);
 143     Matrix_Free(M);
 144     return R;
 145 }
 146
 147 /* Add optimal and suboptimal points found by Polyhedron_Integer_Minimum
 148  * to Q, replacing (and destroying) the old Q.
 149  * rays is a temporary array for use in AddRays that can be reused
 150  * over multiple calls and that is extended as needed.
 151  */
 152 Polyhedron *add_points(Polyhedron *Q, Vector *opt, Matrix *subopt,
 153                         int n_subopt, Matrix *rays,
 154                         struct barvinok_options *options)
 155 {
 156     int i;
 157     Polyhedron *R;
 158
 159     if (1 + n_subopt > rays->NbRows) {
 160         int n = rays->NbRows;
 161         Matrix_Extend(rays, 1+n_subopt);
 162         for (i = n; i < rays->NbRows; ++i) {
 163             value_set_si(rays->p[i][0], 1);
 164             value_set_si(rays->p[i][1+Q->Dimension], 1);
 165         }
 166     }
 167
 168     Vector_Copy(opt->p, rays->p[0]+1, Q->Dimension);
 169     for (i = 0; i < n_subopt; ++i)
 170         Vector_Copy(subopt->p[i], rays->p[1+i]+1, Q->Dimension);
 171
 172     R = AddRays(rays->p[0], 1+n_subopt, Q, options->MaxRays);
 173     Polyhedron_Free(Q);
 174
 175     return R;
 176 }
 177
 178 /* Extends an initial approximation hull->init of the integer hull
 179  * of hull->P using generalized basis reduction.
 180  * In particular, it considers each facet of the current
 181  * approximation add optimizes in the direction of the normal,
 182  * adding the optimal point to the approximation as well as
 183  * all other points that may have been found along the way,
 184  * until no facet yields any more new points.
 185  * Returns a list of the vertices of this integer hull.
 186  */
 187 static Matrix *gbr_hull_extend(struct integer_hull *hull,
 188                                 struct barvinok_options *options)
 189 {
 190     Value min, max;
 191     Polyhedron *Q = hull->init;
 192     Matrix *rays = Matrix_Alloc(1, 2+Q->Dimension);
 193     Matrix *candidates;
 194     int n_candidates = 0;
 195     Matrix *vertices;
 196     int i_min;
 197
 198     hull->F= Matrix_Alloc(Q->NbConstraints, 2+Q->Dimension);
 199     hull->n_F = 0;
 200
 201     candidates = Matrix_Alloc(0, Q->Dimension);
 202
 203     value_init(min);
 204     value_init(max);
 205
 206     value_set_si(rays->p[0][0], 1);
 207     value_set_si(rays->p[0][1+Q->Dimension], 1);
 208
 209     while ((i_min = select_best(hull, Q, &min, &max, options)) != -1) {
 210         Polyhedron *R;
 211         int i;
 212         Vector *opt;
 213
 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;
 218
 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);
 226
 227         hull->n_F = matrix_add(hull->F, hull->n_F, Q->Constraint[i_min]);
 228
 229         if (!opt)
 230             continue;
 231
 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]);
 236
 237         Q = add_points(Q, opt, candidates, n_candidates, rays, options);
 238
 239         Vector_Free(opt);
 240     }
 241
 242     vertices = Polyhedron_Vertices(Q);
 243
 244     Polyhedron_Free(Q);
 245     hull->init = NULL;
 246
 247     value_clear(min);
 248     value_clear(max);
 249
 250     Matrix_Free(rays);
 251     Matrix_Free(hull->F);
 252     hull->F = NULL;
 253     Matrix_Free(candidates);
 254
 255     return vertices;
 256 }
 257
 258 /* Returns the Minkowski sum of the cone and the polytope
 259  * that is the convex hull of its (integer) extremal rays.
 260  */
 261 static Polyhedron *truncate_cone(Polyhedron *C,
 262                                  struct barvinok_options *options)
 263 {
 264     int i;
 265     Matrix *V;
 266     Polyhedron *P;
 267
 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);
 275     }
 276     P = Rays2Polyhedron(V, options->MaxRays);
 277     Matrix_Free(V);
 278     return P;
 279 }
 280
 281 /* Frees original list and transformation matrix CV */
 282 static Matrix *transform_list_of_vertices(Matrix *list, Matrix *CV)
 283 {
 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);
 291
 292     return M;
 293 }
 294
 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.
 303  */
 304 static void set_to_one(struct integer_hull *hull, Value *obj,
 305                        Value *lower, struct barvinok_options *options)
 306 {
 307     if (INCLUDE_KNOWN_FACETS_IN_LP) {
 308         int i;
 309         enum lp_result res;
 310         Value one;
 311
 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);
 317
 318         value_init(one);
 319         value_set_si(one, 1);
 320
 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);
 324
 325         if (value_notzero_p(*lower))
 326             return;
 327     }
 328
 329     if (value_zero_p(obj[hull->P->Dimension]))
 330         value_set_si(*lower, 0);
 331     else
 332         value_set_si(*lower, 1);
 333 }
 334
 335 /* Add to Q those points that the calling function already knows about,
 336  * replacing (and destroying) the original Q.
 337  */
 338 static Polyhedron *add_known_points(Polyhedron *Q,
 339                                         Matrix *c, int n_c, unsigned MaxRays)
 340 {
 341     Polyhedron *R;
 342
 343     R = AddRays(c->p[0], n_c, Q, MaxRays);
 344     Polyhedron_Free(Q);
 345     return R;
 346 }
 347
 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.
 351  */
 352 static Matrix *gbr_cone_hull(Polyhedron *C, Matrix *c, int n_c,
 353                                 struct barvinok_options *options)
 354 {
 355     Matrix *vertices;
 356     Matrix *CV = NULL;
 357     struct integer_hull hull;
 358
 359     POL_ENSURE_VERTICES(C);
 360     remove_all_equalities(&C, NULL, NULL, &CV, 0, options->MaxRays);
 361
 362     POL_ENSURE_VERTICES(C);
 363
 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);
 370
 371     if (CV) {
 372         vertices = transform_list_of_vertices(vertices, CV);
 373         Polyhedron_Free(C);
 374     }
 375
 376     return vertices;
 377 }
 378
 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)
 382 {
 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);
 390     }
 391 }
 392
 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.
 402  */
 403 static Polyhedron *internal_polytope(Polyhedron *P,
 404                                      struct barvinok_options *options)
 405 {
 406     int i, j;
 407     Polyhedron *init;
 408     Matrix *vertices = Matrix_Alloc(P->NbRays, P->Dimension+2);
 409
 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);
 414     }
 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;
 423         }
 424     }
 425     init = Rays2Polyhedron(vertices, options->MaxRays);
 426     Matrix_Free(vertices);
 427
 428     while (!emptyQ(init) && init->NbEq) {
 429         Vector *sample;
 430         Vector *v = Vector_Alloc(init->Dimension+2);
 431         Polyhedron *Q;
 432
 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]);
 436
 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]);
 443
 444             Q = AddConstraints(v->p, 1, P, options->MaxRays);
 445             sample = Polyhedron_Sample(Q, options);
 446             Polyhedron_Free(Q);
 447         }
 448         if (!sample) {
 449             P = AddConstraints(init->Constraint[0], 1, P, options->MaxRays);
 450             Polyhedron_Free(init);
 451             Vector_Free(v);
 452             return P;
 453         }
 454         assert(sample);
 455
 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;
 460
 461         Vector_Free(sample);
 462         Vector_Free(v);
 463     }
 464
 465     return init;
 466 }
 467
 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.
 474  */
 475 static void set_lower(struct integer_hull *hull, Value *obj,
 476                       Value *lower, struct barvinok_options *options)
 477 {
 478     enum lp_result res;
 479     Value one;
 480
 481     value_init(one);
 482     value_set_si(one, 1);
 483
 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);
 487
 488     value_clear(one);
 489 }
 490
 491 /* Computes the vertices of the integer hull of P
 492  */
 493 Matrix *Polyhedron_Integer_Hull(Polyhedron *P,
 494                                 struct barvinok_options *options)
 495 {
 496     Matrix *vertices;
 497     Matrix *CV = NULL;
 498     Polyhedron *init;
 499
 500     POL_ENSURE_VERTICES(P);
 501     remove_all_equalities(&P, NULL, NULL, &CV, 0, options->MaxRays);
 502
 503     POL_ENSURE_VERTICES(P);
 504
 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);
 522         }
 523     }
 524
 525     if (CV) {
 526         vertices = transform_list_of_vertices(vertices, CV);
 527         Polyhedron_Free(P);
 528     }
 529
 530     return vertices;
 531 }