sample.c

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