From 849526838d2a2626a713f943dfbdad6661864045 Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Thu, 1 Nov 2007 10:22:00 +0100 Subject: [PATCH] hilbert.c: add Cone_Integer_Hull for computing vertices of integer hull --- Makefile.am | 2 +- cone_integer_hull.c | 24 ++++++++++ hilbert.c | 128 ++++++++++++++++++++++++++++++++++++++++++++++++++++ hilbert.h | 3 ++ testlib.cc | 8 +++- 5 files changed, 163 insertions(+), 2 deletions(-) create mode 100644 cone_integer_hull.c diff --git a/Makefile.am b/Makefile.am index 6f2a446..a1a8cb9 100644 --- a/Makefile.am +++ b/Makefile.am @@ -33,7 +33,7 @@ noinst_PROGRAMS = test testlib randomtest \ barvinok_union polytope_volume test_approx \ barvinok_summate verify_lexsmaller piptest \ polyhedron_sample 4coins lexmin @bv_barvinok_maximize@ \ - cone_hilbert_basis + cone_hilbert_basis cone_integer_hull EXTRA_PROGRAMS = barvinok_maximize pkginclude_HEADERS = \ barvinok/NTL_QQ.h \ diff --git a/cone_integer_hull.c b/cone_integer_hull.c new file mode 100644 index 0000000..f9f04da --- /dev/null +++ b/cone_integer_hull.c @@ -0,0 +1,24 @@ +#include +#include +#include "hilbert.h" + +int main(int argc, char **argv) +{ + Matrix *M; + Polyhedron *C; + struct barvinok_options *options = barvinok_options_new_with_defaults(); + + M = Matrix_Read(); + C = Constraints2Polyhedron(M, options->MaxRays); + Matrix_Free(M); + + M = Cone_Integer_Hull(C, options); + + Polyhedron_Free(C); + + Matrix_Print(stdout, P_VALUE_FMT, M); + Matrix_Free(M); + + barvinok_options_free(options); + return 0; +} diff --git a/hilbert.c b/hilbert.c index 96db7ce..be3a8d3 100644 --- a/hilbert.c +++ b/hilbert.c @@ -1,9 +1,13 @@ +#include #define Vector ZSolveVector #define Matrix ZSolveMatrix #include "zsolve/libzsolve.h" #undef Vector #undef Matrix +#include +#include #include "hilbert.h" +#include "polysign.h" #include "topcom.h" static ZSolveMatrix Matrix2zsolve(Matrix *M) @@ -165,3 +169,127 @@ Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays) return M3; } + +/* Assumes no two arrays of values are the same, so we can just + * look for the first different elements, without knowing the + * length of the arrays. + */ +static int lex_cmp(const void *va, const void *vb) +{ + int i; + const Value *a = *(const Value **)va; + const Value *b = *(const Value **)vb; + + for (i = 0; ; ++i) { + int sign = mpz_cmp(a[i], b[i]); + if (sign) + return sign; + } +} + +/* Compute integer hull of truncated linear cone C, i.e., of C with + * the origin removed. + * Here, we do this by first computing the Hilbert basis of C + * and then discarding elements from this basis that are rational + * overconvex combinations of other elements in the basis. + */ +Matrix *Cone_Integer_Hull(Polyhedron *C, struct barvinok_options *options) +{ + int i, j, k; + Matrix *hilbert = Cone_Hilbert_Basis(C, options->MaxRays); + Matrix *rays, *hull; + unsigned dim = C->Dimension; + Value tmp; + unsigned MaxRays = options->MaxRays; + + /* When checking for redundant points below, we want to + * check if there are any _rational_ solutions. + */ + POL_UNSET(options->MaxRays, POL_INTEGER); + + POL_ENSURE_VERTICES(C); + rays = Matrix_Alloc(C->NbRays-1, C->Dimension); + for (i = 0, j = 0; i < C->NbRays; ++i) { + if (value_notzero_p(C->Ray[i][1+C->Dimension])) + continue; + Vector_Copy(C->Ray[i]+1, rays->p[j++], C->Dimension); + } + + /* We only sort the pointers into the big Value array */ + qsort(rays->p, rays->NbRows, sizeof(Value *), lex_cmp); + qsort(hilbert->p, hilbert->NbRows, sizeof(Value *), lex_cmp); + + /* Remove rays from Hilbert basis */ + for (i = 0, j = 0, k = 0; i < hilbert->NbRows && j < rays->NbRows; ++i) { + if (Vector_Equal(hilbert->p[i], rays->p[j], C->Dimension)) + ++j; + else + hilbert->p[k++] = hilbert->p[i]; + } + hilbert->NbRows = k; + + /* Now remove points that are overconvex combinations of other points */ + value_init(tmp); + for (i = 0; hilbert->NbRows > 1 && i < hilbert->NbRows; ++i) { + Matrix *LP; + Vector *obj; + int nray = rays->NbRows; + int npoint = hilbert->NbRows; + enum lp_result result; + + LP = Matrix_Alloc(dim + 1 + nray + (npoint-1), 2 + nray + (npoint-1)); + for (j = 0; j < dim; ++j) { + for (k = 0; k < nray; ++k) + value_assign(LP->p[j][k+1], rays->p[k][j]); + for (k = 0; k < npoint; ++k) { + if (k == i) + value_oppose(LP->p[j][1+nray+npoint-1], hilbert->p[k][j]); + else + value_assign(LP->p[j][1+nray+k-(k>i)], hilbert->p[k][j]); + } + } + value_set_si(LP->p[dim][0], 1); + for (k = 0; k < nray+npoint-1; ++k) + value_set_si(LP->p[dim][1+k], 1); + value_set_si(LP->p[dim][LP->NbColumns-1], -1); + for (k = 0; k < LP->NbColumns-2; ++k) { + value_set_si(LP->p[dim+1+k][0], 1); + value_set_si(LP->p[dim+1+k][1+k], 1); + } + + /* Somewhat arbitrary objective function. */ + obj = Vector_Alloc(LP->NbColumns-1); + value_set_si(obj->p[0], 1); + value_set_si(obj->p[obj->Size-1], 1); + + result = constraints_opt(LP, obj->p, obj->p[0], lp_min, &tmp, + options); + + /* If the LP is not empty, the point can be discarded */ + if (result != lp_empty) { + hilbert->NbRows--; + if (i < hilbert->NbRows) + hilbert->p[i] = hilbert->p[hilbert->NbRows]; + --i; + } + + Matrix_Free(LP); + Vector_Free(obj); + } + value_clear(tmp); + + hull = Matrix_Alloc(rays->NbRows + hilbert->NbRows, dim+1); + for (i = 0; i < rays->NbRows; ++i) { + Vector_Copy(rays->p[i], hull->p[i], dim); + value_set_si(hull->p[i][dim], 1); + } + for (i = 0; i < hilbert->NbRows; ++i) { + Vector_Copy(hilbert->p[i], hull->p[rays->NbRows+i], dim); + value_set_si(hull->p[rays->NbRows+i][dim], 1); + } + Matrix_Free(rays); + Matrix_Free(hilbert); + + options->MaxRays = MaxRays; + return hull; +} diff --git a/hilbert.h b/hilbert.h index aff1571..f3bad17 100644 --- a/hilbert.h +++ b/hilbert.h @@ -4,7 +4,10 @@ extern "C" { #endif +struct barvinok_options; + Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays); +Matrix *Cone_Integer_Hull(Polyhedron *C, struct barvinok_options *options); #if defined(__cplusplus) } diff --git a/testlib.cc b/testlib.cc index 4b19838..12f4a7c 100644 --- a/testlib.cc +++ b/testlib.cc @@ -480,10 +480,16 @@ int test_hilbert(struct barvinok_options *options) Matrix_Free(M); M = Cone_Hilbert_Basis(P, options->MaxRays); - Polyhedron_Free(P); assert(M->NbRows = 5); assert(M->NbColumns = 3); Matrix_Free(M); + + M = Cone_Integer_Hull(P, options); + assert(M->NbRows = 4); + assert(M->NbColumns = 3); + Matrix_Free(M); + + Polyhedron_Free(P); } int main(int argc, char **argv) -- 2.11.4.GIT