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 wd_dir_lex_cmp(const void *va, const void *vb)
 181 {
 182     const struct width_direction *a = (const struct width_direction *)va;
 183     const struct width_direction *b = (const struct width_direction *)vb;
 184
 185     return Vector_Compare(a->dir->p, b->dir->p, a->dir->Size);
 186 }
 187
 188 static int add_vertex(Matrix *M, int n, Value *v)
 189 {
 190     if (n >= M->NbRows)
 191         Matrix_Extend(M, 3*(M->NbRows+10)/2);
 192     value_set_si(M->p[n][0], 1);
 193     Vector_Copy(v, M->p[n]+1, M->NbColumns-2);
 194     value_set_si(M->p[n][M->NbColumns-1], 1);
 195     return n+1;
 196 }
 197
 198 /* Puts the points in v that lie in P in front of the list
 199  * and returns their number.
 200  */
 201 static int valid_vertices(Polyhedron *P, Matrix *v, int n_v)
 202 {
 203     int i, j, k;
 204     Value tmp;
 205
 206     assert(v->NbColumns == P->Dimension+2);
 207     value_init(tmp);
 208
 209     for (j = 0, k = 0; j < n_v; ++j) {
 210         for (i = 0; i < P->NbConstraints; ++i) {
 211             Inner_Product(v->p[j]+1, P->Constraint[i]+1, P->Dimension+1, &tmp);
 212             if (value_neg_p(tmp))
 213                 break;
 214         }
 215         if (i < P->NbConstraints)
 216             continue;
 217         if (j != k)
 218             Vector_Exchange(v->p[j]+1, v->p[k]+1, P->Dimension);
 219         ++k;
 220     }
 221
 222     value_clear(tmp);
 223     return k;
 224 }
 225
 226 static struct width_direction_array *
 227 compute_width_directions(Param_Polyhedron *PP, struct barvinok_options *options)
 228 {
 229     Matrix **vertex_dirs;
 230     Param_Vertices *V_max, *V_min;
 231     int i, V_max_i, V_min_i;
 232     unsigned nvar = PP->V->Vertex->NbRows;
 233     struct width_direction_array *width_dirs = new_width_direction_array();
 234     Matrix *all_vertices = Matrix_Alloc(nvar, 1+nvar+1);
 235     int n_vertices = 0;
 236
 237     vertex_dirs = compute_vertex_dirs(PP);
 238
 239     for (V_max = PP->V; V_max; V_max = V_max->next)
 240         Param_Vertex_Common_Denominator(V_max);
 241
 242     for (V_max = PP->V, V_max_i = 0; V_max; V_max = V_max->next, V_max_i++) {
 243         for (V_min = V_max->next, V_min_i = V_max_i+1;
 244                 V_min;
 245                 V_min = V_min->next, V_min_i++) {
 246             Matrix *M;
 247             Matrix *basis;
 248             Polyhedron *C;
 249             unsigned V_max_n = vertex_dirs[V_max_i]->NbRows;
 250             unsigned V_min_n = vertex_dirs[V_min_i]->NbRows;
 251             int sorted_n;
 252             int n_valid;
 253
 254             if (options->verbose)
 255                 fprintf(stderr, "%d/%d %d/%d %d \r",
 256                                     V_max_i, PP->nbV,
 257                                     V_min_i, PP->nbV,
 258                                     width_dirs->n);
 259
 260             M = Matrix_Alloc(V_max_n+V_min_n, 1+nvar+1);
 261             for (i = 0; i < V_max_n; ++i) {
 262                 value_set_si(M->p[i][0], 1);
 263                 Vector_Oppose(vertex_dirs[V_max_i]->p[i], M->p[i]+1, nvar);
 264             }
 265             for (i = 0; i < V_min_n; ++i) {
 266                 value_set_si(M->p[V_max_n+i][0], 1);
 267                 Vector_Copy(vertex_dirs[V_min_i]->p[i], M->p[V_max_n+i]+1, nvar);
 268             }
 269             C = Constraints2Polyhedron(M, options->MaxRays);
 270             Matrix_Free(M);
 271             n_valid = valid_vertices(C, all_vertices, n_vertices);
 272             basis = Cone_Integer_Hull(C, all_vertices, n_valid, options);
 273             grow_width_direction_array(width_dirs, basis->NbRows);
 274             qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
 275                     wd_dir_lex_cmp);
 276             sorted_n = width_dirs->n;
 277             for (i = 0; i < basis->NbRows; ++i) {
 278                 Vector v;
 279                 struct width_direction wd;
 280                 Matrix *VM_min, *VM_max;
 281                 int pos;
 282
 283                 VM_min = V_min->Vertex;
 284                 VM_max = V_max->Vertex;
 285                 pos = First_Non_Zero(basis->p[i], nvar);
 286                 if (value_neg_p(basis->p[i][pos])) {
 287                     Vector_Oppose(basis->p[i], basis->p[i], nvar);
 288                     VM_min = V_max->Vertex;
 289                     VM_max = V_min->Vertex;
 290                 }
 291
 292                 v.Size = nvar;
 293                 v.p = basis->p[i];
 294                 wd.dir = &v;
 295                 if (bsearch(&wd, width_dirs->wd, sorted_n,
 296                             sizeof(struct width_direction),
 297                             wd_dir_lex_cmp))
 298                     continue;
 299
 300                 n_vertices = add_vertex(all_vertices, n_vertices, basis->p[i]);
 301                 compute_width_direction(VM_min, VM_max, basis->p[i],
 302                                         &width_dirs->wd[width_dirs->n++]);
 303             }
 304             Matrix_Free(basis);
 305             Polyhedron_Free(C);
 306         }
 307     }
 308     Matrix_Free(all_vertices);
 309
 310     for (i = 0; i < PP->nbV; ++i)
 311         Matrix_Free(vertex_dirs[i]);
 312     free(vertex_dirs);
 313
 314     return width_dirs;
 315 }
 316
 317 /* Computes the lattice width direction of a parametric polytope.
 318  * The parameter space is allowed to be unbounded.
 319  * Currently, the parametric polytope and the parameter space
 320  * are assumed to be full-dimensional.
 321  *
 322  * First, we compute the parametric vertices.
 323  * Then, for each pair of vertices, we construct a (rational) cone
 324  * of directions for which one vertex attains the minimal value
 325  * and the other vertex attians the maximal value.
 326  * The candidate directions are the elements of the integer hulls
 327  * of these cones.
 328  * The minimal direction is then obtained by computing the
 329  * region in the parameter space where each direction yields
 330  * a smaller (or equal) width than all the other directions.
 331  *
 332  * In principle, we can avoid computing candidate directions
 333  * for vertices with no overlapping activity domains (possibly
 334  * after opening some facets of the activity domains in the
 335  * familiar way).
 336  *
 337  * The output is a list of triples, consisting of a direction,
 338  * the corresponding width and the chamber in the parameter
 339  * space where this direction leads to the minimal width.
 340  *
 341  * The algorithm is described in "Integer points in a parameterised
 342  * polyhedron" by Friedrich Eisenbrand and Gennady Shmonin.
 343  */
 344 struct width_direction_array *
 345 Polyhedron_Lattice_Width_Directions(Polyhedron *P, Polyhedron *C,
 346                                     struct barvinok_options *options)
 347 {
 348     Param_Polyhedron *PP;
 349     unsigned nparam = C->Dimension;
 350     int i, j, k;
 351     struct width_direction_array *width_dirs;
 352     Polyhedron *TC;
 353     Vector *inner;
 354
 355     assert(P->NbEq == 0);
 356     assert(C->NbEq == 0);
 357
 358     /* Use true context since the algorithm assumes P is non-empty
 359      * for every point in the context.
 360      */
 361     TC = true_context(P, C, options->MaxRays);
 362     inner = inner_point(TC);
 363
 364     /* This is overkill, as we discard the computed chambers. */
 365     PP = Polyhedron2Param_Polyhedron(P, TC, options);
 366
 367     width_dirs = compute_width_directions(PP, options);
 368     Param_Polyhedron_Free(PP);
 369
 370     qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
 371             wd_width_lex_cmp);
 372
 373     for (i = 1, j = 1; i < width_dirs->n; ++i) {
 374         /* We could also weed out width_directions that differ by a
 375          * (positive) constant from another width_direction, but then
 376          * we'd have to put the two width_directions on a common
 377          * denominator first.
 378          */
 379         if (Vector_Equal(width_dirs->wd[j-1].width->p,
 380                          width_dirs->wd[i].width->p, nparam+2))
 381             clear_width_direction(&width_dirs->wd[i]);
 382         else
 383             width_dirs->wd[j++] = width_dirs->wd[i];
 384     }
 385     width_dirs->n = j;
 386
 387     for (i = 0, k = 0; i < width_dirs->n; ++i) {
 388         Matrix *M = Matrix_Alloc(TC->NbConstraints+width_dirs->n-(i-k)-1, nparam+2);
 389         for (j = 0; j < TC->NbConstraints; ++j)
 390             Vector_Copy(TC->Constraint[j], M->p[j], nparam+2);
 391         for (j = 0; j < width_dirs->n; ++j) {
 392             int pos;
 393             if (k <= j && j <= i)
 394                 continue;
 395             if (j < k)
 396                 pos = TC->NbConstraints + j;
 397             else
 398                 pos = TC->NbConstraints + j - (i-k) - 1;
 399             Vector_Subtract(width_dirs->wd[j].width->p,
 400                             width_dirs->wd[j].width->p[nparam+1],
 401                             width_dirs->wd[i].width->p,
 402                             width_dirs->wd[i].width->p[nparam+1],
 403                             M->p[pos]+1, M->p[pos], nparam+1);
 404             value_set_si(M->p[pos][0], 1);
 405             Vector_Normalize(M->p[pos]+1, nparam+1);
 406             if (!is_internal(inner, M->p[pos]))
 407                 value_decrement(M->p[pos][nparam+1], M->p[pos][nparam+1]);
 408         }
 409         width_dirs->wd[i].domain = Constraints2Polyhedron(M, options->MaxRays);
 410         if (emptyQ(width_dirs->wd[i].domain))
 411             clear_width_direction(&width_dirs->wd[i]);
 412         else
 413             width_dirs->wd[k++] = width_dirs->wd[i];
 414         Matrix_Free(M);
 415     }
 416     width_dirs->n = k;
 417     Vector_Free(inner);
 418     Polyhedron_Free(TC);
 419
 420     return width_dirs;
 421 }
 422
 423 /* Construct evalue of chambers with their associated widths */
 424 evalue *Polyhedron_Lattice_Width(Polyhedron *P, Polyhedron *C,
 425                                  struct barvinok_options *options)
 426 {
 427     evalue *width;
 428     struct evalue_section *s;
 429     struct width_direction_array *width_dirs;
 430     int i;
 431     unsigned nparam = C->Dimension;
 432
 433     width_dirs = Polyhedron_Lattice_Width_Directions(P, C, options);
 434     s = ALLOCN(struct evalue_section, width_dirs->n);
 435     for (i = 0; i < width_dirs->n; ++i) {
 436         s[i].D = width_dirs->wd[i].domain;
 437         width_dirs->wd[i].domain = NULL;
 438         s[i].E = affine2evalue(width_dirs->wd[i].width->p,
 439                                width_dirs->wd[i].width->p[nparam+1],
 440                                nparam);
 441     }
 442     free_width_direction_array(width_dirs);
 443
 444     width = evalue_from_section_array(s, i);
 445     free(s);
 446
 447     return width;
 448 }
 449
 450 __isl_give isl_pw_qpolynomial *isl_basic_set_lattice_width(
 451         __isl_take isl_basic_set *bset)
 452 {
 453         isl_ctx *ctx;
 454         isl_space *dim;
 455         isl_pw_qpolynomial *pwqp;
 456         unsigned nparam;
 457         Polyhedron *U;
 458         Polyhedron *P;
 459         evalue *E;
 460         struct barvinok_options *options;
 461         int options_allocated = 0;
 462
 463         if (!bset)
 464                 return NULL;
 465
 466         ctx = isl_basic_set_get_ctx(bset);
 467         options = isl_ctx_peek_barvinok_options(ctx);
 468         if (!options) {
 469                 options = barvinok_options_new_with_defaults();
 470                 options_allocated = 1;
 471         }
 472
 473         nparam = isl_basic_set_dim(bset, isl_dim_param);
 474         dim = isl_basic_set_get_space(bset);
 475         dim = isl_space_params(dim);
 476
 477         U = Universe_Polyhedron(nparam);
 478         P = isl_basic_set_to_polylib(bset);
 479
 480         E = Polyhedron_Lattice_Width(P, U, options);
 481
 482         pwqp = isl_pw_qpolynomial_from_evalue(dim, E);
 483         isl_basic_set_free(bset);
 484
 485         evalue_free(E);
 486         Polyhedron_Free(P);
 487         Polyhedron_Free(U);
 488         if (options_allocated)
 489                 barvinok_options_free(options);
 490
 491         return pwqp;
 492 }
 493
 494 __isl_give isl_pw_qpolynomial *isl_set_lattice_width(__isl_take isl_set *set)
 495 {
 496         if (!set)
 497                 return NULL;
 498
 499         if (isl_set_plain_is_empty(set)) {
 500                 isl_space *dim;
 501                 dim = isl_set_get_space(set);
 502                 dim = isl_space_domain(isl_space_from_range(dim));
 503                 isl_set_free(set);
 504                 return isl_pw_qpolynomial_zero(dim);
 505         }
 506
 507         if (isl_set_n_basic_set(set) != 1)
 508                 isl_die(isl_set_get_ctx(set), isl_error_unsupported,
 509                         "unions not supported (yet)", goto error);
 510
 511         return isl_basic_set_lattice_width(isl_set_simple_hull(set));
 512 error:
 513         isl_set_free(set);
 514         return NULL;
 515 }
 516
 517 static int set_lw(__isl_take isl_set *set, void *user)
 518 {
 519         isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
 520         isl_pw_qpolynomial *pwqp;
 521
 522         pwqp = isl_set_lattice_width(set);
 523         *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, pwqp);
 524
 525         return 0;
 526 }
 527
 528 __isl_give isl_union_pw_qpolynomial *isl_union_set_lattice_width(
 529         __isl_take isl_union_set *uset)
 530 {
 531         isl_space *dim;
 532         isl_union_pw_qpolynomial *res;
 533
 534         dim = isl_union_set_get_space(uset);
 535         res = isl_union_pw_qpolynomial_zero(dim);
 536         if (isl_union_set_foreach_set(uset, &set_lw, &res) < 0)
 537                 goto error;
 538         isl_union_set_free(uset);
 539
 540         return res;
 541 error:
 542         isl_union_set_free(uset);
 543         isl_union_pw_qpolynomial_free(res);
 544         return NULL;
 545 }