Started implementing the model checker for MDPs. Added reduce functionality to vector utility. Moved min/max capability to NoBoundOperator.

src/formula/NoBoundOperator.h

@@ -74,8 +74,8 @@ public:
 	/*!
 	 * Empty constructor
 	 */
-	NoBoundOperator() {
-		this->pathFormula = NULL;
+	NoBoundOperator() : optimalityOperator(false), minimumOperator(false) {
+		this->pathFormula = nullptr;
 	}
 
 	/*!
@@ -83,7 +83,19 @@ public:
 	 *
 	 * @param pathFormula The child node.
 	 */
-	NoBoundOperator(AbstractPathFormula<T>* pathFormula) {
+	NoBoundOperator(AbstractPathFormula<T>* pathFormula) : optimalityOperator(false), minimumOperator(false) {
+		this->pathFormula = pathFormula;
+	}
+
+	/*!
+	 * Constructor
+	 *
+	 * @param pathFormula The child node.
+	 * @param minimumOperator A flag indicating whether this operator is a minimizing or a
+	 * maximizing operator.
+	 */
+	NoBoundOperator(AbstractPathFormula<T>* pathFormula, bool minimumOperator)
+		: optimalityOperator(true), minimumOperator(minimumOperator) {
 		this->pathFormula = pathFormula;
 	}
 
@@ -142,8 +154,34 @@ public:
 		return checker.conforms(this->pathFormula);
 	}
 
+	/*!
+	 * Retrieves whether the operator is to be interpreted as an optimizing (i.e. min/max) operator.
+	 * @returns True if the operator is an optimizing operator.
+	 */
+	bool isOptimalityOperator() const {
+		return optimalityOperator;
+	}
+
+	/*!
+	 * Retrieves whether the operator is a minimizing operator given that it is an optimality
+	 * operator.
+	 * @returns True if the operator is an optimizing operator and it is a minimizing operator and
+	 * false otherwise, i.e. if it is either not an optimizing operator or not a minimizing operator.
+	 */
+	bool isMinimumOperator() const {
+		return optimalityOperator && minimumOperator;
+	}
+
 private:
 	AbstractPathFormula<T>* pathFormula;
+
+	// A flag that indicates whether this operator is meant as an optimizing (i.e. min/max) operator
+	// over a nondeterministic model.
+	bool optimalityOperator;
+
+	// In the case this operator is an optimizing operator, this flag indicates whether it is
+	// looking for the minimum or the maximum value.
+	bool minimumOperator;
 };
 
 } /* namespace formula */

src/formula/ProbabilisticNoBoundOperator.h

@@ -51,7 +51,7 @@ public:
 	/*!
 	 * Empty constructor
 	 */
