diff --git a/src/modelchecker/MdpPrctlModelChecker.h b/src/modelchecker/MdpPrctlModelChecker.h index 8b48f11a3..acc75b2fb 100644 --- a/src/modelchecker/MdpPrctlModelChecker.h +++ b/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) { diff --git a/src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h b/src/modelchecker/TopologicalValueIterationMdpPrctlModelChecker.h new file mode 100644 index 000000000..23b1ac6fc --- /dev/null +++ b/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 + +#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 TopologicalValueIterationMdpPrctlModelChecker : public MdpPrctlModelChecker { + +public: + explicit TopologicalValueIterationMdpPrctlModelChecker(storm::models::Mdp& mdp) : MdpPrctlModelChecker(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 const& matrix, std::vector& x, std::vector const& b, std::vector 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("precision"); + unsigned maxIterations = s->get("maxiter"); + bool relative = s->get("relative"); + + std::vector> stronglyConnectedComponents; + storm::models::GraphTransitions stronglyConnectedComponentsDependencyGraph; + storm::utility::GraphAnalyzer::performSccDecomposition(matrix, nondeterministicChoiceIndices, stronglyConnectedComponents, stronglyConnectedComponentsDependencyGraph); + + std::vector topologicalSort; + storm::utility::GraphAnalyzer::getTopologicalSort(stronglyConnectedComponentsDependencyGraph, topologicalSort); + + // Set up the environment for the power method. + std::vector multiplyResult(matrix.getRowCount()); + std::vector* currentX = &x; + std::vector* newX = new std::vector(x.size()); + std::vector* 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 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_ */ diff --git a/src/models/GraphTransitions.h b/src/models/GraphTransitions.h index 359fcfcf2..ec8bb2bfd 100644 --- a/src/models/GraphTransitions.h +++ b/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 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; } /*! diff --git a/src/parser/NondeterministicSparseTransitionParser.cpp b/src/parser/NondeterministicSparseTransitionParser.cpp index 4d9491ce1..c0bf308b9 100644 --- a/src/parser/NondeterministicSparseTransitionParser.cpp +++ b/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."; diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h index 0080fa2af..570283c4d 100644 --- a/src/storage/BitVector.h +++ b/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); } diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 5cab2deb3..5259a8052 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -830,17 +830,22 @@ public: } } + void multiplyWithVectorInPlace(std::vector& vector) const { + typename std::vector::iterator resultIt = vector.begin(); + typename std::vector::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 const& states, std::vector const& nondeterministicChoiceIndices, std::vector& vector, std::vector& 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 const& states, std::vector const& nondeterministicChoiceIndices, std::vector& 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 const& nondeterministicChoiceIndices, std::vector& vector, std::vector 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. diff --git a/src/utility/GraphAnalyzer.h b/src/utility/GraphAnalyzer.h index 7321926cf..f99c3e319 100644 --- a/src/utility/GraphAnalyzer.h +++ b/src/utility/GraphAnalyzer.h @@ -424,11 +424,72 @@ public: storm::models::GraphTransitions 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 + static void getTopologicalSort(storm::models::GraphTransitions const& transitions, std::vector& topologicalSort) { + topologicalSort.reserve(transitions.getNumberOfStates()); + + std::vector recursionStack; + recursionStack.reserve(transitions.getNumberOfStates()); + + std::vector::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::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 static void performSccDecomposition(uint_fast64_t numberOfStates, storm::models::GraphTransitions const& forwardTransitions, std::vector>& stronglyConnectedComponents, storm::models::GraphTransitions& 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]); diff --git a/src/utility/Vector.h b/src/utility/Vector.h index cdc7cd318..ab1fc461a 100644 --- a/src/utility/Vector.h +++ b/src/utility/Vector.h @@ -107,6 +107,21 @@ void reduceVectorMin(std::vector const& source, std::vector* target, std:: } } +template +void reduceVectorMin(std::vector const& source, std::vector* target, std::vector const& scc, std::vector 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 void reduceVectorMax(std::vector const& source, std::vector* target, std::vector const& filter) { uint_fast64_t currentSourceRow = 0; @@ -127,6 +142,21 @@ void reduceVectorMax(std::vector const& source, std::vector* target, std:: } } +template +void reduceVectorMax(std::vector const& source, std::vector* target, std::vector const& scc, std::vector 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 bool equalModuloPrecision(std::vector const& vectorLeft, std::vector const& vectorRight, T precision, bool relativeError) { if (vectorLeft.size() != vectorRight.size()) { @@ -145,6 +175,24 @@ bool equalModuloPrecision(std::vector const& vectorLeft, std::vector const return true; } +template +bool equalModuloPrecision(std::vector const& vectorLeft, std::vector const& vectorRight, std::vector 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