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 "remove_equalities.h"
  10 #include "config.h"
  11
  12 #ifndef USE_ZSOLVE
  13 Matrix *Cone_Hilbert_Integer_Hull(Polyhedron *C,
  14                                   struct barvinok_options *options)
  15 {
  16     assert(0);
  17 }
  18 #endif
  19
  20 struct integer_hull {
  21     Polyhedron  *P;     /* Original polyhedron or cone */
  22     Polyhedron  *init;  /* Initial (under)approximation of integer hull */
  23
  24     Matrix      *F;     /* Set of already computed facets */
  25     int         n_F;    /* The number of computed facets  */
  26
  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.
  30      */
  31     void (*set_lower_bound)(struct integer_hull *hull, Value *obj,
  32                             Value *lower, struct barvinok_options *options);
  33 };
  34
  35 static int matrix_add(Matrix *M, int n, Value *v)
  36 {
  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;
  41 }
  42
  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.
  48  */
  49 static int select_best(struct integer_hull *hull,
  50                        Polyhedron *P,
  51                        Value *min, Value *max,
  52                        struct barvinok_options *options)
  53 {
  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;
  60
  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++;
  78     }
  79
  80     if (i_min != -1) {
  81         value_decrement(*max, upper->p[i_min]);
  82         value_assign(*min, lower->p[i_min]);
  83     }
  84
  85     Vector_Free(lower);
  86     Vector_Free(upper);
  87     Vector_Free(distances);
  88
  89     return i_min;
  90 }
  91
  92 /* Return the (integer) vertices of P */
  93 static Matrix *Polyhedron_Vertices(Polyhedron *P)
  94 {
  95     int i, j;
  96     Matrix *M = Matrix_Alloc(P->NbRays, P->Dimension+1);
  97
  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);
 104     }
 105     M->NbRows = j;
 106
 107     return M;
 108 }
 109
 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.
 119  */
 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).
 128  */
 129 #define INCLUDE_KNOWN_FACETS_IN_ILP     0
 130
 131 Polyhedron *add_facets(Polyhedron *P, Matrix *facets, int n,
 132                         struct barvinok_options *options)
 133 {
 134     int i;
 135     Polyhedron *R;
 136     Matrix *M;
 137
 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;
 146 }
 147
 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.
 152  */
 153 Polyhedron *add_points(Polyhedron *Q, Vector *opt, Matrix *subopt,
 154                         int n_subopt, Matrix *rays,
 155                         struct barvinok_options *options)
 156 {
 157     int i;
 158     Polyhedron *R;
 159
 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);
 166         }
 167     }
 168
 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);
 172
 173     R = AddRays(rays->p[0], 1+n_subopt, Q, options->MaxRays);
 174     Polyhedron_Free(Q);
 175
 176     return R;
 177 }
 178
 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.
 187  */
 188 static Matrix *gbr_hull_extend(struct integer_hull *hull,
 189                                 struct barvinok_options *options)
 190 {
 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;
 198
 199     hull->F= Matrix_Alloc(Q->NbConstraints, 2+Q->Dimension);
 200     hull->n_F = 0;
 201
 202     candidates = Matrix_Alloc(0, Q->Dimension);
 203
 204     value_init(min);
 205     value_init(max);
 206
 207     value_set_si(rays->p[0][0], 1);
 208     value_set_si(rays->p[0][1+Q->Dimension], 1);
 209
 210     while ((i_min = select_best(hull, Q, &min, &max, options)) != -1) {
 211         Polyhedron *R;
 212         int i;
 213         Vector *opt;
 214
 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;
 219
 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);
 227
 228         hull->n_F = matrix_add(hull->F, hull->n_F, Q->Constraint[i_min]);
 229
 230         if (!opt)
 231             continue;
 232
 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]);
 237
 238         Q = add_points(Q, opt, candidates, n_candidates, rays, options);
 239
 240         Vector_Free(opt);
 241     }
 242
 243     vertices = Polyhedron_Vertices(Q);
 244
 245     Polyhedron_Free(Q);
 246     hull->init = NULL;
 247
 248     value_clear(min);
 249     value_clear(max);
 250
 251     Matrix_Free(rays);
 252     Matrix_Free(hull->F);
 253     hull->F = NULL;
 254     Matrix_Free(candidates);
 255
 256     return vertices;
 257 }
 258
 259 /* Returns the Minkowski sum of the cone and the polytope
 260  * that is the convex hull of its (integer) extremal rays.
 261  */
 262 static Polyhedron *truncate_cone(Polyhedron *C,
 263                                  struct barvinok_options *options)
 264 {
 265     int i;
 266     Matrix *V;
 267     Polyhedron *P;
 268
 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);
 276     }
 277     P = Rays2Polyhedron(V, options->MaxRays);
 278     Matrix_Free(V);
 279     return P;
 280 }
 281
 282 /* Frees original list and transformation matrix CV */
 283 static Matrix *transform_list_of_vertices(Matrix *list, Matrix *CV)
 284 {
 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);
 292
 293     return M;
 294 }
 295
 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.
 304  */
 305 static void set_to_one(struct integer_hull *hull, Value *obj,
 306                        Value *lower, struct barvinok_options *options)
 307 {
 308     if (INCLUDE_KNOWN_FACETS_IN_LP) {
 309         int i;
 310         enum lp_result res;
 311         Value one;
 312
 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);
 318
 319         value_init(one);
 320         value_set_si(one, 1);
 321
 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);
 325
 326         if (value_notzero_p(*lower))
 327             return;
 328     }
 329
 330     if (value_zero_p(obj[hull->P->Dimension]))
 331         value_set_si(*lower, 0);
 332     else
 333         value_set_si(*lower, 1);
 334 }
 335
 336 /* Add to Q those points that the calling function already knows about,
 337  * replacing (and destroying) the original Q.
 338  */
 339 static Polyhedron *add_known_points(Polyhedron *Q,
 340                                         Matrix *c, int n_c, unsigned MaxRays)
 341 {
 342     Polyhedron *R;
 343
 344     R = AddRays(c->p[0], n_c, Q, MaxRays);
 345     Polyhedron_Free(Q);
 346     return R;
 347 }
 348
 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.
 352  */
 353 static Matrix *gbr_cone_hull(Polyhedron *C, Matrix *c, int n_c,
 354                                 struct barvinok_options *options)
 355 {
 356     Matrix *vertices;
 357     Matrix *CV = NULL;
 358     struct integer_hull hull;
 359
 360     POL_ENSURE_VERTICES(C);
 361     remove_all_equalities(&C, NULL, NULL, &CV, 0, options->MaxRays);
 362
 363     POL_ENSURE_VERTICES(C);
 364
 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);
 371
 372     if (CV) {
 373         vertices = transform_list_of_vertices(vertices, CV);
 374         Polyhedron_Free(C);
 375     }
 376
 377     return vertices;
 378 }
 379
 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)
 383 {
 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);
 391     }
 392 }
 393
 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.
 403  */
 404 static Polyhedron *internal_polytope(Polyhedron *P,
 405                                      struct barvinok_options *options)
 406 {
 407     int i, j;
 408     Polyhedron *init;
 409     Matrix *vertices = Matrix_Alloc(P->NbRays, P->Dimension+2);
 410
 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);
 415     }
 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;
 424         }
 425     }
 426     init = Rays2Polyhedron(vertices, options->MaxRays);
 427     Matrix_Free(vertices);
 428
 429     while (!emptyQ(init) && init->NbEq) {
 430         Vector *sample;
 431         Vector *v = Vector_Alloc(init->Dimension+2);
 432         Polyhedron *Q;
 433
 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]);
 437
 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]);
 444
 445             Q = AddConstraints(v->p, 1, P, options->MaxRays);
 446             sample = Polyhedron_Sample(Q, options);
 447             Polyhedron_Free(Q);
 448         }
 449         if (!sample) {
 450             P = AddConstraints(init->Constraint[0], 1, P, options->MaxRays);
 451             Polyhedron_Free(init);
 452             Vector_Free(v);
 453             return P;
 454         }
 455         assert(sample);
 456
 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;
 461
 462         Vector_Free(sample);
 463         Vector_Free(v);
 464     }
 465
 466     return init;
 467 }
 468
 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.
 475  */
 476 static void set_lower(struct integer_hull *hull, Value *obj,
 477                       Value *lower, struct barvinok_options *options)
 478 {
 479     enum lp_result res;
 480     Value one;
 481
 482     value_init(one);
 483     value_set_si(one, 1);
 484
 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);
 488
 489     value_clear(one);
 490 }
 491
 492 /* Computes the vertices of the integer hull of P
 493  */
 494 Matrix *Polyhedron_Integer_Hull(Polyhedron *P,
 495                                 struct barvinok_options *options)
 496 {
 497     Matrix *vertices;
 498     Matrix *CV = NULL;
 499     Polyhedron *init;
 500
 501     POL_ENSURE_VERTICES(P);
 502     remove_all_equalities(&P, NULL, NULL, &CV, 0, options->MaxRays);
 503
 504     POL_ENSURE_VERTICES(P);
 505
 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);
 523         }
 524     }
 525
 526     if (CV) {
 527         vertices = transform_list_of_vertices(vertices, CV);
 528         Polyhedron_Free(P);
 529     }
 530
 531     return vertices;
 532 }