From af1aa4e1e53e59768cc113f5785e157441e37267 Mon Sep 17 00:00:00 2001 From: dehnert Date: Thu, 14 Mar 2013 17:49:14 +0100 Subject: [PATCH] Added native matrix-vector multiplication for our matrix format (as fast as gmm++). Fixed bug in bit vector. Fixed some issues in SCC decomposition. MDP model checkers now have the solving methods by default (native ones) and may override them with their own ones, if desired. Added some aux stuff, like vector helper methods. --- src/modelchecker/MdpPrctlModelChecker.h | 79 ++++++++++++++++++++++- src/models/GraphTransitions.h | 42 ++++++------ src/storage/BitVector.h | 4 +- src/storage/SparseMatrix.h | 86 ++++++++++++++++++++++--- src/utility/GraphAnalyzer.h | 23 ++++--- src/utility/Vector.h | 9 +++ 6 files changed, 199 insertions(+), 44 deletions(-) diff --git a/src/modelchecker/MdpPrctlModelChecker.h b/src/modelchecker/MdpPrctlModelChecker.h index 0f50dfae9..8b48f11a3 100644 --- a/src/modelchecker/MdpPrctlModelChecker.h +++ b/src/modelchecker/MdpPrctlModelChecker.h @@ -500,10 +500,85 @@ protected: mutable std::stack minimumOperatorStack; private: - virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& matrix, std::vector& vector, std::vector* summand = nullptr, uint_fast64_t repetitions = 1) const = 0; + virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix const& matrix, std::vector& vector, std::vector* summand = nullptr, uint_fast64_t repetitions = 1) const { + // Get the starting row indices for the non-deterministic choices to reduce the resulting + // vector properly. + std::shared_ptr> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices(); + + // Create vector for result of multiplication, which is reduced to the result vector after + // each multiplication. + std::vector multiplyResult(matrix.getRowCount()); + + // Now perform matrix-vector multiplication as long as we meet the bound of the formula. + for (uint_fast64_t i = 0; i < repetitions; ++i) { + matrix.multiplyWithVector(vector, multiplyResult); + if (summand != nullptr) { + gmm::add(*summand, multiplyResult); + } + if (this->minimumOperatorStack.top()) { + storm::utility::reduceVectorMin(multiplyResult, &vector, *nondeterministicChoiceIndices); + } else { + storm::utility::reduceVectorMax(multiplyResult, &vector, *nondeterministicChoiceIndices); + } + } + } - virtual void solveEquationSystem(storm::storage::SparseMatrix const& matrix, std::vector& vector, std::vector const& b, std::vector const& nondeterministicChoiceIndices) const = 0; + virtual 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"); + + // 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 iterations = 0; + bool converged = false; + + // Proceed with the iterations as long as the method did not converge or reach the + // user-specified maximum number of iterations. + while (!converged && iterations < maxIterations) { + // Compute x' = A*x + b. + matrix.multiplyWithVector(*currentX, multiplyResult); + + gmm::add(b, multiplyResult); + + // Reduce the vector x' by applying min/max for all non-deterministic choices. + if (this->minimumOperatorStack.top()) { + storm::utility::reduceVectorMin(multiplyResult, newX, nondeterministicChoiceIndices); + } else { + storm::utility::reduceVectorMax(multiplyResult, newX, nondeterministicChoiceIndices); + } + // Determine whether the method converged. + converged = storm::utility::equalModuloPrecision(*currentX, *newX, precision, relative); + + // Update environment variables. + swap = currentX; + currentX = newX; + newX = swap; + ++iterations; + } + + if (iterations % 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 " << iterations << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge."); + } + } std::shared_ptr> computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector constraint) const { std::shared_ptr> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices(); diff --git a/src/models/GraphTransitions.h b/src/models/GraphTransitions.h index 7629a58f8..5e54d577d 100644 --- a/src/models/GraphTransitions.h +++ b/src/models/GraphTransitions.h @@ -39,8 +39,8 @@ public: * @param forward If set to true, this objects will store the graph structure * of the backwards transition relation. */ - GraphTransitions(std::shared_ptr> transitionMatrix, bool forward) - : successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix->getColumnCount()), numberOfTransitions(transitionMatrix->getNonZeroEntryCount()) { + GraphTransitions(storm::storage::SparseMatrix const& transitionMatrix, bool forward) + : successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix.getColumnCount()), numberOfTransitions(transitionMatrix.getNonZeroEntryCount()) { if (forward) { this->initializeForward(transitionMatrix); } else { @@ -48,12 +48,12 @@ public: } } - GraphTransitions(std::shared_ptr> transitionMatrix, std::shared_ptr> choiceIndices, bool forward) - : successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix->getColumnCount()), numberOfTransitions(transitionMatrix->getNonZeroEntryCount()) { + GraphTransitions(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& nondeterministicChoiceIndices, bool forward) + : successorList(nullptr), stateIndications(nullptr), numberOfStates(transitionMatrix.getColumnCount()), numberOfTransitions(transitionMatrix.getNonZeroEntryCount()) { if (forward) { - this->initializeForward(transitionMatrix, choiceIndices); + this->initializeForward(transitionMatrix, nondeterministicChoiceIndices); } else { - this->initializeBackward(transitionMatrix, choiceIndices); + this->initializeBackward(transitionMatrix, nondeterministicChoiceIndices); } } @@ -126,37 +126,37 @@ private: * Initializes this graph transitions object using the forward transition * relation given by means of a sparse matrix. */ - void initializeForward(std::shared_ptr> transitionMatrix) { + void initializeForward(storm::storage::SparseMatrix const& transitionMatrix) { this->successorList = new uint_fast64_t[numberOfTransitions]; this->stateIndications = new uint_fast64_t[numberOfStates + 1]; // First, we copy the index list from the sparse matrix as this will // stay the same. - std::copy(transitionMatrix->getRowIndications().begin(), transitionMatrix->getRowIndications().end(), this->stateIndications); + std::copy(transitionMatrix.getRowIndications().begin(), transitionMatrix.getRowIndications().end(), this->stateIndications); // Now we can iterate over all rows of the transition matrix and record // the target state. for (uint_fast64_t i = 0, currentNonZeroElement = 0; i < numberOfStates; i++) { - for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) { + for (auto rowIt = transitionMatrix.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) { this->successorList[currentNonZeroElement++] = *rowIt; } } } - void initializeForward(std::shared_ptr> transitionMatrix, std::shared_ptr> choiceIndices) { + void initializeForward(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& nondeterministicChoiceIndices) { this->successorList = new uint_fast64_t[numberOfTransitions]; this->stateIndications = new uint_fast64_t[numberOfStates + 1]; for (uint_fast64_t i = 0; i < numberOfStates; ++i) { - this->stateIndications[i] = transitionMatrix->getRowIndications().at(choiceIndices->at(i)); + this->stateIndications[i] = transitionMatrix.getRowIndications().at(nondeterministicChoiceIndices[i]); } this->stateIndications[numberOfStates] = numberOfTransitions; // Now we can iterate over all rows of the transition matrix and record // the target state. for (uint_fast64_t i = 0, currentNonZeroElement = 0; i < numberOfStates; i++) { - for (uint_fast64_t j = choiceIndices->at(i); j < choiceIndices->at(i + 1); ++j) { - for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) { + for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) { + for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) { this->successorList[currentNonZeroElement++] = *rowIt; } } @@ -168,13 +168,13 @@ private: * relation, whose forward transition relation is given by means of a sparse * matrix. */ - void initializeBackward(std::shared_ptr> transitionMatrix) { + void initializeBackward(storm::storage::SparseMatrix const& transitionMatrix) { this->successorList = new uint_fast64_t[numberOfTransitions]; this->stateIndications = new uint_fast64_t[numberOfStates + 1](); // First, we need to count how many backward transitions each state has. for (uint_fast64_t i = 0; i < numberOfStates; ++i) { - for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) { + for (auto rowIt = transitionMatrix.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) { this->stateIndications[*rowIt + 1]++; } } @@ -197,7 +197,7 @@ private: // Now we are ready to actually fill in the list of predecessors for // every state. Again, we start by considering all but the last row. for (uint_fast64_t i = 0; i < numberOfStates; ++i) { - for (auto rowIt = transitionMatrix->beginConstColumnIterator(i); rowIt != transitionMatrix->endConstColumnIterator(i); ++rowIt) { + for (auto rowIt = transitionMatrix.beginConstColumnIterator(i); rowIt != transitionMatrix.endConstColumnIterator(i); ++rowIt) { this->successorList[nextIndicesList[*rowIt]++] = i; } } @@ -206,14 +206,14 @@ private: delete[] nextIndicesList; } - void initializeBackward(std::shared_ptr> transitionMatrix, std::shared_ptr> choiceIndices) { + void initializeBackward(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& nondeterministicChoiceIndices) { this->successorList = new uint_fast64_t[numberOfTransitions]; this->stateIndications = new uint_fast64_t[numberOfStates + 1](); // First, we need to count how many backward transitions each state has. for (uint_fast64_t i = 0; i < numberOfStates; ++i) { - for (uint_fast64_t j = choiceIndices->at(i); j < choiceIndices->at(i + 1); ++j) { - for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) { + for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) { + for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) { this->stateIndications[*rowIt + 1]++; } } @@ -237,8 +237,8 @@ private: // Now we are ready to actually fill in the list of predecessors for // every state. Again, we start by considering all but the last row. for (uint_fast64_t i = 0; i < numberOfStates; i++) { - for (uint_fast64_t j = (*choiceIndices)[i]; j < (*choiceIndices)[i + 1]; ++j) { - for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) { + for (uint_fast64_t j = nondeterministicChoiceIndices[i]; j < nondeterministicChoiceIndices[i + 1]; ++j) { + for (auto rowIt = transitionMatrix.beginConstColumnIterator(j); rowIt != transitionMatrix.endConstColumnIterator(j); ++rowIt) { this->successorList[nextIndicesList[*rowIt]++] = i; } } diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h index 65ace6e3c..0080fa2af 100644 --- a/src/storage/BitVector.h +++ b/src/storage/BitVector.h @@ -547,7 +547,7 @@ private: startingIndex >>= 6; uint64_t* bucketPtr = this->bucketArray + startingIndex; - do { + while ((startingIndex << 6) < endIndex) { // Compute the remaining bucket content by a right shift // to the current bit. uint_fast64_t remainingInBucket = *bucketPtr >> currentBitInByte; @@ -570,7 +570,7 @@ private: // Advance to the next bucket. ++startingIndex; ++bucketPtr; currentBitInByte = 0; - } while ((startingIndex << 6) < endIndex); + } return endIndex; } diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index f5c4c954f..a96c21923 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -62,6 +62,44 @@ public: */ typedef const T* const constIterator; + class constRowsIterator { + public: + constRowsIterator(SparseMatrix const& matrix) : matrix(matrix), offset(0) { + // Intentionally left empty. + } + + constRowsIterator& operator++() { + ++offset; + return *this; + } + + bool operator==(constRowsIterator const& other) const { + return offset == other.offset; + } + + bool operator!=(uint_fast64_t offset) const { + return this->offset != offset; + } + + uint_fast64_t column() const { + return matrix.columnIndications[offset]; + } + + T const& value() const { + return matrix.valueStorage[offset]; + } + + void setOffset(uint_fast64_t offset) { + this->offset = offset; + } + + private: + SparseMatrix const& matrix; + uint_fast64_t offset; + }; + + friend class constRowsIterator; + /*! * An enum representing the internal state of the Matrix. * After creating the Matrix using the Constructor, the Object is in state UnInitialized. After calling initialize(), that state changes to Initialized and after all entries have been entered and finalize() has been called, to ReadReady. @@ -773,18 +811,44 @@ public: return result; } - T getRowVectorProduct(uint_fast64_t row, std::vector& vector) { - T result = storm::utility::constGetZero();; - auto valueIterator = valueStorage.begin() + rowIndications[row]; - const auto valueIteratorEnd = valueStorage.begin() + rowIndications[row + 1]; - auto columnIterator = columnIndications.begin() + rowIndications[row]; - const auto columnIteratorEnd = columnIndications.begin() + rowIndications[row + 1]; - for (; valueIterator != valueIteratorEnd; ++valueIterator, ++columnIterator) { - result += *valueIterator * vector[*columnIterator]; + T multiplyRowWithVector(constRowsIterator& rowsIt, uint_fast64_t rowsIte, std::vector& vector) const { + T result = storm::utility::constGetZero(); + for (; rowsIt != rowsIte; ++rowsIt) { + result += (rowsIt.value()) * vector[rowsIt.column()]; } return result; } + void multiplyWithVector(std::vector& vector, std::vector& result) const { + typename std::vector::iterator resultIt = result.begin(); + typename std::vector::iterator resultIte = result.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) { + result[rowIt] = multiplyRowWithVector(rowsIt, this->rowIndications[nextRow], vector); + } + } + } + /*! * Returns the size of the matrix in memory measured in bytes. * @return The size of the matrix in memory measured in bytes. @@ -800,6 +864,10 @@ public: return size; } + constRowsIterator constRowsIteratorBegin() const { + return constRowsIterator(*this); + } + /*! * Returns an iterator to the columns of the non-zero entries of the given * row. @@ -910,7 +978,7 @@ public: * to separate groups will be separated by a dashed line. * @return a (non-compressed) string representation of the matrix. */ - std::string toString(std::shared_ptr> nondeterministicChoiceIndices) const { + std::string toString(std::vector const* nondeterministicChoiceIndices) const { std::stringstream result; uint_fast64_t currentNondeterministicChoiceIndex = 0; diff --git a/src/utility/GraphAnalyzer.h b/src/utility/GraphAnalyzer.h index eeb619a6c..06f6742db 100644 --- a/src/utility/GraphAnalyzer.h +++ b/src/utility/GraphAnalyzer.h @@ -71,7 +71,7 @@ public: } // Get the backwards transition relation from the model to ease the search. - storm::models::GraphTransitions backwardTransitions(model.getTransitionMatrix(), false); + storm::models::GraphTransitions backwardTransitions(*model.getTransitionMatrix(), false); // Add all psi states as the already satisfy the condition. *statesWithProbabilityGreater0 |= psiStates; @@ -174,7 +174,7 @@ public: } // Get the backwards transition relation from the model to ease the search. - storm::models::GraphTransitions backwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), false); + storm::models::GraphTransitions backwardTransitions(*model.getTransitionMatrix(), *model.getNondeterministicChoiceIndices(), false); // Add all psi states as the already satisfy the condition. *statesWithProbability0 |= psiStates; @@ -213,7 +213,7 @@ public: std::shared_ptr> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices(); // Get the backwards transition relation from the model to ease the search. - storm::models::GraphTransitions backwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), false); + storm::models::GraphTransitions backwardTransitions(*model.getTransitionMatrix(), *model.getNondeterministicChoiceIndices(), false); storm::storage::BitVector* currentStates = new storm::storage::BitVector(model.getNumberOfStates(), true); @@ -300,7 +300,7 @@ public: std::shared_ptr> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices(); // Get the backwards transition relation from the model to ease the search. - storm::models::GraphTransitions backwardTransitions(transitionMatrix, nondeterministicChoiceIndices, false); + storm::models::GraphTransitions backwardTransitions(*transitionMatrix, *nondeterministicChoiceIndices, false); // Add all psi states as the already satisfy the condition. *statesWithProbability0 |= psiStates; @@ -360,7 +360,7 @@ public: std::shared_ptr> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices(); // Get the backwards transition relation from the model to ease the search. - storm::models::GraphTransitions backwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), false); + storm::models::GraphTransitions backwardTransitions(*model.getTransitionMatrix(), *model.getNondeterministicChoiceIndices(), false); storm::storage::BitVector* currentStates = new storm::storage::BitVector(model.getNumberOfStates(), true); @@ -417,14 +417,17 @@ public: } template - static uint_fast64_t performSccDecomposition(storm::models::AbstractNondeterministicModel& model, std::vector*>* stronglyConnectedComponents) { + static uint_fast64_t performSccDecomposition(storm::storage::SparseMatrix const& matrix, std::vector const& nondeterministicChoiceIndices, std::vector*>* stronglyConnectedComponents) { LOG4CPLUS_INFO(logger, "Computing SCC decomposition."); // Get the forward transition relation from the model to ease the search. - storm::models::GraphTransitions forwardTransitions(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), true); + storm::models::GraphTransitions forwardTransitions(matrix, nondeterministicChoiceIndices, true); + + std::cout << matrix.toString(&nondeterministicChoiceIndices) << std::endl; + std::cout << forwardTransitions.toString() << std::endl; // Perform the actual SCC decomposition based on the graph-transitions of the system. - uint_fast64_t result = performSccDecomposition(model.getNumberOfStates(), *model.getLabeledStates("init"), forwardTransitions, stronglyConnectedComponents); + uint_fast64_t result = performSccDecomposition(nondeterministicChoiceIndices.size(), forwardTransitions, stronglyConnectedComponents); LOG4CPLUS_INFO(logger, "Done computing SCC decomposition."); return result; @@ -432,7 +435,7 @@ public: private: template - static uint_fast64_t performSccDecomposition(uint_fast64_t numberOfStates, storm::storage::BitVector const& initialStates, storm::models::GraphTransitions const& forwardTransitions, std::vector*>* stronglyConnectedComponents) { + static uint_fast64_t performSccDecomposition(uint_fast64_t numberOfStates, storm::models::GraphTransitions const& forwardTransitions, std::vector*>* stronglyConnectedComponents) { std::vector tarjanStack; tarjanStack.reserve(numberOfStates); storm::storage::BitVector tarjanStackStates(numberOfStates); @@ -442,7 +445,7 @@ private: storm::storage::BitVector visitedStates(numberOfStates); uint_fast64_t currentIndex = 0; - for (auto state : initialStates) { + for (uint_fast64_t state = 0; state < numberOfStates; ++state) { if (!visitedStates.get(state)) { performSccDecompositionHelper(state, currentIndex, stateIndices, lowlinks, tarjanStack, tarjanStackStates, visitedStates, forwardTransitions, stronglyConnectedComponents); } diff --git a/src/utility/Vector.h b/src/utility/Vector.h index fce8c18e1..cdc7cd318 100644 --- a/src/utility/Vector.h +++ b/src/utility/Vector.h @@ -78,6 +78,15 @@ void subtractFromConstantOneVector(std::vector* vector) { } } +template +void addVectors(std::vector const& states, std::vector const& nondeterministicChoiceIndices, std::vector& original, std::vector const& summand) { + for (auto stateIt = states.cbegin(), stateIte = states.cend(); stateIt != stateIte; ++stateIt) { + for (auto rowIt = nondeterministicChoiceIndices[*stateIt], rowIte = nondeterministicChoiceIndices[*stateIt + 1]; rowIt != rowIte; ++rowIt) { + original[rowIt] += summand[rowIt]; + } + } +} + template void reduceVectorMin(std::vector const& source, std::vector* target, std::vector const& filter) { uint_fast64_t currentSourceRow = 0;