hilbert.c

   1 #include <stdlib.h>
   2 #define Vector ZSolveVector
   3 #define Matrix ZSolveMatrix
   4 #include "zsolve/libzsolve.h"
   5 #undef Vector
   6 #undef Matrix
   7 #include <barvinok/options.h>
   8 #include <barvinok/util.h>
   9 #include "hilbert.h"
  10 #include "polysign.h"
  11
  12 static ZSolveMatrix Matrix2zsolve(Matrix *M)
  13 {
  14     int i, j;
  15     ZSolveMatrix zmatrix;
  16
  17     zmatrix = createMatrix(M->NbColumns-2, M->NbRows);
  18     for (i = 0; i < M->NbRows; ++i)
  19         for (j = 0; j < M->NbColumns-2; ++j) {
  20             assert(mpz_cmp_si(M->p[i][1+j], -MAXINT) > 0);
  21             assert(mpz_cmp_si(M->p[i][1+j], MAXINT) < 0);
  22             zmatrix->Data[i*zmatrix->Width+j] = mpz_get_si(M->p[i][1+j]);
  23         }
  24
  25     return zmatrix;
  26 }
  27
  28 static Matrix *VectorArray2Matrix(VectorArray array, unsigned cols)
  29 {
  30     int i, j;
  31     Matrix *M = Matrix_Alloc(array->Size, cols+1);
  32
  33     for (i = 0; i < array->Size; ++i) {
  34         for (j = 0; j < cols; ++j)
  35             value_set_si(M->p[i][j], array->Data[i][j]);
  36         value_set_si(M->p[i][cols], 1);
  37     }
  38     return M;
  39 }
  40
  41 static void Polyhedron_Remove_Positivity_Constraint(Polyhedron *P)
  42 {
  43     int i;
  44
  45     for (i = 0; i < P->NbConstraints; ++i) {
  46         if (First_Non_Zero(P->Constraint[i]+1, P->Dimension) != -1)
  47             continue;
  48         if (i < P->NbConstraints-1)
  49             Vector_Exchange(P->Constraint[i],
  50                             P->Constraint[P->NbConstraints-1],
  51                             P->Dimension+2);
  52         P->NbConstraints--;
  53         --i;
  54     }
  55 }
  56
  57 /* Return
  58  *             T 0
  59  *             0 1
  60  */
  61 static Matrix *homogenize(Matrix *T)
  62 {
  63     int i;
  64     Matrix *H = Matrix_Alloc(T->NbRows+1, T->NbColumns+1);
  65
  66     for (i = 0; i < T->NbRows; ++i)
  67         Vector_Copy(T->p[i], H->p[i], T->NbColumns);
  68     value_set_si(H->p[T->NbRows][T->NbColumns], 1);
  69     return H;
  70 }
  71
  72 static Matrix *Polyhedron2standard_form(Polyhedron *P, Matrix **T)
  73 {
  74     int i, j;
  75     int rows;
  76     unsigned dim = P->Dimension;
  77     Matrix *M2;
  78     Matrix *H, *U;
  79
  80     assert(P->NbEq == 0);
  81     Polyhedron_Remove_Positivity_Constraint(P);
  82     for (i = 0; i < P->NbConstraints; ++i)
  83         assert(value_zero_p(P->Constraint[i][1+dim]));
  84
  85     H = standard_constraints(P, 0, &rows, &U);
  86     *T = homogenize(U);
  87     Matrix_Free(U);
  88
  89     M2 = Matrix_Alloc(rows, 2+dim+rows);
  90
  91     for (i = dim; i < H->NbRows; ++i) {
  92         Vector_Copy(H->p[i], M2->p[i-dim]+1, dim);
  93         value_set_si(M2->p[i-dim][1+i], -1);
  94     }
  95     for (i = 0, j = H->NbRows-dim; i < dim; ++i) {
  96         if (First_Non_Zero(H->p[i], i) == -1)
  97             continue;
  98         Vector_Oppose(H->p[i], M2->p[j]+1, dim);
  99         value_set_si(M2->p[j][1+j+dim], 1);
 100         ++j;
 101     }
 102     Matrix_Free(H);
 103     return M2;
 104 }
 105
 106 /* Assumes C is a linear cone (i.e. with apex zero).
 107  * All equalities are removed first to speed up the computation
 108  * in zsolve.
 109  */
 110 Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays)
 111 {
 112     unsigned dim;
 113     int i;
 114     Matrix *M2, *M3, *T;
 115     Matrix *CV = NULL;
 116     LinearSystem initialsystem;
 117     ZSolveMatrix matrix;
 118     ZSolveVector rhs;
 119     ZSolveContext ctx;
 120
 121     remove_all_equalities(&C, NULL, NULL, &CV, 0, MaxRays);
 122     dim = C->Dimension;
 123
 124     for (i = 0; i < C->NbConstraints; ++i)
 125         assert(value_zero_p(C->Constraint[i][1+dim]) ||
 126                First_Non_Zero(C->Constraint[i]+1, dim) == -1);
 127
 128     M2 = Polyhedron2standard_form(C, &T);
 129     matrix = Matrix2zsolve(M2);
 130     Matrix_Free(M2);
 131
 132     rhs = createVector(matrix->Height);
 133     for (i = 0; i < matrix->Height; i++)
 134         rhs[i] = 0;
 135
 136     initialsystem = createLinearSystem();
 137     setLinearSystemMatrix(initialsystem, matrix);
 138     deleteMatrix(matrix);
 139
 140     setLinearSystemRHS(initialsystem, rhs);
 141     deleteVector(rhs);
 142
 143     setLinearSystemLimit(initialsystem, -1, 0, MAXINT, 0);
 144     setLinearSystemEquationType(initialsystem, -1, EQUATION_EQUAL, 0);
 145
 146     ctx = createZSolveContextFromSystem(initialsystem, NULL, 0, 0, NULL, NULL);
 147     zsolveSystem(ctx, 0);
 148
 149     M2 = VectorArray2Matrix(ctx->Homs, C->Dimension);
 150     deleteZSolveContext(ctx, 1);
 151     Matrix_Transposition(T);
 152     M3 = Matrix_Alloc(M2->NbRows, M2->NbColumns);
 153     Matrix_Product(M2, T, M3);
 154     Matrix_Free(M2);
 155     Matrix_Free(T);
 156
 157     if (CV) {
 158         Matrix *T, *M;
 159         T = Transpose(CV);
 160         M = Matrix_Alloc(M3->NbRows, T->NbColumns);
 161         Matrix_Product(M3, T, M);
 162         Matrix_Free(M3);
 163         Matrix_Free(CV);
 164         Matrix_Free(T);
 165         Polyhedron_Free(C);
 166         M3 = M;
 167     }
 168
 169     return M3;
 170 }
 171
 172 /* Assumes no two arrays of values are the same, so we can just
 173  * look for the first different elements, without knowing the
 174  * length of the arrays.
 175  */
 176 static int lex_cmp(const void *va, const void *vb)
 177 {
 178     int i;
 179     const Value *a = *(const Value **)va;
 180     const Value *b = *(const Value **)vb;
 181
 182     for (i = 0; ; ++i) {
 183         int sign = mpz_cmp(a[i], b[i]);
 184         if (sign)
 185             return sign;
 186     }
 187 }
 188
 189 /* Compute integer hull of truncated linear cone C, i.e., of C with
 190  * the origin removed.
 191  * Here, we do this by first computing the Hilbert basis of C
 192  * and then discarding elements from this basis that are rational
 193  * overconvex combinations of other elements in the basis.
 194  */
 195 Matrix *Cone_Integer_Hull(Polyhedron *C, struct barvinok_options *options)
 196 {
 197     int i, j, k;
 198     Matrix *hilbert = Cone_Hilbert_Basis(C, options->MaxRays);
 199     Matrix *rays, *hull;
 200     unsigned dim = C->Dimension;
 201     Value tmp;
 202     unsigned MaxRays = options->MaxRays;
 203
 204     /* When checking for redundant points below, we want to
 205      * check if there are any _rational_ solutions.
 206      */
 207     POL_UNSET(options->MaxRays, POL_INTEGER);
 208
 209     POL_ENSURE_VERTICES(C);
 210     rays = Matrix_Alloc(C->NbRays-1, C->Dimension);
 211     for (i = 0, j = 0; i < C->NbRays; ++i) {
 212         if (value_notzero_p(C->Ray[i][1+C->Dimension]))
 213             continue;
 214         Vector_Copy(C->Ray[i]+1, rays->p[j++], C->Dimension);
 215     }
 216
 217     /* We only sort the pointers into the big Value array */
 218     qsort(rays->p, rays->NbRows, sizeof(Value *), lex_cmp);
 219     qsort(hilbert->p, hilbert->NbRows, sizeof(Value *), lex_cmp);
 220
 221     /* Remove rays from Hilbert basis */
 222     for (i = 0, j = 0, k = 0; i < hilbert->NbRows && j < rays->NbRows; ++i) {
 223         if (Vector_Equal(hilbert->p[i], rays->p[j], C->Dimension))
 224             ++j;
 225         else
 226             hilbert->p[k++] = hilbert->p[i];
 227     }
 228     hilbert->NbRows = k;
 229
 230     /* Now remove points that are overconvex combinations of other points */
 231     value_init(tmp);
 232     for (i = 0; hilbert->NbRows > 1 && i < hilbert->NbRows; ++i) {
 233         Matrix *LP;
 234         Vector *obj;
 235         int nray = rays->NbRows;
 236         int npoint = hilbert->NbRows;
 237         enum lp_result result;
 238
 239         LP = Matrix_Alloc(dim + 1 + nray + (npoint-1), 2 + nray + (npoint-1));
 240         for (j = 0; j < dim; ++j) {
 241             for (k = 0; k < nray; ++k)
 242                 value_assign(LP->p[j][k+1], rays->p[k][j]);
 243             for (k = 0; k < npoint; ++k) {
 244                 if (k == i)
 245                     value_oppose(LP->p[j][1+nray+npoint-1], hilbert->p[k][j]);
 246                 else
 247                     value_assign(LP->p[j][1+nray+k-(k>i)], hilbert->p[k][j]);
 248             }
 249         }
 250         value_set_si(LP->p[dim][0], 1);
 251         for (k = 0; k < nray+npoint-1; ++k)
 252             value_set_si(LP->p[dim][1+k], 1);
 253         value_set_si(LP->p[dim][LP->NbColumns-1], -1);
 254         for (k = 0; k < LP->NbColumns-2; ++k) {
 255             value_set_si(LP->p[dim+1+k][0], 1);
 256             value_set_si(LP->p[dim+1+k][1+k], 1);
 257         }
 258
 259         /* Somewhat arbitrary objective function. */
 260         obj = Vector_Alloc(LP->NbColumns-1);
 261         value_set_si(obj->p[0], 1);
 262         value_set_si(obj->p[obj->Size-1], 1);
 263
 264         result = constraints_opt(LP, obj->p, obj->p[0], lp_min, &tmp,
 265                                  options);
 266
 267         /* If the LP is not empty, the point can be discarded */
 268         if (result != lp_empty) {
 269             hilbert->NbRows--;
 270             if (i < hilbert->NbRows)
 271                 hilbert->p[i] = hilbert->p[hilbert->NbRows];
 272             --i;
 273         }
 274
 275         Matrix_Free(LP);
 276         Vector_Free(obj);
 277     }
 278     value_clear(tmp);
 279
 280     hull = Matrix_Alloc(rays->NbRows + hilbert->NbRows, dim+1);
 281     for (i = 0; i < rays->NbRows; ++i) {
 282         Vector_Copy(rays->p[i], hull->p[i], dim);
 283         value_set_si(hull->p[i][dim], 1);
 284     }
 285     for (i = 0; i < hilbert->NbRows; ++i) {
 286         Vector_Copy(hilbert->p[i], hull->p[rays->NbRows+i], dim);
 287         value_set_si(hull->p[rays->NbRows+i][dim], 1);
 288     }
 289     Matrix_Free(rays);
 290     Matrix_Free(hilbert);
 291
 292     options->MaxRays = MaxRays;
 293     return hull;
 294 }