hilbert.c

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