compile parker
[barvinok.git] / hilbert.c
blob3ea20bd24edc759eabf59a8378a41722f9f2e338
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 "polysign.h"
13 static ZSolveMatrix Matrix2zsolve(Matrix *M)
15 int i, j;
16 ZSolveMatrix zmatrix;
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]);
26 return zmatrix;
29 static Matrix *VectorArray2Matrix(VectorArray array, unsigned cols)
31 int i, j;
32 Matrix *M = Matrix_Alloc(array->Size, cols+1);
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);
39 return M;
42 static void Polyhedron_Remove_Positivity_Constraint(Polyhedron *P)
44 int i;
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;
58 /* Return
59 * T 0
60 * 0 1
62 static Matrix *homogenize(Matrix *T)
64 int i;
65 Matrix *H = Matrix_Alloc(T->NbRows+1, T->NbColumns+1);
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;
73 static Matrix *Polyhedron2standard_form(Polyhedron *P, Matrix **T)
75 int i, j;
76 int rows;
77 unsigned dim = P->Dimension;
78 Matrix *M2;
79 Matrix *H, *U;
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]));
86 H = standard_constraints(P, 0, &rows, &U);
87 *T = homogenize(U);
88 Matrix_Free(U);
90 M2 = Matrix_Alloc(rows, 2+dim+rows);
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);
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;
103 Matrix_Free(H);
104 return M2;
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.
111 Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays)
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;
122 remove_all_equalities(&C, NULL, NULL, &CV, 0, MaxRays);
123 dim = C->Dimension;
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);
129 M2 = Polyhedron2standard_form(C, &T);
130 matrix = Matrix2zsolve(M2);
131 Matrix_Free(M2);
133 rhs = createVector(matrix->Height);
134 for (i = 0; i < matrix->Height; i++)
135 rhs[i] = 0;
137 initialsystem = createLinearSystem();
138 setLinearSystemMatrix(initialsystem, matrix);
139 deleteMatrix(matrix);
141 setLinearSystemRHS(initialsystem, rhs);
142 deleteVector(rhs);
144 setLinearSystemLimit(initialsystem, -1, 0, MAXINT, 0);
145 setLinearSystemEquationType(initialsystem, -1, EQUATION_EQUAL, 0);
147 ctx = createZSolveContextFromSystem(initialsystem, NULL, 0, 0, NULL, NULL);
148 zsolveSystem(ctx, 0);
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);
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;
170 return M3;
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.
177 static int lex_cmp(const void *va, const void *vb)
179 int i;
180 const Value *a = *(const Value **)va;
181 const Value *b = *(const Value **)vb;
183 for (i = 0; ; ++i) {
184 int sign = mpz_cmp(a[i], b[i]);
185 if (sign)
186 return sign;
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.
196 Matrix *Cone_Hilbert_Integer_Hull(Polyhedron *C,
197 struct barvinok_options *options)
199 int i, j, k;
200 Matrix *hilbert = Cone_Hilbert_Basis(C, options->MaxRays);
201 Matrix *rays, *hull;
202 unsigned dim = C->Dimension;
203 Value tmp;
204 unsigned MaxRays = options->MaxRays;
206 /* When checking for redundant points below, we want to
207 * check if there are any _rational_ solutions.
209 POL_UNSET(options->MaxRays, POL_INTEGER);
211 POL_ENSURE_VERTICES(C);
212 rays = Matrix_Alloc(C->NbRays-1, C->Dimension);
213 for (i = 0, j = 0; i < C->NbRays; ++i) {
214 if (value_notzero_p(C->Ray[i][1+C->Dimension]))
215 continue;
216 Vector_Copy(C->Ray[i]+1, rays->p[j++], C->Dimension);
219 /* We only sort the pointers into the big Value array */
220 qsort(rays->p, rays->NbRows, sizeof(Value *), lex_cmp);
221 qsort(hilbert->p, hilbert->NbRows, sizeof(Value *), lex_cmp);
223 /* Remove rays from Hilbert basis */
224 for (i = 0, j = 0, k = 0; i < hilbert->NbRows && j < rays->NbRows; ++i) {
225 if (Vector_Equal(hilbert->p[i], rays->p[j], C->Dimension))
226 ++j;
227 else
228 hilbert->p[k++] = hilbert->p[i];
230 hilbert->NbRows = k;
232 /* Now remove points that are overconvex combinations of other points */
233 value_init(tmp);
234 for (i = 0; hilbert->NbRows > 1 && i < hilbert->NbRows; ++i) {
235 Matrix *LP;
236 Vector *obj;
237 int nray = rays->NbRows;
238 int npoint = hilbert->NbRows;
239 enum lp_result result;
241 LP = Matrix_Alloc(dim + 1 + nray + (npoint-1), 2 + nray + (npoint-1));
242 for (j = 0; j < dim; ++j) {
243 for (k = 0; k < nray; ++k)
244 value_assign(LP->p[j][k+1], rays->p[k][j]);
245 for (k = 0; k < npoint; ++k) {
246 if (k == i)
247 value_oppose(LP->p[j][1+nray+npoint-1], hilbert->p[k][j]);
248 else
249 value_assign(LP->p[j][1+nray+k-(k>i)], hilbert->p[k][j]);
252 value_set_si(LP->p[dim][0], 1);
253 for (k = 0; k < nray+npoint-1; ++k)
254 value_set_si(LP->p[dim][1+k], 1);
255 value_set_si(LP->p[dim][LP->NbColumns-1], -1);
256 for (k = 0; k < LP->NbColumns-2; ++k) {
257 value_set_si(LP->p[dim+1+k][0], 1);
258 value_set_si(LP->p[dim+1+k][1+k], 1);
261 /* Somewhat arbitrary objective function. */
262 obj = Vector_Alloc(LP->NbColumns-1);
263 value_set_si(obj->p[0], 1);
264 value_set_si(obj->p[obj->Size-1], 1);
266 result = constraints_opt(LP, obj->p, obj->p[0], lp_min, &tmp,
267 options);
269 /* If the LP is not empty, the point can be discarded */
270 if (result != lp_empty) {
271 hilbert->NbRows--;
272 if (i < hilbert->NbRows)
273 hilbert->p[i] = hilbert->p[hilbert->NbRows];
274 --i;
277 Matrix_Free(LP);
278 Vector_Free(obj);
280 value_clear(tmp);
282 hull = Matrix_Alloc(rays->NbRows + hilbert->NbRows, dim+1);
283 for (i = 0; i < rays->NbRows; ++i) {
284 Vector_Copy(rays->p[i], hull->p[i], dim);
285 value_set_si(hull->p[i][dim], 1);
287 for (i = 0; i < hilbert->NbRows; ++i) {
288 Vector_Copy(hilbert->p[i], hull->p[rays->NbRows+i], dim);
289 value_set_si(hull->p[rays->NbRows+i][dim], 1);
291 Matrix_Free(rays);
292 Matrix_Free(hilbert);
294 options->MaxRays = MaxRays;
295 return hull;