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