lattice_width.c

   1 #include <stdlib.h>
   2 #include <barvinok/options.h>
   3 #include <barvinok/util.h>
   4 #include "hilbert.h"
   5 #include "lattice_width.h"
   6 #include "param_util.h"
   7 #include "reduce_domain.h"
   8
   9 #define ALLOC(type) (type*)malloc(sizeof(type))
  10 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  11
  12 static void clear_width_direction(struct width_direction *wd)
  13 {
  14     Vector_Free(wd->width);
  15     Vector_Free(wd->dir);
  16     if (wd->domain)
  17         Polyhedron_Free(wd->domain);
  18 }
  19
  20 static struct width_direction_array *new_width_direction_array(void)
  21 {
  22     struct width_direction_array *dirs = ALLOC(struct width_direction_array);
  23     assert(dirs);
  24     dirs->n = 0;
  25     dirs->alloc = 32;
  26     dirs->wd = ALLOCN(struct width_direction, dirs->alloc);
  27     assert(dirs->wd);
  28     return dirs;
  29 }
  30
  31 static void grow_width_direction_array(struct width_direction_array *dirs,
  32                                        int extra)
  33 {
  34     if (dirs->n + extra <= dirs->alloc)
  35         return;
  36     dirs->alloc = (5*(dirs->n+extra))/4;
  37     dirs->wd = (struct width_direction*)realloc(dirs->wd,
  38                             dirs->alloc * sizeof(struct width_direction));
  39     assert(dirs->wd);
  40 }
  41
  42 void free_width_direction_array(struct width_direction_array *dirs)
  43 {
  44     int i;
  45
  46     for (i = 0; i < dirs->n; ++i)
  47         clear_width_direction(&dirs->wd[i]);
  48     free(dirs->wd);
  49     free(dirs);
  50 }
  51
  52 #define INT_BITS (sizeof(unsigned) * 8)
  53
  54 /* For each parametric vertex, compute cone of directions
  55  * for which this vertex attains the minimal value.
  56  */
  57 static Matrix **compute_vertex_dirs(Param_Polyhedron *PP)
  58 {
  59     int i;
  60     unsigned nvar = PP->V->Vertex->NbRows;
  61     Param_Vertices *V;
  62     Matrix **vertex_dirs = ALLOCN(Matrix *, PP->nbV);
  63
  64     for (i = 0, V = PP->V; V; ++i, V = V->next) {
  65         int kx;
  66         unsigned bx;
  67         int j, k;
  68         unsigned *facets;
  69         int n;
  70         Matrix *M;
  71         Polyhedron *P;
  72
  73         if (V->Facets) {
  74             int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
  75             facets = V->Facets;
  76             n = bit_vector_count(facets, len);
  77         } else
  78             facets = supporting_constraints(PP->Constraints, V, &n);
  79         M = Matrix_Alloc(n, 1+nvar+1);
  80         for (k = 0, j = 0, kx = 0, bx = MSB; j < n; ++k) {
  81             if (facets[kx] & bx) {
  82                 value_set_si(M->p[j][0], 1);
  83                 Vector_Copy(PP->Constraints->p[k]+1, M->p[j++]+1, nvar);
  84             }
  85             NEXT(kx, bx);
  86         }
  87         P = Constraints2Polyhedron(M, 0);
  88         Matrix_Free(M);
  89         vertex_dirs[i] = Matrix_Alloc(P->NbRays-1, nvar);
  90         for (k = 0, j = 0; k < P->NbRays; ++k) {
  91             if (value_notzero_p(P->Ray[k][1+nvar]))
  92                 continue;
  93             Vector_Copy(P->Ray[k]+1, vertex_dirs[i]->p[j++], nvar);
  94         }
  95         Polyhedron_Free(P);
  96         if (!V->Facets)
  97             free(facets);
  98     }
  99     return vertex_dirs;
 100 }
 101
 102 /* Compute
 103  *
 104  *       c     a     b
 105  *      --- = --- - ---
 106  *      c_d   a_d   b_d
 107  */
 108 static void Vector_Subtract(Value *a, Value a_d,
 109                             Value *b, Value b_d,
 110                             Value *c, Value *c_d, int len)
 111 {
 112     Value ma, mb;
 113     value_init(ma);
 114     value_init(mb);
 115     value_lcm(*c_d, a_d, b_d);
 116     value_divexact(ma, *c_d, a_d);
 117     value_divexact(mb, *c_d, b_d);
 118     value_oppose(mb, mb);
 119     Vector_Combine(a, b, c, ma, mb, len);
 120     value_clear(ma);
 121     value_clear(mb);
 122 }
 123
 124 /* Compute width for a given direction dir and initialize width_direction
 125  * structure.
 126  */
 127 static void compute_width_direction(Matrix *V_min, Matrix *V_max,
 128                                     Value *dir, struct width_direction *wd)
 129 {
 130     Vector *max = Vector_Alloc(V_min->NbColumns);
 131     unsigned nvar = V_min->NbRows;
 132     unsigned nparam = V_min->NbColumns-2;
 133
 134     wd->width = Vector_Alloc(V_min->NbColumns);
 135     wd->dir = Vector_Alloc(nvar);
 136     Vector_Copy(dir, wd->dir->p, nvar);
 137     wd->domain = NULL;
 138
 139     V_min->NbColumns--;
 140     V_max->NbColumns--;
 141
 142     Vector_Matrix_Product(dir, V_max, max->p);
 143     Vector_Matrix_Product(dir, V_min, wd->width->p);
 144     Vector_Subtract(max->p, V_max->p[0][V_max->NbColumns],
 145                     wd->width->p, V_min->p[0][V_min->NbColumns],
 146                     wd->width->p, &wd->width->p[nparam+1],
 147                     nparam+1);
 148
 149     V_min->NbColumns++;
 150     V_max->NbColumns++;
 151
 152     Vector_Normalize(wd->width->p, nparam+2);
 153
 154     Vector_Free(max);
 155 }
 156
 157 static int Vector_Compare(Value *p1, Value *p2, unsigned len)
 158 {
 159     int i;
 160
 161     for (i = 0; i < len; ++i) {
 162         int sign = mpz_cmp(p1[i], p2[i]);
 163         if (sign)
 164             return sign;
 165     }
 166     return 0;
 167 }
 168
 169 static int wd_width_lex_cmp(const void *va, const void *vb)
 170 {
 171     const struct width_direction *a = (const struct width_direction *)va;
 172     const struct width_direction *b = (const struct width_direction *)vb;
 173
 174     return Vector_Compare(a->width->p, b->width->p, a->width->Size);
 175 }
 176
 177 static int wd_dir_lex_cmp(const void *va, const void *vb)
 178 {
 179     const struct width_direction *a = (const struct width_direction *)va;
 180     const struct width_direction *b = (const struct width_direction *)vb;
 181
 182     return Vector_Compare(a->dir->p, b->dir->p, a->dir->Size);
 183 }
 184
 185 static struct width_direction_array *
 186 compute_width_directions(Param_Polyhedron *PP, struct barvinok_options *options)
 187 {
 188     Matrix **vertex_dirs;
 189     Param_Vertices *V_max, *V_min;
 190     int i, V_max_i, V_min_i;
 191     unsigned nvar = PP->V->Vertex->NbRows;
 192     struct width_direction_array *width_dirs = new_width_direction_array();
 193
 194     vertex_dirs = compute_vertex_dirs(PP);
 195
 196     for (V_max = PP->V; V_max; V_max = V_max->next)
 197         Param_Vertex_Common_Denominator(V_max);
 198
 199     for (V_max = PP->V, V_max_i = 0; V_max; V_max = V_max->next, V_max_i++) {
 200         for (V_min = V_max->next, V_min_i = V_max_i+1;
 201                 V_min;
 202                 V_min = V_min->next, V_min_i++) {
 203             Matrix *M;
 204             Matrix *basis;
 205             Polyhedron *C;
 206             unsigned V_max_n = vertex_dirs[V_max_i]->NbRows;
 207             unsigned V_min_n = vertex_dirs[V_min_i]->NbRows;
 208             int sorted_n;
 209
 210             if (options->verbose)
 211                 fprintf(stderr, "%d/%d %d/%d %d \r",
 212                                     V_max_i, PP->nbV,
 213                                     V_min_i, PP->nbV,
 214                                     width_dirs->n);
 215
 216             M = Matrix_Alloc(V_max_n+V_min_n, 1+nvar+1);
 217             for (i = 0; i < V_max_n; ++i) {
 218                 value_set_si(M->p[i][0], 1);
 219                 Vector_Oppose(vertex_dirs[V_max_i]->p[i], M->p[i]+1, nvar);
 220             }
 221             for (i = 0; i < V_min_n; ++i) {
 222                 value_set_si(M->p[V_max_n+i][0], 1);
 223                 Vector_Copy(vertex_dirs[V_min_i]->p[i], M->p[V_max_n+i]+1, nvar);
 224             }
 225             C = Constraints2Polyhedron(M, options->MaxRays);
 226             Matrix_Free(M);
 227             basis = Cone_Integer_Hull(C, options);
 228             grow_width_direction_array(width_dirs, basis->NbRows);
 229             qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
 230                     wd_dir_lex_cmp);
 231             sorted_n = width_dirs->n;
 232             for (i = 0; i < basis->NbRows; ++i) {
 233                 Vector v;
 234                 struct width_direction wd;
 235
 236                 v.Size = nvar;
 237                 v.p = basis->p[i];
 238                 wd.dir = &v;
 239                 if (bsearch(&wd, width_dirs->wd, sorted_n,
 240                             sizeof(struct width_direction),
 241                             wd_dir_lex_cmp))
 242                     continue;
 243
 244                 compute_width_direction(V_min->Vertex, V_max->Vertex,
 245                                         basis->p[i],
 246                                         &width_dirs->wd[width_dirs->n++]);
 247             }
 248             Matrix_Free(basis);
 249             Polyhedron_Free(C);
 250         }
 251     }
 252
 253     for (i = 0; i < PP->nbV; ++i)
 254         Matrix_Free(vertex_dirs[i]);
 255     free(vertex_dirs);
 256
 257     return width_dirs;
 258 }
 259
 260 /* Computes the lattice width direction of a parametric polytope.
 261  * The parameter space is allowed to be unbounded.
 262  * Currently, the parametric polytope and the parameter space
 263  * are assumed to be full-dimensional.
 264  *
 265  * First, we compute the parametric vertices.
 266  * Then, for each pair of vertices, we construct a (rational) cone
 267  * of directions for which one vertex attains the minimal value
 268  * and the other vertex attians the maximal value.
 269  * The candidate directions are the elements of the integer hulls
 270  * of these cones.
 271  * The minimal direction is then obtained by computing the
 272  * region in the parameter space where each direction yields
 273  * a smaller (or equal) width than all the other directions.
 274  *
 275  * In principle, we can avoid computing candidate directions
 276  * for vertices with no overlapping activity domains (possibly
 277  * after opening some facets of the activity domains in the
 278  * familiar way).
 279  *
 280  * The output is a list of triples, consisting of a direction,
 281  * the corresponding width and the chamber in the parameter
 282  * space where this direction leads to the minimal width.
 283  *
 284  * The algorithm is described in "Integer points in a parameterised
 285  * polyhedron" by Friedrich Eisenbrand and Gennady Shmonin.
 286  */
 287 struct width_direction_array *
 288 Polyhedron_Lattice_Width_Directions(Polyhedron *P, Polyhedron *C,
 289                                     struct barvinok_options *options)
 290 {
 291     Param_Polyhedron *PP;
 292     unsigned nparam = C->Dimension;
 293     int i, j, k;
 294     struct width_direction_array *width_dirs;
 295     Polyhedron *TC;
 296     Vector *inner;
 297
 298     assert(P->NbEq == 0);
 299     assert(C->NbEq == 0);
 300
 301     /* Use true context since the algorithm assumes P is non-empty
 302      * for every point in the context.
 303      */
 304     TC = true_context(P, C, options->MaxRays);
 305     inner = inner_point(TC);
 306
 307     /* This is overkill, as we discard the computed chambers. */
 308     PP = Polyhedron2Param_Polyhedron(P, TC, options);
 309
 310     width_dirs = compute_width_directions(PP, options);
 311     Param_Polyhedron_Free(PP);
 312
 313     qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
 314             wd_width_lex_cmp);
 315
 316     for (i = 1, j = 1; i < width_dirs->n; ++i) {
 317         /* We could also weed out width_directions that differ by a
 318          * (positive) constant from another width_direction, but then
 319          * we'd have to put the two width_directions on a common
 320          * denominator first.
 321          */
 322         if (Vector_Equal(width_dirs->wd[j-1].width->p,
 323                          width_dirs->wd[i].width->p, nparam+2))
 324             clear_width_direction(&width_dirs->wd[i]);
 325         else
 326             width_dirs->wd[j++] = width_dirs->wd[i];
 327     }
 328     width_dirs->n = j;
 329
 330     for (i = 0, k = 0; i < width_dirs->n; ++i) {
 331         Matrix *M = Matrix_Alloc(TC->NbConstraints+width_dirs->n-(i-k)-1, nparam+2);
 332         for (j = 0; j < TC->NbConstraints; ++j)
 333             Vector_Copy(TC->Constraint[j], M->p[j], nparam+2);
 334         for (j = 0; j < width_dirs->n; ++j) {
 335             int pos;
 336             if (k <= j && j <= i)
 337                 continue;
 338             if (j < k)
 339                 pos = TC->NbConstraints + j;
 340             else
 341                 pos = TC->NbConstraints + j - (i-k) - 1;
 342             Vector_Subtract(width_dirs->wd[j].width->p,
 343                             width_dirs->wd[j].width->p[nparam+1],
 344                             width_dirs->wd[i].width->p,
 345                             width_dirs->wd[i].width->p[nparam+1],
 346                             M->p[pos]+1, M->p[pos], nparam+1);
 347             value_set_si(M->p[pos][0], 1);
 348             Vector_Normalize(M->p[pos]+1, nparam+1);
 349             if (!is_internal(inner, M->p[pos]))
 350                 value_decrement(M->p[pos][nparam+1], M->p[pos][nparam+1]);
 351         }
 352         width_dirs->wd[i].domain = Constraints2Polyhedron(M, options->MaxRays);
 353         if (emptyQ(width_dirs->wd[i].domain))
 354             clear_width_direction(&width_dirs->wd[i]);
 355         else
 356             width_dirs->wd[k++] = width_dirs->wd[i];
 357         Matrix_Free(M);
 358     }
 359     width_dirs->n = k;
 360     Vector_Free(inner);
 361     Polyhedron_Free(TC);
 362
 363     return width_dirs;
 364 }
 365
 366 /* Construct evalue of chambers with their associated widths */
 367 evalue *Polyhedron_Lattice_Width(Polyhedron *P, Polyhedron *C,
 368                                  struct barvinok_options *options)
 369 {
 370     evalue *width;
 371     struct evalue_section *s;
 372     struct width_direction_array *width_dirs;
 373     int i;
 374     unsigned nparam = C->Dimension;
 375
 376     width_dirs = Polyhedron_Lattice_Width_Directions(P, C, options);
 377     s = ALLOCN(struct evalue_section, width_dirs->n);
 378     for (i = 0; i < width_dirs->n; ++i) {
 379         s[i].D = width_dirs->wd[i].domain;
 380         width_dirs->wd[i].domain = NULL;
 381         s[i].E = affine2evalue(width_dirs->wd[i].width->p,
 382                                width_dirs->wd[i].width->p[nparam+1],
 383                                nparam);
 384     }
 385     free_width_direction_array(width_dirs);
 386
 387     width = evalue_from_section_array(s, i);
 388     free(s);
 389
 390     return width;
 391 }