-	ProbabilisticNoBoundOperator() : NoBoundOperator<T>(nullptr), optimalityOperator(false), minimumOperator(false) {
+	ProbabilisticNoBoundOperator() : NoBoundOperator<T>(nullptr) {
 		// Intentionally left empty
 	}
 
@@ -60,18 +60,7 @@ public:
 	 *
 	 * @param pathFormula The child node.
 	 */
-	ProbabilisticNoBoundOperator(AbstractPathFormula<T>* pathFormula) : NoBoundOperator<T>(pathFormula),
-			optimalityOperator(false), minimumOperator(false) {
-		// Intentionally left empty
-	}
-
-	/*!
-	 * Constructor
-	 *
-	 * @param pathFormula The child node.
-	 */
-	ProbabilisticNoBoundOperator(AbstractPathFormula<T>* pathFormula, bool minimumOperator) : NoBoundOperator<T>(pathFormula),
-			optimalityOperator(true), minimumOperator(minimumOperator) {
+	ProbabilisticNoBoundOperator(AbstractPathFormula<T>* pathFormula) : NoBoundOperator<T>(pathFormula) {
 		// Intentionally left empty
 	}
 
@@ -80,8 +69,8 @@ public:
 	 */
 	virtual std::string toString() const {
 		std::string result = "P";
-		if (optimalityOperator) {
-			if (minimumOperator) {
+		if (this->isOptimalityOperator()) {
+			if (this->isMinimumOperator()) {
 				result += "min";
 			} else {
 				result += "max";
@@ -92,15 +81,6 @@ public:
 		result += "]";
 		return result;
 	}
-
-private:
-	// A flag that indicates whether this operator is meant as an optimizing (i.e. min/max) operator
-	// over a nondeterministic model.
-	bool optimalityOperator;
-
-	// In the case this operator is an optimizing operator, this flag indicates whether it is
-	// looking for the minimum or the maximum value.
-	bool minimumOperator;
 };
 
 } /* namespace formula */

src/modelchecker/DtmcPrctlModelChecker.h

@@ -178,6 +178,7 @@ public:
 		}
 
 		if (!this->getModel().hasAtomicProposition(formula.getAp())) {
+			LOG4CPLUS_ERROR(logger, "Atomic proposition '" << formula.getAp() << "' is invalid.");
 			throw storm::exceptions::InvalidPropertyException() << "Atomic proposition '" << formula.getAp() << "' is invalid.";
 			return nullptr;
 		}
@@ -246,6 +247,10 @@ public:
 	 * @returns The set of states satisfying the formula, represented by a bit vector
 	 */
 	std::vector<Type>* checkNoBoundOperator(const storm::formula::NoBoundOperator<Type>& formula) const {
+		// Check if the operator was an optimality operator and report a warning in that case.
+		if (formula.isOptimalityOperator()) {
+			LOG4CPLUS_WARN(logger, "Formula contains min/max operator which is not meaningful over deterministic models.");
+		}
 		return formula.getPathFormula().check(*this);
 	}
 

src/modelchecker/GmmxxDtmcPrctlModelChecker.h

@@ -329,7 +329,7 @@ public:
 	 * @return The name of this module.
 	 */
 	static std::string getModuleName() {
-		return "gmm++";
+		return "gmm++det";
 	}
 
 	/*!

src/modelchecker/GmmxxMdpPrctlModelChecker.h

@@ -0,0 +1,409 @@
+/*
+ * GmmxxDtmcPrctlModelChecker.h
+ *
+ *  Created on: 06.12.2012
+ *      Author: Christian Dehnert
+ */
+
+#ifndef STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_
+
+#include <cmath>
+
+#include "src/models/Mdp.h"
+#include "src/modelchecker/MdpPrctlModelChecker.h"
+#include "src/utility/GraphAnalyzer.h"
+#include "src/utility/Vector.h"
+#include "src/utility/ConstTemplates.h"
+#include "src/utility/Settings.h"
+#include "src/adapters/GmmxxAdapter.h"
+#include "src/exceptions/InvalidPropertyException.h"
+#include "src/storage/JacobiDecomposition.h"
+
+#include "gmm/gmm_matrix.h"
+#include "gmm/gmm_iter_solvers.h"
+
+#include "log4cplus/logger.h"
+#include "log4cplus/loggingmacros.h"
+
+extern log4cplus::Logger logger;
+
+namespace storm {
+
+namespace modelChecker {
+
+/*
+ * A model checking engine that makes use of the gmm++ backend.
+ */
+template <class Type>
+class GmmxxDtmcPrctlModelChecker : public MdpPrctlModelChecker<Type> {
+
+public:
+	explicit GmmxxDtmcPrctlModelChecker(storm::models::Mdp<Type>& mdp) : MdpPrctlModelChecker<Type>(mdp) { }
+
+	virtual ~GmmxxDtmcPrctlModelChecker() { }
+
+	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.
+		storm::storage::BitVector* leftStates = this->checkStateFormula(formula.getLeft());
+		storm::storage::BitVector* rightStates = this->checkStateFormula(formula.getRight());
+
+		// Copy the matrix before we make any changes.
+		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionMatrix());
+
+		// Get the starting row indices for the non-deterministic choices to reduce the resulting
+		// vector properly.
+		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+
+		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
+		tmpMatrix.makeRowsAbsorbing(~(*leftStates | *rightStates) | *rightStates);
+
+		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(tmpMatrix);
+
+		// Create the vector with which to multiply.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		storm::utility::setVectorValues(result, *rightStates, storm::utility::constGetOne<Type>());
+
+		// 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());
+
+		// 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()) {
+				storm::utility::reduceVectorMin(*multiplyResult, result, *nondeterministicChoiceIndices);
+			} else {
+				storm::utility::reduceVectorMax(*multiplyResult, result, *nondeterministicChoiceIndices);
+			}
+		}
+		delete multiplyResult;
+
+		// Delete intermediate results and return result.
+		delete gmmxxMatrix;
+		delete leftStates;
+		delete rightStates;
+		return result;
+	}
+
+	virtual std::vector<Type>* checkNext(const storm::formula::Next<Type>& formula) const {
+		// First, we need to compute the states that satisfy the sub-formula of the next-formula.
+		storm::storage::BitVector* nextStates = this->checkStateFormula(formula.getChild());
+
+		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*this->getModel().getTransitionMatrix());
+
+		// Create the vector with which to multiply and initialize it correctly.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		storm::utility::setVectorValues(result, *nextStates, storm::utility::constGetOne<Type>());
+
+		// Delete obsolete sub-result.
+		delete nextStates;
+
+		// Create resulting vector.
+		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);
+
+		// Get the starting row indices for the non-deterministic choices to reduce the resulting
+		// vector properly.
+		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+
+		if (minimumOperatorStack.top()) {
+			storm::utility::reduceVectorMin(*temporaryResult, result, *nondeterministicChoiceIndices);
+		} else {
+			storm::utility::reduceVectorMax(*temporaryResult, result, *nondeterministicChoiceIndices);
+		}
+
+		// Delete temporary matrix plus temporary result and return result.
+		delete gmmxxMatrix;
+		delete temporaryResult;
+		return result;
+	}
+
+	virtual std::vector<Type>* checkUntil(const storm::formula::Until<Type>& formula) const {
+		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
+		storm::storage::BitVector* leftStates = this->checkStateFormula(formula.getLeft());
+		storm::storage::BitVector* rightStates = this->checkStateFormula(formula.getRight());
+
+		// Then, we need to identify the states which have to be taken out of the matrix, i.e.
+		// 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()) {
+			storm::utility::GraphAnalyzer::performProb01Min(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
+		} else {
+			storm::utility::GraphAnalyzer::performProb01Max(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
+		}
+
+		// Delete sub-results that are obsolete now.
+		delete leftStates;
+		delete rightStates;
+
+		LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
+		LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
+		storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
+		LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+
+		// Create resulting vector.
+		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);
+
+			// Transform the submatrix to the gmm++ format to use its solvers.
+			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));
+
+			// 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);
+
+			// Solve the corresponding system of linear equations.
+			this->solveEquationSystem(*gmmxxMatrix, x, b);
+
+			// Set values of resulting vector according to result.
+			storm::utility::setVectorValues<Type>(result, maybeStates, x);
+
+			// Delete temporary matrix.
+			delete gmmxxMatrix;
+		}
+
+		// Set values of resulting vector that are known exactly.
+		storm::utility::setVectorValues<Type>(result, statesWithProbability0, storm::utility::constGetZero<Type>());
+		storm::utility::setVectorValues<Type>(result, statesWithProbability1, storm::utility::constGetOne<Type>());
+
+		return result;
+	}
+
+	virtual std::vector<Type>* checkInstantaneousReward(const storm::formula::InstantaneousReward<Type>& formula) const {
+		// Only compute the result if the model has a state-based reward model.
+		if (!this->getModel().hasStateRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing (state-based) reward model for formula.");
+			throw storm::exceptions::InvalidPropertyException() << "Missing (state-based) reward model for formula.";
+		}
+
+		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*this->getModel().getTransitionMatrix());
+
+		// Initialize result to state rewards of the model.
+		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+
+		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
+		std::vector<Type>* swap = nullptr;
+		std::vector<Type>* tmpResult = new std::vector<Type>(this->getModel().getNumberOfStates());
+		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
+			gmm::mult(*gmmxxMatrix, *result, *tmpResult);
+			swap = tmpResult;
+			tmpResult = result;
+			result = swap;
+		}
+
+		// Delete temporary variables and return result.
+		delete tmpResult;
+		delete gmmxxMatrix;
+		return result;
+	}
+
+	virtual std::vector<Type>* checkCumulativeReward(const storm::formula::CumulativeReward<Type>& formula) const {
+		// Only compute the result if the model has at least one reward model.
+		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing reward model for formula.");
+			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
+		}
+
+		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*this->getModel().getTransitionMatrix());
+
+		// Compute the reward vector to add in each step based on the available reward models.
+		std::vector<Type>* totalRewardVector = nullptr;
+		if (this->getModel().hasTransitionRewards()) {
+			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+			if (this->getModel().hasStateRewards()) {
+				gmm::add(*this->getModel().getStateRewardVector(), *totalRewardVector);
+			}
+		} else {
+			totalRewardVector = new std::vector<Type>(*this->getModel().getStateRewardVector());
+		}
+
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+
+		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
+		std::vector<Type>* swap = nullptr;
+		std::vector<Type>* tmpResult = new std::vector<Type>(this->getModel().getNumberOfStates());
+		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
+			gmm::mult(*gmmxxMatrix, *result, *tmpResult);
+			swap = tmpResult;
+			tmpResult = result;
+			result = swap;
+
+			// Add the reward vector to the result.
+			gmm::add(*totalRewardVector, *result);
+		}
+
+		// Delete temporary variables and return result.
+		delete tmpResult;
+		delete gmmxxMatrix;
+		delete totalRewardVector;
+		return result;
+	}
+
+	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula) const {
+		// Only compute the result if the model has at least one reward model.
+		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula");
+			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
+		}
+
+		// Determine the states for which the target predicate holds.
+		storm::storage::BitVector* targetStates = this->checkStateFormula(formula.getChild());
+
+		// 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);
+		infinityStates.complement();
+
+		// Create resulting vector.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+
+		// Check whether there are states for which we have to compute the result.
+		storm::storage::BitVector maybeStates = ~(*targetStates) & ~infinityStates;
+		const int 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);
+			// Converting the matrix from the fixpoint notation to the form needed for the equation
+			// system. That is, we go from x = A*x + b to (I-A)x = b.
+			submatrix->convertToEquationSystem();
+
+			// Transform the submatrix to the gmm++ format to use its solvers.
+			gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*submatrix);
+			delete submatrix;
+
+			// Initialize the x vector with 1 for each element. This is the initial guess for
+			// the iterative solvers.
+			std::vector<Type> x(maybeStatesSetBitCount, storm::utility::constGetOne<Type>());
+
+			// Prepare the right-hand side of the equation system.
+			std::vector<Type>* b = new std::vector<Type>(maybeStatesSetBitCount);
+			if (this->getModel().hasTransitionRewards()) {
+				// If a transition-based reward model is available, we initialize the right-hand
+				// side to the vector resulting from summing the rows of the pointwise product
+				// of the transition probability matrix and the transition reward matrix.
+				std::vector<Type>* pointwiseProductRowSumVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+				storm::utility::selectVectorValues(b, maybeStates, *pointwiseProductRowSumVector);
+				delete pointwiseProductRowSumVector;
+
+				if (this->getModel().hasStateRewards()) {
+					// If a state-based reward model is also available, we need to add this vector
+					// as well. As the state reward vector contains entries not just for the states
+					// that we still consider (i.e. maybeStates), we need to extract these values
+					// first.
+					std::vector<Type>* subStateRewards = new std::vector<Type>(maybeStatesSetBitCount);
+					storm::utility::setVectorValues(subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
+					gmm::add(*subStateRewards, *b);
+					delete subStateRewards;
+				}
+			} else {
+				// If only a state-based reward model is  available, we take this vector as the
+				// right-hand side. As the state reward vector contains entries not just for the
+				// states that we still consider (i.e. maybeStates), we need to extract these values
+				// first.
+				storm::utility::setVectorValues(b, maybeStates, *this->getModel().getStateRewardVector());
+			}
+
+			// Solve the corresponding system of linear equations.
+			this->solveLinearEquationSystem(*gmmxxMatrix, x, *b);
+
+			// Set values of resulting vector according to result.
+			storm::utility::setVectorValues<Type>(result, maybeStates, x);
+
+			// Delete temporary matrix and right-hand side.
+			delete gmmxxMatrix;
+			delete b;
+		}
+
+		// Set values of resulting vector that are known exactly.
+		storm::utility::setVectorValues(result, *targetStates, storm::utility::constGetZero<Type>());
+		storm::utility::setVectorValues(result, infinityStates, storm::utility::constGetInfinity<Type>());
+
+		// Delete temporary storages and return result.
+		delete targetStates;
+		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.
+	 *
+	 * @param A The matrix A specifying the coefficients of the equations.
+	 * @param x The vector x for which to solve the equations. The initial value of the elements of
+	 * this vector are used as the initial guess and might thus influence performance and convergence.
+	 * @param b The vector b specifying the values on the right-hand-sides of the equations.
+	 * @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.
+		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");
+
+
+
+
+
+		// Check if the solver converged and issue a warning otherwise.
+		if (iter.converged()) {
+			LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations.");
+		} else {
+			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+		}
+
+	}
+};
+
+} //namespace modelChecker
+
+} //namespace storm
+
+#endif /* STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_ */

