lattice_width.c

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