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, look for the coordinate axis with the smallest maximal non-zero
  10  * coefficient over all rays and a constraint that bounds the values on
  11  * this axis to the maximal value over the vertices plus the above maximal
  12  * non-zero coefficient minus 1.
  13  * Any integer point outside this region should be the sum of a point inside
  14  * the region and an integer multiple of the rays.
  15  */
  16 static Polyhedron *remove_ray(Polyhedron *P, unsigned MaxRays)
  17 {
  18     int r = 0;
  19     Vector *min, *max, *c;
  20     int i;
  21     Value s, v, tmp;
  22     int pos;
  23     Polyhedron *R;
  24
  25     if (P->NbBid == 0)
  26         for (; r < P->NbRays; ++r)
  27             if (value_zero_p(P->Ray[r][P->Dimension+1]))
  28                 break;
  29     if (P->NbBid == 0 && r == P->NbRays)
  30         return NULL;
  31
  32     max = Vector_Alloc(P->Dimension);
  33     min = Vector_Alloc(P->Dimension);
  34     for (r = 0; r < P->NbBid; ++r)
  35         for (i = 0 ; i < P->Dimension; ++i)
  36             if (value_abs_gt(P->Ray[r][1+i], max->p[i]))
  37                 value_absolute(max->p[i], P->Ray[r][1+i]);
  38
  39     for (i = 0 ; i < P->Dimension; ++i)
  40         value_oppose(min->p[i], max->p[i]);
  41
  42     for (r = P->NbBid; r < P->NbRays; ++r) {
  43         if (value_notzero_p(P->Ray[r][P->Dimension+1]))
  44             continue;
  45         for (i = 0 ; i < P->Dimension; ++i) {
  46             if (value_gt(P->Ray[r][1+i], max->p[i]))
  47                 value_assign(max->p[i], P->Ray[r][1+i]);
  48             if (value_lt(P->Ray[r][1+i], min->p[i]))
  49                 value_assign(min->p[i], P->Ray[r][1+i]);
  50         }
  51     }
  52
  53     value_init(s);
  54     value_init(v);
  55     value_init(tmp);
  56
  57     for (i = 0 ; i < P->Dimension; ++i) {
  58         if (value_notzero_p(min->p[i]) &&
  59             (value_zero_p(s) || value_abs_lt(min->p[i], s))) {
  60             value_assign(s, min->p[i]);
  61             pos = i;
  62         }
  63         if (value_notzero_p(max->p[i]) &&
  64             (value_zero_p(s) || value_abs_lt(max->p[i], s))) {
  65             value_assign(s, max->p[i]);
  66             pos = i;
  67         }
  68     }
  69
  70     for (r = P->NbBid; r < P->NbRays; ++r)
  71         if (value_notzero_p(P->Ray[r][P->Dimension+1]))
  72             break;
  73
  74     if (value_pos_p(s))
  75         mpz_cdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  76     else
  77         mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  78
  79     for ( ; r < P->NbRays; ++r) {
  80         if (value_zero_p(P->Ray[r][P->Dimension+1]))
  81             continue;
  82
  83         if (value_pos_p(s)) {
  84             mpz_cdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  85             if (value_gt(tmp, v))
  86                 value_assign(v, tmp);
  87         } else {
  88             mpz_fdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  89             if (value_lt(tmp, v))
  90                 value_assign(v, tmp);
  91         }
  92     }
  93
  94     c = Vector_Alloc(1+P->Dimension+1);
  95
  96     value_addto(v, v, s);
  97     value_set_si(c->p[0], 1);
  98     if (value_pos_p(s)) {
  99         value_set_si(c->p[1+pos], -1);
 100         value_assign(c->p[1+P->Dimension], v);
 101     } else {
 102         value_set_si(c->p[1+pos], 1);
 103         value_oppose(c->p[1+P->Dimension], v);
 104     }
 105     value_decrement(c->p[1+P->Dimension], c->p[1+P->Dimension]);
 106
 107     R = AddConstraints(c->p, 1, P, MaxRays);
 108
 109     Vector_Free(c);
 110
 111     Vector_Free(min);
 112     Vector_Free(max);
 113
 114     value_clear(tmp);
 115     value_clear(s);
 116     value_clear(v);
 117
 118     return R;
 119 }
 120
 121 static void print_minmax(Polyhedron *P)
 122 {
 123     int i, j;
 124     POL_ENSURE_VERTICES(P);
 125     Polyhedron_Print(stderr, P_VALUE_FMT, P);
 126     for (i = 0; i < P->Dimension; ++i) {
 127         Value min, max, tmp;
 128         value_init(min);
 129         value_init(max);
 130         value_init(tmp);
 131
 132         mpz_cdiv_q(min, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
 133         mpz_fdiv_q(max, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
 134
 135         for (j = 1; j < P->NbRays; ++j) {
 136             mpz_cdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
 137             if (value_lt(tmp, min))
 138                 value_assign(min, tmp);
 139             mpz_fdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
 140             if (value_gt(tmp, max))
 141                 value_assign(max, tmp);
 142         }
 143         fprintf(stderr, "i: %d,  min: ", i);
 144         value_print(stderr, VALUE_FMT, min);
 145         fprintf(stderr, ", max: ");
 146         value_print(stderr, VALUE_FMT, max);
 147         fprintf(stderr, "\n");
 148
 149         value_clear(min);
 150         value_clear(max);
 151         value_clear(tmp);
 152     }
 153 }
 154
 155 /* Remove coordinates that have a fixed value and return the matrix
 156  * that adds these fixed coordinates again through T.
 157  */
 158 static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T)
 159 {
 160     int i, j, k, n;
 161     int dim = P->Dimension;
 162     int *remove = ALLOCN(int, dim);
 163     Polyhedron *Q;
 164     int NbEq;
 165
 166     assert(POL_HAS(P, POL_INEQUALITIES));
 167     for (i = 0; i < dim; ++i)
 168         remove[i] = 0;
 169     NbEq = 0;
 170     for (i = 0; i < P->NbEq; ++i) {
 171         int pos = First_Non_Zero(P->Constraint[i]+1, dim);
 172         if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) != -1)
 173             continue;
 174         remove[pos] = 1;
 175         ++NbEq;
 176     }
 177     assert(NbEq > 0);
 178     Q = Polyhedron_Alloc(P->Dimension-NbEq, P->NbConstraints-NbEq, P->NbRays);
 179     for (i = 0, k = 0; i < P->NbConstraints; ++i) {
 180         if (i < P->NbEq) {
 181             int pos = First_Non_Zero(P->Constraint[i]+1, dim);
 182             if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) == -1)
 183                 continue;
 184         }
 185         value_assign(Q->Constraint[k][0], P->Constraint[i][0]);
 186         for (j = 0, n = 0; j < P->Dimension; ++j) {
 187             if (remove[j])
 188                 ++n;
 189             else
 190                 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
 191         }
 192         value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
 193         ++k;
 194     }
 195     for (i = 0; i < Q->NbRays; ++i) {
 196         value_assign(Q->Ray[i][0], P->Ray[i][0]);
 197         for (j = 0, n = 0; j < P->Dimension; ++j) {
 198             if (remove[j])
 199                 ++n;
 200             else
 201                 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
 202         }
 203         value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
 204     }
 205     *T = Matrix_Alloc(P->Dimension+1, Q->Dimension+1);
 206     for (i = 0, n = 0; i < P->Dimension; ++i) {
 207         if (remove[i]) {
 208             value_oppose((*T)->p[i][Q->Dimension], P->Constraint[n][1+P->Dimension]);
 209             ++n;
 210         } else
 211             value_set_si((*T)->p[i][i-n], 1);
 212     }
 213     value_set_si((*T)->p[i][i-n], 1);
 214     POL_SET(Q, POL_VALID);
 215     if (POL_HAS(P, POL_INEQUALITIES))
 216         POL_SET(Q, POL_INEQUALITIES);
 217     if (POL_HAS(P, POL_FACETS))
 218         POL_SET(Q, POL_FACETS);
 219     if (POL_HAS(P, POL_POINTS))
 220         POL_SET(Q, POL_POINTS);
 221     if (POL_HAS(P, POL_VERTICES))
 222         POL_SET(Q, POL_VERTICES);
 223     free(remove);
 224     return Q;
 225 }
 226
 227 /* This function implements the algorithm described in
 228  * "An Implementation of the Generalized Basis Reduction Algorithm
 229  *  for Integer Programming" of Cook el al. to find an integer point
 230  * in a polyhedron.
 231  * If the polyhedron is unbounded, we first remove its rays.
 232  */
 233 Vector *Polyhedron_Sample(Polyhedron *P, unsigned MaxRays)
 234 {
 235     int i, j;
 236     Vector *sample = NULL;
 237     Polyhedron *Q;
 238     Matrix *T, *inv, *M;
 239     Value min, max, tmp;
 240     Vector *v;
 241     int ok;
 242
 243     POL_ENSURE_VERTICES(P);
 244     if (emptyQ(P))
 245         return NULL;
 246
 247     if (P->Dimension == 0) {
 248         sample = Vector_Alloc(1);
 249         value_set_si(sample->p[0], 1);
 250         return sample;
 251     }
 252
 253     for (i = 0; i < P->NbRays; ++i)
 254         if (value_one_p(P->Ray[i][1+P->Dimension])) {
 255             sample = Vector_Alloc(P->Dimension+1);
 256             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
 257             return sample;
 258         }
 259
 260     Q = remove_ray(P, MaxRays);
 261     if (Q) {
 262         sample = Polyhedron_Sample(Q, MaxRays);
 263         Polyhedron_Free(Q);
 264         return sample;
 265     }
 266
 267     Matrix *basis = reduced_basis(P);
 268
 269     T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 270     inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 271     for (i = 0; i < P->Dimension; ++i)
 272         for (j = 0; j < P->Dimension; ++j)
 273             value_assign(T->p[i][j], basis->p[i][j]);
 274     value_set_si(T->p[P->Dimension][P->Dimension], 1);
 275     Matrix_Free(basis);
 276
 277     M = Matrix_Copy(T);
 278     ok = Matrix_Inverse(M, inv);
 279     assert(ok);
 280     Matrix_Free(M);
 281
 282     Q = Polyhedron_Image(P, T, MaxRays);
 283
 284     POL_ENSURE_VERTICES(Q);
 285
 286     value_init(min);
 287     value_init(max);
 288     value_init(tmp);
 289
 290     mpz_cdiv_q(min, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
 291     mpz_fdiv_q(max, Q->Ray[0][1], Q->Ray[0][1+Q->Dimension]);
 292
 293     for (j = 1; j < Q->NbRays; ++j) {
 294         mpz_cdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
 295         if (value_lt(tmp, min))
 296             value_assign(min, tmp);
 297         mpz_fdiv_q(tmp, Q->Ray[j][1], Q->Ray[j][1+Q->Dimension]);
 298         if (value_gt(tmp, max))
 299             value_assign(max, tmp);
 300     }
 301
 302     v = Vector_Alloc(1+Q->Dimension+1);
 303     value_set_si(v->p[1], -1);
 304
 305     for (value_assign(tmp, min); value_le(tmp, max); value_increment(tmp, tmp)) {
 306         Polyhedron *R, *S;
 307         Matrix *T;
 308         Vector *S_sample;
 309         value_assign(v->p[1+Q->Dimension], tmp);
 310
 311         R = AddConstraints(v->p, 1, Q, MaxRays);
 312         R = DomainConstraintSimplify(R, MaxRays);
 313         if (emptyQ(R)) {
 314             Polyhedron_Free(R);
 315             continue;
 316         }
 317
 318         S = Polyhedron_RemoveFixedColumns(R, &T);
 319         Polyhedron_Free(R);
 320         S_sample = Polyhedron_Sample(S, MaxRays);
 321         Polyhedron_Free(S);
 322         if (S_sample) {
 323             Vector *Q_sample = Vector_Alloc(Q->Dimension + 1);
 324             Matrix_Vector_Product(T, S_sample->p, Q_sample->p);
 325             Matrix_Free(T);
 326             Vector_Free(S_sample);
 327             sample = Vector_Alloc(P->Dimension + 1);
 328             Matrix_Vector_Product(inv, Q_sample->p, sample->p);
 329             Vector_Free(Q_sample);
 330             break;
 331         }
 332         Matrix_Free(T);
 333     }
 334
 335     Matrix_Free(T);
 336     Matrix_Free(inv);
 337     Polyhedron_Free(Q);
 338     Vector_Free(v);
 339
 340     value_clear(min);
 341     value_clear(max);
 342     value_clear(tmp);
 343
 344     return sample;
 345 }