sample.c

   1 #include <polylib/polylibgmp.h>
   2 #include <barvinok/util.h>
   3 #include "basis_reduction.h"
   4 #include "sample.h"
   5
   6 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   7
   8 /* If P has no rays, then we return NULL.
   9  * Otherwise, transform the polyhedron such that one of the rays
  10  * is the first unit vector and cut it off at a height that ensures
  11  * that if the whole polyhedron has any points, then the remaining part
  12  * has integer points.  In particular we add the largest coefficient
  13  * of a ray to the highest vertex (rounded up).
  14  *
  15  * The matrix that transforms the resulting polytope to part of the
  16  * original polyhedron is returned through T.
  17  */
  18 static Polyhedron *remove_ray(Polyhedron *P, Matrix **T, unsigned MaxRays)
  19 {
  20     int r = 0;
  21     Matrix *M, *M2;
  22     Value c, tmp;
  23     Value g;
  24     int first;
  25     Vector *v;
  26     Value offset, size;
  27     Polyhedron *R;
  28     int i;
  29
  30     if (P->NbBid == 0)
  31         for (; r < P->NbRays; ++r)
  32             if (value_zero_p(P->Ray[r][P->Dimension+1]))
  33                 break;
  34     if (P->NbBid == 0 && r == P->NbRays)
  35         return NULL;
  36
  37     value_init(g);
  38     v = Vector_Alloc(P->Dimension+1);
  39     Vector_Gcd(P->Ray[r]+1, P->Dimension, &g);
  40     Vector_AntiScale(P->Ray[r]+1, v->p, g, P->Dimension+1);
  41     M = unimodular_complete(v);
  42     value_set_si(M->p[P->Dimension][P->Dimension], 1);
  43     M2 = Transpose(M);
  44     Matrix_Free(M);
  45     P = Polyhedron_Preimage(P, M2, 0);
  46     *T = M2;
  47     value_clear(g);
  48     Vector_Free(v);
  49
  50     first = 1;
  51     value_init(offset);
  52     value_init(size);
  53     value_init(tmp);
  54     value_set_si(size, 0);
  55
  56     for (i = 0; i < P->NbBid; ++i) {
  57         value_absolute(tmp, P->Ray[i][1]);
  58         if (value_gt(tmp, size))
  59             value_assign(size, tmp);
  60     }
  61     for (i = P->NbBid; i < P->NbRays; ++i) {
  62         if (value_zero_p(P->Ray[i][P->Dimension+1])) {
  63             if (value_gt(P->Ray[i][1], size))
  64                 value_assign(size, P->Ray[i][1]);
  65             continue;
  66         }
  67         mpz_cdiv_q(tmp, P->Ray[i][1], P->Ray[i][P->Dimension+1]);
  68         if (first || value_gt(tmp, offset)) {
  69             value_assign(offset, tmp);
  70             first = 0;
  71         }
  72     }
  73     value_addto(offset, offset, size);
  74     value_clear(size);
  75     value_clear(tmp);
  76
  77     v = Vector_Alloc(P->Dimension+2);
  78     value_set_si(v->p[0], 1);
  79     value_set_si(v->p[1], -1);
  80     value_assign(v->p[1+P->Dimension], offset);
  81     R = AddConstraints(v->p, 1, P, MaxRays);
  82     value_clear(offset);
  83     Vector_Free(v);
  84     Polyhedron_Free(P);
  85
  86     return R;
  87 }
  88
  89 static void print_minmax(Polyhedron *P)
  90 {
  91     int i, j;
  92     POL_ENSURE_VERTICES(P);
  93     Polyhedron_Print(stderr, P_VALUE_FMT, P);
  94     for (i = 0; i < P->Dimension; ++i) {
  95         Value min, max, tmp;
  96         value_init(min);
  97         value_init(max);
  98         value_init(tmp);
  99
 100         mpz_cdiv_q(min, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
 101         mpz_fdiv_q(max, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
 102
 103         for (j = 1; j < P->NbRays; ++j) {
 104             mpz_cdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
 105             if (value_lt(tmp, min))
 106                 value_assign(min, tmp);
 107             mpz_fdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
 108             if (value_gt(tmp, max))
 109                 value_assign(max, tmp);
 110         }
 111         fprintf(stderr, "i: %d,  min: ", i);
 112         value_print(stderr, VALUE_FMT, min);
 113         fprintf(stderr, ", max: ");
 114         value_print(stderr, VALUE_FMT, max);
 115         fprintf(stderr, "\n");
 116
 117         value_clear(min);
 118         value_clear(max);
 119         value_clear(tmp);
 120     }
 121 }
 122
 123 /* Remove coordinates that have a fixed value and return the matrix
 124  * that adds these fixed coordinates again through T.
 125  */
 126 static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T)
 127 {
 128     int i, j, n;
 129     int dim = P->Dimension;
 130     int *remove = ALLOCN(int, dim);
 131     Polyhedron *Q;
 132
 133     assert(POL_HAS(P, POL_INEQUALITIES));
 134     for (i = 0; i < dim; ++i)
 135         remove[i] = 0;
 136     for (i = 0; i < P->NbEq; ++i) {
 137         int pos = First_Non_Zero(P->Constraint[i]+1, dim);
 138         assert(First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) == -1);
 139         remove[pos] = 1;
 140     }
 141     Q = Polyhedron_Alloc(P->Dimension-P->NbEq, P->NbConstraints-P->NbEq, P->NbRays);
 142     for (i = 0; i < Q->NbConstraints; ++i) {
 143         value_assign(Q->Constraint[i][0], P->Constraint[P->NbEq+i][0]);
 144         for (j = 0, n = 0; j < P->Dimension; ++j) {
 145             if (remove[j])
 146                 ++n;
 147             else
 148                 value_assign(Q->Constraint[i][1+j-n], P->Constraint[P->NbEq+i][1+j]);
 149         }
 150         value_assign(Q->Constraint[i][1+j-n], P->Constraint[P->NbEq+i][1+j]);
 151     }
 152     for (i = 0; i < Q->NbRays; ++i) {
 153         value_assign(Q->Ray[i][0], P->Ray[i][0]);
 154         for (j = 0, n = 0; j < P->Dimension; ++j) {
 155             if (remove[j])
 156                 ++n;
 157             else
 158                 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
 159         }
 160         value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
 161     }
 162     *T = Matrix_Alloc(P->Dimension+1, Q->Dimension+1);
 163     for (i = 0, n = 0; i < P->Dimension; ++i) {
 164         if (remove[i]) {
 165             value_oppose((*T)->p[i][Q->Dimension], P->Constraint[n][1+P->Dimension]);
 166             ++n;
 167         } else
 168             value_set_si((*T)->p[i][i-n], 1);
 169     }
 170     value_set_si((*T)->p[i][i-n], 1);
 171     POL_SET(Q, POL_VALID);
 172     if (POL_HAS(P, POL_INEQUALITIES))
 173         POL_SET(Q, POL_INEQUALITIES);
 174     if (POL_HAS(P, POL_FACETS))
 175         POL_SET(Q, POL_FACETS);
 176     if (POL_HAS(P, POL_POINTS))
 177         POL_SET(Q, POL_POINTS);
 178     if (POL_HAS(P, POL_VERTICES))
 179         POL_SET(Q, POL_VERTICES);
 180     free(remove);
 181     return Q;
 182 }
 183
 184 /* This function implements the algorithm described in
 185  * "An Implementation of the Generalized Basis Reduction Algorithm
 186  *  for Integer Programming" of Cook el al. to find an integer point
 187  * in a polyhedron.
 188  * If the polyhedron is unbounded, we first remove its rays.
 189  */
 190 Vector *Polyhedron_Sample(Polyhedron *P, unsigned MaxRays)
 191 {
 192     int i, j;
 193     Vector *sample = NULL;
 194     Polyhedron *Q;
 195     Matrix *T, *inv, *M;
 196     Value min, max, tmp;
 197     Vector *v;
 198     int ok;
 199
 200     POL_ENSURE_VERTICES(P);
 201
 202     for (i = 0; i < P->NbRays; ++i)
 203         if (value_one_p(P->Ray[i][1+P->Dimension])) {
 204             sample = Vector_Alloc(P->Dimension+1);
 205             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
 206             return sample;
 207         }
 208
 209     /* for now */
 210     assert(P->NbEq == 0);
 211
 212     Q = remove_ray(P, &T, MaxRays);
 213     if (Q) {
 214         Vector *Q_sample;
 215
 216         Q_sample = Polyhedron_Sample(Q, MaxRays);
 217         Polyhedron_Free(Q);
 218
 219         if (Q_sample) {
 220             sample = Vector_Alloc(P->Dimension + 1);
 221             Matrix_Vector_Product(T, Q_sample->p, sample->p);
 222             Vector_Free(Q_sample);
 223         }
 224
 225         Matrix_Free(T);
 226         return sample;
 227     }
 228
 229     Matrix *basis = reduced_basis(P);
 230
 231     T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 232     inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 233     for (i = 0; i < P->Dimension; ++i)
 234         for (j = 0; j < P->Dimension; ++j)
 235             value_assign(T->p[i][j], basis->p[i][j]);
 236     value_set_si(T->p[P->Dimension][P->Dimension], 1);
 237     Matrix_Free(basis);
 238
 239     M = Matrix_Copy(T);
 240     ok = Matrix_Inverse(M, inv);
 241     assert(ok);
 242     Matrix_Free(M);
 243
 244     Q = Polyhedron_Image(P, T, MaxRays);
 245
 246     value_init(min);
 247     value_init(max);
 248     value_init(tmp);
 249
 250     mpz_cdiv_q(min, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
 251     mpz_fdiv_q(max, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
 252
 253     for (j = 1; j < Q->NbRays; ++j) {
 254         mpz_cdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
 255         if (value_lt(tmp, min))
 256             value_assign(min, tmp);
 257         mpz_fdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
 258         if (value_gt(tmp, max))
 259             value_assign(max, tmp);
 260     }
 261
 262     v = Vector_Alloc(1+Q->Dimension+1);
 263     value_set_si(v->p[1], -1);
 264
 265     for (value_assign(tmp, min); value_le(tmp, max); value_increment(tmp, tmp)) {
 266         Polyhedron *R, *S;
 267         Matrix *T;
 268         Vector *S_sample;
 269         value_assign(v->p[1+Q->Dimension], tmp);
 270
 271         R = AddConstraints(v->p, 1, Q, MaxRays);
 272         S = Polyhedron_RemoveFixedColumns(R, &T);
 273         Polyhedron_Free(R);
 274         S_sample = Polyhedron_Sample(S, MaxRays);
 275         Polyhedron_Free(S);
 276         if (S_sample) {
 277             Vector *Q_sample = Vector_Alloc(Q->Dimension + 1);
 278             Matrix_Vector_Product(T, S_sample->p, Q_sample->p);
 279             Matrix_Free(T);
 280             Vector_Free(S_sample);
 281             sample = Vector_Alloc(P->Dimension + 1);
 282             Matrix_Vector_Product(inv, Q_sample->p, sample->p);
 283             Vector_Free(Q_sample);
 284             break;
 285         }
 286         Matrix_Free(T);
 287     }
 288
 289     Matrix_Free(T);
 290     Matrix_Free(inv);
 291     Polyhedron_Free(Q);
 292     Vector_Free(v);
 293
 294     value_clear(min);
 295     value_clear(max);
 296     value_clear(tmp);
 297
 298     return sample;
 299 }