sample.c

   1 #include <barvinok/util.h>
   2 #include <barvinok/basis_reduction.h>
   3 #include <barvinok/sample.h>
   4 #include <barvinok/options.h>
   5 #include "polysign.h"
   6
   7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   8
   9 /* If P has no rays, then we return NULL.
  10  * Otherwise, look for the coordinate axis with the smallest maximal non-zero
  11  * coefficient over all rays and a constraint that bounds the values on
  12  * this axis to the maximal value over the vertices plus the above maximal
  13  * non-zero coefficient times the number of rays minus 1.
  14  * Any integer point x outside this region is the sum of a point inside
  15  * the region and an integer multiple of the rays.
  16  * Write x = \sum_i a_i v_i + \sum_j b_j r_j
  17  * with \sum_i a_i = 1.
  18  * Then x = \sum_i a_i v_i + \sum_j {b_j} r_j + \sum_j [b_j] r_j
  19  * with y = \sum_i a_i v_i + \sum_j {b_j} r_j a point inside the region.
  20  */
  21 static Polyhedron *remove_ray(Polyhedron *P, unsigned MaxRays)
  22 {
  23     int r = 0;
  24     Vector *min, *max, *c;
  25     int i;
  26     Value s, v, tmp;
  27     int pos;
  28     Polyhedron *R;
  29     int rays;
  30
  31     POL_ENSURE_VERTICES(P);
  32     if (P->NbBid == 0)
  33         for (; r < P->NbRays; ++r)
  34             if (value_zero_p(P->Ray[r][P->Dimension+1]))
  35                 break;
  36     if (P->NbBid == 0 && r == P->NbRays)
  37         return NULL;
  38
  39     max = Vector_Alloc(P->Dimension);
  40     min = Vector_Alloc(P->Dimension);
  41     for (r = 0; r < P->NbBid; ++r)
  42         for (i = 0 ; i < P->Dimension; ++i)
  43             if (value_abs_gt(P->Ray[r][1+i], max->p[i]))
  44                 value_absolute(max->p[i], P->Ray[r][1+i]);
  45
  46     for (i = 0 ; i < P->Dimension; ++i)
  47         value_oppose(min->p[i], max->p[i]);
  48
  49     rays = P->NbBid;
  50     for (r = P->NbBid; r < P->NbRays; ++r) {
  51         if (value_notzero_p(P->Ray[r][P->Dimension+1]))
  52             continue;
  53         for (i = 0 ; i < P->Dimension; ++i) {
  54             if (value_gt(P->Ray[r][1+i], max->p[i]))
  55                 value_assign(max->p[i], P->Ray[r][1+i]);
  56             if (value_lt(P->Ray[r][1+i], min->p[i]))
  57                 value_assign(min->p[i], P->Ray[r][1+i]);
  58         }
  59         ++rays;
  60     }
  61
  62     value_init(s);
  63     value_init(v);
  64     value_init(tmp);
  65
  66     for (i = 0 ; i < P->Dimension; ++i) {
  67         if (value_notzero_p(min->p[i]) &&
  68             (value_zero_p(s) || value_abs_lt(min->p[i], s))) {
  69             value_assign(s, min->p[i]);
  70             pos = i;
  71         }
  72         if (value_notzero_p(max->p[i]) &&
  73             (value_zero_p(s) || value_abs_lt(max->p[i], s))) {
  74             value_assign(s, max->p[i]);
  75             pos = i;
  76         }
  77     }
  78
  79     for (r = P->NbBid; r < P->NbRays; ++r)
  80         if (value_notzero_p(P->Ray[r][P->Dimension+1]))
  81             break;
  82
  83     if (value_pos_p(s))
  84         mpz_cdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  85     else
  86         mpz_fdiv_q(v, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  87
  88     for ( ; r < P->NbRays; ++r) {
  89         if (value_zero_p(P->Ray[r][P->Dimension+1]))
  90             continue;
  91
  92         if (value_pos_p(s)) {
  93             mpz_cdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  94             if (value_gt(tmp, v))
  95                 value_assign(v, tmp);
  96         } else {
  97             mpz_fdiv_q(tmp, P->Ray[r][1+pos], P->Ray[r][P->Dimension+1]);
  98             if (value_lt(tmp, v))
  99                 value_assign(v, tmp);
 100         }
 101     }
 102
 103     c = Vector_Alloc(1+P->Dimension+1);
 104
 105     value_set_si(tmp, rays);
 106     value_addmul(v, tmp, s);
 107     value_set_si(c->p[0], 1);
 108     if (value_pos_p(s)) {
 109         value_set_si(c->p[1+pos], -1);
 110         value_assign(c->p[1+P->Dimension], v);
 111     } else {
 112         value_set_si(c->p[1+pos], 1);
 113         value_oppose(c->p[1+P->Dimension], v);
 114     }
 115     value_decrement(c->p[1+P->Dimension], c->p[1+P->Dimension]);
 116
 117     R = AddConstraints(c->p, 1, P, MaxRays);
 118
 119     Vector_Free(c);
 120
 121     Vector_Free(min);
 122     Vector_Free(max);
 123
 124     value_clear(tmp);
 125     value_clear(s);
 126     value_clear(v);
 127
 128     return R;
 129 }
 130
 131 static void print_minmax(Polyhedron *P)
 132 {
 133     int i, j;
 134     POL_ENSURE_VERTICES(P);
 135     Polyhedron_Print(stderr, P_VALUE_FMT, P);
 136     for (i = 0; i < P->Dimension; ++i) {
 137         Value min, max, tmp;
 138         value_init(min);
 139         value_init(max);
 140         value_init(tmp);
 141
 142         mpz_cdiv_q(min, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
 143         mpz_fdiv_q(max, P->Ray[0][1+i], P->Ray[0][1+P->Dimension]);
 144
 145         for (j = 1; j < P->NbRays; ++j) {
 146             mpz_cdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
 147             if (value_lt(tmp, min))
 148                 value_assign(min, tmp);
 149             mpz_fdiv_q(tmp, P->Ray[j][1+i], P->Ray[j][1+P->Dimension]);
 150             if (value_gt(tmp, max))
 151                 value_assign(max, tmp);
 152         }
 153         fprintf(stderr, "i: %d,  min: ", i);
 154         value_print(stderr, VALUE_FMT, min);
 155         fprintf(stderr, ", max: ");
 156         value_print(stderr, VALUE_FMT, max);
 157         fprintf(stderr, "\n");
 158
 159         value_clear(min);
 160         value_clear(max);
 161         value_clear(tmp);
 162     }
 163 }
 164
 165 /* Remove coordinates that have a fixed value and return the matrix
 166  * that adds these fixed coordinates again through T.
 167  */
 168 static Polyhedron *Polyhedron_RemoveFixedColumns(Polyhedron *P, Matrix **T)
 169 {
 170     int i, j, k, n;
 171     int dim = P->Dimension;
 172     int *remove = ALLOCN(int, dim);
 173     Polyhedron *Q;
 174     int NbEq;
 175
 176     assert(POL_HAS(P, POL_INEQUALITIES));
 177     for (i = 0; i < dim; ++i)
 178         remove[i] = -1;
 179     NbEq = 0;
 180     for (i = 0; i < P->NbEq; ++i) {
 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         remove[pos] = i;
 185         ++NbEq;
 186     }
 187     assert(NbEq > 0);
 188     Q = Polyhedron_Alloc(P->Dimension-NbEq, P->NbConstraints-NbEq, P->NbRays);
 189     for (i = 0, k = 0; i < P->NbConstraints; ++i) {
 190         if (i < P->NbEq) {
 191             int pos = First_Non_Zero(P->Constraint[i]+1, dim);
 192             if (First_Non_Zero(P->Constraint[i]+1+pos+1, dim-pos-1) == -1)
 193                 continue;
 194         }
 195         value_assign(Q->Constraint[k][0], P->Constraint[i][0]);
 196         for (j = 0, n = 0; j < P->Dimension; ++j) {
 197             if (remove[j] != -1)
 198                 ++n;
 199             else
 200                 value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
 201         }
 202         value_assign(Q->Constraint[k][1+j-n], P->Constraint[i][1+j]);
 203         ++k;
 204     }
 205     for (i = 0; i < Q->NbRays; ++i) {
 206         value_assign(Q->Ray[i][0], P->Ray[i][0]);
 207         for (j = 0, n = 0; j < P->Dimension; ++j) {
 208             if (remove[j] != -1)
 209                 ++n;
 210             else
 211                 value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
 212         }
 213         value_assign(Q->Ray[i][1+j-n], P->Ray[i][1+j]);
 214     }
 215     *T = Matrix_Alloc(P->Dimension+1, Q->Dimension+1);
 216     for (i = 0, n = 0; i < P->Dimension; ++i) {
 217         if (remove[i] != -1) {
 218             value_oppose((*T)->p[i][Q->Dimension],
 219                          P->Constraint[remove[i]][1+P->Dimension]);
 220             ++n;
 221         } else
 222             value_set_si((*T)->p[i][i-n], 1);
 223     }
 224     value_set_si((*T)->p[i][i-n], 1);
 225     POL_SET(Q, POL_VALID);
 226     if (POL_HAS(P, POL_INEQUALITIES))
 227         POL_SET(Q, POL_INEQUALITIES);
 228     if (POL_HAS(P, POL_FACETS))
 229         POL_SET(Q, POL_FACETS);
 230     if (POL_HAS(P, POL_POINTS))
 231         POL_SET(Q, POL_POINTS);
 232     if (POL_HAS(P, POL_VERTICES))
 233         POL_SET(Q, POL_VERTICES);
 234     free(remove);
 235     return Q;
 236 }
 237
 238 static Polyhedron *remove_all_equalities(Polyhedron *P, Matrix **T,
 239                                          unsigned MaxRays)
 240 {
 241     /* Matrix "view" of equalities */
 242     Matrix M;
 243     M.NbRows = P->NbEq;
 244     M.NbColumns = P->Dimension+2;
 245     M.p_Init = P->p_Init;
 246     M.p = P->Constraint;
 247
 248     *T = compress_variables(&M, 0);
 249
 250     if (!*T)
 251         return NULL;
 252     P = Polyhedron_Preimage(P, *T, MaxRays);
 253
 254     return P;
 255 }
 256
 257 static Vector *product_sample(Polyhedron *P, Matrix *T,
 258                               struct barvinok_options *options)
 259 {
 260     int i;
 261     Vector *sample = NULL;
 262     Vector *tmp = Vector_Alloc(T->NbRows);
 263     i = 0;
 264     for (; P; P = P->next) {
 265         Vector *P_sample;
 266         Polyhedron *next = P->next;
 267         P->next = NULL;
 268         P_sample = Polyhedron_Sample(P, options);
 269         if (!P_sample) {
 270             Vector_Free(tmp);
 271             tmp = NULL;
 272             break;
 273         }
 274         Vector_Copy(P_sample->p, tmp->p+i, P->Dimension);
 275         Vector_Free(P_sample);
 276         i += P->Dimension;
 277         P->next = next;
 278     }
 279     if (tmp) {
 280         sample = Vector_Alloc(T->NbRows + 1);
 281         Matrix_Vector_Product(T, tmp->p, sample->p);
 282         value_set_si(sample->p[T->NbRows], 1);
 283         Vector_Free(tmp);
 284     }
 285     return sample;
 286 }
 287
 288 /* This function implements the algorithm described in
 289  * "An Implementation of the Generalized Basis Reduction Algorithm
 290  *  for Integer Programming" of Cook el al. to find an integer point
 291  * in a polyhedron.
 292  * If the polyhedron is unbounded, we first remove its rays.
 293  */
 294 Vector *Polyhedron_Sample(Polyhedron *P, struct barvinok_options *options)
 295 {
 296     int i, j;
 297     Vector *sample = NULL, *obj = NULL;
 298     Polyhedron *Q;
 299     Matrix *inv = NULL;
 300     Value min, max, tmp;
 301     Vector *v;
 302     int ok;
 303     enum lp_result res;
 304     Matrix *T;
 305
 306     if (emptyQ2(P))
 307         return NULL;
 308
 309     if (P->Dimension == 0) {
 310         sample = Vector_Alloc(1);
 311         value_set_si(sample->p[0], 1);
 312         return sample;
 313     }
 314
 315     if (P->Dimension == 1)
 316         POL_ENSURE_VERTICES(P);
 317
 318     for (i = 0; i < P->NbRays; ++i)
 319         if (value_one_p(P->Ray[i][1+P->Dimension])) {
 320             sample = Vector_Alloc(P->Dimension+1);
 321             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
 322             return sample;
 323         }
 324
 325     if (P->NbEq > 0) {
 326         Vector *Q_sample;
 327         Polyhedron *Q = remove_all_equalities(P, &T, options->MaxRays);
 328         if (!Q)
 329             return NULL;
 330         Q_sample = Polyhedron_Sample(Q, options);
 331         Polyhedron_Free(Q);
 332         if (Q_sample) {
 333             sample = Vector_Alloc(P->Dimension + 1);
 334             Matrix_Vector_Product(T, Q_sample->p, sample->p);
 335             Vector_Free(Q_sample);
 336         }
 337         Matrix_Free(T);
 338         return sample;
 339     }
 340
 341     Q = Polyhedron_Factor(P, 0, &T, options->MaxRays);
 342     if (Q) {
 343         sample = product_sample(Q, T, options);
 344         Domain_Free(Q);
 345         Matrix_Free(T);
 346         return sample;
 347     }
 348
 349     value_init(min);
 350     value_init(max);
 351
 352     obj = Vector_Alloc(P->Dimension+1);
 353     value_set_si(obj->p[0], 1);
 354     res = polyhedron_range(P, obj->p, obj->p[0], &min, &max, options);
 355     if (res == lp_unbounded) {
 356 unbounded:
 357         value_clear(min);
 358         value_clear(max);
 359         Vector_Free(obj);
 360
 361         Q = remove_ray(P, options->MaxRays);
 362         assert(Q);
 363         sample = Polyhedron_Sample(Q, options);
 364         Polyhedron_Free(Q);
 365         return sample;
 366     }
 367     if (res == lp_empty) {
 368         value_clear(min);
 369         value_clear(max);
 370         Vector_Free(obj);
 371         return NULL;
 372     }
 373     assert(res == lp_ok);
 374
 375     if (value_eq(min, max)) {
 376         Q = P;
 377     } else {
 378         Matrix *M, *T;
 379         Matrix *basis = Polyhedron_Reduced_Basis(P, options);
 380
 381         if (!basis)
 382             goto unbounded;
 383         T = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 384         inv = Matrix_Alloc(P->Dimension+1, P->Dimension+1);
 385         for (i = 0; i < P->Dimension; ++i)
 386             for (j = 0; j < P->Dimension; ++j)
 387                 value_assign(T->p[i][j], basis->p[i][j]);
 388         value_set_si(T->p[P->Dimension][P->Dimension], 1);
 389         Matrix_Free(basis);
 390
 391         M = Matrix_Copy(T);
 392         ok = Matrix_Inverse(M, inv);
 393         assert(ok);
 394         Matrix_Free(M);
 395
 396         Q = Polyhedron_Image(P, T, options->MaxRays);
 397         res = polyhedron_range(Q, obj->p, obj->p[0], &min, &max, options);
 398         assert(res == lp_ok);
 399
 400         Matrix_Free(T);
 401     }
 402
 403     value_init(tmp);
 404
 405     v = Vector_Alloc(1+Q->Dimension+1);
 406     value_set_si(v->p[1], -1);
 407
 408     for (value_assign(tmp, min); value_le(tmp, max); value_increment(tmp, tmp)) {
 409         Polyhedron *R, *S;
 410         Matrix *T;
 411         Vector *S_sample;
 412         value_assign(v->p[1+Q->Dimension], tmp);
 413
 414         R = AddConstraints(v->p, 1, Q, options->MaxRays);
 415         R = DomainConstraintSimplify(R, options->MaxRays);
 416         if (emptyQ(R)) {
 417             Polyhedron_Free(R);
 418             continue;
 419         }
 420
 421         S = Polyhedron_RemoveFixedColumns(R, &T);
 422         Polyhedron_Free(R);
 423         S_sample = Polyhedron_Sample(S, options);
 424         Polyhedron_Free(S);
 425         if (S_sample) {
 426             Vector *Q_sample = obj;
 427             obj = NULL;
 428             Matrix_Vector_Product(T, S_sample->p, Q_sample->p);
 429             Matrix_Free(T);
 430             Vector_Free(S_sample);
 431             if (!inv)
 432                 sample = Q_sample;
 433             else {
 434                 sample = Vector_Alloc(P->Dimension + 1);
 435                 Matrix_Vector_Product(inv, Q_sample->p, sample->p);
 436                 Vector_Free(Q_sample);
 437             }
 438             break;
 439         }
 440         Matrix_Free(T);
 441     }
 442
 443     if (obj)
 444         Vector_Free(obj);
 445     if (inv)
 446         Matrix_Free(inv);
 447     if (P != Q)
 448         Polyhedron_Free(Q);
 449     Vector_Free(v);
 450
 451     value_clear(min);
 452     value_clear(max);
 453     value_clear(tmp);
 454
 455     return sample;
 456 }