Merge branch 'master' of https://sselab.de/lab9/private/git/storm

src/modelchecker/MdpPrctlModelChecker.h

@@ -545,6 +545,7 @@ private:
 		while (!converged && iterations < maxIterations) {
 			// Compute x' = A*x + b.
 			matrix.multiplyWithVector(*currentX, multiplyResult);
+			// matrix.multiplyAddAndReduceInPlace(nondeterministicChoiceIndices, *currentX, b, this->minimumOperatorStack.top());
 
 			gmm::add(b, multiplyResult);
 
@@ -558,11 +559,15 @@ private:
 			// Determine whether the method converged.
 			converged = storm::utility::equalModuloPrecision(*currentX, *newX, precision, relative);
 
+
 			// Update environment variables.
 			swap = currentX;
 			currentX = newX;
 			newX = swap;
 			++iterations;
+
+			// *newX = *currentX,
+
 		}
 
 		if (iterations % 2 == 1) {

src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h

@@ -0,0 +1,139 @@
+/*
+ * GmmxxDtmcPrctlModelChecker.h
+ *
+ *  Created on: 06.12.2012
+ *      Author: Christian Dehnert
+ */
+
+#ifndef STORM_MODELCHECKER_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_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 TopologicalValueIterationMdpPrctlModelChecker : public MdpPrctlModelChecker<Type> {
+
+public:
+	explicit TopologicalValueIterationMdpPrctlModelChecker(storm::models::Mdp<Type>& mdp) : MdpPrctlModelChecker<Type>(mdp) { }
+
+	virtual ~TopologicalValueIterationMdpPrctlModelChecker() { }
+
+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 solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, 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 precision = s->get<double>("precision");
+		unsigned maxIterations = s->get<unsigned>("maxiter");
+		bool relative = s->get<bool>("relative");
+
+		std::vector<std::vector<uint_fast64_t>> stronglyConnectedComponents;
+		storm::models::GraphTransitions<Type> stronglyConnectedComponentsDependencyGraph;
+		storm::utility::GraphAnalyzer::performSccDecomposition(matrix, nondeterministicChoiceIndices, stronglyConnectedComponents, stronglyConnectedComponentsDependencyGraph);
+
+		std::vector<uint_fast64_t> topologicalSort;
+		storm::utility::GraphAnalyzer::getTopologicalSort(stronglyConnectedComponentsDependencyGraph, topologicalSort);
+
+		// Set up the environment for the power method.
+		std::vector<Type> multiplyResult(matrix.getRowCount());
+		std::vector<Type>* currentX = &x;
+		std::vector<Type>* newX = new std::vector<Type>(x.size());
+		std::vector<Type>* swap = nullptr;
+		uint_fast64_t currentMaxLocalIterations = 0;
+		uint_fast64_t localIterations = 0;
+		uint_fast64_t globalIterations = 0;
+		bool converged = true;
+
+		for (auto sccIndexIt = topologicalSort.begin(); sccIndexIt != topologicalSort.end() && converged; ++sccIndexIt) {
+			std::vector<uint_fast64_t> const& scc = stronglyConnectedComponents[*sccIndexIt];
+
+			localIterations = 0;
+			converged = false;
+			while (!converged && localIterations < maxIterations) {
+				// Compute x' = A*x + b.
+				matrix.multiplyWithVector(scc, nondeterministicChoiceIndices, *currentX, multiplyResult);
+				storm::utility::addVectors(scc, nondeterministicChoiceIndices, multiplyResult, b);
+
+				// Reduce the vector x' by applying min/max for all non-deterministic choices.
+				if (this->minimumOperatorStack.top()) {
+					storm::utility::reduceVectorMin(multiplyResult, newX, scc, nondeterministicChoiceIndices);
+				} else {
+					storm::utility::reduceVectorMax(multiplyResult, newX, scc, nondeterministicChoiceIndices);
+				}
+
+				// Determine whether the method converged.
+				// converged = storm::utility::equalModuloPrecision(*currentX, *newX, scc, precision, relative);
+				converged = storm::utility::equalModuloPrecision(*currentX, *newX, precision, relative);
+
+				// Update environment variables.
+				swap = currentX;
+				currentX = newX;
+				newX = swap;
+				++localIterations;
+				++globalIterations;
+			}
+
+			std::cout << "converged locally for scc of size " << scc.size() << std::endl;
+
+			if (localIterations > currentMaxLocalIterations) {
+				currentMaxLocalIterations = localIterations;
+			}
+		}
+
+		if (globalIterations % 2 == 1) {
+			std::swap(x, *currentX);
+			delete currentX;
+		} else {
+			delete newX;
+		}
+
+		// Check if the solver converged and issue a warning otherwise.
+		if (converged) {
+			LOG4CPLUS_INFO(logger, "Iterative solver converged after " << currentMaxLocalIterations << " iterations.");
+		} else {
+			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
+		}
+	}
+};
+
+} //namespace modelChecker
+
+} //namespace storm
+
+#endif /* STORM_MODELCHECKER_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_ */

