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