From f705694a358bbc6e1a7b37b3b80537c236118f98 Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Tue, 1 Apr 2008 18:16:34 +0200 Subject: [PATCH] barvinok_enumerate_e: optionally use parker's method --- Makefile.am | 15 +++- barvinok_enumerate_e.cc | 57 +++++++++++++++- configure.in | 1 + normalization.c | 177 ++++++++++++++++++++++++++++++++++++++++++++++-- normalization.h | 2 + version.c | 5 ++ 6 files changed, 249 insertions(+), 8 deletions(-) diff --git a/Makefile.am b/Makefile.am index d997fc7..fe3939e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -27,6 +27,8 @@ piplib/libpiplibMP.la: FORCE $(MAKE) $(AM_MAKEFLAGS) -C piplib libpiplibMP.la zsolve/libzsolve.la: FORCE $(MAKE) $(AM_MAKEFLAGS) -C zsolve libzsolve.la +parker/libparker.la: FORCE + $(MAKE) $(AM_MAKEFLAGS) -C parker libparker.la AM_CPPFLAGS = -I$(srcdir)/bernstein/include -I$(srcdir)/lib \ @POLYLIB_CPPFLAGS@ @PIPLIB_CPPFLAGS@ @ginac_CFLAGS@ @@ -235,10 +237,17 @@ barvinok_enumerate_e_SOURCES = \ verify_series.cc \ $(BEEO_SOURCES) EXTRA_barvinok_enumerate_e_SOURCES = fdstream.cc fdstream.h -barvinok_enumerate_e_CPPFLAGS = @OMEGA_CPPFLAGS@ $(AM_CPPFLAGS) -barvinok_enumerate_e_LDFLAGS = @OMEGA_LDFLAGS@ @LDFLAGS@ +if USE_PARKER + MONA_LIBS = -ldfa -lbdd -lmem + PARKER_LA = parker/libparker.la + PARKERLDFLAGS = -L$(MONAPATH)/DFA -L$(MONAPATH)/BDD -L$(MONAPATH)/Mem + PARKERCPPFLAGS = -I$(MONAPATH)/Mem -I$(MONAPATH)/DFA -I$(MONAPATH)/BDD +endif +barvinok_enumerate_e_CPPFLAGS = \ + $(PARKERCPPFLAGS) @OMEGA_CPPFLAGS@ $(AM_CPPFLAGS) +barvinok_enumerate_e_LDFLAGS = $(PARKERLDFLAGS) @OMEGA_LDFLAGS@ @LDFLAGS@ barvinok_enumerate_e_LDADD = \ - @OMEGA_LIBS@ libbarvinok-core.la \ + $(PARKER_LA) $(MONA_LIBS) @OMEGA_LIBS@ libbarvinok-core.la \ @POLYLIB_LIBS@ @PIPLIB_LIBS@ barvinok_bound_SOURCES = \ evalue_convert.cc \ diff --git a/barvinok_enumerate_e.cc b/barvinok_enumerate_e.cc index e90b343..afe80ef 100644 --- a/barvinok_enumerate_e.cc +++ b/barvinok_enumerate_e.cc @@ -10,11 +10,15 @@ #ifdef HAVE_OMEGA #include "omega/convert.h" #endif +#ifdef USE_PARKER +#include "parker/count_solutions.h" +#endif #include "skewed_genfun.h" #include "verify.h" #include "verif_ehrhart.h" #include "verify_series.h" #include "evalue_convert.h" +#include "normalization.h" /* The input of this example program is a polytope in combined * data and parameter space followed by two lines indicating @@ -25,8 +29,11 @@ * The polytope is in PolyLib notation. */ +#define ALLOC(t,p) p = (t*)malloc(sizeof(*p)) + struct argp_option argp_options[] = { { "omega", 'o', 0, 0 }, + { "parker", 'P', 0, 0 }, { "pip", 'p', 0, 0 }, { "series", 's', 0, 0 }, { "series", 's', 0, 0, "compute rational generating function" }, @@ -39,6 +46,7 @@ struct arguments { struct verify_options verify; struct convert_options convert; int omega; + int parker; int pip; int scarf; int series; @@ -71,6 +79,13 @@ error_t parse_opt(int key, char *arg, struct argp_state *state) error(0, 0, "--omega option not supported"); #endif break; + case 'P': +#ifdef USE_PARKER + arguments->parker = 1; +#else + error(0, 0, "--parker option not supported"); +#endif + break; case 'p': arguments->pip = 1; break; @@ -99,6 +114,43 @@ Polyhedron *Omega_simplify(Polyhedron *P, } #endif +#ifdef USE_PARKER +/* + * Use parker's method to compute the number of integer points in P. + * Since this method assumes all variables are non-negative, + * we have to transform the input polytope first. + */ +evalue *barvinok_enumerate_parker(Polyhedron *P, + unsigned exist, unsigned nparam, + unsigned MaxRays) +{ + Polyhedron *R; + evalue *res; + + assert(nparam == 0); + R = skew_to_positive_orthant(P, P->Dimension-exist, MaxRays); + Relation r = Polyhedron2relation(R, exist, 0, NULL); + Polyhedron_Free(R); + double d = count_solutions(r); + ALLOC(evalue, res); + value_init(res->d); + value_set_si(res->d, 0); + res->x.p = new_enode(partition, 2, 0); + EVALUE_SET_DOMAIN(res->x.p->arr[0], Universe_Polyhedron(0)); + value_set_si(res->x.p->arr[1].d, 1); + value_init(res->x.p->arr[1].x.n); + value_set_double(res->x.p->arr[1].x.n, d); + return res; +} +#else +evalue *barvinok_enumerate_parker(Polyhedron *P, + unsigned exist, unsigned nparam, + unsigned MaxRays) +{ + assert(0); +} +#endif + static void verify_results(Polyhedron *P, evalue *EP, gen_fun *gf, int exist, int nparam, arguments *options); @@ -125,6 +177,7 @@ int main(int argc, char **argv) arguments.verify.barvinok = options; arguments.omega = 0; + arguments.parker = 0; arguments.pip = 0; arguments.scarf = 0; arguments.series = 0; @@ -179,7 +232,9 @@ int main(int argc, char **argv) print_evalue(stdout, EP, param_name); } } else { - if (arguments.scarf) + if (arguments.parker) + EP = barvinok_enumerate_parker(A, exist, nparam, options->MaxRays); + else if (arguments.scarf) EP = barvinok_enumerate_scarf(A, exist, nparam, options); else if (arguments.pip && exist > 0) EP = barvinok_enumerate_pip_with_options(A, exist, nparam, options); diff --git a/configure.in b/configure.in index c943bf0..d772c0b 100644 --- a/configure.in +++ b/configure.in @@ -379,6 +379,7 @@ if test "$with_parker" != "no"; then MONAPATH=$with_mona AC_MSG_RESULT($MONAPATH) use_parker=true + AC_DEFINE(USE_PARKER,[],[use parker]) fi if test "x$with_mona" = "x"; then diff --git a/normalization.c b/normalization.c index a16b96e..ebb8cfe 100644 --- a/normalization.c +++ b/normalization.c @@ -32,9 +32,6 @@ Matrix *standard_constraints(Matrix *C, unsigned nparam, int *rows_p, Matrix *M; Matrix *H, *U, *Q; - for (i = 0; i < C->NbRows; ++i) - assert(value_one_p(C->p[i][0])); - /* move constraints only involving parameters down * and move unit vectors (if there are any) to the right place. */ @@ -83,7 +80,179 @@ Matrix *standard_constraints(Matrix *C, unsigned nparam, int *rows_p, if (First_Non_Zero(H->p[i], i) != -1) rows++; } - *rows_p = rows; + if (rows_p) + *rows_p = rows; return H; } + +/* + * Transform the constraints of P, which may have existentially + * quantified variables to "standard form". + * + * We need to make sure that the new variables are affine + * combinations of only the original variables (and not + * of any of the existentially quantified variables). + * We therefore first apply standard_constraints on + * the constraints involving only the variables and + * then apply the same procedure on the remaining + * constraints with the variables as parameters. + * + * Perhaps we are overly careful here. Moving the constraints + * involving existentially quantified variables down should prevent + * the regular variables from being replaced by linear combinations + * of existential variables. + */ +static Matrix *exist_standard_constraints(Polyhedron *P, unsigned nvar) +{ + int i, i1, i2; + unsigned exist = P->Dimension - nvar; + Matrix *M = Matrix_Alloc(P->NbConstraints, 2+P->Dimension); + Matrix *M1 = Matrix_Alloc(P->NbConstraints, 2+nvar); + Matrix *M2 = Matrix_Alloc(P->NbConstraints, 2+P->Dimension); + Matrix *S1, *S2; + Value c; + + for (i = 0; i < P->NbConstraints; ++i) { + value_assign(M->p[i][0], P->Constraint[i][0]); + Vector_Copy(P->Constraint[i]+1, M->p[i]+1+exist, nvar); + Vector_Copy(P->Constraint[i]+1+nvar, M->p[i]+1, exist); + value_assign(M->p[i][1+P->Dimension], P->Constraint[i][1+P->Dimension]); + } + Gauss(M, P->NbEq, P->Dimension+1); + for (i = i1 = i2 = 0; i < M->NbRows; ++i) { + if (First_Non_Zero(M->p[i]+1, exist) == -1) { + value_assign(M1->p[i1][0], M->p[i][0]); + Vector_Copy(M->p[i]+1+exist, M1->p[i1]+1, nvar+1); + ++i1; + } else { + Vector_Copy(M->p[i], M2->p[i2], M->NbColumns); + ++i2; + } + } + M1->NbRows = i1; + M2->NbRows = i2; + Matrix_Free(M); + + S1 = standard_constraints(M1, 0, NULL, NULL); + S2 = standard_constraints(M2, nvar, NULL, NULL); + + M = Matrix_Alloc(P->NbConstraints, 2+P->Dimension); + for (i = 0; i < M1->NbRows; ++i) { + value_assign(M->p[i][0], M1->p[i][0]); + Vector_Copy(S1->p[i], M->p[i]+1, nvar); + value_assign(M->p[i][1+P->Dimension], M1->p[i][1+nvar]); + } + for (i = 0; i < M2->NbRows; ++i) { + value_assign(M->p[M1->NbRows+i][0], M2->p[i][0]); + Vector_Copy(M2->p[i]+1+exist, M->p[M1->NbRows+i]+1, nvar); + Vector_Copy(S2->p[i], M->p[M1->NbRows+i]+1+nvar, exist); + value_assign(M->p[M1->NbRows+i][1+P->Dimension], + M2->p[i][1+P->Dimension]); + } + + /* Make sure "standard" constraints for existential variables + * have only non-positive coefficients for the variables. + */ + assert(M2->NbRows >= exist); + value_init(c); + for (i = 0; i < exist; ++i) { + assert(value_pos_p(S2->p[i][i])); + for (i1 = 0; i1 < nvar; ++i1) { + if (value_negz_p(M->p[M1->NbRows+i][1+i1])) + continue; + mpz_cdiv_q(c, M->p[M1->NbRows+i][1+i1], S2->p[i][i]); + value_oppose(c, c); + for (i2 = i; i2 < M2->NbRows; ++i2) + value_addmul(M->p[M1->NbRows+i2][1+i1], c, S2->p[i2][i]); + } + } + value_clear(c); + + Matrix_Free(M1); + Matrix_Free(M2); + + return M; +} + +/* + * For those variables x that may attain negative values, + * compute a translations t such that x' = x + t is sure + * to attain only non-negative values, i.e., if M is + * + * A x + b >= 0 + * + * then the x' in + * + * A x' + b - A t >= 0 + * + * are sure to attain only non-negative values. + */ +static Vector *compute_shifts(Matrix *M) +{ + int c, r; + unsigned dim = M->NbColumns-2; + Vector *shifts = Vector_Alloc(M->NbColumns); + for (c = 0, r = 0; c < dim; ++c) { + for ( ; r < M->NbRows; ++r) { + if (value_notzero_p(M->p[r][1+c])) + break; + } + assert(r < M->NbRows); + assert(value_pos_p(M->p[r][1+c])); + if (c > 0) + Inner_Product(M->p[r]+1, shifts->p+1, c, &shifts->p[1+c]); + value_subtract(shifts->p[1+c], M->p[r][1+dim], shifts->p[1+c]); + if (value_pos_p(shifts->p[1+c])) + mpz_cdiv_q(shifts->p[1+c], shifts->p[1+c], M->p[r][1+c]); + else + value_set_si(shifts->p[1+c], 0); + } + return shifts; +} + +/* + * Given a union of polyhedra, with possibly existentially + * quantified variables, but no parameters, apply a unimodular + * transformation with constant translation to the variables, + * and a unimodular transformation with a translation possibly + * depending on the variables to the existentially quantified variables, + * such that all variables and all existentially quantified variables + * attain only non-negative values. + */ +Polyhedron *skew_to_positive_orthant(Polyhedron *D, unsigned nvar, + unsigned MaxRays) +{ + Polyhedron *P; + Polyhedron *R = NULL; + Polyhedron **next = &R; + + for (P = D; P; P = P->next) { + Matrix *M, *C; + Vector *shifts, *offsets; + int i; + + M = exist_standard_constraints(P, nvar); + shifts = compute_shifts(M); + + offsets = Vector_Alloc(P->NbConstraints); + Matrix_Vector_Product(M, shifts->p, offsets->p); + Vector_Free(shifts); + + C = Matrix_Alloc(P->NbConstraints, P->Dimension+2); + for (i = 0; i < P->NbConstraints; ++i) { + value_assign(C->p[i][0], M->p[i][0]); + Vector_Copy(M->p[i]+1, C->p[i]+1, P->Dimension); + value_subtract(C->p[i][1+P->Dimension], + M->p[i][1+P->Dimension], offsets->p[i]); + } + Vector_Free(offsets); + Matrix_Free(M); + + *next = Constraints2Polyhedron(C, MaxRays); + Matrix_Free(C); + + next = &(*next)->next; + } + return R; +} diff --git a/normalization.h b/normalization.h index b656c93..64c2ccc 100644 --- a/normalization.h +++ b/normalization.h @@ -6,6 +6,8 @@ extern "C" { Matrix *standard_constraints(Matrix *C, unsigned nparam, int *rows_p, Matrix **T); +Polyhedron *skew_to_positive_orthant(Polyhedron *D, unsigned nvar, + unsigned MaxRays); #if defined(__cplusplus) } diff --git a/version.c b/version.c index 670fda9..7f655a3 100644 --- a/version.c +++ b/version.c @@ -49,6 +49,11 @@ const char *barvinok_version(void) #else " -ZSOLVE" #endif +#ifdef USE_PARKER + " +PARKER" +#else + " -PARKER" +#endif "\n" ; } -- 2.11.4.GIT