diff --git a/prism/.classpath b/prism/.classpath
index 1bef3aa5a..07828a2a9 100644
--- a/prism/.classpath
+++ b/prism/.classpath
@@ -12,5 +12,6 @@
+
diff --git a/prism/Makefile b/prism/Makefile
index 4a8523861..a2836e8ba 100644
--- a/prism/Makefile
+++ b/prism/Makefile
@@ -298,11 +298,46 @@ MAKE_DIRS = dd jdd odd dv prism mtbdd sparse hybrid parser settings userinterfac
EXT_PACKAGES = lpsolve55 lp_solve_5.5_java
+EXACT_LP_PACKAGES = glpk glpk_java
+
+EXACT_LP_REQUIREMENTS_INSTALLED = 1
+
+# Find glpk_java requirements
+DETECT_SWIG = $(shell which swig 2> /dev/null)
+DETECT_MVN = $(shell which mvn 2> /dev/null)
+DETECT_LIBTOOL = $(shell which libtool 2> /dev/null)
+DETECT_ACLOCAL = $(shell which aclocal 2> /dev/null)
+
+ifeq ($(DETECT_SWIG),)
+ EXACT_LP_REQUIREMENTS_INSTALLED = 0
+endif
+ifeq ($(DETECT_MVN),)
+ EXACT_LP_REQUIREMENTS_INSTALLED = 0
+endif
+ifeq ($(DETECT_LIBTOOL),)
+ EXACT_LP_REQUIREMENTS_INSTALLED = 0
+endif
+ifeq ($(DETECT_ACLOCAL),)
+ EXACT_LP_REQUIREMENTS_INSTALLED = 0
+endif
+
.PHONY: clean javadoc tests
default: all
-all: cuddpackage extpackages prism
+all: cuddpackage extpackages exact_lp_requirements exact_lp prism
+
+exact_lp_requirements: checks
+ cp ext/glpk_java/jar/glpk-java.jar lib
+ @if [ "$(DETECT_SWIG)" != "" ]; then echo "Found SWIG in $(DETECT_SWIG)"; \
+ else echo "Warning: SWIG not found. Required for GLPK java bindings."; fi
+ @if [ "$(DETECT_MVN)" != "" ]; then echo "Found MAVEN in $(DETECT_MVN)"; \
+ else echo "Warning: MAVEN not found. Required for GLPK java bindings."; fi
+ @if [ "$(DETECT_LIBTOOL)" != "" ]; then echo "Found LIBTOOL in $(DETECT_LIBTOOL)"; \
+ else echo "Warning: LIBTOOL not found. Required for GLPK java bindings."; fi
+ @if [ "$(DETECT_ACLOCAL)" != "" ]; then echo "Found ACLOCAL in $(DETECT_ACLOCAL)"; \
+ else echo "Warning: ACLOCAL not found. Required for both GLPK and GLPK java bindings."; fi
+ @sed -i '/EXACT_LP_REQUIREMENTS_INSTALLED/c\ public final int EXACT_LP_REQUIREMENTS_INSTALLED = $(EXACT_LP_REQUIREMENTS_INSTALLED);' src/prism/PrismSettings.java;
cuddpackage: checks
@if [ "$(CUDD_DIR)" = "" ]; then echo "Error: Cannot find CUDD"; exit 1; fi
@@ -360,6 +395,25 @@ extpackages: checks
) || exit 1; \
done
+exact_lp: checks
+ @(if [ "$(EXACT_LP_REQUIREMENTS_INSTALLED)" = 1 ]; then \
+ for ext in $(EXACT_LP_PACKAGES); do \
+ echo Making $$ext; \
+ (cd ext/$$ext && \
+ $(MAKE) \
+ OSTYPE="$(OSTYPE)" \
+ ARCH="$(ARCH)" \
+ LIBPREFIX="$(LIBPREFIX)" \
+ LIBSUFFIX="$(LIBSUFFIX)" \
+ LIBMATH="$(LIBMATH)" \
+ BINDISTSUFFIX="$(BINDISTSUFFIX)" \
+ JAVA_DIR="$(JAVA_DIR)" \
+ JAVA_JNI_H_DIR="$(JAVA_JNI_H_DIR)" \
+ JAVA_JNI_MD_H_DIR="$(JAVA_JNI_MD_H_DIR)" \
+ ) || exit 1; \
+ done \
+ fi)
+
prism: checks make_dirs bin_scripts
make_dirs:
@@ -543,7 +597,7 @@ clean_cudd:
@(cd $(CUDD_DIR) && $(MAKE) distclean)
clean_ext:
- @(for ext in $(EXT_PACKAGES); do \
+ @(for ext in $(EXT_PACKAGES) $(EXACT_LP_PACKAGES); do \
echo Cleaning $$ext ...; \
(cd ext/$$ext && \
$(MAKE) -s LIBPREFIX="$(LIBPREFIX)" LIBSUFFIX="$(LIBSUFFIX)" clean) \
diff --git a/prism/ext/glpk/Makefile b/prism/ext/glpk/Makefile
new file mode 100644
index 000000000..9325fda35
--- /dev/null
+++ b/prism/ext/glpk/Makefile
@@ -0,0 +1,57 @@
+################################################
+# NB: This Makefile is designed to be called #
+# from the main PRISM Makefile. It won't #
+# work on its own because it needs #
+# various options to be passed in #
+################################################
+
+default: all
+
+all: checks ../../lib/glpk
+
+glpkvers = "4.60"
+
+IS_SRC_CREATED = $(shell find . -name "src" 2> /dev/null)
+ALREADY_DOWNLOADED = $(shell ls src | grep glpk-$(glpkvers).tar.gz 2> /dev/null)
+ALREADY_INSTALLED_LIB = $(shell find src/glpk/src/.libs -name $(LIBPREFIX)glpk$(LIBSUFFIX).40 2> /dev/null)
+ALREADY_INSTALLED = 1
+
+ifeq ($(ALREADY_INSTALLED_LIB),)
+ ALREADY_INSTALLED = 0
+endif
+
+# Try and prevent accidental makes (i.e. called manually, not from top-level Makefile)
+checks:
+ @if [ "$(LIBSUFFIX)" = "" ]; then \
+ (echo "Error: This Makefile is designed to be called from the main PRISM Makefile"; exit 1) \
+ fi;
+
+../../lib/glpk:
+ @(if [ "$(OSTYPE)" != "cygwin" ]; then \
+ (if [ "$(IS_SRC_CREATED)" = "" ]; then \
+ mkdir src; fi) && \
+ (cd src/ && \
+ (if [ "$(ALREADY_DOWNLOADED)" = "" ]; then \
+ wget http://ftp.gnu.org/gnu/glpk/glpk-$(glpkvers).tar.gz && \
+ tar -xzf glpk-$(glpkvers).tar.gz && \
+ mv glpk-$(glpkvers) glpk; \
+ fi) && \
+ cd glpk/ && \
+ (if [ $(ALREADY_INSTALLED) = 0 ]; then \
+ patch -p1 < ../../glpk_exact_solutions.patch && \
+ ./autogen.sh > /dev/null && \
+ ./configure --with-gmp > /dev/null && \
+ echo Building GLPK. This may take several minutes... && \
+ make > /dev/null; \
+ fi) \
+ ) && cp --remove-destination src/glpk/src/.libs/$(LIBPREFIX)glpk$(LIBSUFFIX).40 ../../lib/; \
+ fi)
+
+clean: checks
+ rm -f ../../lib/$(LIBPREFIX)glpk$(LIBSUFFIX)
+ rm -f ../../lib/$(LIBPREFIX)glpk$(LIBSUFFIX).40
+ rm -rf src/*
+
+celan: clean
+
+#################################################
diff --git a/prism/ext/glpk/glpk_exact_solutions.patch b/prism/ext/glpk/glpk_exact_solutions.patch
new file mode 100644
index 000000000..7da40ddd3
--- /dev/null
+++ b/prism/ext/glpk/glpk_exact_solutions.patch
@@ -0,0 +1,320 @@
+diff --git a/src/glpapi06.c b/src/glpapi06.c
+index f526e72..e67f5e1 100644
+--- a/src/glpapi06.c
++++ b/src/glpapi06.c
+@@ -209,7 +209,8 @@ up: col->stat = GLP_NU, col->prim = col->ub;
+ }
+ if (parm->msg_lev >= GLP_MSG_ALL && parm->out_dly == 0)
+ { if (P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)
+- xprintf("OPTIMAL SOLUTION FOUND\n");
++ return;
++ //xprintf("OPTIMAL SOLUTION FOUND\n");
+ else if (P->pbs_stat == GLP_NOFEAS)
+ xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
+ else if (parm->meth == GLP_PRIMAL)
+@@ -501,6 +502,8 @@ void glp_init_smcp(glp_smcp *parm)
+ parm->out_frq = 500;
+ parm->out_dly = 0;
+ parm->presolve = GLP_OFF;
++ parm->precision = 64;
++ parm->num_states = 10;
+ return;
+ }
+
+@@ -783,6 +786,90 @@ double glp_get_col_dual(glp_prob *lp, int j)
+ /***********************************************************************
+ * NAME
+ *
++* glp_get_probs - retrieve fixed precision solutions array
++*
++* SYNOPSIS
++*
++* unsigned char **glp_get_probs(glp_prob *lp);
++*
++* RETURNS
++*
++* The routine glp_get_probs returns the array of fixed precision
++* solutions. */
++
++unsigned char **glp_get_probs(glp_prob *lp)
++{
++ unsigned char **probs;
++ probs = lp->probs;
++ return probs;
++}
++
++/***********************************************************************
++* NAME
++*
++* glp_get_exps - retrieve exponent positions for each solution
++*
++* SYNOPSIS
++*
++* int *glp_get_exps(glp_prob *lp);
++*
++* RETURNS
++*
++* The routine glp_get_exps returns the positions of the exponents
++* for each solution in probs array. */
++
++int *glp_get_exps(glp_prob *lp)
++{
++ int *exps;
++ exps = lp->exps;
++ return exps;
++}
++
++/***********************************************************************
++* NAME
++*
++* glp_get_nums - retrieve numerators of exact rational solutions
++*
++* SYNOPSIS
++*
++* unsigned char *glp_get_nums(glp_prob *lp);
++*
++* RETURNS
++*
++* The routine glp_get_nums returns the numerators for
++* each exact rational arithmetic solution. */
++
++unsigned char **glp_get_nums(glp_prob *lp)
++{
++ unsigned char **nums;
++ nums = lp->nums;
++ return nums;
++}
++
++/***********************************************************************
++* NAME
++*
++* glp_get_dens - retrieve denominators of exact rational solutions
++*
++* SYNOPSIS
++*
++* unsigned char **glp_get_dens(glp_prob *lp);
++*
++* RETURNS
++*
++* The routine glp_get_dens returns the denominators for
++* each exact rational arithmetic solution. */
++
++unsigned char **glp_get_dens(glp_prob *lp)
++{
++ unsigned char **dens;
++ dens = lp->dens;
++ return dens;
++}
++
++/***********************************************************************
++* NAME
++*
+ * glp_get_unbnd_ray - determine variable causing unboundedness
+ *
+ * SYNOPSIS
+diff --git a/src/glpapi07.c b/src/glpapi07.c
+index 89a06bf..3fe33a2 100644
+--- a/src/glpapi07.c
++++ b/src/glpapi07.c
+@@ -125,7 +125,8 @@ static void set_d_eps(mpq_t x, double val)
+ if (s < 0) mpq_neg(x, x);
+ /* check that the desired tolerance has been attained */
+ xassert(fabs(val - mpq_get_d(x)) <= eps * (1.0 + fabs(val)));
+-done: return;
++done:
++ return;
+ }
+
+ static void load_data(SSX *ssx, glp_prob *lp)
+@@ -159,8 +160,10 @@ static void load_data(SSX *ssx, glp_prob *lp)
+ default: xassert(type != type);
+ }
+ ssx->type[k] = type;
+- set_d_eps(ssx->lb[k], lb);
+- set_d_eps(ssx->ub[k], ub);
++ //set_d_eps(ssx->lb[k], lb);
++ //set_d_eps(ssx->ub[k], ub);
++ mpq_set_d(ssx->lb[k], lb);
++ mpq_set_d(ssx->ub[k], ub);
+ }
+ /* optimization direction */
+ switch (lp->dir)
+@@ -176,7 +179,8 @@ static void load_data(SSX *ssx, glp_prob *lp)
+ coef = 0.0;
+ else
+ coef = lp->col[k-m]->coef;
+- set_d_eps(ssx->coef[k], coef);
++ //set_d_eps(ssx->coef[k], coef);
++ mpq_set_d(ssx->coef[k], coef);
+ }
+ /* constraint coefficients */
+ ind = xcalloc(1+m, sizeof(int));
+@@ -188,7 +192,8 @@ static void load_data(SSX *ssx, glp_prob *lp)
+ for (k = 1; k <= len; k++)
+ { loc++;
+ ssx->A_ind[loc] = ind[k];
+- set_d_eps(ssx->A_val[loc], val[k]);
++ //set_d_eps(ssx->A_val[loc], val[k]);
++ mpq_set_d(ssx->A_val[loc], val[k]);
+ }
+ }
+ xassert(loc == nnz);
+@@ -262,8 +267,14 @@ int glp_exact(glp_prob *lp, const glp_smcp *parm)
+ int m = lp->m;
+ int n = lp->n;
+ int nnz = lp->nnz;
+- int i, j, k, type, pst, dst, ret, stat;
++ int i, j, k, type, pst, dst, ret, stat, actual_state, index;
+ double lb, ub, prim, dual, sum;
++ char *exact_prob;
++ /* fixed precision initial parameters */
++ mpf_t fix_prob;
++ mpf_init(fix_prob);
++ mpf_set_prec(fix_prob, parm->precision);
++ mp_exp_t exponent;
+ if (parm == NULL)
+ parm = &_parm, glp_init_smcp((glp_smcp *)parm);
+ /* check control parameters */
+@@ -273,11 +284,26 @@ int glp_exact(glp_prob *lp, const glp_smcp *parm)
+ if (parm->tm_lim < 0)
+ xerror("glp_exact: tm_lim = %d; invalid parameter\n",
+ parm->tm_lim);
++ if (parm->num_states <= 0)
++ xerror("glp_exact: num_states = %d; invalid parameter\n",
++ parm->num_states);
++ if (parm->precision <= 0)
++ xerror("glp_exact: precision = %d; invalid parameter\n",
++ parm->precision);
+ /* the problem must have at least one row and one column */
+ if (!(m > 0 && n > 0))
+ { xprintf("glp_exact: problem has no rows/columns\n");
+ return GLP_EFAIL;
+ }
++ actual_state = 0;
++ /* alloc memory for storing exact probabilities */
++ lp->nums = (unsigned char **) malloc(parm->num_states * sizeof(unsigned char*));
++ lp->dens = (unsigned char **) malloc(parm->num_states * sizeof(unsigned char*));
++ /* alloc memory for storing fixed probabilities */
++ lp->exps = (int *) malloc(parm->num_states * sizeof(int));
++ lp->probs = (unsigned char **) malloc(parm->num_states * sizeof(unsigned char*));
++ for (index = 0; index < parm->num_states; index++)
++ lp->probs[index] = (unsigned char *) malloc(mpf_get_prec(fix_prob) * sizeof(unsigned char));
+ #if 1
+ /* basic solution is currently undefined */
+ lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
+@@ -303,14 +329,14 @@ int glp_exact(glp_prob *lp, const glp_smcp *parm)
+ }
+ }
+ /* create the simplex solver workspace */
+- xprintf("glp_exact: %d rows, %d columns, %d non-zeros\n",
+- m, n, nnz);
++ //xprintf("glp_exact: %d rows, %d columns, %d non-zeros\n",
++ // m, n, nnz);
+ #ifdef HAVE_GMP
+- xprintf("GNU MP bignum library is being used\n");
++ //xprintf("GNU MP bignum library is being used\n");
+ #else
+- xprintf("GLPK bignum module is being used\n");
+- xprintf("(Consider installing GNU MP to attain a much better perf"
+- "ormance.)\n");
++ //xprintf("GLPK bignum module is being used\n");
++ //xprintf("(Consider installing GNU MP to attain a much better perf"
++ //"ormance.)\n");
+ #endif
+ ssx = ssx_create(m, n, nnz);
+ /* load LP problem data into the workspace */
+@@ -433,6 +459,14 @@ int glp_exact(glp_prob *lp, const glp_smcp *parm)
+ xassert(ssx != ssx);
+ }
+ dual = mpq_get_d(ssx->cbar[j]);
++ /* exact rational results */
++ lp->nums[actual_state] = mpz_get_str(NULL, 10, mpq_numref(ssx->cbar[j]));
++ lp->dens[actual_state] = mpz_get_str(NULL, 10, mpq_denref(ssx->cbar[j]));
++ /* fixed results */
++ mpf_set_q(fix_prob, ssx->cbar[j]);
++ lp->probs[actual_state] = mpf_get_str(NULL, &exponent, 10, 0, fix_prob);
++ lp->exps[actual_state] = exponent;
++ actual_state++;
+ }
+ if (k <= m)
+ { glp_set_row_stat(lp, k, stat);
+@@ -456,5 +490,3 @@ done: /* delete the simplex solver workspace */
+ /* return to the application program */
+ return ret;
+ }
+-
+-/* eof */
+diff --git a/src/glpk.h b/src/glpk.h
+index 137801c..654a232 100644
+--- a/src/glpk.h
++++ b/src/glpk.h
+@@ -137,6 +137,8 @@ typedef struct
+ int out_frq; /* spx.out_frq */
+ int out_dly; /* spx.out_dly (milliseconds) */
+ int presolve; /* enable/disable using LP presolver */
++ int precision; /* current precision */
++ int num_states; /* number of states */
+ double foo_bar[36]; /* (reserved) */
+ } glp_smcp;
+
+@@ -496,6 +498,18 @@ double glp_get_col_prim(glp_prob *P, int j);
+ double glp_get_col_dual(glp_prob *P, int j);
+ /* retrieve column dual value (basic solution) */
+
++unsigned char **glp_get_probs(glp_prob *lp);
++/* retrieve fixed precision solutions array */
++
++int *glp_get_exps(glp_prob *lp);
++/* retrieve exponent positions for each solution */
++
++unsigned char **glp_get_nums(glp_prob *lp);
++/* retrieve the numerators for each exact rational solution */
++
++unsigned char **glp_get_dens(glp_prob *lp);
++/* retrieve the denominators for each exact rational solution */
++
+ int glp_get_unbnd_ray(glp_prob *P);
+ /* determine variable causing unboundedness */
+
+diff --git a/src/glpssx02.c b/src/glpssx02.c
+index 4b3ea97..4a810a4 100644
+--- a/src/glpssx02.c
++++ b/src/glpssx02.c
+@@ -31,9 +31,9 @@ static void show_progress(SSX *ssx, int phase)
+ int i, def = 0;
+ for (i = 1; i <= ssx->m; i++)
+ if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
+- xprintf("%s%6d: %s = %22.15g (%d)\n", phase == 1 ? " " : "*",
+- ssx->it_cnt, phase == 1 ? "infsum" : "objval",
+- mpq_get_d(ssx->bbar[0]), def);
++ //xprintf("%s%6d: %s = %22.15g (%d)\n", phase == 1 ? " " : "*",
++ //ssx->it_cnt, phase == 1 ? "infsum" : "objval",
++ //mpq_get_d(ssx->bbar[0]), def);
+ #if 0
+ ssx->tm_lag = utime();
+ #else
+@@ -446,7 +446,7 @@ skip: /* compute simplex multipliers */
+ ret = ssx_phase_II(ssx);
+ switch (ret)
+ { case 0:
+- xprintf("OPTIMAL SOLUTION FOUND\n");
++ //xprintf("OPTIMAL SOLUTION FOUND\n");
+ ret = 0;
+ break;
+ case 1:
+diff --git a/src/prob.h b/src/prob.h
+index 4c1d34e..fc675d9 100644
+--- a/src/prob.h
++++ b/src/prob.h
+@@ -143,6 +143,11 @@ struct glp_prob
+ GLP_NOFEAS - no integer solution exists */
+ double mip_obj;
+ /* objective function value */
++ /* for exact solver values storage */
++ unsigned char **probs; /* values of the probabilities */
++ int *exps; /* exponents of the probabilities */
++ unsigned char **nums; /* rational numerators */
++ unsigned char **dens; /* rational denominators */
+ };
+
+ struct GLPROW
diff --git a/prism/ext/glpk_java/Makefile b/prism/ext/glpk_java/Makefile
new file mode 100644
index 000000000..cbeac64a5
--- /dev/null
+++ b/prism/ext/glpk_java/Makefile
@@ -0,0 +1,68 @@
+################################################
+# NB: This Makefile is designed to be called #
+# from the main PRISM Makefile. It won't #
+# work on its own because it needs #
+# various options to be passed in #
+################################################
+
+default: all
+
+all: checks ../../lib/glpk_java
+
+LDPATH = $(shell cd ../glpk/src/glpk/src/.libs; pwd)
+CFLAGSPATH = $(shell cd ../glpk/src/glpk/src; pwd)
+SWIGFLAGSPATH = $(shell cd ../glpk/src/glpk/src; pwd)
+
+glpkjavavers = "1.7.0"
+
+IS_SRC_CREATED = $(shell find . -name "src" 2> /dev/null)
+ALREADY_DOWNLOADED = $(shell ls src | grep libglpk-java-$(glpkjavavers).tar.gz 2> /dev/null)
+ALREADY_INSTALLED_JAR = $(shell find src/glpk_java/swig -name glpk-java.jar 2> /dev/null)
+ALREADY_INSTALLED_LIB = $(shell find src/glpk_java/swig/.libs -name $(LIBPREFIX)glpk_java$(LIBSUFFIX) 2> /dev/null)
+ALREADY_INSTALLED = 1
+
+ifeq ($(ALREADY_INSTALLED_JAR),)
+ ALREADY_INSTALLED = 0
+endif
+ifeq ($(ALREADY_INSTALLED_LIB),)
+ ALREADY_INSTALLED = 0
+endif
+
+# Try and prevent accidental makes (i.e. called manually, not from top-level Makefile)
+checks:
+ @if [ "$(LIBSUFFIX)" = "" ]; then \
+ (echo "Error: This Makefile is designed to be called from the main PRISM Makefile"; exit 1) \
+ fi;
+
+../../lib/glpk_java:
+ @(if [ "$(OSTYPE)" != "cygwin" ]; then \
+ (if [ "$(IS_SRC_CREATED)" = "" ]; then \
+ mkdir src ;fi) && \
+ export JAVA_HOME=$(JAVA_DIR) && \
+ (cd src/ && \
+ (if [ "$(ALREADY_DOWNLOADED)" = "" ]; then \
+ wget http://download.sourceforge.net/project/glpk-java/glpk-java/glpk-java-$(glpkjavavers)/libglpk-java-$(glpkjavavers).tar.gz && \
+ tar -xzf libglpk-java-$(glpkjavavers).tar.gz && \
+ mv libglpk-java-$(glpkjavavers) glpk_java; \
+ fi) && \
+ cd glpk_java/ && \
+ (if [ $(ALREADY_INSTALLED) = 0 ]; then \
+ patch -p1 < ../../libglpk_java_exact.patch && \
+ export LD_LIBRARY_PATH=$(LD_LIBRARY_PATH):$(LDPATH) && \
+ export CFLAGS=-I$(CFLAGSPATH) && \
+ export SWIGFLAGS=-I$(SWIGFLAGSPATH) && \
+ ./autogen.sh > /dev/null && \
+ ./configure > /dev/null && \
+ make > /dev/null; \
+ fi) \
+ ) && cp --remove-destination src/glpk_java/swig/glpk-java.jar ../../lib/; \
+ cp --remove-destination src/glpk_java/swig/.libs/$(LIBPREFIX)glpk_java$(LIBSUFFIX) ../../lib/; \
+ fi)
+
+clean: checks
+ rm -f ../../lib/$(LIBPREFIX)glpk_java$(LIBSUFFIX)
+ rm -rf src/*
+
+celan: clean
+
+#################################################
diff --git a/prism/ext/glpk_java/jar/glpk-java.jar b/prism/ext/glpk_java/jar/glpk-java.jar
new file mode 100644
index 000000000..1dbb9e97b
Binary files /dev/null and b/prism/ext/glpk_java/jar/glpk-java.jar differ
diff --git a/prism/ext/glpk_java/libglpk_java_exact.patch b/prism/ext/glpk_java/libglpk_java_exact.patch
new file mode 100644
index 000000000..eace2c298
--- /dev/null
+++ b/prism/ext/glpk_java/libglpk_java_exact.patch
@@ -0,0 +1,54 @@
+diff --git a/swig/glpk.i b/swig/glpk.i
+index bf7d91f..e5ba65f 100644
+--- a/swig/glpk.i
++++ b/swig/glpk.i
+@@ -334,6 +334,8 @@ glp_java_vertex_data *glp_java_vertex_get_data(
+ %include "glpk_java_arrays.i"
+ %glp_array_functions(int, intArray)
+ %glp_array_functions(double, doubleArray)
++%glp_array_functions(unsigned char, uint8Array)
++%glp_array_functions(unsigned char *, uint8ArrayArray)
+
+ // Add handling for String arrays in glp_main
+ %include "various.i"
+diff --git a/swig/glpk_java_arrays.i b/swig/glpk_java_arrays.i
+index 214cec7..e4ecaaa 100644
+--- a/swig/glpk_java_arrays.i
++++ b/swig/glpk_java_arrays.i
+@@ -29,11 +29,19 @@ static TYPE NAME##_getitem(TYPE *ary, int index) {
+ return (TYPE) 0;
+ }
+ }
++
+ static void NAME##_setitem(TYPE *ary, int index, TYPE value) {
+ if (ary != NULL) {
+ ary[index] = value;
+ }
+ }
++
++static char *NAME##_toString(TYPE *ary) {
++ char *str;
++ str = (char *)ary;
++ return str;
++}
++
+ %}
+
+ %javamethodmodifiers new_##NAME(int nelements) "
+@@ -85,6 +93,16 @@ TYPE NAME##_getitem(TYPE *ary, int index);
+ public";
+ void NAME##_setitem(TYPE *ary, int index, TYPE value);
+
++%javamethodmodifiers NAME##_toString(TYPE* ary) "
++/**
++ * Sets the value of an element of an array of TYPE.
++ * BEWARE: The validity of the index is not checked.
++ *
++ * @param ary array
++ */
++public";
++char *NAME##_toString(TYPE *ary);
++
+ %enddef
+ /* Old Swig versions require a LF after enddef */
+
diff --git a/prism/src/explicit/DTMCModelChecker.java b/prism/src/explicit/DTMCModelChecker.java
index e6c5b3059..6dc532708 100644
--- a/prism/src/explicit/DTMCModelChecker.java
+++ b/prism/src/explicit/DTMCModelChecker.java
@@ -29,9 +29,18 @@
import java.io.File;
import java.util.Arrays;
import java.util.BitSet;
+import java.util.HashMap;
+import java.util.Iterator;
import java.util.List;
import java.util.Map;
+import org.gnu.glpk.GLPK;
+import org.gnu.glpk.GLPKConstants;
+import org.gnu.glpk.SWIGTYPE_p_double;
+import org.gnu.glpk.SWIGTYPE_p_int;
+import org.gnu.glpk.glp_prob;
+import org.gnu.glpk.glp_smcp;
+
import parser.VarList;
import parser.ast.Declaration;
import parser.ast.DeclarationIntUnbounded;
@@ -519,7 +528,8 @@ public ModelCheckerResult computeReachProbs(DTMC dtmc, BitSet remain, BitSet tar
LinEqMethod linEqMethod = this.linEqMethod;
// Switch to a supported method, if necessary
- if (!(linEqMethod == LinEqMethod.POWER || linEqMethod == LinEqMethod.GAUSS_SEIDEL)) {
+ if (!(linEqMethod == LinEqMethod.POWER || linEqMethod == LinEqMethod.GAUSS_SEIDEL
+ || linEqMethod == LinEqMethod.EXACT_LINEAR_PROGRAMMING || settings.EXACT_LP_REQUIREMENTS_INSTALLED == 0)) {
linEqMethod = LinEqMethod.GAUSS_SEIDEL;
mainLog.printWarning("Switching to linear equation solution method \"" + linEqMethod.fullName() + "\"");
}
@@ -598,6 +608,9 @@ public ModelCheckerResult computeReachProbs(DTMC dtmc, BitSet remain, BitSet tar
case GAUSS_SEIDEL:
res = computeReachProbsGaussSeidel(dtmc, no, yes, init, known);
break;
+ case EXACT_LINEAR_PROGRAMMING:
+ res = computeReachProbsExactLinearProg(dtmc, no, yes, init, known);
+ break;
default:
throw new PrismException("Unknown linear equation solution method " + linEqMethod.fullName());
}
@@ -1040,6 +1053,265 @@ protected ModelCheckerResult computeReachProbsGaussSeidel(DTMC dtmc, BitSet no,
return res;
}
+ /**
+ * Compute arbitrary precision reachability probabilities using linear programming.
+ * @param dtmc The DTMC
+ * @param no Probability 0 states
+ * @param yes Probability 1 states
+ * @param init Optionally, an initial solution vector (will be overwritten)
+ * @param known Optionally, a set of states for which the exact answer is known
+ * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values.
+ */
+ protected ModelCheckerResult computeReachProbsExactLinearProg(DTMC dtmc, BitSet no, BitSet yes, double init[], BitSet known) throws PrismException
+ {
+ ModelCheckerResult res;
+ BitSet unknown;
+ int i, nStates, maybeSize, countMaybeTrans, ret;
+ double sol, soln[], initVal, matAEntry, bEntry;
+ long timer, matATimer, exactSolverTimer;
+ SWIGTYPE_p_int ind;
+ SWIGTYPE_p_double val;
+ glp_smcp parm;
+ glp_prob lp;
+
+ // Start value iteration
+ timer = System.currentTimeMillis();
+ mainLog.println("Starting exact linear programming...");
+
+ // Store num states
+ nStates = dtmc.getNumStates();
+
+ // Determine the set of maybe states (S?)
+ unknown = new BitSet();
+ unknown.set(0, nStates);
+ unknown.andNot(yes);
+ unknown.andNot(no);
+ if (known != null)
+ unknown.andNot(known);
+
+ // Create solution vector
+ soln = (init == null) ? new double[nStates] : init;
+
+ // Initialise solution vector. Use (where available) the following in order of preference:
+ // (1) exact answer, if already known; (2) 1.0/0.0 if in yes/no; (3) passed in initial value; (4) initVal
+ // where initVal is 0.0 or 1.0, depending on whether we converge from below/above.
+ initVal = (valIterDir == ValIterDir.BELOW) ? 0.0 : 1.0;
+ if (init != null) {
+ if (known != null) {
+ for (i = 0; i < nStates; i++)
+ soln[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : init[i];
+ } else {
+ for (i = 0; i < nStates; i++)
+ soln[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 : init[i];
+ }
+ } else {
+ for (i = 0; i < nStates; i++)
+ soln[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 : initVal;
+ }
+
+ // Store the cardinality of the maybe state
+ maybeSize = unknown.cardinality();
+
+ // If we have states to compute values for
+ if (maybeSize > 0) {
+ // A mapping between maybe states and its positions in the maybe bitset.
+ // This is useful for calculating the `A` matrix later on.
+ // Also, calculate the number of the transitions of all maybe states `countMaybeTrans`
+ countMaybeTrans = 0;
+ HashMap maybeToState = new HashMap();
+ HashMap stateToMaybe = new HashMap();
+ int maybeStateCount = 1;
+ for(int state = unknown.nextSetBit(0); state >= 0; state = unknown.nextSetBit(state + 1)) {
+ maybeToState.put(maybeStateCount, state);
+ stateToMaybe.put(state, maybeStateCount);
+ countMaybeTrans++;
+ maybeStateCount++;
+ }
+
+ // Memory allocation for index and value arrays.
+ // For each transition, GLPK stores the values for each state in the maybe set
+ // and uses these arrays as temporal storage
+ ind = GLPK.new_intArray(maybeSize + 1);
+ val = GLPK.new_doubleArray(maybeSize + 1);
+
+ // Construct max (min) linear programming problem.
+ // First, we put the problem in the canonical form
+ // and then apply the duality definition
+ lp = GLPK.glp_create_prob();
+ mainLog.println("Linear programming problem created");
+
+ // Set the lp constants for the dual simplex method.
+ // columns -> one for each transition
+ // rows -> maybe states
+ GLPK.glp_add_cols(lp, countMaybeTrans);
+ GLPK.glp_add_rows(lp, maybeSize);
+ GLPK.intArray_setitem(ind, 0, 0);
+ GLPK.doubleArray_setitem(val, 0, 0);
+
+ // Constraints for each equation expressed as inequalities.
+ // For the min (max) problem, each equation has -1 as a upper (lower) bound and we compute the max (min).
+ // This was obtained rewriting the problems (6) (7) in their canonical form, and then construct the dual.
+ // The `-1` bounds correspond to the values of the cost vector (remember that we put the problems
+ // in canonical form first).
+ // Also, we set each row (representing a maybe state) as part of the basis defined in (10)
+ for (int row = 1; row <= maybeSize; row++) {
+ // The name of the state
+ GLPK.glp_set_row_name(lp, row, "s" + unknown.nextSetBit(row));
+ GLPK.glp_set_row_bnds(lp, row, GLPKConstants.GLP_UP, 0, -1);
+ GLPK.glp_set_obj_dir(lp, GLPKConstants.GLP_MAX);
+ GLPK.glp_set_row_stat(lp, row, GLPKConstants.GLP_BS);
+ GLPK.intArray_setitem(ind, row, 0);
+ GLPK.doubleArray_setitem(val, row, 0);
+ }
+
+ // The current column of the matrix (transition)
+ int currentTransition = 1;
+ matATimer = System.currentTimeMillis();
+ // Fill out the (A | I) matrix (lp problem equations) of the reachability problem.
+ // In this step, we calculate the values for each entry of the A matrix,
+ // going through all maybe states and their transitions for each choice
+ for (int state = unknown.nextSetBit(0); state >= 0; state = unknown.nextSetBit(state + 1)) {
+ // Index of the i-th equation in the lp's matrix
+ int eqIndex = 1;
+ Iterator> iter = dtmc.getTransitionsIterator(state);
+ bEntry = 0.;
+ boolean pointsToItself = false;
+ // Iterate over the transitions between st and choice
+ while (iter.hasNext()) {
+ matAEntry = 0.;
+ Map.Entry entry = iter.next();
+ int actualKey = entry.getKey();
+ // Skip if the probability is 0
+ if (entry.getValue() == 0)
+ continue;
+ // Calculate the entries (summation) for the b vector (constraint vector)
+ if (yes.get(actualKey))
+ bEntry += entry.getValue();
+ // If it is a maybe state, we calculate the matrix entry as defined on
+ // page 8 of the paper
+ if (unknown.get(actualKey)) {
+ if (actualKey == state) {
+ pointsToItself = true;
+ matAEntry = entry.getValue() - 1;
+ } else {
+ matAEntry = entry.getValue();
+ }
+ // Write the entry of each equation in the array
+ GLPK.intArray_setitem(ind, eqIndex, stateToMaybe.get(actualKey));
+ GLPK.doubleArray_setitem(val, eqIndex, matAEntry);
+ eqIndex++;
+ }
+ }
+ // If the state does not point to itself, then fill out the missing entry
+ if (!pointsToItself) {
+ GLPK.intArray_setitem(ind, eqIndex, stateToMaybe.get(state));
+ GLPK.doubleArray_setitem(val, eqIndex, -1);
+ eqIndex++;
+ }
+
+ // Set the values of the arrays in the full linear programming equation matrix.
+ // Negate the summation calculated for the constraint vector entries
+ GLPK.glp_set_col_bnds(lp, currentTransition, GLPKConstants.GLP_LO, 0, 0);
+ GLPK.glp_set_mat_col(lp, currentTransition, eqIndex - 1, ind, val);
+ GLPK.glp_set_obj_coef(lp, currentTransition, -bEntry);
+
+ currentTransition++;
+ }
+ mainLog.println("Linear programming equations construction time: " +
+ ((double)System.currentTimeMillis() - matATimer)/1000 + " seconds");
+
+ // Solving the linear problem in the linear programming library GLPK.
+ // probs[]: an array for storing the exact probabilities as string in the form `p/q`
+ // `num` and `den` store the temporal `p` and `q` respectively (this could be avoided)
+ String probs[] = new String[nStates];
+
+ // Set up exact parameters
+ parm = new glp_smcp();
+ GLPK.glp_init_smcp(parm);
+ parm.setPrecision(this.getPrecision());
+ parm.setNum_states(nStates);
+
+ // Solve lp problem
+ exactSolverTimer = System.currentTimeMillis();
+ ret = GLPK.glp_exact(lp, parm);
+
+ // Check return value: if `ret` is not zero there has been an error
+ if (ret != 0) {
+ String msg;
+ if (ret == GLPKConstants.GLP_EBADB) {
+ msg = "The initial basis specified in the problem object is invalid: wrong number of basic variables";
+ } else if (ret == GLPKConstants.GLP_ESING) {
+ msg = "The basis matrix corresponding to the initial basis is singular";
+ } else if (ret == GLPKConstants.GLP_EBOUND) {
+ msg = "Some double-bounded variables have incorrect bound";
+ } else {
+ msg = "The linear programming problem could not be solved";
+ }
+ mainLog.println(msg);
+ throw new PrismException(msg);
+ }
+
+ // Calculate and store exact probabilities.
+ // For each maybe state, we extract the fixed precision results from GLPK and store them in `probs[]`.
+ // `sol` stores the exact values as doubles. That implies that the precision is at most 64-bit
+ // However, the exact value is expressed as a string in `probs[]`.
+ sol = GLPK.glp_get_obj_val(lp);
+ String currentResult;
+ int exp;
+ mainLog.println("\n" + this.getPrecision() + "-bit results (non-zero only) for each state of the DTMC:");
+ for (int dtmcState = 0; dtmcState < nStates; dtmcState++) {
+ if (unknown.get(dtmcState)) {
+ sol = GLPK.glp_get_row_dual(lp, stateToMaybe.get(dtmcState));
+ soln[dtmcState] = sol;
+ currentResult = GLPK.uint8Array_toString(GLPK.uint8ArrayArray_getitem(GLPK.glp_get_probs(lp), stateToMaybe.get(dtmcState) - 1));
+ exp = GLPK.intArray_getitem(GLPK.glp_get_exps(lp), stateToMaybe.get(dtmcState) - 1);
+ // We need to format the output since we get the exponent as a position in the string
+ String zeroes = "";
+ for (int nZero = 0; nZero < Math.abs(exp); nZero++)
+ zeroes += "0";
+ currentResult = exp > 0 ? new StringBuilder(currentResult).insert(exp, ".").toString() : "0." + zeroes + currentResult;
+ probs[dtmcState] = currentResult;
+ }
+ // If it is a `yes` state, its probability is 1
+ else if (yes.get(dtmcState)) {
+ soln[dtmcState] = 1.;
+ probs[dtmcState] = "1";
+ }
+ // If it is a `no` state, its probability is 0
+ else if (no.get(dtmcState)) {
+ soln[dtmcState] = 0.;
+ probs[dtmcState] = "0";
+ }
+ // Avoid zero entries
+ if (probs[dtmcState] != "0")
+ mainLog.println("[" + dtmcState + "]: " + probs[dtmcState]);
+ }
+ mainLog.println("\nExact probabilities calculation time: " +
+ ((double)System.currentTimeMillis() - exactSolverTimer)/1000 + " seconds");
+
+ // Free memory
+ GLPK.delete_intArray(ind);
+ GLPK.delete_doubleArray(val);
+ GLPK.delete_intArray(GLPK.glp_get_exps(lp));
+ GLPK.delete_uint8ArrayArray(GLPK.glp_get_probs(lp));
+ GLPK.delete_uint8ArrayArray(GLPK.glp_get_nums(lp));
+ GLPK.delete_uint8ArrayArray(GLPK.glp_get_dens(lp));
+ GLPK.glp_delete_prob(lp);
+
+ }
+
+ // Finished linear programming
+ timer = System.currentTimeMillis() - timer;
+ mainLog.print("Exact linear programming");
+ mainLog.println(" took " + timer / 1000.0 + " seconds.");
+
+ // Return results
+ res = new ModelCheckerResult();
+ res.soln = soln;
+ res.timeTaken = timer / 1000.0;
+ return res;
+ }
+
/**
* Compute bounded reachability probabilities.
* i.e. compute the probability of reaching a state in {@code target} within k steps.
diff --git a/prism/src/explicit/MDPModelChecker.java b/prism/src/explicit/MDPModelChecker.java
index 30879e735..48b1bd354 100644
--- a/prism/src/explicit/MDPModelChecker.java
+++ b/prism/src/explicit/MDPModelChecker.java
@@ -28,9 +28,12 @@
import java.util.Arrays;
import java.util.BitSet;
+import java.util.HashMap;
+import java.util.HashSet;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
+import java.util.Set;
import parser.VarList;
import parser.ast.Declaration;
@@ -42,6 +45,7 @@
import prism.PrismException;
import prism.PrismFileLog;
import prism.PrismLog;
+import prism.PrismSettings;
import prism.PrismUtils;
import strat.MDStrategyArray;
import acceptance.AcceptanceReach;
@@ -52,6 +56,13 @@
import explicit.rewards.MDPRewards;
import explicit.rewards.Rewards;
+import org.gnu.glpk.GLPK;
+import org.gnu.glpk.GLPKConstants;
+import org.gnu.glpk.SWIGTYPE_p_double;
+import org.gnu.glpk.SWIGTYPE_p_int;
+import org.gnu.glpk.glp_prob;
+import org.gnu.glpk.glp_smcp;
+
/**
* Explicit-state model checker for Markov decision processes (MDPs).
*/
@@ -335,8 +346,9 @@ public ModelCheckerResult computeReachProbs(MDP mdp, BitSet remain, BitSet targe
MDPSolnMethod mdpSolnMethod = this.mdpSolnMethod;
// Switch to a supported method, if necessary
- if (mdpSolnMethod == MDPSolnMethod.LINEAR_PROGRAMMING) {
+ if (settings.EXACT_LP_REQUIREMENTS_INSTALLED == 0) {
mdpSolnMethod = MDPSolnMethod.GAUSS_SEIDEL;
+ mainLog.printWarning("Cannot compute exact linear programming. GLPK libraries are not properly installed");
mainLog.printWarning("Switching to MDP solution method \"" + mdpSolnMethod.fullName() + "\"");
}
@@ -447,6 +459,9 @@ public ModelCheckerResult computeReachProbs(MDP mdp, BitSet remain, BitSet targe
case MODIFIED_POLICY_ITERATION:
res = computeReachProbsModPolIter(mdp, no, yes, min, strat);
break;
+ case EXACT_LINEAR_PROGRAMMING:
+ res = computeReachProbsExactLinearProg(mdp, no, yes, target, min, strat);
+ break;
default:
throw new PrismException("Unknown MDP solution method " + mdpSolnMethod.fullName());
}
@@ -1049,6 +1064,277 @@ protected ModelCheckerResult computeReachProbsModPolIter(MDP mdp, BitSet no, Bit
res.timeTaken = timer / 1000.0;
return res;
}
+
+ /**
+ * Compute arbitrary precision reachability probabilities using linear programming.
+ * @param mdp: The MDP
+ * @param no: Probability 0 states
+ * @param yes: Probability 1 states
+ * @param target: Target states
+ * @param min: Min or max probabilities (true=min, false=max)
+ * @param strat Storage for (memoryless) strategy choice indices (ignored if null)
+ */
+ protected ModelCheckerResult computeReachProbsExactLinearProg(MDP mdp, BitSet no, BitSet yes, BitSet target, boolean min, int strat[]) throws PrismException
+ {
+ ModelCheckerResult res;
+ int nStates, maybeSize, countMaybeTrans, ret;
+ double sol, soln[], matAEntry, bEntry;
+ long timer, basisTimer, matATimer, exactSolverTimer;
+ BitSet maybe;
+ SWIGTYPE_p_int ind;
+ SWIGTYPE_p_double val;
+ glp_smcp parm;
+ glp_prob lp;
+
+ // Start linear programming
+ timer = System.currentTimeMillis();
+ mainLog.println("Starting exact linear programming (" + (min ? "min" : "max") + ")...");
+
+ // Store num states
+ nStates = mdp.getNumStates();
+
+ // Solution vector where we store the exact
+ // probabilities as doubles
+ soln = new double[nStates];
+ BitSet states = new BitSet();
+ states.set(0, nStates);
+
+ // Determine the set of maybe states (S?)
+ maybe = new BitSet();
+ maybe.set(0, nStates);
+ maybe.andNot(yes);
+ maybe.andNot(no);
+
+ // Store the cardinality of the maybe state
+ maybeSize = maybe.cardinality();
+
+ // The set of basic columns.
+ // This set will be filled as defined in (10)
+ Set basicCols = new HashSet();
+
+ // A mapping between maybe states and its positions in the maybe bitset.
+ // This is useful for calculating the `A` matrix later on.
+ // Also, it is calculated the number of the transitions of all maybe states `countMaybeTrans`
+ countMaybeTrans = 0;
+ HashMap maybeToState = new HashMap();
+ HashMap stateToMaybe = new HashMap();
+ int maybeStateCount = 1;
+ for(int state = maybe.nextSetBit(0); state >= 0; state = maybe.nextSetBit(state + 1)) {
+ maybeToState.put(maybeStateCount, state);
+ stateToMaybe.put(state, maybeStateCount);
+ countMaybeTrans += mdp.getNumChoices(state);
+ maybeStateCount++;
+ }
+
+ // Memory allocation for index and value arrays.
+ // For each transition, GLPK stores the values for each state in the maybe set
+ // and uses these arrays as temporal storage
+ ind = GLPK.new_intArray(maybeSize + 1);
+ val = GLPK.new_doubleArray(maybeSize + 1);
+
+ // Construct max (min) linear programming problem.
+ // First, we put the problem in the canonical form
+ // and then apply the duality definition
+ lp = GLPK.glp_create_prob();
+ mainLog.println("Linear programming problem created");
+
+ // Set the lp constants for the dual simplex method.
+ // columns -> one for each transition
+ // rows -> maybe states
+ GLPK.glp_add_cols(lp, countMaybeTrans);
+ GLPK.glp_add_rows(lp, maybeSize);
+ GLPK.intArray_setitem(ind, 0, 0);
+ GLPK.doubleArray_setitem(val, 0, 0);
+
+ // Constraints for each equation expressed as inequalities.
+ // For the min (max) problem, each equation has -1 as a upper (lower) bound and we compute the max (min).
+ // This was obtained rewriting the problems (6) (7) in their canonical form, and then construct the dual.
+ // The `-1` bounds correspond to the values of the cost vector (remember that we put the problems
+ // in canonical form first).
+ // Also, we set each row (representing a maybe state) as part of the basis defined in (10)
+ for (int row = 1; row <= maybeSize; row++) {
+ // The name of the state
+ GLPK.glp_set_row_name(lp, row, "s" + maybe.nextSetBit(row));
+ if (min) {
+ GLPK.glp_set_row_bnds(lp, row, GLPKConstants.GLP_UP, 0, -1);
+ GLPK.glp_set_obj_dir(lp, GLPKConstants.GLP_MAX);
+ } else {
+ GLPK.glp_set_row_bnds(lp, row, GLPKConstants.GLP_LO, -1, 0);
+ GLPK.glp_set_obj_dir(lp, GLPKConstants.GLP_MIN);
+ }
+ GLPK.glp_set_row_stat(lp, row, GLPKConstants.GLP_BS);
+ GLPK.intArray_setitem(ind, row, 0);
+ GLPK.doubleArray_setitem(val, row, 0);
+ }
+
+ // The current column of the matrix (transition)
+ int currentTransition = 1;
+ matATimer = System.currentTimeMillis();
+ // Fill out the (A | I) matrix (lp problem equations) of the reachability problem.
+ // In this step, we calculate the values for each entry of the A matrix,
+ // going through all maybe states and their transitions for each choice
+ for (int state = maybe.nextSetBit(0); state >= 0; state = maybe.nextSetBit(state + 1)) {
+ for (int choice = 0; choice < mdp.getNumChoices(state); choice++) {
+ // Index of the i-th equation in the lp's matrix
+ int eqIndex = 1;
+ Iterator> iter = mdp.getTransitionsIterator(state, choice);
+ bEntry = 0.;
+ boolean pointsToItself = false;
+ // Iterate over the transitions between st and choice
+ while (iter.hasNext()) {
+ matAEntry = 0.;
+ Map.Entry entry = iter.next();
+ int actualKey = entry.getKey();
+ // Skip if the probability is 0
+ if (entry.getValue() == 0)
+ continue;
+ // Calculate the entries (summation) for the b vector (constraint vector)
+ if (yes.get(actualKey))
+ bEntry += entry.getValue();
+ // If it is a maybe state, we calculate the matrix entry as defined on
+ // page 8 of the paper
+ if (maybe.get(actualKey)) {
+ if (actualKey == state) {
+ pointsToItself = true;
+ matAEntry = entry.getValue() - 1;
+ } else {
+ matAEntry = entry.getValue();
+ }
+ // Write the entry of each equation in the array
+ GLPK.intArray_setitem(ind, eqIndex, stateToMaybe.get(actualKey));
+ GLPK.doubleArray_setitem(val, eqIndex, matAEntry);
+ eqIndex++;
+ }
+ }
+ // If the state does not point to itself, then fill out the missing entry
+ if (!pointsToItself) {
+ GLPK.intArray_setitem(ind, eqIndex, stateToMaybe.get(state));
+ GLPK.doubleArray_setitem(val, eqIndex, -1);
+ eqIndex++;
+ }
+ // This is described in (10) on the paper
+ // If strat is null, it would mean that the strategy is not apt, so we should ignore
+ // the basis construction for this case
+ if (strat != null && choice == strat[state])
+ basicCols.add(currentTransition);
+
+ // Set the values of the arrays in the full linear programming equation matrix.
+ // Negate the summation calculated for the constraint vector entries
+ GLPK.glp_set_col_bnds(lp, currentTransition, GLPKConstants.GLP_LO, 0, 0);
+ GLPK.glp_set_mat_col(lp, currentTransition, eqIndex - 1, ind, val);
+ GLPK.glp_set_obj_coef(lp, currentTransition, -bEntry);
+
+ currentTransition++;
+ }
+ }
+ mainLog.println("Linear programming equations construction time: " +
+ ((double)System.currentTimeMillis() - matATimer)/1000 + " seconds");
+
+ // Determine the valid initial basis defined in (10).
+ // Set every basic column as basic (GLP_BS constant), otherwise
+ // it is non-basic having an active lower bound (GLP_NL)
+ basisTimer = System.currentTimeMillis();
+ for (int col = 1; col < currentTransition; col++) {
+ if (basicCols.contains(col)) {
+ GLPK.glp_set_col_stat(lp, col, GLPKConstants.GLP_BS);
+ } else {
+ GLPK.glp_set_col_stat(lp, col, GLPKConstants.GLP_NL);
+ }
+ }
+ mainLog.println("Basis construction time: " +
+ ((double)System.currentTimeMillis() - basisTimer)/1000 + " seconds");
+
+ // Solving the linear problem in the linear programming library GLPK.
+ // probs[]: an array for storing the exact probabilities as string in the form `p/q`
+ // `num` and `den` store the temporal `p` and `q` respectively (this could be avoided)
+ String probs[] = new String[nStates];
+
+ // Set up fixed parameters
+ parm = new glp_smcp();
+ GLPK.glp_init_smcp(parm);
+ parm.setPrecision(this.getPrecision());
+ parm.setNum_states(nStates);
+
+ // Solve lp problem
+ exactSolverTimer = System.currentTimeMillis();
+ ret = GLPK.glp_exact(lp, parm);
+
+ // Check return value: if `ret` is not zero there has been an error
+ if (ret != 0) {
+ String msg;
+ if (ret == GLPKConstants.GLP_EBADB) {
+ msg = "The initial basis specified in the problem object is invalid: wrong number of basic variables";
+ } else if (ret == GLPKConstants.GLP_ESING) {
+ msg = "The basis matrix corresponding to the initial basis is singular";
+ } else if (ret == GLPKConstants.GLP_EBOUND) {
+ msg = "Some double-bounded variables have incorrect bound";
+ } else {
+ msg = "The linear programming problem could not be solved";
+ }
+ mainLog.println(msg);
+ throw new PrismException(msg);
+ }
+
+ // Calculate and store exact probabilities.
+ // For each maybe state, we extract the fixed precision results from GLPK and store them in `probs[]`.
+ // `sol` stores the exact values as doubles. That implies that the precision is at most 64-bit
+ // However, the exact value is expressed as a string in `probs[]`.
+ sol = GLPK.glp_get_obj_val(lp);
+ String currentResult;
+ int exp;
+ mainLog.println("\n" + this.getPrecision() + "-bit results (non-zero only) for each state of the MDP:");
+ for (int mdpState = 0; mdpState < nStates; mdpState++) {
+ if (maybe.get(mdpState)) {
+ sol = GLPK.glp_get_row_dual(lp, stateToMaybe.get(mdpState));
+ soln[mdpState] = sol;
+ currentResult = GLPK.uint8Array_toString(GLPK.uint8ArrayArray_getitem(GLPK.glp_get_probs(lp), stateToMaybe.get(mdpState) - 1));
+ exp = GLPK.intArray_getitem(GLPK.glp_get_exps(lp), stateToMaybe.get(mdpState) - 1);
+ // We need to format the output since we get the exponent as a position in the string
+ String zeroes = "";
+ for (int nZero = 0; nZero < Math.abs(exp); nZero++)
+ zeroes += "0";
+ currentResult = exp > 0 ? new StringBuilder(currentResult).insert(exp, ".").toString() : "0." + zeroes + currentResult;
+ probs[mdpState] = currentResult;
+ }
+ // If it is a `yes` state, its probability is 1
+ else if (yes.get(mdpState)) {
+ soln[mdpState] = 1.;
+ probs[mdpState] = "1";
+ }
+ // If it is a `no` state, its probability is 0
+ else if (no.get(mdpState)) {
+ soln[mdpState] = 0.;
+ probs[mdpState] = "0";
+ }
+ // Avoid zero entries
+ if (probs[mdpState] != "0")
+ mainLog.println("[" + mdpState + "]: " + probs[mdpState]);
+ }
+ mainLog.println("\nExact probabilities calculation time: " +
+ ((double)System.currentTimeMillis() - exactSolverTimer)/1000 + " seconds");
+
+ // Free memory
+ GLPK.delete_intArray(ind);
+ GLPK.delete_doubleArray(val);
+ GLPK.delete_intArray(GLPK.glp_get_exps(lp));
+ GLPK.delete_uint8ArrayArray(GLPK.glp_get_probs(lp));
+ GLPK.delete_uint8ArrayArray(GLPK.glp_get_nums(lp));
+ GLPK.delete_uint8ArrayArray(GLPK.glp_get_dens(lp));
+ GLPK.glp_delete_prob(lp);
+
+ // Finished linear programming
+ timer = System.currentTimeMillis() - timer;
+ mainLog.print("Exact linear programming");
+ mainLog.println(" took " + timer / 1000.0 + " seconds.");
+
+ // Return results
+ // (Note we don't add the strategy - the one passed in is already there
+ // and might have some existing choices stored for other states).
+ res = new ModelCheckerResult();
+ res.timeTaken = timer / 1000.0;
+ res.soln = soln;
+ return res;
+ }
/**
* Construct strategy information for min/max reachability probabilities.
diff --git a/prism/src/explicit/ProbModelChecker.java b/prism/src/explicit/ProbModelChecker.java
index 40f9d1a57..13a0d9ef3 100644
--- a/prism/src/explicit/ProbModelChecker.java
+++ b/prism/src/explicit/ProbModelChecker.java
@@ -73,6 +73,8 @@ public class ProbModelChecker extends NonProbModelChecker
protected double termCritParam = 1e-8;
// Max iterations for numerical solution
protected int maxIters = 100000;
+ // Floating point precision
+ protected int precision = 64;
// Use precomputation algorithms in model checking?
protected boolean precomp = true;
protected boolean prob0 = true;
@@ -93,7 +95,7 @@ public class ProbModelChecker extends NonProbModelChecker
// Method used for numerical solution
public enum LinEqMethod {
- POWER, JACOBI, GAUSS_SEIDEL, BACKWARDS_GAUSS_SEIDEL, JOR, SOR, BACKWARDS_SOR;
+ POWER, JACOBI, GAUSS_SEIDEL, BACKWARDS_GAUSS_SEIDEL, EXACT_LINEAR_PROGRAMMING, JOR, SOR, BACKWARDS_SOR;
public String fullName()
{
switch (this) {
@@ -105,6 +107,8 @@ public String fullName()
return "Gauss-Seidel";
case BACKWARDS_GAUSS_SEIDEL:
return "Backwards Gauss-Seidel";
+ case EXACT_LINEAR_PROGRAMMING:
+ return "Exact linear programming";
case JOR:
return "JOR";
case SOR:
@@ -119,7 +123,7 @@ public String fullName()
// Method used for solving MDPs
public enum MDPSolnMethod {
- VALUE_ITERATION, GAUSS_SEIDEL, POLICY_ITERATION, MODIFIED_POLICY_ITERATION, LINEAR_PROGRAMMING;
+ VALUE_ITERATION, GAUSS_SEIDEL, POLICY_ITERATION, MODIFIED_POLICY_ITERATION, LINEAR_PROGRAMMING, EXACT_LINEAR_PROGRAMMING;
public String fullName()
{
switch (this) {
@@ -133,6 +137,8 @@ public String fullName()
return "Modified policy iteration";
case LINEAR_PROGRAMMING:
return "Linear programming";
+ case EXACT_LINEAR_PROGRAMMING:
+ return "Exact linear programming";
default:
return this.toString();
}
@@ -151,7 +157,7 @@ public enum ValIterDir {
// Method used for numerical solution
public enum SolnMethod {
- VALUE_ITERATION, GAUSS_SEIDEL, POLICY_ITERATION, MODIFIED_POLICY_ITERATION, LINEAR_PROGRAMMING
+ VALUE_ITERATION, GAUSS_SEIDEL, POLICY_ITERATION, MODIFIED_POLICY_ITERATION, LINEAR_PROGRAMMING, EXACT_LINEAR_PROGRAMMING
};
/**
@@ -180,6 +186,8 @@ public ProbModelChecker(PrismComponent parent) throws PrismException
setLinEqMethod(LinEqMethod.SOR);
} else if (s.equals("Backwards SOR")) {
setLinEqMethod(LinEqMethod.BACKWARDS_SOR);
+ } else if (s.equals("Exact linear programming")) {
+ setLinEqMethod(LinEqMethod.EXACT_LINEAR_PROGRAMMING);
} else {
throw new PrismNotSupportedException("Explicit engine does not support linear equation solution method \"" + s + "\"");
}
@@ -195,6 +203,8 @@ public ProbModelChecker(PrismComponent parent) throws PrismException
setMDPSolnMethod(MDPSolnMethod.MODIFIED_POLICY_ITERATION);
} else if (s.equals("Linear programming")) {
setMDPSolnMethod(MDPSolnMethod.LINEAR_PROGRAMMING);
+ } else if (s.equals("Exact linear programming")) {
+ setMDPSolnMethod(MDPSolnMethod.EXACT_LINEAR_PROGRAMMING);
} else {
throw new PrismNotSupportedException("Explicit engine does not support MDP solution method \"" + s + "\"");
}
@@ -211,6 +221,8 @@ public ProbModelChecker(PrismComponent parent) throws PrismException
setTermCritParam(settings.getDouble(PrismSettings.PRISM_TERM_CRIT_PARAM));
// PRISM_MAX_ITERS
setMaxIters(settings.getInteger(PrismSettings.PRISM_MAX_ITERS));
+ // PRISM_PRECISION
+ setPrecision(settings.getInteger(PrismSettings.PRISM_PRECISION));
// PRISM_PRECOMPUTATION
setPrecomp(settings.getBoolean(PrismSettings.PRISM_PRECOMPUTATION));
// PRISM_PROB0
@@ -248,6 +260,7 @@ public void inheritSettings(ProbModelChecker other)
setTermCrit(other.getTermCrit());
setTermCritParam(other.getTermCritParam());
setMaxIters(other.getMaxIters());
+ setPrecision(other.getPrecision());
setPrecomp(other.getPrecomp());
setProb0(other.getProb0());
setProb1(other.getProb1());
@@ -267,6 +280,7 @@ public void printSettings()
mainLog.print("termCrit = " + termCrit + " ");
mainLog.print("termCritParam = " + termCritParam + " ");
mainLog.print("maxIters = " + maxIters + " ");
+ mainLog.print("precision = " + precision + " ");
mainLog.print("precomp = " + precomp + " ");
mainLog.print("prob0 = " + prob0 + " ");
mainLog.print("prob1 = " + prob1 + " ");
@@ -325,6 +339,14 @@ public void setMaxIters(int maxIters)
this.maxIters = maxIters;
}
+ /**
+ * Set accuracy in bits for arbitrary precision results.
+ */
+ public void setPrecision(int precision)
+ {
+ this.precision = precision;
+ }
+
/**
* Set whether or not to use precomputation (Prob0, Prob1, etc.).
*/
@@ -423,6 +445,11 @@ public int getMaxIters()
return maxIters;
}
+ public int getPrecision()
+ {
+ return precision;
+ }
+
public boolean getPrecomp()
{
return precomp;
diff --git a/prism/src/prism/Makefile b/prism/src/prism/Makefile
index 62e786b97..1c9e087ca 100644
--- a/prism/src/prism/Makefile
+++ b/prism/src/prism/Makefile
@@ -64,6 +64,7 @@ $(PRISM_DIR_REL)/$(OBJ_DIR)/$(THIS_DIR)/ngprism$(EXE): ngprism.c
clean: checks
@rm -f $(CLASS_FILES) $(PRISM_DIR_REL)/$(LIB_DIR)/$(LIBPREFIX)prism$(LIBSUFFIX) $(O_FILES) $(PRISM_DIR_REL)/$(OBJ_DIR)/$(THIS_DIR)/ngprism$(EXE)
+ @sed -i '/EXACT_LP_REQUIREMENTS_INSTALLED/c\ public final int EXACT_LP_REQUIREMENTS_INSTALLED = 1;' PrismSettings.java
celan: clean
diff --git a/prism/src/prism/Prism.java b/prism/src/prism/Prism.java
index 82d2b5e03..574eee7df 100644
--- a/prism/src/prism/Prism.java
+++ b/prism/src/prism/Prism.java
@@ -2787,7 +2787,7 @@ public Result modelCheck(PropertiesFile propertiesFile, Property prop) throws Pr
}
// For exact model checking
- if (settings.getBoolean(PrismSettings.PRISM_EXACT_ENABLED)) {
+ if (settings.getBoolean(PrismSettings.PRISM_EXACT_ENABLED) && !settings.getBoolean(PrismSettings.PRISM_EXACT_LP_ENABLED)) {
return modelCheckExact(propertiesFile, prop);
}
// For fast adaptive uniformisation
diff --git a/prism/src/prism/PrismSettings.java b/prism/src/prism/PrismSettings.java
index a3e9f4004..3a91ddf1d 100644
--- a/prism/src/prism/PrismSettings.java
+++ b/prism/src/prism/PrismSettings.java
@@ -89,7 +89,8 @@ public class PrismSettings implements Observer
public static final String PRISM_TERM_CRIT = "prism.termCrit";//"prism.termination";
public static final String PRISM_TERM_CRIT_PARAM = "prism.termCritParam";//"prism.terminationEpsilon";
public static final String PRISM_MAX_ITERS = "prism.maxIters";//"prism.maxIterations";
-
+ public static final String PRISM_PRECISION = "prism.precision";//"prism.precision";
+
public static final String PRISM_CUDD_MAX_MEM = "prism.cuddMaxMem";
public static final String PRISM_CUDD_EPSILON = "prism.cuddEpsilon";
public static final String PRISM_NUM_SB_LEVELS = "prism.numSBLevels";//"prism.hybridNumLevels";
@@ -102,6 +103,7 @@ public class PrismSettings implements Observer
public static final String PRISM_SCC_METHOD = "prism.sccMethod";
public static final String PRISM_SYMM_RED_PARAMS = "prism.symmRedParams";
public static final String PRISM_EXACT_ENABLED = "prism.exact.enabled";
+ public static final String PRISM_EXACT_LP_ENABLED = "prism.exact.enabled";
public static final String PRISM_PTA_METHOD = "prism.ptaMethod";
public static final String PRISM_TRANSIENT_METHOD = "prism.transientMethod";
public static final String PRISM_AR_OPTIONS = "prism.arOptions";
@@ -185,8 +187,8 @@ public class PrismSettings implements Observer
public static final String LOG_SELECTION_COLOUR = "log.selectionColour";
public static final String LOG_BG_COLOUR = "log.bgColour";
public static final String LOG_BUFFER_LENGTH = "log.bufferLength";
-
-
+
+
//Defaults, types and constaints
public static final String[] propertyOwnerNames =
@@ -228,17 +230,19 @@ public class PrismSettings implements Observer
"Which engine (hybrid, sparse, MTBDD, explicit) should be used for model checking." },
{ BOOLEAN_TYPE, PRISM_EXACT_ENABLED, "Do exact model checking", "4.2.1", new Boolean(false), "",
"Perform exact model checking." },
+ { BOOLEAN_TYPE, PRISM_EXACT_LP_ENABLED, "Do exact model checking using linear programming", "4.2.1", new Boolean(false), "",
+ "Perform exact model checking using linear programming." },
{ CHOICE_TYPE, PRISM_PTA_METHOD, "PTA model checking method", "3.3", "Stochastic games", "Digital clocks,Stochastic games,Backwards reachability",
"Which method to use for model checking of PTAs." },
{ CHOICE_TYPE, PRISM_TRANSIENT_METHOD, "Transient probability computation method", "3.3", "Uniformisation", "Uniformisation,Fast adaptive uniformisation",
"Which method to use for computing transient probabilities in CTMCs." },
// NUMERICAL SOLUTION OPTIONS:
- { CHOICE_TYPE, PRISM_LIN_EQ_METHOD, "Linear equations method", "2.1", "Jacobi", "Power,Jacobi,Gauss-Seidel,Backwards Gauss-Seidel,Pseudo-Gauss-Seidel,Backwards Pseudo-Gauss-Seidel,JOR,SOR,Backwards SOR,Pseudo-SOR,Backwards Pseudo-SOR",
+ { CHOICE_TYPE, PRISM_LIN_EQ_METHOD, "Linear equations method", "2.1", "Jacobi", "Power,Jacobi,Gauss-Seidel,Backwards Gauss-Seidel,Pseudo-Gauss-Seidel,Backwards Pseudo-Gauss-Seidel,JOR,SOR,Backwards SOR,Pseudo-SOR,Backwards Pseudo-SOR,Exact linear programming",
"Which iterative method to use when solving linear equation systems." },
{ DOUBLE_TYPE, PRISM_LIN_EQ_METHOD_PARAM, "Over-relaxation parameter", "2.1", new Double(0.9), "",
"Over-relaxation parameter for iterative numerical methods such as JOR/SOR." },
- { CHOICE_TYPE, PRISM_MDP_SOLN_METHOD, "MDP solution method", "4.0", "Value iteration", "Value iteration,Gauss-Seidel,Policy iteration,Modified policy iteration,Linear programming",
+ { CHOICE_TYPE, PRISM_MDP_SOLN_METHOD, "MDP solution method", "4.0", "Value iteration", "Value iteration,Gauss-Seidel,Policy iteration,Modified policy iteration,Linear programming,Exact linear programming",
"Which method to use when solving Markov decision processes." },
{ CHOICE_TYPE, PRISM_MDP_MULTI_SOLN_METHOD, "MDP multi-objective solution method", "4.0.3", "Value iteration", "Value iteration,Gauss-Seidel,Linear programming",
"Which method to use when solving multi-objective queries on Markov decision processes." },
@@ -248,6 +252,8 @@ public class PrismSettings implements Observer
"Epsilon value to use for checking termination of iterative numerical methods." },
{ INTEGER_TYPE, PRISM_MAX_ITERS, "Termination max. iterations", "2.1", new Integer(10000), "0,",
"Maximum number of iterations to perform if iterative methods do not converge." },
+ { INTEGER_TYPE, PRISM_PRECISION, "Floating point precision", "2.1", new Integer(64), "1,",
+ "The accuracy of the mantissa (in bits) is determined by this variable. This user-selectable precision is a minimum value: it is rounded up to a whole limb." },
// MODEL CHECKING OPTIONS:
{ BOOLEAN_TYPE, PRISM_PRECOMPUTATION, "Use precomputation", "2.1", new Boolean(true), "",
"Whether to use model checking precomputation algorithms (Prob0, Prob1, etc.), where optional." },
@@ -416,13 +422,16 @@ public class PrismSettings implements Observer
{ INTEGER_TYPE, LOG_BUFFER_LENGTH, "Buffer length", "2.1", new Integer(10000), "1,", "Length of the buffer for the log display." }
}
};
-
+
public static final String[] oldPropertyNames = {"simulator.apmcStrategy", "simulator.engine", "simulator.newPathAskDefault"};
-
+
+ //Exact LP requirements
+ public final int EXACT_LP_REQUIREMENTS_INSTALLED = 1;
+
public DefaultSettingOwner[] optionOwners;
private Hashtable data;
private boolean modified;
-
+
private ArrayList settingsListeners;
/**
@@ -906,6 +915,12 @@ else if (sw.equals("explicit") || sw.equals("ex")) {
else if (sw.equals("exact")) {
set(PRISM_EXACT_ENABLED, true);
}
+ // Exact linear programming model checking
+ else if (sw.equals("exact-lp")) {
+ set(PRISM_EXACT_LP_ENABLED, true);
+ set(PRISM_MDP_SOLN_METHOD, "Exact linear programming");
+ set(PRISM_LIN_EQ_METHOD, "Exact linear programming");
+ }
// PTA model checking methods
else if (sw.equals("ptamethod")) {
if (i < args.length - 1) {
@@ -1025,7 +1040,22 @@ else if (sw.equals("maxiters")) {
throw new PrismException("No value specified for -" + sw + " switch");
}
}
-
+ // Precision
+ else if (sw.equals("precision")) {
+ if (i < args.length - 1) {
+ try {
+ j = Integer.parseInt(args[++i]);
+ if (j <= 0)
+ throw new NumberFormatException("");
+ set(PRISM_PRECISION, j);
+ } catch (NumberFormatException e) {
+ throw new PrismException("Invalid value for -" + sw + " switch");
+ }
+ } else {
+ throw new PrismException("No value specified for -" + sw + " switch");
+ }
+ }
+
// MODEL CHECKING OPTIONS:
// Precomputation algs off
@@ -1585,6 +1615,7 @@ public static void printHelp(PrismLog mainLog)
mainLog.println("-hybrid (or -h) ................ Use the Hybrid engine [default]");
mainLog.println("-explicit (or -ex) ............. Use the explicit engine");
mainLog.println("-exact ......................... Perform exact (arbitrary precision) model checking");
+ mainLog.println("-exact-lp ...................... Perform exact (arbitrary precision) model checking using linear programming");
mainLog.println("-ptamethod .............. Specify PTA engine (games, digital, backwards) [default: games]");
mainLog.println("-transientmethod ........ CTMC transient analysis methof (unif, fau) [default: unif]");
mainLog.println();