sample.c

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