Added until-model checking for MDPs. Implemented Prob1A algorithm. Added asynchronous leader example.

src/formula/ProbabilisticNoBoundOperator.h

@@ -64,6 +64,15 @@ public:
 		// Intentionally left empty
 	}
 
+	/*!
+	 * Constructor
+	 *
+	 * @param pathFormula The child node.
+	 */
+	ProbabilisticNoBoundOperator(AbstractPathFormula<T>* pathFormula, bool minimumOperator) : NoBoundOperator<T>(pathFormula, minimumOperator) {
+		// Intentionally left empty
+	}
+
 	/*!
 	 * @returns a string representation of the formula
 	 */

src/formula/RewardNoBoundOperator.h

@@ -64,11 +64,28 @@ public:
 		// Intentionally left empty
 	}
 
+	/*!
+	 * Constructor
+	 *
+	 * @param pathFormula The child node.
+	 */
+	RewardNoBoundOperator(AbstractPathFormula<T>* pathFormula, bool minimumOperator) : NoBoundOperator<T>(pathFormula, minimumOperator) {
+		// Intentionally left empty
+	}
+
 	/*!
 	 * @returns a string representation of the formula
 	 */
 	virtual std::string toString() const {
-		std::string result = "R = ? [";
+		std::string result = "R";
+		if (this->isOptimalityOperator()) {
+			if (this->isMinimumOperator()) {
+				result += "min";
+			} else {
+				result += "max";
+			}
+		}
+		result += "=? [";
 		result += this->getPathFormula().toString();
 		result += "]";
 		return result;

src/modelchecker/GmmxxDtmcPrctlModelChecker.h

@@ -148,7 +148,7 @@ public:
 			// Prepare the right-hand side of the equation system. For entry i this corresponds to
 			// the accumulated probability of going from state i to some 'yes' state.
 			std::vector<Type> b(mayBeStatesSetBitCount);
-			this->getModel().getTransitionMatrix()->getConstrainedRowCountVector(maybeStates, statesWithProbability1, &b);
+			this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, statesWithProbability1, &b);
 
 			// Solve the corresponding system of linear equations.
 			this->solveLinearEquationSystem(*gmmxxMatrix, x, b);
