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