lattice_width.c

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