lattice_width.c

   1 #include <assert.h>
   2 #include <stdlib.h>
   3 #include <isl_set_polylib.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/util.h>
   6 #include "hilbert.h"
   7 #include "hull.h"
   8 #include "lattice_width.h"
   9 #include "param_util.h"
  10 #include "reduce_domain.h"
  11
  12 #define ALLOC(type) (type*)malloc(sizeof(type))
  13 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  14
  15 static void clear_width_direction(struct width_direction *wd)
  16 {
  17     Vector_Free(wd->width);
  18     Vector_Free(wd->dir);
  19     if (wd->domain)
  20         Polyhedron_Free(wd->domain);
  21 }
  22
  23 static struct width_direction_array *new_width_direction_array(void)
  24 {
  25     struct width_direction_array *dirs = ALLOC(struct width_direction_array);
  26     assert(dirs);
  27     dirs->n = 0;
  28     dirs->alloc = 32;
  29     dirs->wd = ALLOCN(struct width_direction, dirs->alloc);
  30     assert(dirs->wd);
  31     return dirs;
  32 }
  33
  34 static void grow_width_direction_array(struct width_direction_array *dirs,
  35                                        int extra)
  36 {
  37     if (dirs->n + extra <= dirs->alloc)
  38         return;
  39     dirs->alloc = (5*(dirs->n+extra))/4;
  40     dirs->wd = (struct width_direction*)realloc(dirs->wd,
  41                             dirs->alloc * sizeof(struct width_direction));
  42     assert(dirs->wd);
  43 }
  44
  45 void free_width_direction_array(struct width_direction_array *dirs)
  46 {
  47     int i;
  48
  49     for (i = 0; i < dirs->n; ++i)
  50         clear_width_direction(&dirs->wd[i]);
  51     free(dirs->wd);
  52     free(dirs);
  53 }
  54
  55 #define INT_BITS (sizeof(unsigned) * 8)
  56
  57 /* For each parametric vertex, compute cone of directions
  58  * for which this vertex attains the minimal value.
  59  */
  60 static Matrix **compute_vertex_dirs(Param_Polyhedron *PP)
  61 {
  62     int i;
  63     unsigned nvar = PP->V->Vertex->NbRows;
  64     Param_Vertices *V;
  65     Matrix **vertex_dirs = ALLOCN(Matrix *, PP->nbV);
  66
  67     for (i = 0, V = PP->V; V; ++i, V = V->next) {
  68         int kx;
  69         unsigned bx;
  70         int j, k;
  71         unsigned *facets;
  72         int n;
  73         Matrix *M;
  74         Polyhedron *P;
  75
  76         if (V->Facets) {
  77             int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
  78             facets = V->Facets;
  79             n = bit_vector_count(facets, len);
  80         } else
  81             facets = supporting_constraints(PP->Constraints, V, &n);
  82         M = Matrix_Alloc(n, 1+nvar+1);
  83         for (k = 0, j = 0, kx = 0, bx = MSB; j < n; ++k) {
  84             if (facets[kx] & bx) {
  85                 value_set_si(M->p[j][0], 1);
  86                 Vector_Copy(PP->Constraints->p[k]+1, M->p[j++]+1, nvar);
  87             }
  88             NEXT(kx, bx);
  89         }
  90         P = Constraints2Polyhedron(M, 0);
  91         Matrix_Free(M);
  92         vertex_dirs[i] = Matrix_Alloc(P->NbRays-1, nvar);
  93         for (k = 0, j = 0; k < P->NbRays; ++k) {
  94             if (value_notzero_p(P->Ray[k][1+nvar]))
  95                 continue;
  96             Vector_Copy(P->Ray[k]+1, vertex_dirs[i]->p[j++], nvar);
  97         }
  98         Polyhedron_Free(P);
  99         if (!V->Facets)
 100             free(facets);
 101     }
 102     return vertex_dirs;
 103 }
 104
 105 /* Compute
 106  *
 107  *       c     a     b
 108  *      --- = --- - ---
 109  *      c_d   a_d   b_d
 110  */
 111 static void Vector_Subtract(Value *a, Value a_d,
 112                             Value *b, Value b_d,
 113                             Value *c, Value *c_d, int len)
 114 {
 115     Value ma, mb;
 116     value_init(ma);
 117     value_init(mb);
 118     value_lcm(*c_d, a_d, b_d);
 119     value_divexact(ma, *c_d, a_d);
 120     value_divexact(mb, *c_d, b_d);
 121     value_oppose(mb, mb);
 122     Vector_Combine(a, b, c, ma, mb, len);
 123     value_clear(ma);
 124     value_clear(mb);
 125 }
 126
 127 /* Compute width for a given direction dir and initialize width_direction
 128  * structure.
 129  */
 130 static void compute_width_direction(Matrix *V_min, Matrix *V_max,
 131                                     Value *dir, struct width_direction *wd)
 132 {
 133     Vector *max = Vector_Alloc(V_min->NbColumns);
 134     unsigned nvar = V_min->NbRows;
 135     unsigned nparam = V_min->NbColumns-2;
 136
 137     wd->width = Vector_Alloc(V_min->NbColumns);
 138     wd->dir = Vector_Alloc(nvar);
 139     Vector_Copy(dir, wd->dir->p, nvar);
 140     wd->domain = NULL;
 141
 142     V_min->NbColumns--;
 143     V_max->NbColumns--;
 144
 145     Vector_Matrix_Product(dir, V_max, max->p);
 146     Vector_Matrix_Product(dir, V_min, wd->width->p);
 147     Vector_Subtract(max->p, V_max->p[0][V_max->NbColumns],
 148                     wd->width->p, V_min->p[0][V_min->NbColumns],
 149                     wd->width->p, &wd->width->p[nparam+1],
 150                     nparam+1);
 151
 152     V_min->NbColumns++;
 153     V_max->NbColumns++;
 154
 155     Vector_Normalize(wd->width->p, nparam+2);
 156
 157     Vector_Free(max);
 158 }
 159
 160 static int Vector_Compare(Value *p1, Value *p2, unsigned len)
 161 {
 162     int i;
 163
 164     for (i = 0; i < len; ++i) {
 165         int sign = mpz_cmp(p1[i], p2[i]);
 166         if (sign)
 167             return sign;
 168     }
 169     return 0;
 170 }
 171
 172 static int wd_width_lex_cmp(const void *va, const void *vb)
 173 {
 174     const struct width_direction *a = (const struct width_direction *)va;
 175     const struct width_direction *b = (const struct width_direction *)vb;
 176
 177     return Vector_Compare(a->width->p, b->width->p, a->width->Size);
 178 }
 179
 180 static int add_vertex(Matrix *M, int n, Value *v)
 181 {
 182     if (n >= M->NbRows)
 183         Matrix_Extend(M, 3*(M->NbRows+10)/2);
 184     value_set_si(M->p[n][0], 1);
 185     Vector_Copy(v, M->p[n]+1, M->NbColumns-2);
 186     value_set_si(M->p[n][M->NbColumns-1], 1);
 187     return n+1;
 188 }
 189
 190 /* Puts the points in v that lie in P in front of the list
 191  * and returns their number.
 192  */
 193 static int valid_vertices(Polyhedron *P, Matrix *v, int n_v)
 194 {
 195     int i, j, k;
 196     Value tmp;
 197
 198     assert(v->NbColumns == P->Dimension+2);
 199     value_init(tmp);
 200
 201     for (j = 0, k = 0; j < n_v; ++j) {
 202         for (i = 0; i < P->NbConstraints; ++i) {
 203             Inner_Product(v->p[j]+1, P->Constraint[i]+1, P->Dimension+1, &tmp);
 204             if (value_neg_p(tmp))
 205                 break;
 206         }
 207         if (i < P->NbConstraints)
 208             continue;
 209         if (j != k)
 210             Vector_Exchange(v->p[j]+1, v->p[k]+1, P->Dimension);
 211         ++k;
 212     }
 213
 214     value_clear(tmp);
 215     return k;
 216 }
 217
 218 static struct width_direction_array *
 219 compute_width_directions(Param_Polyhedron *PP, struct barvinok_options *options)
 220 {
 221     Matrix **vertex_dirs;
 222     Param_Vertices *V_max, *V_min;
 223     int i, V_max_i, V_min_i;
 224     unsigned nvar = PP->V->Vertex->NbRows;
 225     struct width_direction_array *width_dirs = new_width_direction_array();
 226     Matrix *all_vertices = Matrix_Alloc(nvar, 1+nvar+1);
 227     int n_vertices = 0;
 228
 229     vertex_dirs = compute_vertex_dirs(PP);
 230
 231     for (V_max = PP->V; V_max; V_max = V_max->next)
 232         Param_Vertex_Common_Denominator(V_max);
 233
 234     for (V_max = PP->V, V_max_i = 0; V_max; V_max = V_max->next, V_max_i++) {
 235         for (V_min = V_max->next, V_min_i = V_max_i+1;
 236                 V_min;
 237                 V_min = V_min->next, V_min_i++) {
 238             Matrix *M;
 239             Matrix *basis;
 240             Polyhedron *C;
 241             unsigned V_max_n = vertex_dirs[V_max_i]->NbRows;
 242             unsigned V_min_n = vertex_dirs[V_min_i]->NbRows;
 243             int n_valid;
 244
 245             if (options->verbose)
 246                 fprintf(stderr, "%d/%d %d/%d %d \r",
 247                                     V_max_i, PP->nbV,
 248                                     V_min_i, PP->nbV,
 249                                     width_dirs->n);
 250
 251             M = Matrix_Alloc(V_max_n+V_min_n, 1+nvar+1);
 252             for (i = 0; i < V_max_n; ++i) {
 253                 value_set_si(M->p[i][0], 1);
 254                 Vector_Oppose(vertex_dirs[V_max_i]->p[i], M->p[i]+1, nvar);
 255             }
 256             for (i = 0; i < V_min_n; ++i) {
 257                 value_set_si(M->p[V_max_n+i][0], 1);
 258                 Vector_Copy(vertex_dirs[V_min_i]->p[i], M->p[V_max_n+i]+1, nvar);
 259             }
 260             C = Constraints2Polyhedron(M, options->MaxRays);
 261             Matrix_Free(M);
 262             n_valid = valid_vertices(C, all_vertices, n_vertices);
 263             basis = Cone_Integer_Hull(C, all_vertices, n_valid, options);
 264             grow_width_direction_array(width_dirs, basis->NbRows);
 265             for (i = 0; i < basis->NbRows; ++i) {
 266                 Matrix *VM_min, *VM_max;
 267                 int pos;
 268
 269                 VM_min = V_min->Vertex;
 270                 VM_max = V_max->Vertex;
 271                 pos = First_Non_Zero(basis->p[i], nvar);
 272                 if (value_neg_p(basis->p[i][pos])) {
 273                     Vector_Oppose(basis->p[i], basis->p[i], nvar);
 274                     VM_min = V_max->Vertex;
 275                     VM_max = V_min->Vertex;
 276                 }
 277
 278                 n_vertices = add_vertex(all_vertices, n_vertices, basis->p[i]);
 279                 compute_width_direction(VM_min, VM_max, basis->p[i],
 280                                         &width_dirs->wd[width_dirs->n++]);
 281             }
 282             Matrix_Free(basis);
 283             Polyhedron_Free(C);
 284         }
 285     }
 286     Matrix_Free(all_vertices);
 287
 288     for (i = 0; i < PP->nbV; ++i)
 289         Matrix_Free(vertex_dirs[i]);
 290     free(vertex_dirs);
 291
 292     return width_dirs;
 293 }
 294
 295 /* Computes the lattice width direction of a parametric polytope.
 296  * The parameter space is allowed to be unbounded.
 297  * Currently, the parametric polytope and the parameter space
 298  * are assumed to be full-dimensional.
 299  *
 300  * First, we compute the parametric vertices.
 301  * Then, for each pair of vertices, we construct a (rational) cone
 302  * of directions for which one vertex attains the minimal value
 303  * and the other vertex attains the maximal value.
 304  * The candidate directions are the elements of the integer hulls
 305  * of these cones.
 306  * The minimal direction is then obtained by computing the
 307  * region in the parameter space where each direction yields
 308  * a smaller (or equal) width than all the other directions.
 309  *
 310  * In principle, we can avoid computing candidate directions
 311  * for vertices with no overlapping activity domains (possibly
 312  * after opening some facets of the activity domains in the
 313  * familiar way).
 314  *
 315  * The output is a list of triples, consisting of a direction,
 316  * the corresponding width and the chamber in the parameter
 317  * space where this direction leads to the minimal width.
 318  *
 319  * The algorithm is described in "Integer points in a parameterised
 320  * polyhedron" by Friedrich Eisenbrand and Gennady Shmonin.
 321  */
 322 struct width_direction_array *
 323 Polyhedron_Lattice_Width_Directions(Polyhedron *P, Polyhedron *C,
 324                                     struct barvinok_options *options)
 325 {
 326     Param_Polyhedron *PP;
 327     unsigned nparam = C->Dimension;
 328     int i, j, k;
 329     struct width_direction_array *width_dirs;
 330     Polyhedron *TC;
 331     Vector *inner;
 332
 333     assert(P->NbEq == 0);
 334     assert(C->NbEq == 0);
 335
 336     /* Use true context since the algorithm assumes P is non-empty
 337      * for every point in the context.
 338      */
 339     TC = true_context(P, C, options->MaxRays);
 340     inner = inner_point(TC);
 341
 342     /* This is overkill, as we discard the computed chambers. */
 343     PP = Polyhedron2Param_Polyhedron(P, TC, options);
 344
 345     width_dirs = compute_width_directions(PP, options);
 346     Param_Polyhedron_Free(PP);
 347
 348     qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
 349             wd_width_lex_cmp);
 350
 351     for (i = 1, j = 1; i < width_dirs->n; ++i) {
 352         /* We could also weed out width_directions that differ by a
 353          * (positive) constant from another width_direction, but then
 354          * we'd have to put the two width_directions on a common
 355          * denominator first.
 356          */
 357         if (Vector_Equal(width_dirs->wd[j-1].width->p,
 358                          width_dirs->wd[i].width->p, nparam+2))
 359             clear_width_direction(&width_dirs->wd[i]);
 360         else
 361             width_dirs->wd[j++] = width_dirs->wd[i];
 362     }
 363     width_dirs->n = j;
 364
 365     for (i = 0, k = 0; i < width_dirs->n; ++i) {
 366         Matrix *M = Matrix_Alloc(TC->NbConstraints+width_dirs->n-(i-k)-1, nparam+2);
 367         for (j = 0; j < TC->NbConstraints; ++j)
 368             Vector_Copy(TC->Constraint[j], M->p[j], nparam+2);
 369         for (j = 0; j < width_dirs->n; ++j) {
 370             int pos;
 371             if (k <= j && j <= i)
 372                 continue;
 373             if (j < k)
 374                 pos = TC->NbConstraints + j;
 375             else
 376                 pos = TC->NbConstraints + j - (i-k) - 1;
 377             Vector_Subtract(width_dirs->wd[j].width->p,
 378                             width_dirs->wd[j].width->p[nparam+1],
 379                             width_dirs->wd[i].width->p,
 380                             width_dirs->wd[i].width->p[nparam+1],
 381                             M->p[pos]+1, M->p[pos], nparam+1);
 382             value_set_si(M->p[pos][0], 1);
 383             Vector_Normalize(M->p[pos]+1, nparam+1);
 384             if (!is_internal(inner, M->p[pos]))
 385                 value_decrement(M->p[pos][nparam+1], M->p[pos][nparam+1]);
 386         }
 387         width_dirs->wd[i].domain = Constraints2Polyhedron(M, options->MaxRays);
 388         if (emptyQ(width_dirs->wd[i].domain))
 389             clear_width_direction(&width_dirs->wd[i]);
 390         else
 391             width_dirs->wd[k++] = width_dirs->wd[i];
 392         Matrix_Free(M);
 393     }
 394     width_dirs->n = k;
 395     Vector_Free(inner);
 396     Polyhedron_Free(TC);
 397
 398     return width_dirs;
 399 }
 400
 401 /* Construct evalue of chambers with their associated widths */
 402 evalue *Polyhedron_Lattice_Width(Polyhedron *P, Polyhedron *C,
 403                                  struct barvinok_options *options)
 404 {
 405     evalue *width;
 406     struct evalue_section *s;
 407     struct width_direction_array *width_dirs;
 408     int i;
 409     unsigned nparam = C->Dimension;
 410
 411     width_dirs = Polyhedron_Lattice_Width_Directions(P, C, options);
 412     s = ALLOCN(struct evalue_section, width_dirs->n);
 413     for (i = 0; i < width_dirs->n; ++i) {
 414         s[i].D = width_dirs->wd[i].domain;
 415         width_dirs->wd[i].domain = NULL;
 416         s[i].E = affine2evalue(width_dirs->wd[i].width->p,
 417                                width_dirs->wd[i].width->p[nparam+1],
 418                                nparam);
 419     }
 420     free_width_direction_array(width_dirs);
 421
 422     width = evalue_from_section_array(s, i);
 423     free(s);
 424
 425     return width;
 426 }
 427
 428 __isl_give isl_pw_qpolynomial *isl_basic_set_lattice_width(
 429         __isl_take isl_basic_set *bset)
 430 {
 431         isl_ctx *ctx;
 432         isl_space *dim;
 433         isl_pw_qpolynomial *pwqp;
 434         unsigned nparam;
 435         Polyhedron *U;
 436         Polyhedron *P;
 437         evalue *E;
 438         struct barvinok_options *options;
 439         int options_allocated = 0;
 440
 441         if (!bset)
 442                 return NULL;
 443
 444         ctx = isl_basic_set_get_ctx(bset);
 445         options = isl_ctx_peek_barvinok_options(ctx);
 446         if (!options) {
 447                 options = barvinok_options_new_with_defaults();
 448                 options_allocated = 1;
 449         }
 450
 451         nparam = isl_basic_set_dim(bset, isl_dim_param);
 452         dim = isl_basic_set_get_space(bset);
 453         dim = isl_space_params(dim);
 454
 455         U = Universe_Polyhedron(nparam);
 456         P = isl_basic_set_to_polylib(bset);
 457
 458         E = Polyhedron_Lattice_Width(P, U, options);
 459
 460         pwqp = isl_pw_qpolynomial_from_evalue(dim, E);
 461         isl_basic_set_free(bset);
 462
 463         evalue_free(E);
 464         Polyhedron_Free(P);
 465         Polyhedron_Free(U);
 466         if (options_allocated)
 467                 barvinok_options_free(options);
 468
 469         return pwqp;
 470 }
 471
 472 __isl_give isl_pw_qpolynomial *isl_set_lattice_width(__isl_take isl_set *set)
 473 {
 474         if (!set)
 475                 return NULL;
 476
 477         if (isl_set_plain_is_empty(set)) {
 478                 isl_space *dim;
 479                 dim = isl_set_get_space(set);
 480                 dim = isl_space_domain(isl_space_from_range(dim));
 481                 isl_set_free(set);
 482                 return isl_pw_qpolynomial_zero(dim);
 483         }
 484
 485         if (isl_set_n_basic_set(set) != 1)
 486                 isl_die(isl_set_get_ctx(set), isl_error_unsupported,
 487                         "unions not supported (yet)", goto error);
 488
 489         return isl_basic_set_lattice_width(isl_set_simple_hull(set));
 490 error:
 491         isl_set_free(set);
 492         return NULL;
 493 }
 494
 495 static isl_stat set_lw(__isl_take isl_set *set, void *user)
 496 {
 497         isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
 498         isl_pw_qpolynomial *pwqp;
 499
 500         pwqp = isl_set_lattice_width(set);
 501         *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, pwqp);
 502
 503         return isl_stat_ok;
 504 }
 505
 506 __isl_give isl_union_pw_qpolynomial *isl_union_set_lattice_width(
 507         __isl_take isl_union_set *uset)
 508 {
 509         isl_space *dim;
 510         isl_union_pw_qpolynomial *res;
 511
 512         dim = isl_union_set_get_space(uset);
 513         res = isl_union_pw_qpolynomial_zero(dim);
 514         if (isl_union_set_n_set(uset) > 1)
 515                 isl_die(isl_union_set_get_ctx(uset), isl_error_unsupported,
 516                         "unions not supported (yet)", goto error);
 517         if (isl_union_set_foreach_set(uset, &set_lw, &res) < 0)
 518                 goto error;
 519         isl_union_set_free(uset);
 520
 521         return res;
 522 error:
 523         isl_union_set_free(uset);
 524         isl_union_pw_qpolynomial_free(res);
 525         return NULL;
 526 }