src/models/GraphTransitions.h

@@ -99,13 +99,21 @@ public:
 		return result;
 	}
 
+	uint_fast64_t getNumberOfStates() const {
+		return numberOfStates;
+	}
+
+	uint_fast64_t getNumberOfTransitions() const {
+		return numberOfTransitions;
+	}
+
 	/*!
 	 * Returns an iterator to the successors of the given state.
 	 * @param state The state for which to get the successor iterator.
 	 * @return An iterator to the predecessors of the given states.
 	 */
 	stateSuccessorIterator beginStateSuccessorsIterator(uint_fast64_t state) const {
-		return &this->successorList[this->stateIndications[state]];
+		return &(this->successorList[0]) + this->stateIndications[state];
 	}
 
 	/*!
@@ -116,12 +124,7 @@ public:
 	 * the given state.
 	 */
 	stateSuccessorIterator endStateSuccessorsIterator(uint_fast64_t state) const {
-		std::cout << "size: " << this->stateIndications.size() << std::endl;
-		std::cout << "idx accessed " << state << std::endl;
-		std::cout << this->stateIndications[state + 1] << std::endl;
-		std::cout << "other size: " << this->successorList.size() << std::endl;
-		std::cout << this->successorList[this->stateIndications[state + 1]] << std::endl;
-		return &this->successorList[this->stateIndications[state + 1]];
+		return &(this->successorList[0]) + this->stateIndications[state + 1];
 	}
 
 	/*!
@@ -130,8 +133,11 @@ public:
 	 */
 	std::string toString() const {
 		std::stringstream stream;
-		stream << successorList << std::endl;
-		stream << stateIndications << std::endl;
+		for (uint_fast64_t state = 0; state < numberOfStates; ++state) {
+			for (auto succIt = this->beginStateSuccessorsIterator(state), succIte = this->endStateSuccessorsIterator(state); succIt != succIte; ++succIt) {
+				stream << state << " -> " << *succIt << std::endl;
+			}
+		}
 		return stream.str();
 	}
 
@@ -147,7 +153,7 @@ private:
 			// Now, we determine the SCCs which are reachable (in one step) from the current SCC.
 			std::set<uint_fast64_t> allTargetSccs;
 			for (auto state : scc) {
-				for (stateSuccessorIterator succIt = beginStateSuccessorsIterator(state), succIte = endStateSuccessorsIterator(state); succIt != succIte; ++succIt) {
+				for (stateSuccessorIterator succIt = transitions.beginStateSuccessorsIterator(state), succIte = transitions.endStateSuccessorsIterator(state); succIt != succIte; ++succIt) {
 					uint_fast64_t targetScc = stateToSccMap.find(*succIt)->second;
 
 					// We only need to consider transitions that are actually leaving the SCC.
@@ -166,8 +172,6 @@ private:
 		// Put the sentinel element at the end and initialize the number of transitions.
 		stateIndications[numberOfStates] = successorList.size();
 		numberOfTransitions = successorList.size();
-
-		std::cout << this->toString() << std::endl;
 	}
 
 	/*!

src/parser/NondeterministicSparseTransitionParser.cpp

@@ -337,7 +337,9 @@ NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(s
 		buf = trimWhitespaces(buf);
 	}
 
-	this->rowMapping->at(maxnode+1) = curRow + 1;
+	for (int_fast64_t node = lastsource + 1; node <= maxnode + 1; node++) {
+		this->rowMapping->at(node) = curRow + 1;
+	}
 
 	if (!fixDeadlocks && hadDeadlocks && !isRewardFile) throw storm::exceptions::WrongFileFormatException() << "Some of the nodes had deadlocks. You can use --fix-deadlocks to insert self-loops on the fly.";
 

src/storage/BitVector.h

@@ -137,7 +137,7 @@ public:
 	 * @param bv A reference to the bit vector to be copied.
 	 */
 	BitVector(BitVector const& bv) : bucketCount(bv.bucketCount), bitCount(bv.bitCount), endIterator(*this, bitCount, bitCount, false), truncateMask((1ll << (bitCount & mod64mask)) - 1ll) {
-		LOG4CPLUS_WARN(logger, "Invoking copy constructor.");
+		LOG4CPLUS_DEBUG(logger, "Invoking copy constructor.");
 		bucketArray = new uint64_t[bucketCount];
 		std::copy(bv.bucketArray, bv.bucketArray + this->bucketCount, this->bucketArray);
 	}

src/storage/SparseMatrix.h

@@ -830,17 +830,22 @@ public:
 		}
 	}
 