src/modelchecker/MdpPrctlModelChecker.h

@@ -0,0 +1,375 @@
+/*
+ * MdpPrctlModelChecker.h
+ *
+ *  Created on: 15.02.2013
+ *      Author: Christian Dehnert
+ */
+
+#ifndef STORM_MODELCHECKER_MDPPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_MDPPRCTLMODELCHECKER_H_
+
+#include "src/formula/Formulas.h"
+#include "src/utility/Vector.h"
+#include "src/storage/SparseMatrix.h"
+
+#include "src/models/Mdp.h"
+#include "src/storage/BitVector.h"
+#include "src/exceptions/InvalidPropertyException.h"
+#include "src/utility/Vector.h"
+#include "src/modelchecker/AbstractModelChecker.h"
+#include <vector>
+
+#include "log4cplus/logger.h"
+#include "log4cplus/loggingmacros.h"
+
+extern log4cplus::Logger logger;
+
+namespace storm {
+
+namespace modelChecker {
+
+/*!
+ * @brief
+ * Interface for model checker classes.
+ *
+ * This class provides basic functions that are the same for all subclasses, but mainly only declares
+ * abstract methods that are to be implemented in concrete instances.
+ *
+ * @attention This class is abstract.
+ */
+template<class Type>
+class MdpPrctlModelChecker :
+	public virtual AbstractModelChecker<Type> {
+public:
+	/*!
+	 * Constructor
+	 *
+	 * @param model The dtmc model which is checked.
+	 */
+	explicit MdpPrctlModelChecker(storm::models::Mdp<Type>& model) : model(model), minimumOperatorStack() {
+
+	}
+
+	/*!
+	 * Copy constructor
+	 *
+	 * @param modelChecker The model checker that is copied.
+	 */
+	explicit MdpPrctlModelChecker(const storm::modelChecker::MdpPrctlModelChecker<Type>* modelChecker) : model(new storm::models::Mdp<Type>(modelChecker->getModel())),  minimumOperatorStack() {
+
+	}
+
+	/*!
+	 * Destructor
+	 */
+	virtual ~MdpPrctlModelChecker() {
+		// Intentionally left empty.
+	}
+
+	/*!
+	 * @returns A reference to the dtmc of the model checker.
+	 */
+	storm::models::Mdp<Type>& getModel() const {
+		return this->model;
+	}
+
+	/*!
+	 * Sets the DTMC model which is checked
+	 * @param model
+	 */
+	void setModel(storm::models::Mdp<Type>& model) {
+		this->model = &model;
+	}
+
+	/*!
+	 * Checks the given state formula on the DTMC and prints the result (true/false) for all initial
+	 * states.
+	 * @param stateFormula The formula to be checked.
+	 */
+	void check(const storm::formula::AbstractStateFormula<Type>& stateFormula) const {
+		std::cout << std::endl;
+		LOG4CPLUS_INFO(logger, "Model checking formula\t" << stateFormula.toString());
+		std::cout << "Model checking formula:\t" << stateFormula.toString() << std::endl;
+		storm::storage::BitVector* result = nullptr;
+		try {
+			result = stateFormula.check(*this);
+			LOG4CPLUS_INFO(logger, "Result for initial states:");
+			std::cout << "Result for initial states:" << std::endl;
+			for (auto initialState : *this->getModel().getLabeledStates("init")) {
+				LOG4CPLUS_INFO(logger, "\t" << initialState << ": " << (result->get(initialState) ? "satisfied" : "not satisfied"));
+				std::cout << "\t" << initialState << ": " << (*result)[initialState] << std::endl;
+			}
+			delete result;
+		} catch (std::exception& e) {
+			std::cout << "Error during computation: " << e.what() << "Skipping property." << std::endl;
+			if (result != nullptr) {
+				delete result;
+			}
+		}
+		std::cout << std::endl;
+		storm::utility::printSeparationLine(std::cout);
+	}
+
+	/*!
+	 * Checks the given operator (with no bound) on the DTMC and prints the result
+	 * (probability/rewards) for all initial states.
+	 * @param noBoundFormula The formula to be checked.
+	 */
+	void check(const storm::formula::NoBoundOperator<Type>& noBoundFormula) const {
+		std::cout << std::endl;
+		LOG4CPLUS_INFO(logger, "Model checking formula\t" << noBoundFormula.toString());
+		std::cout << "Model checking formula:\t" << noBoundFormula.toString() << std::endl;
+		std::vector<Type>* result = nullptr;
+		try {
+			result = noBoundFormula.check(*this);
+			LOG4CPLUS_INFO(logger, "Result for initial states:");
+			std::cout << "Result for initial states:" << std::endl;
+			for (auto initialState : *this->getModel().getLabeledStates("init")) {
+				LOG4CPLUS_INFO(logger, "\t" << initialState << ": " << (*result)[initialState]);
+				std::cout << "\t" << initialState << ": " << (*result)[initialState] << std::endl;
+			}
+			delete result;
+		} catch (std::exception& e) {
+			std::cout << "Error during computation: " << e.what() << " Skipping property." << std::endl;
+			if (result != nullptr) {
+				delete result;
+			}
+		}
+		std::cout << std::endl;
+		storm::utility::printSeparationLine(std::cout);
+	}
+
+	/*!
+	 * The check method for a state formula; Will infer the actual type of formula and delegate it
+	 * to the specialized method
+	 *
+	 * @param formula The state formula to check
+	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 */
+	storm::storage::BitVector* checkStateFormula(const storm::formula::AbstractStateFormula<Type>& formula) const {
+		return formula.check(*this);
+	}
+
+	/*!
+	 * The check method for a state formula with an And node as root in its formula tree
+	 *
+	 * @param formula The And formula to check
+	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 */
+	storm::storage::BitVector* checkAnd(const storm::formula::And<Type>& formula) const {
+		storm::storage::BitVector* result = checkStateFormula(formula.getLeft());
+		storm::storage::BitVector* right = checkStateFormula(formula.getRight());
+		(*result) &= (*right);
+		delete right;
+		return result;
+	}
+
+	/*!
+	 * The check method for a formula with an AP node as root in its formula tree
+	 *
+	 * @param formula The Ap state formula to check
+	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 */
+	storm::storage::BitVector* checkAp(const storm::formula::Ap<Type>& formula) const {
+		if (formula.getAp().compare("true") == 0) {
+			return new storm::storage::BitVector(this->getModel().getNumberOfStates(), true);
+		} else if (formula.getAp().compare("false") == 0) {
+			return new storm::storage::BitVector(this->getModel().getNumberOfStates());
+		}
+
+		if (!this->getModel().hasAtomicProposition(formula.getAp())) {
+			LOG4CPLUS_ERROR(logger, "Atomic proposition '" << formula.getAp() << "' is invalid.");
+			throw storm::exceptions::InvalidPropertyException() << "Atomic proposition '" << formula.getAp() << "' is invalid.";
+			return nullptr;
+		}
+
+		return new storm::storage::BitVector(*this->getModel().getLabeledStates(formula.getAp()));
+	}
+
+	/*!
+	 * The check method for a formula with a Not node as root in its formula tree
+	 *
+	 * @param formula The Not state formula to check
+	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 */
+	storm::storage::BitVector* checkNot(const storm::formula::Not<Type>& formula) const {
+		storm::storage::BitVector* result = checkStateFormula(formula.getChild());
+		result->complement();
+		return result;
+	}
+
+	/*!
+	 * The check method for a state formula with an Or node as root in its formula tree
+	 *
+	 * @param formula The Or state formula to check
+	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 */
+	virtual storm::storage::BitVector* checkOr(const storm::formula::Or<Type>& formula) const {
+		storm::storage::BitVector* result = checkStateFormula(formula.getLeft());
+		storm::storage::BitVector* right = checkStateFormula(formula.getRight());
+		(*result) |= (*right);
+		delete right;
+		return result;
+	}
+
+	/*!
+	 * The check method for a state formula with a bound operator node as root in
+	 * its formula tree
+	 *
+	 * @param formula The state formula to check
+	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 */
+	storm::storage::BitVector* checkPathBoundOperator(const storm::formula::PathBoundOperator<Type>& formula) const {
+		// First, we need to compute the probability for satisfying the path formula for each state.
+		std::vector<Type>* quantitativeResult = this->checkPathFormula(formula.getPathFormula());
+
+		// Create resulting bit vector, which will hold the yes/no-answer for every state.
+		storm::storage::BitVector* result = new storm::storage::BitVector(this->getModel().getNumberOfStates());
+
+		// Now, we can compute which states meet the bound specified in this operator and set the
+		// corresponding bits to true in the resulting vector.
+		for (uint_fast64_t i = 0; i < this->getModel().getNumberOfStates(); ++i) {
+			if (formula.meetsBound((*quantitativeResult)[i])) {
+				result->set(i, true);
+			}
+		}
+
+		// Delete the probabilities computed for the states and return result.
+		delete quantitativeResult;
+		return result;
+	}
+
+	/*!
+	 * The check method for a state formula with a probabilistic operator node without bounds as root
+	 * in its formula tree
+	 *
+	 * @param formula The state formula to check
+	 * @returns The set of states satisfying the formula, represented by a bit vector
+	 */
+	std::vector<Type>* checkNoBoundOperator(const storm::formula::NoBoundOperator<Type>& formula) const {
+		// Check if the operator was an non-optimality operator and report an error in that case.
+		if (!formula.isOptimalityOperator()) {
+			LOG4CPLUS_ERROR(logger, "Formula does not specify neither min nor max optimality, which is not meaningful over nondeterministic models.");
+			throw storm::exceptions::InvalidArgumentException() << "Formula does not specify neither min nor max optimality, which is not meaningful over nondeterministic models.";
+		}
+		minimumOperatorStack.push(formula.isMinimumOperator());
+		std::vector<Type>* result = formula.getPathFormula().check(*this);
+		minimumOperatorStack.pop();
+		return result;
+	}
+
+	/*!
+	 * The check method for a path formula; Will infer the actual type of formula and delegate it
+	 * to the specialized method
+	 *
+	 * @param formula The path formula to check
+	 * @returns for each state the probability that the path formula holds.
+	 */
+	std::vector<Type>* checkPathFormula(const storm::formula::AbstractPathFormula<Type>& formula) const {
+		return formula.check(*this);
+	}
+
+	/*!
+	 * The check method for a path formula with a Bounded Until operator node as root in its formula tree
+	 *
+	 * @param formula The Bounded Until path formula to check
+	 * @returns for each state the probability that the path formula holds.
+	 */
+	virtual std::vector<Type>* checkBoundedUntil(const storm::formula::BoundedUntil<Type>& formula) const = 0;
+
+	/*!
+	 * The check method for a path formula with a Next operator node as root in its formula tree
+	 *
+	 * @param formula The Next path formula to check
+	 * @returns for each state the probability that the path formula holds.
+	 */
+	virtual std::vector<Type>* checkNext(const storm::formula::Next<Type>& formula) const = 0;
+
+	/*!
+	 * The check method for a path formula with a Bounded Eventually operator node as root in its
+	 * formula tree
+	 *
+	 * @param formula The Bounded Eventually path formula to check
+	 * @returns for each state the probability that the path formula holds
+	 */
+	virtual std::vector<Type>* checkBoundedEventually(const storm::formula::BoundedEventually<Type>& formula) const {
+		// Create equivalent temporary bounded until formula and check it.
+		storm::formula::BoundedUntil<Type> temporaryBoundedUntilFormula(new storm::formula::Ap<Type>("true"), formula.getChild().clone(), formula.getBound());
+		return this->checkBoundedUntil(temporaryBoundedUntilFormula);
+	}
+
+	/*!
+	 * The check method for a path formula with an Eventually operator node as root in its formula tree
+	 *
+	 * @param formula The Eventually path formula to check
+	 * @returns for each state the probability that the path formula holds
+	 */
+	virtual std::vector<Type>* checkEventually(const storm::formula::Eventually<Type>& formula) const {
+		// Create equivalent temporary until formula and check it.
+		storm::formula::Until<Type> temporaryUntilFormula(new storm::formula::Ap<Type>("true"), formula.getChild().clone());
+		return this->checkUntil(temporaryUntilFormula);
+	}
+
+	/*!
+	 * The check method for a path formula with a Globally operator node as root in its formula tree
+	 *
+	 * @param formula The Globally path formula to check
+	 * @returns for each state the probability that the path formula holds
+	 */
+	virtual std::vector<Type>* checkGlobally(const storm::formula::Globally<Type>& formula) const {
+		// Create "equivalent" temporary eventually formula and check it.
+		storm::formula::Eventually<Type> temporaryEventuallyFormula(new storm::formula::Not<Type>(formula.getChild().clone()));
+		std::vector<Type>* result = this->checkEventually(temporaryEventuallyFormula);
+
+		// Now subtract the resulting vector from the constant one vector to obtain final result.
+		storm::utility::subtractFromConstantOneVector(result);
+		return result;
+	}
+
+	/*!
+	 * The check method for a path formula with an Until operator node as root in its formula tree
+	 *
+	 * @param formula The Until path formula to check
+	 * @returns for each state the probability that the path formula holds.
+	 */
+	virtual std::vector<Type>* checkUntil(const storm::formula::Until<Type>& formula) const = 0;
+
+	/*!
+	 * The check method for a path formula with an Instantaneous Reward operator node as root in its
+	 * formula tree
+	 *
+	 * @param formula The Instantaneous Reward formula to check
+	 * @returns for each state the reward that the instantaneous reward yields
+	 */
+	virtual std::vector<Type>* checkInstantaneousReward(const storm::formula::InstantaneousReward<Type>& formula) const = 0;
+
+	/*!
+	 * The check method for a path formula with a Cumulative Reward operator node as root in its
+	 * formula tree
+	 *
+	 * @param formula The Cumulative Reward formula to check
+	 * @returns for each state the reward that the cumulative reward yields
+	 */
+	virtual std::vector<Type>* checkCumulativeReward(const storm::formula::CumulativeReward<Type>& formula) const = 0;
+
+	/*!
+	 * The check method for a path formula with a Reachability Reward operator node as root in its
+	 * formula tree
+	 *
+	 * @param formula The Reachbility Reward formula to check
+	 * @returns for each state the reward that the reachability reward yields
+	 */
+	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula) const = 0;
+
+protected:
+	std::stack<bool> minimumOperatorStack;
+
+private:
+	storm::models::Mdp<Type>& model;
+};
+
+} //namespace modelChecker
+
+} //namespace storm
+
+#endif /* STORM_MODELCHECKER_DTMCPRCTLMODELCHECKER_H_ */

src/utility/Vector.h

@@ -53,6 +53,48 @@ 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) {
+			++currentTargetRow;
+			(*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;
+		}
+	}
+}
+
+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) {
+			++currentTargetRow;
+			(*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;
+		}
+	}
+}
+
 } //namespace utility
 
 } //namespace storm