lattice_width.c

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