Fixed several bugs in several modules (bit vector, parser, etc.). Topological value iteration now works for the consensus protocol and the two dice example.

src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h

@@ -67,47 +67,56 @@ private:
 		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 iterations = 0;
-		bool converged = false;
+		uint_fast64_t currentMaxLocalIterations = 0;
+		uint_fast64_t localIterations = 0;
+		uint_fast64_t globalIterations = 0;
+		bool converged = true;
 
-		/*
-		for (auto sccIt = stronglyConnectedComponents.rbegin(), sccIte = stronglyConnectedComponents.rend(); sccIt != sccIte; ++sccIt) {
-			converged = false;
+		for (auto sccIndexIt = topologicalSort.begin(); sccIndexIt != topologicalSort.end() && converged; ++sccIndexIt) {
+			std::vector<uint_fast64_t> const& scc = stronglyConnectedComponents[*sccIndexIt];
 
-			std::cout << "for scc! " << *sccIt << std::endl;
-			while (!converged && iterations < maxIterations) {
-				std::cout << "first line " << std::endl;
+			localIterations = 0;
+			converged = false;
+			while (!converged && localIterations < maxIterations) {
 				// Compute x' = A*x + b.
-				matrix.multiplyWithVector(**sccIt, nondeterministicChoiceIndices, *currentX, multiplyResult);
-				std::cout << "after mult" << std::endl;
-				storm::utility::addVectors(**sccIt, nondeterministicChoiceIndices, multiplyResult, b);
-				std::cout << "after add " << std::endl;
+				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, nondeterministicChoiceIndices);
+					storm::utility::reduceVectorMin(multiplyResult, newX, scc, nondeterministicChoiceIndices);
 				} else {
-					storm::utility::reduceVectorMax(multiplyResult, newX, nondeterministicChoiceIndices);
+					storm::utility::reduceVectorMax(multiplyResult, newX, scc, nondeterministicChoiceIndices);
 				}
 
-				std::cout << "comp" << std::endl;
 				// 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;
-				++iterations;
+				++localIterations;
+				++globalIterations;
+			}
+
+			std::cout << "converged locally for scc of size " << scc.size() << std::endl;
+
+			if (localIterations > currentMaxLocalIterations) {
+				currentMaxLocalIterations = localIterations;
 			}
 		}
 
-		if (iterations % 2 == 1) {
+		if (globalIterations % 2 == 1) {
 			std::swap(x, *currentX);
 			delete currentX;
 		} else {
@@ -116,12 +125,10 @@ private:
 
 		// Check if the solver converged and issue a warning otherwise.
 		if (converged) {
-			LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations.");
+			LOG4CPLUS_INFO(logger, "Iterative solver converged after " << currentMaxLocalIterations << " iterations.");
 		} else {
 			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
 		}
-
-		*/
 	}
 };
 

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

@@ -834,13 +834,7 @@ public:
 		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) {

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