+	void multiplyWithVectorInPlace(std::vector<T>& vector) const {
+		typename std::vector<T>::iterator resultIt = vector.begin();
+		typename std::vector<T>::iterator resultIte = vector.end();
+		constRowsIterator rowIt = this->constRowsIteratorBegin();
+		uint_fast64_t nextRow = 1;
+
+		for (; resultIt != resultIte; ++resultIt, ++nextRow) {
+			*resultIt = multiplyRowWithVector(rowIt, this->rowIndications[nextRow], vector);
+		}
+	}
+
 	void multiplyWithVector(std::vector<uint_fast64_t> const& states, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<T>& vector, std::vector<T>& result) const {
 		constRowsIterator rowsIt = this->constRowsIteratorBegin();
 		uint_fast64_t nextRow = 1;
 
-		std::cout << nondeterministicChoiceIndices << std::endl;
-		std::cout << rowCount << std::endl;
-
 		for (auto stateIt = states.cbegin(), stateIte = states.cend(); stateIt != stateIte; ++stateIt) {
-			std::cout << "stateIt " << *stateIt << std::endl;
-			std::cout << nondeterministicChoiceIndices[*stateIt] << std::endl;
-			std::cout << this->rowIndications[nondeterministicChoiceIndices[*stateIt]] << std::endl;
 			rowsIt.setOffset(this->rowIndications[nondeterministicChoiceIndices[*stateIt]]);
 			nextRow = nondeterministicChoiceIndices[*stateIt] + 1;
 			for (auto rowIt = nondeterministicChoiceIndices[*stateIt], rowIte = nondeterministicChoiceIndices[*stateIt + 1]; rowIt != rowIte; ++rowIt, ++nextRow) {
@@ -849,6 +854,37 @@ public:
 		}
 	}
 
+	void multiplyWithVectorInPlace(std::vector<uint_fast64_t> const& states, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<T>& vector) const {
+		constRowsIterator rowsIt = this->constRowsIteratorBegin();
+		uint_fast64_t nextRow = 1;
+
+		for (auto stateIt = states.cbegin(), stateIte = states.cend(); stateIt != stateIte; ++stateIt) {
+			rowsIt.setOffset(this->rowIndications[nondeterministicChoiceIndices[*stateIt]]);
+			nextRow = nondeterministicChoiceIndices[*stateIt] + 1;
+			for (auto rowIt = nondeterministicChoiceIndices[*stateIt], rowIte = nondeterministicChoiceIndices[*stateIt + 1]; rowIt != rowIte; ++rowIt, ++nextRow) {
+				vector[rowIt] = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector);
+			}
+		}
+	}
+
+	void multiplyAddAndReduceInPlace(std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<T>& vector, std::vector<T> const& summand, bool minimize) const {
+		constRowsIterator rowsIt = this->constRowsIteratorBegin();
+		uint_fast64_t nextRow = 1;
+
+		for (uint_fast64_t stateIt = 0, stateIte = nondeterministicChoiceIndices.size() - 1; stateIt < stateIte; ++stateIt) {
+			vector[stateIt] = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector) + summand[nondeterministicChoiceIndices[stateIt]];
+			++nextRow;
+			for (uint_fast64_t rowIt = nondeterministicChoiceIndices[stateIt] + 1, rowIte = nondeterministicChoiceIndices[stateIt + 1]; rowIt != rowIte; ++rowIt, ++nextRow) {
+				T value = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector) + summand[rowIt];
+				if (minimize && value < vector[stateIt]) {
+					vector[stateIt] = value;
+				} else if (!minimize && value > vector[stateIt]) {
+					vector[stateIt] = value;
+				}
+			}
+		}
+	}
+
 	/*!
 	 * Returns the size of the matrix in memory measured in bytes.
 	 * @return The size of the matrix in memory measured in bytes.

src/utility/GraphAnalyzer.h

@@ -424,11 +424,72 @@ public:
 		storm::models::GraphTransitions<T> forwardTransitions(matrix, nondeterministicChoiceIndices, true);
 
 		// Perform the actual SCC decomposition based on the graph-transitions of the system.
-		performSccDecomposition(nondeterministicChoiceIndices.size(), forwardTransitions, stronglyConnectedComponents, stronglyConnectedComponentsDependencyGraph);
+		performSccDecomposition(nondeterministicChoiceIndices.size() - 1, forwardTransitions, stronglyConnectedComponents, stronglyConnectedComponentsDependencyGraph);
 
 		LOG4CPLUS_INFO(logger, "Done computing SCC decomposition.");
 	}
 
+	template <typename T>
+	static void getTopologicalSort(storm::models::GraphTransitions<T> const& transitions, std::vector<uint_fast64_t>& topologicalSort) {
+		topologicalSort.reserve(transitions.getNumberOfStates());
+
+		std::vector<uint_fast64_t> recursionStack;
+		recursionStack.reserve(transitions.getNumberOfStates());
+
+		std::vector<typename storm::models::GraphTransitions<T>::stateSuccessorIterator> iteratorRecursionStack;
+		iteratorRecursionStack.reserve(transitions.getNumberOfStates());
+
+		storm::storage::BitVector visitedStates(transitions.getNumberOfStates());
+		for (uint_fast64_t state = 0; state < transitions.getNumberOfStates(); ++state) {
+			if (!visitedStates.get(state)) {
+				recursionStack.push_back(state);
+				iteratorRecursionStack.push_back(transitions.beginStateSuccessorsIterator(state));
+
+				recursionStepForward:
+				while (!recursionStack.empty()) {
+					uint_fast64_t currentState = recursionStack.back();
+					typename storm::models::GraphTransitions<T>::stateSuccessorIterator currentIt = iteratorRecursionStack.back();
+
+					visitedStates.set(currentState, true);
+
+					recursionStepBackward:
+					for (; currentIt != transitions.endStateSuccessorsIterator(currentState); ++currentIt) {
+						if (!visitedStates.get(*currentIt)) {
+							// Put unvisited successor on top of our recursion stack and remember that.
+							recursionStack.push_back(*currentIt);
+
+							// Save current iterator position so we can continue where we left off later.
+							iteratorRecursionStack.pop_back();
+							iteratorRecursionStack.push_back(currentIt + 1);
+
+							// Also, put initial value for iterator on corresponding recursion stack.
+							iteratorRecursionStack.push_back(transitions.beginStateSuccessorsIterator(*currentIt));
+
+							goto recursionStepForward;
+						}
+					}
+
+					topologicalSort.push_back(currentState);
+
+					// If we reach this point, we have completed the recursive descent for the current state.
+					// That is, we need to pop it from the recursion stacks.
+					recursionStack.pop_back();
+					iteratorRecursionStack.pop_back();
+
+					// If there is at least one state under the current one in our recursion stack, we need
+					// to restore the topmost state as the current state and jump to the part after the
+					// original recursive call.
+					if (recursionStack.size() > 0) {
+						currentState = recursionStack.back();
+						currentIt = iteratorRecursionStack.back();
+
+						goto recursionStepBackward;
+					}
+				}
+			}
+		}
+	}
+
 private:
 	template <typename T>
 	static void performSccDecomposition(uint_fast64_t numberOfStates, storm::models::GraphTransitions<T> const& forwardTransitions, std::vector<std::vector<uint_fast64_t>>& stronglyConnectedComponents, storm::models::GraphTransitions<T>& stronglyConnectedComponentsDependencyGraph) {
@@ -482,24 +543,26 @@ private:
 			tarjanStackStates.set(currentState, true);
 
 			// Now, traverse all successors of the current state.
-			recursionStepBackward:
 			for(; currentIt != forwardTransitions.endStateSuccessorsIterator(currentState); ++currentIt) {
 				// If we have not visited the successor already, we need to perform the procedure
 				// recursively on the newly found state.
 				if (!visitedStates.get(*currentIt)) {
+					// Save current iterator position so we can continue where we left off later.
+					recursionIteratorStack.pop_back();
+					recursionIteratorStack.push_back(currentIt);
+
 					// Put unvisited successor on top of our recursion stack and remember that.
 					recursionStateStack.push_back(*currentIt);
 					statesInStack[*currentIt] = true;
 
-					// Save current iterator position so we can continue where we left off later.
-					recursionIteratorStack.pop_back();
-					recursionIteratorStack.push_back(currentIt + 1);
-
 					// Also, put initial value for iterator on corresponding recursion stack.
 					recursionIteratorStack.push_back(forwardTransitions.beginStateSuccessorsIterator(*currentIt));
 
 					// Perform the actual recursion step in an iterative way.
 					goto recursionStepForward;
+
+					recursionStepBackward:
+					lowlinks[currentState] = std::min(lowlinks[currentState], lowlinks[*currentIt]);
 				} else if (tarjanStackStates.get(*currentIt)) {
 					// Update the lowlink of the current state.
 					lowlinks[currentState] = std::min(lowlinks[currentState], stateIndices[*currentIt]);

src/utility/Vector.h

@@ -107,6 +107,21 @@ void reduceVectorMin(std::vector<T> const& source, std::vector<T>* target, std::
 	}
 }
 
+template<class T>
+void reduceVectorMin(std::vector<T> const& source, std::vector<T>* target, std::vector<uint_fast64_t> const& scc, std::vector<uint_fast64_t> const& filter) {
+	for (auto stateIt = scc.cbegin(); stateIt != scc.cend(); ++stateIt) {
+		(*target)[*stateIt] = source[filter[*stateIt]];
+
+		for (auto row = filter[*stateIt] + 1; row < filter[*stateIt + 1]; ++row) {
+			// We have to minimize the value, so only overwrite the current value if the
+			// value is actually lower.
+			if (source[row] < (*target)[*stateIt]) {
+				(*target)[*stateIt] = source[row];
+			}
+		}
+	}
+}
+
 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;
@@ -127,6 +142,21 @@ void reduceVectorMax(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& scc, std::vector<uint_fast64_t> const& filter) {
+	for (auto stateIt = scc.cbegin(); stateIt != scc.cend(); ++stateIt) {
+		(*target)[*stateIt] = source[filter[*stateIt]];
+
+		for (auto row = filter[*stateIt] + 1; row < filter[*stateIt + 1]; ++row) {
+			// We have to maximize the value, so only overwrite the current value if the
+			// value is actually lower.
+			if (source[row] > (*target)[*stateIt]) {
+				(*target)[*stateIt] = source[row];
+			}
+		}
+	}
+}
+
 template<class T>
 bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight, T precision, bool relativeError) {
 	if (vectorLeft.size() != vectorRight.size()) {
@@ -145,6 +175,24 @@ bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const
 	return true;
 }
 
+template<class T>
+bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight, std::vector<uint_fast64_t> const& scc, T precision, bool relativeError) {
+	if (vectorLeft.size() != vectorRight.size()) {
+		LOG4CPLUS_ERROR(logger, "Lengths 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 state : scc) {
+		if (relativeError) {
+			if (std::abs(vectorLeft[state] - vectorRight[state])/vectorRight[state] > precision) return false;
+		} else {
+			if (std::abs(vectorLeft[state] - vectorRight[state]) > precision) return false;
+		}
+	}
+
+	return true;
+}
+
 } //namespace utility
 
 } //namespace storm