@@ -346,9 +346,10 @@ public:
 	 */
 	static void putOptions(boost::program_options::options_description* desc) {
 		desc->add_options()("lemethod", boost::program_options::value<std::string>()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr}.");
-		desc->add_options()("lemaxiter", boost::program_options::value<unsigned>()->default_value(10000), "Sets the maximal number of iterations for iterative linear equation solving.");
-		desc->add_options()("precision", boost::program_options::value<double>()->default_value(10e-6), "Sets the precision for iterative linear equation solving.");
+		desc->add_options()("maxiter", boost::program_options::value<unsigned>()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving.");
+		desc->add_options()("precision", boost::program_options::value<double>()->default_value(1e-6), "Sets the precision for iterative equation solving.");
 		desc->add_options()("precond", boost::program_options::value<std::string>()->default_value("ilu")->notifier(&validatePreconditioner), "Sets the preconditioning technique for linear equation solving. Must be in {ilu, diagonal, ildlt, none}.");
+		desc->add_options()("relative", boost::program_options::value<bool>()->default_value(true), "Sets whether the relative or absolute error is considered for detecting convergence.");
 	}
 
 	/*!
@@ -388,7 +389,7 @@ private:
 
 		// Prepare an iteration object that determines the accuracy, maximum number of iterations
 		// and the like.
-		gmm::iteration iter(s->get<double>("precision"), 0, s->get<unsigned>("lemaxiter"));
+		gmm::iteration iter(s->get<double>("precision"), 0, s->get<unsigned>("maxiter"));
 
 		// Now do the actual solving.
 		LOG4CPLUS_INFO(logger, "Starting iterative solver.");

src/modelchecker/GmmxxMdpPrctlModelChecker.h

@@ -5,8 +5,8 @@
  *      Author: Christian Dehnert
  */
 
-#ifndef STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_
-#define STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_
+#ifndef STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_
 
 #include <cmath>
 
@@ -36,12 +36,12 @@ namespace modelChecker {
  * A model checking engine that makes use of the gmm++ backend.
  */
 template <class Type>
-class GmmxxDtmcPrctlModelChecker : public MdpPrctlModelChecker<Type> {
+class GmmxxMdpPrctlModelChecker : public MdpPrctlModelChecker<Type> {
 
 public:
-	explicit GmmxxDtmcPrctlModelChecker(storm::models::Mdp<Type>& mdp) : MdpPrctlModelChecker<Type>(mdp) { }
+	explicit GmmxxMdpPrctlModelChecker(storm::models::Mdp<Type>& mdp) : MdpPrctlModelChecker<Type>(mdp) { }
 
-	virtual ~GmmxxDtmcPrctlModelChecker() { }
+	virtual ~GmmxxMdpPrctlModelChecker() { }
 
 	virtual std::vector<Type>* checkBoundedUntil(const storm::formula::BoundedUntil<Type>& formula) const {
 		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
@@ -67,13 +67,13 @@ public:
 
 		// Create vector for result of multiplication, which is reduced to the result vector after
 		// each multiplication.
-		std::vector<Type>* multiplyResult = new std::vector<Type>(this->getModel().getTransitionMatrix().getRowCount());
+		std::vector<Type>* multiplyResult = new std::vector<Type>(this->getModel().getTransitionMatrix()->getRowCount());
 
 		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
 		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
 			gmm::mult(*gmmxxMatrix, *result, *multiplyResult);
 
-			if (minimumOperatorStack.top()) {
+			if (this->minimumOperatorStack.top()) {
 				storm::utility::reduceVectorMin(*multiplyResult, result, *nondeterministicChoiceIndices);
 			} else {
 				storm::utility::reduceVectorMax(*multiplyResult, result, *nondeterministicChoiceIndices);
@@ -103,7 +103,7 @@ public:
 		delete nextStates;
 
 		// Create resulting vector.
-		std::vector<Type>* temporaryResult = new std::vector<Type>(this->getModel().getTransitionMatrix().getRowCount());
+		std::vector<Type>* temporaryResult = new std::vector<Type>(this->getModel().getTransitionMatrix()->getRowCount());
 
 		// Perform the actual computation, namely matrix-vector multiplication.
 		gmm::mult(*gmmxxMatrix, *result, *temporaryResult);
@@ -112,7 +112,7 @@ public:
 		// vector properly.
 		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
 
-		if (minimumOperatorStack.top()) {
+		if (this->minimumOperatorStack.top()) {
 			storm::utility::reduceVectorMin(*temporaryResult, result, *nondeterministicChoiceIndices);
 		} else {
 			storm::utility::reduceVectorMax(*temporaryResult, result, *nondeterministicChoiceIndices);
@@ -133,7 +133,7 @@ public:
 		// all states that have probability 0 and 1 of satisfying the until-formula.
 		storm::storage::BitVector statesWithProbability0(this->getModel().getNumberOfStates());
 		storm::storage::BitVector statesWithProbability1(this->getModel().getNumberOfStates());
-		if (minimumOperatorStack.top()) {
+		if (this->minimumOperatorStack.top()) {
 			storm::utility::GraphAnalyzer::performProb01Min(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
 		} else {
 			storm::utility::GraphAnalyzer::performProb01Max(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
@@ -152,32 +152,36 @@ public:
 		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
 
 		// Only try to solve system if there are states for which the probability is unknown.
-		uint_fast64_t mayBeStatesSetBitCount = maybeStates.getNumberOfSetBits();
-		if (mayBeStatesSetBitCount > 0) {
-			// Now we can eliminate the rows and columns from the original transition probability matrix.
-			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates);
+		uint_fast64_t maybeStatesSetBitCount = maybeStates.getNumberOfSetBits();
+		if (maybeStatesSetBitCount > 0) {
+			// First, we can eliminate the rows and columns from the original transition probability matrix for states
+			// whose probabilities are already known.
+			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices());
 
-			// Transform the submatrix to the gmm++ format to use its solvers.
+			// Get the "new" nondeterministic choice indices for the submatrix.
+			std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
+
+			// Transform the submatrix to the gmm++ format to use its capabilities.
 			gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*submatrix);
-			delete submatrix;
 
-			// Initialize the x vector with 0.5 for each element. This is the initial guess for
-			// iteratively solving the equations.
-			std::vector<Type> x(mayBeStatesSetBitCount, Type(0.5));
+			// Create vector for results for maybe states.
+			std::vector<Type>* x = new std::vector<Type>(maybeStatesSetBitCount);
 
 			// Prepare the right-hand side of the equation system. For entry i this corresponds to
 			// the accumulated probability of going from state i to some 'yes' state.
-			std::vector<Type> b(mayBeStatesSetBitCount);
-			this->getModel().getTransitionMatrix()->getConstrainedRowCountVector(maybeStates, statesWithProbability1, &b);
+			std::vector<Type> b(submatrix->getRowCount());
+			this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, *this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, &b);
+			delete submatrix;
 
-			// Solve the corresponding system of linear equations.
-			this->solveEquationSystem(*gmmxxMatrix, x, b);
+			// Solve the corresponding system of equations.
+			this->solveEquationSystem(*gmmxxMatrix, x, b, *subNondeterministicChoiceIndices);
 
 			// Set values of resulting vector according to result.
-			storm::utility::setVectorValues<Type>(result, maybeStates, x);
+			storm::utility::setVectorValues<Type>(result, maybeStates, *x);
 
-			// Delete temporary matrix.
+			// Delete temporary matrix and vector.
 			delete gmmxxMatrix;
+			delete x;
 		}
 
 		// Set values of resulting vector that are known exactly.
@@ -272,7 +276,8 @@ public:
 		// Determine which states have a reward of infinity by definition.
 		storm::storage::BitVector infinityStates(this->getModel().getNumberOfStates());
 		storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true);
-		storm::utility::GraphAnalyzer::performProb1(this->getModel(), trueStates, *targetStates, &infinityStates);
+		// TODO: just commented out to make it compile
+		//storm::utility::GraphAnalyzer::performProb1(this->getModel(), trueStates, *targetStates, &infinityStates);
 		infinityStates.complement();
 
 		// Create resulting vector.
@@ -325,7 +330,8 @@ public:
 			}
 
 			// Solve the corresponding system of linear equations.
-			this->solveLinearEquationSystem(*gmmxxMatrix, x, *b);
+			// TODO: just commented out to make it compile
+			// this->solveEquationSystem(*gmmxxMatrix, x, *b);
 
 			// Set values of resulting vector according to result.
 			storm::utility::setVectorValues<Type>(result, maybeStates, x);
@@ -344,31 +350,6 @@ public:
 		return result;
 	}
 
-	/*!
-	 * Returns the name of this module.
-	 * @return The name of this module.
-	 */
-	static std::string getModuleName() {
-		return "gmm++nondet";
-	}
-
-	/*!
-	 * Returns a trigger such that if the option "matrixlib" is set to "gmm++", this model checker
-	 * is to be used.
-	 * @return An option trigger for this module.
-	 */
-	static std::pair<std::string, std::string> getOptionTrigger() {
-		return std::pair<std::string, std::string>("matrixlib", "gmm++");
-	}
-
-	/*!
-	 * Registers all options associated with the gmm++ matrix library.
-	 */
-	static void putOptions(boost::program_options::options_description* desc) {
-		desc->add_options()("lemaxiter", boost::program_options::value<unsigned>()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving.");
-		desc->add_options()("precision", boost::program_options::value<double>()->default_value(10e-6), "Sets the precision for iterative linear equation solving.");
-	}
-
 private:
 	/*!
 	 * Solves the given equation system under the given parameters using the power method.
@@ -380,25 +361,70 @@ private:
 	 * @return The solution of the system of linear equations in form of the elements of the vector
 	 * x.
 	 */
-	void solveLinearEquationSystem(gmm::csr_matrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b) const {
-		// Get the settings object to customize linear solving.
+	void solveEquationSystem(gmm::csr_matrix<Type> const& A, std::vector<Type>* x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
+		// Get the settings object to customize solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
-		// Get relevant user-defined settings for solving the equations
-		double precison = s->get<double>("precision");
-		unsigned maxIterations = s->get<unsigned>("lemaxiter");
+		// Get relevant user-defined settings for solving the equations.
+		double precision = s->get<double>("precision");
+		unsigned maxIterations = s->get<unsigned>("maxiter");
+		bool relative = s->get<bool>("relative");
 
+		// Set up the environment for the power method.
+		std::vector<Type>* temporaryResult = new std::vector<Type>(b.size());
+		std::vector<Type>* newX = new std::vector<Type>(x->size());
+		std::vector<Type>* swap = nullptr;
+		bool converged = false;
+		uint_fast64_t iterations = 0;
+
+		// Proceed with the iterations as long as the method did not converge or reach the
+		// user-specified maximum number of iterations.
+		while (!converged && iterations < maxIterations) {
+			// Compute x' = A*x + b.
+			gmm::mult(A, *x, *temporaryResult);
+			gmm::add(b, *temporaryResult);
+
+			// Reduce the vector x' by applying min/max for all non-deterministic choices.
+			if (this->minimumOperatorStack.top()) {
+				storm::utility::reduceVectorMin(*temporaryResult, newX, nondeterministicChoiceIndices);
+			} else {
+				storm::utility::reduceVectorMax(*temporaryResult, newX, nondeterministicChoiceIndices);
+			}
 
+			// Determine whether the method converged.
+			converged = storm::utility::equalModuloPrecision(*x, *newX, precision, relative);
 
+			// Update environment variables.
+			swap = x;
+			x = newX;
+			newX = swap;
+			++iterations;
+		}
 
+		delete temporaryResult;
 
 		// Check if the solver converged and issue a warning otherwise.
-		if (iter.converged()) {
-			LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
+		if (converged) {
+			LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
 		} else {
 			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
 		}
+	}
+
+	std::shared_ptr<std::vector<uint_fast64_t>> computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector constraint) const {
+		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+		std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices(new std::vector<uint_fast64_t>(constraint.getNumberOfSetBits() + 1));
+		uint_fast64_t currentRowCount = 0;
+		uint_fast64_t currentIndexCount = 1;
+		(*subNondeterministicChoiceIndices)[0] = 0;
+		for (auto index : constraint) {
+			(*subNondeterministicChoiceIndices)[currentIndexCount] = currentRowCount + (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
+			currentRowCount += (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
+			++currentIndexCount;
+		}
+		(*subNondeterministicChoiceIndices)[constraint.getNumberOfSetBits()] = currentRowCount;
 
+		return subNondeterministicChoiceIndices;
 	}
 };
 
@@ -406,4 +432,4 @@ private:
 
 } //namespace storm
 
-#endif /* STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_ */
+#endif /* STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_ */

src/modelchecker/MdpPrctlModelChecker.h

@@ -18,6 +18,7 @@
 #include "src/utility/Vector.h"
 #include "src/modelchecker/AbstractModelChecker.h"
 #include <vector>
+#include <stack>
 
 #include "log4cplus/logger.h"
 #include "log4cplus/loggingmacros.h"
@@ -361,11 +362,11 @@ public:
 	 */
 	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula) const = 0;
 
-protected:
-	std::stack<bool> minimumOperatorStack;
-
 private:
 	storm::models::Mdp<Type>& model;
+
+protected:
+	mutable std::stack<bool> minimumOperatorStack;
 };
 
 } //namespace modelChecker

src/models/GraphTransitions.h

@@ -216,7 +216,7 @@ private:
 		// Now we are ready to actually fill in the list of predecessors for
 		// every state. Again, we start by considering all but the last row.
 		for (uint_fast64_t i = 0; i < numberOfStates; i++) {
-			for (uint_fast64_t j = choiceIndices->at(i); j < choiceIndices->at(i + 1); ++j) {
+			for (uint_fast64_t j = (*choiceIndices)[i]; j < (*choiceIndices)[i + 1]; ++j) {
 				for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) {
 					this->successorList[nextIndicesList[*rowIt]++] = i;
 				}

src/parser/NondeterministicSparseTransitionParser.cpp

@@ -131,9 +131,17 @@ uint_fast64_t NondeterministicSparseTransitionParser::firstPass(char* buf, uint_
 		 */
 		nonzero++;
 
+		// The PRISM output format lists the name of the transition in the fourth column,
+		// but omits the fourth column if it is an internal action. In either case, however, the third column
+		// is followed by a space. We need to skip over that space first (instead of trimming whitespaces),
+		// before we can skip to the line end, because trimming the white spaces will proceed to the next line
+		// in case there is no action label in the fourth column.
+		++buf;
+
 		/*
 		 *	Proceed to beginning of next line.
 		 */
+		buf += strcspn(buf, " \t\n\r");
 		buf = trimWhitespaces(buf);
 	}
 
@@ -277,10 +285,12 @@ NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(s
 		/*
 		 *	Proceed to beginning of next line in file and next row in matrix.
 		 */
+		++buf;
+		buf += strcspn(buf, " \t\n\r");
 		buf = trimWhitespaces(buf);
 	}
 
-	this->rowMapping->at(maxnode+1) = curRow;
+	this->rowMapping->at(maxnode+1) = curRow + 1;
 
 	if (!fixDeadlocks && hadDeadlocks) throw storm::exceptions::WrongFileFormatException() << "Some of the nodes had deadlocks. You can use --fix-deadlocks to insert self-loops on the fly.";
 

src/storage/SparseMatrix.h

@@ -668,7 +668,6 @@ public:
 		uint_fast64_t rowEnd = rowIndications[row + 1];
 
 		if (rowStart >= rowEnd) {
-			this->print();
 			LOG4CPLUS_ERROR(logger, "Cannot make row absorbing, because there is no entry in this row.");
 			throw storm::exceptions::InvalidStateException("Cannot make row absorbing, because there is no entry in this row.");
 		}
@@ -722,13 +721,30 @@ public:
 	 * @param resultVector A pointer to the resulting vector that has at least
 	 * as many elements as there are bits set to true in the constraint.
 	 */
-	void getConstrainedRowCountVector(const storm::storage::BitVector& rowConstraint, const storm::storage::BitVector& columnConstraint, std::vector<T>* resultVector) const {
+	void getConstrainedRowSumVector(const storm::storage::BitVector& rowConstraint, const storm::storage::BitVector& columnConstraint, std::vector<T>* resultVector) const {
 		uint_fast64_t currentRowCount = 0;
 		for (auto row : rowConstraint) {
 			(*resultVector)[currentRowCount++] = getConstrainedRowSum(row, columnConstraint);
 		}
 	}
 
+	/*!
+	 * Computes a vector in which each element is the sum of those elements in the
+	 * corresponding row whose column bits are set to one in the given constraint.
+	 * @param rowConstraint A bit vector that indicates for which rows to perform summation.
+	 * @param columnConstraint A bit vector that indicates which columns to add.
+	 * @param resultVector A pointer to the resulting vector that has at least
+	 * as many elements as there are bits set to true in the constraint.
+	 */
+	void getConstrainedRowSumVector(const storm::storage::BitVector& rowConstraint, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, const storm::storage::BitVector& columnConstraint, std::vector<T>* resultVector) const {
+		uint_fast64_t currentRowCount = 0;
+		for (auto index : rowConstraint) {
+			for (uint_fast64_t row = nondeterministicChoiceIndices[index]; row < nondeterministicChoiceIndices[index + 1]; ++row) {
+				(*resultVector)[currentRowCount++] = getConstrainedRowSum(row, columnConstraint);
+			}
+		}
+	}
+
 	/*!
 	 * Creates a sub-matrix of the current matrix by dropping all rows and
 	 * columns whose bits are not set to one in the given bit vector.
@@ -761,7 +777,7 @@ public:
 
 		// Create a temporary array that stores for each index whose bit is set
 		// to true the number of bits that were set before that particular index.
-		uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[rowCount];
+		uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[colCount];
 		uint_fast64_t lastIndex = 0;
 		uint_fast64_t currentNumberOfSetBits = 0;
 		for (auto index : constraint) {
@@ -792,6 +808,70 @@ public:
 		return result;
 	}
 
+	SparseMatrix* getSubmatrix(storm::storage::BitVector& constraint, std::vector<uint_fast64_t> const& rowIndices) const {
+		LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size).");
+
+		// Check for valid constraint.
+		if (constraint.getNumberOfSetBits() == 0) {
+			LOG4CPLUS_ERROR(logger, "Trying to create a sub-matrix of size 0.");
+			throw storm::exceptions::InvalidArgumentException("Trying to create a sub-matrix of size 0.");
+		}
+
+		// First, we need to determine the number of non-zero entries and the number of rows of the
+		// sub-matrix.
+		uint_fast64_t subNonZeroEntries = 0;
+		uint_fast64_t subRowCount = 0;
+		for (auto index : constraint) {
+			subRowCount += rowIndices[index + 1] - rowIndices[index];
+			for (uint_fast64_t i = rowIndices[index]; i < rowIndices[index + 1]; ++i) {
+				for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
+					if (constraint.get(columnIndications[j])) {
+						++subNonZeroEntries;
+					}
+				}
+			}
+		}
+
+		LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << constraint.getNumberOfSetBits() << ".");
+
+		// Create and initialize resulting matrix.
+		SparseMatrix* result = new SparseMatrix(subRowCount, constraint.getNumberOfSetBits());
+		result->initialize(subNonZeroEntries);
+
+		// Create a temporary array that stores for each index whose bit is set
+		// to true the number of bits that were set before that particular index.
+		uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[colCount];
+		uint_fast64_t lastIndex = 0;
+		uint_fast64_t currentNumberOfSetBits = 0;
+		for (auto index : constraint) {
+			while (lastIndex <= index) {
+				bitsSetBeforeIndex[lastIndex++] = currentNumberOfSetBits;
+			}
+			++currentNumberOfSetBits;
+		}
+
+		// Copy over selected entries.
+		uint_fast64_t rowCount = 0;
+		for (auto index : constraint) {
+			for (uint_fast64_t i = rowIndices[index]; i < rowIndices[index + 1]; ++i) {
+				for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) {
+					if (constraint.get(columnIndications[j])) {
+						result->addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]);
+					}
+				}
+				++rowCount;
+			}
+		}
+
+		// Dispose of the temporary array.
+		delete[] bitsSetBeforeIndex;
+
+		// Finalize sub-matrix and return result.
+		result->finalize();
+		LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix.");
+		return result;
+	}
+
 	void convertToEquationSystem() {
 		invertDiagonal();
 		negateAllNonDiagonalElements();
@@ -1005,11 +1085,82 @@ public:
 		return true;
 	}
 
-	void print() const {
-		std::cout << "entries in (" << rowCount << "x" << colCount << " matrix):" << std::endl;
-		std::cout << rowIndications << std::endl;
-		std::cout << columnIndications << std::endl;
-		std::cout << valueStorage << std::endl;
+	/*!
+	 * Retrieves a compressed string representation of the matrix.
+	 * @return a compressed string representation of the matrix.
+	 */
+	std::string toStringCompressed() const {
+		std::stringstream result;
+		result << rowIndications << std::endl;
+		result << columnIndications << std::endl;
+		result << valueStorage << std::endl;
+		return result.str();
+	}
+
+	/*!
+	 * Retrieves a (non-compressed) string representation of the matrix.
+	 * Note: the matrix is presented densely. That is, all zeros are filled in and are part of the string
+	 * representation.
+	 * @param nondeterministicChoiceIndices A vector indicating which rows belong together. If given, rows belonging
+	 * to separate groups will be separated by a dashed line.
+	 * @return a (non-compressed) string representation of the matrix.
+	 */
+	std::string toString(std::shared_ptr<std::vector<std::uint_fast64_t>> nondeterministicChoiceIndices) const {
+		std::stringstream result;
+		uint_fast64_t currentNondeterministicChoiceIndex = 0;
+
+		// Print column numbers in header.
+		result << "\t\t";
+		for (uint_fast64_t i = 0; i < colCount; ++i) {
+			result << i << "\t";
+		}
+		result << std::endl;
+
+		// Iterate over all rows.
+		for (uint_fast64_t i = 0; i < rowCount; ++i) {
+			uint_fast64_t nextIndex = rowIndications[i];
+
+			// If we need to group of rows, print a dashed line in case we have moved to the next group of rows.
+			if (nondeterministicChoiceIndices != nullptr) {
+				if (i == (*nondeterministicChoiceIndices)[currentNondeterministicChoiceIndex]) {
+					if (i != 0) {
+						result << "\t(\t";
+						for (uint_fast64_t j = 0; j < colCount - 2; ++j) {
+							result << "----";
+							if (j == 1) {
+								result << "\t" << currentNondeterministicChoiceIndex << "\t";
+							}
+						}
+						result << "\t)" << std::endl;
+					}
+					++currentNondeterministicChoiceIndex;
+				}
+			}
+
+			// Print the actual row.
+			result << i << "\t(\t";
+			uint_fast64_t currentRealIndex = 0;
+			while (currentRealIndex < colCount) {
+				if (currentRealIndex == columnIndications[nextIndex] && nextIndex < rowIndications[i + 1]) {
+					result << valueStorage[nextIndex] << "\t";
+					++nextIndex;
+				} else {
+					result << "0\t";
+				}
+				++currentRealIndex;
+			}
+			result << "\t)\t" << i << std::endl;
+		}
+
+		// Print column numbers in footer.
+		result << "\t\t";
+		for (uint_fast64_t i = 0; i < colCount; ++i) {
+			result << i << "\t";
+		}
+		result << std::endl;
+
+		// Return final result.
+		return result.str();
 	}
 
 private:

src/storm.cpp

@@ -24,6 +24,7 @@
 #include "src/models/AtomicPropositionsLabeling.h"
 #include "src/modelchecker/EigenDtmcPrctlModelChecker.h"
 #include "src/modelchecker/GmmxxDtmcPrctlModelChecker.h"
+#include "src/modelchecker/GmmxxMdpPrctlModelChecker.h"
 #include "src/parser/AutoParser.h"
 #include "src/parser/PrctlParser.h"
 #include "src/utility/Settings.h"
@@ -125,6 +126,10 @@ bool parseOptions(const int argc, const char* argv[]) {
 	return true;
 }
 
+void setUp() {
+	std::cout.precision(10);
+}
+
 /*!
  * Function to perform some cleanup.
  */
@@ -206,18 +211,36 @@ void testCheckingSynchronousLeader(storm::models::Dtmc<double>& dtmc, uint_fast6
 }
 
 void testCheckingDice(storm::models::Mdp<double>& mdp) {
-	storm::storage::BitVector* targetStates = mdp.getLabeledStates(std::string("two"));
-	*targetStates |= *mdp.getLabeledStates(std::string("three"));
+	storm::formula::Ap<double>* threeFormula = new storm::formula::Ap<double>("three");
+	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(threeFormula);
+	storm::formula::ProbabilisticNoBoundOperator<double>* probFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
+
+	storm::modelChecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxMdpPrctlModelChecker<double>(mdp);
+
+	mc->check(*probFormula);
+	delete probFormula;
+
+	delete mc;
+}
+
+void testCheckingAsynchLeader(storm::models::Mdp<double>& mdp) {
+	storm::formula::Ap<double>* electedFormula = new storm::formula::Ap<double>("elected");
+	storm::formula::Eventually<double>* eventuallyFormula = new storm::formula::Eventually<double>(electedFormula);
+	storm::formula::ProbabilisticNoBoundOperator<double>* probMaxFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula, false);
 
-	storm::storage::BitVector* trueStates = new storm::storage::BitVector(mdp.getNumberOfStates(), true);
-	storm::storage::BitVector* statesWithProbability0 = new storm::storage::BitVector(mdp.getNumberOfStates());
-	storm::storage::BitVector* statesWithProbability1 = new storm::storage::BitVector(mdp.getNumberOfStates());
+	storm::modelChecker::GmmxxMdpPrctlModelChecker<double>* mc = new storm::modelChecker::GmmxxMdpPrctlModelChecker<double>(mdp);
 
-	storm::utility::GraphAnalyzer::performProb01Max(mdp, *trueStates, *targetStates, statesWithProbability0, statesWithProbability1);
+	mc->check(*probMaxFormula);
+	delete probMaxFormula;
 
-	delete statesWithProbability0;
-	delete statesWithProbability1;
-	delete trueStates;
+	electedFormula = new storm::formula::Ap<double>("elected");
+	eventuallyFormula = new storm::formula::Eventually<double>(electedFormula);
+	storm::formula::ProbabilisticNoBoundOperator<double>* probMinFormula = new storm::formula::ProbabilisticNoBoundOperator<double>(eventuallyFormula, true);
+
+	mc->check(*probMinFormula);
+	delete probMinFormula;
+
+	delete mc;
 }
 
 /*!
@@ -240,6 +263,7 @@ void testChecking() {
 		mdp->printModelInformationToStream(std::cout);
 
 		// testCheckingDice(*mdp);
+		// testCheckingAsynchLeader(*mdp);
 	} else {
 		std::cout << "Input is neither a DTMC nor an MDP." << std::endl;
 	}
@@ -253,6 +277,7 @@ int main(const int argc, const char* argv[]) {
 	if (!parseOptions(argc, argv)) {
 		return 0;
 	}
+	setUp();
 
 	LOG4CPLUS_INFO(logger, "StoRM was invoked.");
 	printHeader(argc, argv);

src/utility/GraphAnalyzer.h

@@ -232,13 +232,12 @@ public:
 
 				for(auto it = backwardTransitions.beginStateSuccessorsIterator(currentState); it != backwardTransitions.endStateSuccessorsIterator(currentState); ++it) {
 					if (phiStates.get(*it) && !nextStates->get(*it)) {
-
-						// Check whether the predecessor has at only successors in the current state set for one of the
+						// Check whether the predecessor has only successors in the current state set for one of the
 						// nondeterminstic choices.
-						for (auto rowIt = nondeterministicChoiceIndices->begin() + *it; rowIt != nondeterministicChoiceIndices->begin() + *it + 1; ++rowIt) {
+						for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) {
 							bool allSuccessorsInCurrentStates = true;
-							for (auto colIt = transitionMatrix->beginConstColumnIterator(*rowIt); colIt != transitionMatrix->endConstColumnIterator(*rowIt); ++colIt) {
-								if (currentStates->get(*colIt)) {
+							for (auto colIt = transitionMatrix->beginConstColumnIterator(row); colIt != transitionMatrix->endConstColumnIterator(row); ++colIt) {
+								if (!currentStates->get(*colIt)) {
 									allSuccessorsInCurrentStates = false;
 									break;
 								}
@@ -316,7 +315,6 @@ public:
 
 			for(auto it = backwardTransitions.beginStateSuccessorsIterator(currentState); it != backwardTransitions.endStateSuccessorsIterator(currentState); ++it) {
 				if (phiStates.get(*it) && !statesWithProbability0->get(*it)) {
-
 					// Check whether the predecessor has at least one successor in the current state
 					// set for every nondeterministic choice.
 					bool addToStatesWithProbability0 = true;
@@ -349,9 +347,70 @@ public:
 
 	template <class T>
 	static void performProb1A(storm::models::AbstractNondeterministicModel<T>& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) {
-		// This result is a rough guess and does not compute all states with probability 1.
-		// TODO: Check whether it makes sense to implement the precise but complicated algorithm here.
-		*statesWithProbability1 = psiStates;
+		// Check for valid parameters.
+		if (statesWithProbability1 == nullptr) {
+			LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability1' must not be null.");
+			throw storm::exceptions::InvalidArgumentException("Parameter 'statesWithProbability1' must not be null.");
+		}
+
+		// Get some temporaries for convenience.
+		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
+		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+
+		// Get the backwards transition relation from the model to ease the search.
+		storm::models::GraphTransitions<T> backwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), false);
+
+		storm::storage::BitVector* currentStates = new storm::storage::BitVector(model.getNumberOfStates(), true);
+
+		std::vector<uint_fast64_t> stack;
+		stack.reserve(model.getNumberOfStates());
+
+		bool done = false;
+		while (!done) {
+			stack.clear();
+			storm::storage::BitVector* nextStates = new storm::storage::BitVector(psiStates);
+			psiStates.addSetIndicesToList(stack);
+
+			while (!stack.empty()) {
+				uint_fast64_t currentState = stack.back();
+				stack.pop_back();
+
+				for(auto it = backwardTransitions.beginStateSuccessorsIterator(currentState); it != backwardTransitions.endStateSuccessorsIterator(currentState); ++it) {
+					if (phiStates.get(*it) && !nextStates->get(*it)) {
+						// Check whether the predecessor has only successors in the current state set for all of the
+						// nondeterminstic choices.
+						bool allSuccessorsInCurrentStatesForAllChoices = true;
+						for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) {
+							for (auto colIt = transitionMatrix->beginConstColumnIterator(row); colIt != transitionMatrix->endConstColumnIterator(row); ++colIt) {
+								if (!currentStates->get(*colIt)) {
+									allSuccessorsInCurrentStatesForAllChoices = false;
+									goto afterCheckLoop;
+								}
+							}
+						}
+
+				afterCheckLoop:
+						// If all successors for all nondeterministic choices are in the current state set, we
+						// add it to the set of states for the next iteration and perform a backward search from
+						// that state.
+						if (allSuccessorsInCurrentStatesForAllChoices) {
+							nextStates->set(*it, true);
+							stack.push_back(*it);
+						}
+					}
+				}
+			}
+
+			// Check whether we need to perform an additional iteration.
+			if (*currentStates == *nextStates) {
+				done = true;
+			} else {
+				*currentStates = *nextStates;
+			}
+		}
+
+		*statesWithProbability1 = *currentStates;
+		delete currentStates;
 	}
 };
 

src/utility/Vector.h

@@ -12,6 +12,11 @@
 #include "ConstTemplates.h"
 #include <iostream>
 
+#include "log4cplus/logger.h"
+#include "log4cplus/loggingmacros.h"
+
+extern log4cplus::Logger logger;
+
 namespace storm {
 
 namespace utility {
@@ -56,20 +61,19 @@ void subtractFromConstantOneVector(std::vector<T>* vector) {
 template<class T>
 void reduceVectorMin(std::vector<T> const& source, std::vector<T>* target, std::vector<uint_fast64_t> const& filter) {
 	uint_fast64_t currentSourceRow = 0;
-	uint_fast64_t currentTargetRow = 0;
-	for (auto it = source->cbegin(); it != source->cend(); ++it, ++currentSourceRow) {
-		// Check whether we have considered all from rows for the current to row.
-		if (filter[currentTargetRow + 1] <= currentSourceRow) {
+	uint_fast64_t currentTargetRow = -1;
+	for (auto it = source.cbegin(); it != source.cend(); ++it, ++currentSourceRow) {
+		// Check whether we have considered all source rows for the current target row.
+		if (filter[currentTargetRow + 1] <= currentSourceRow || currentSourceRow == 0) {
 			++currentTargetRow;
-			(*target)[currentTargetRow] = (*source)[currentSourceRow];
+			(*target)[currentTargetRow] = source[currentSourceRow];
 			continue;
 		}
 
-
 		// We have to minimize the value, so only overwrite the current value if the
 		// value is actually lower.
 		if (*it < (*target)[currentTargetRow]) {
-			(*source)[currentTargetRow] = *it;
+			(*target)[currentTargetRow] = *it;
 		}
 	}
 }
@@ -77,22 +81,39 @@ void reduceVectorMin(std::vector<T> const& source, std::vector<T>* target, std::
 template<class T>
 void reduceVectorMax(std::vector<T> const& source, std::vector<T>* target, std::vector<uint_fast64_t> const& filter) {
 	uint_fast64_t currentSourceRow = 0;
-	uint_fast64_t currentTargetRow = 0;
-	for (auto it = source->cbegin(); it != source->cend(); ++it, ++currentSourceRow) {
-		// Check whether we have considered all from rows for the current to row.
-		if (filter[currentTargetRow + 1] <= currentSourceRow) {
+	uint_fast64_t currentTargetRow = -1;
+	for (auto it = source.cbegin(); it != source.cend(); ++it, ++currentSourceRow) {
+		// Check whether we have considered all source rows for the current target row.
+		if (filter[currentTargetRow + 1] <= currentSourceRow || currentSourceRow == 0) {
 			++currentTargetRow;
-			(*target)[currentTargetRow] = (*source)[currentSourceRow];
+			(*target)[currentTargetRow] = source[currentSourceRow];
 			continue;
 		}
 
-
 		// We have to maximize the value, so only overwrite the current value if the
 		// value is actually greater.
 		if (*it > (*target)[currentTargetRow]) {
-			(*source)[currentTargetRow] = *it;
+			(*target)[currentTargetRow] = *it;
+		}
+	}
+}
+
+template<class T>
+bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight, T precision, bool relativeError) {
+	if (vectorLeft.size() != vectorRight.size()) {
+		LOG4CPLUS_ERROR(logger, "Length of vectors does not match and makes comparison impossible.");
+		throw storm::exceptions::InvalidArgumentException() << "Length of vectors does not match and makes comparison impossible.";
+	}
+
+	for (uint_fast64_t i = 0; i < vectorLeft.size(); ++i) {
+		if (relativeError) {
+			if (std::abs(vectorLeft[i] - vectorRight[i])/vectorRight[i] > precision) return false;
+		} else {
+			if (std::abs(vectorLeft[i] - vectorRight[i]) > precision) return false;
 		}
 	}
+
+	return true;
 }
 
 } //namespace utility