From d4cf812c5ec211b4b73471e53fe8cc8704444553 Mon Sep 17 00:00:00 2001 From: dehnert Date: Sun, 17 Feb 2013 13:43:52 +0100 Subject: [PATCH] Added until-model checking for MDPs. Implemented Prob1A algorithm. Added asynchronous leader example. --- src/formula/ProbabilisticNoBoundOperator.h | 9 + src/formula/RewardNoBoundOperator.h | 19 +- src/modelchecker/GmmxxDtmcPrctlModelChecker.h | 9 +- src/modelchecker/GmmxxMdpPrctlModelChecker.h | 146 ++++++++------- src/modelchecker/MdpPrctlModelChecker.h | 7 +- src/models/GraphTransitions.h | 2 +- ...NondeterministicSparseTransitionParser.cpp | 12 +- src/storage/SparseMatrix.h | 167 +++++++++++++++++- src/storm.cpp | 43 ++++- src/utility/GraphAnalyzer.h | 77 +++++++- src/utility/Vector.h | 49 +++-- 11 files changed, 430 insertions(+), 110 deletions(-) diff --git a/src/formula/ProbabilisticNoBoundOperator.h b/src/formula/ProbabilisticNoBoundOperator.h index cae4df96e..ee7f50f48 100644 --- a/src/formula/ProbabilisticNoBoundOperator.h +++ b/src/formula/ProbabilisticNoBoundOperator.h @@ -64,6 +64,15 @@ public: // Intentionally left empty } + /*! + * Constructor + * + * @param pathFormula The child node. + */ + ProbabilisticNoBoundOperator(AbstractPathFormula* pathFormula, bool minimumOperator) : NoBoundOperator(pathFormula, minimumOperator) { + // Intentionally left empty + } + /*! * @returns a string representation of the formula */ diff --git a/src/formula/RewardNoBoundOperator.h b/src/formula/RewardNoBoundOperator.h index 6765bc854..63aef3413 100644 --- a/src/formula/RewardNoBoundOperator.h +++ b/src/formula/RewardNoBoundOperator.h @@ -64,11 +64,28 @@ public: // Intentionally left empty } + /*! + * Constructor + * + * @param pathFormula The child node. + */ + RewardNoBoundOperator(AbstractPathFormula* pathFormula, bool minimumOperator) : NoBoundOperator(pathFormula, minimumOperator) { + // Intentionally left empty + } + /*! * @returns a string representation of the formula */ virtual std::string toString() const { - std::string result = "R = ? ["; + std::string result = "R"; + if (this->isOptimalityOperator()) { + if (this->isMinimumOperator()) { + result += "min"; + } else { + result += "max"; + } + } + result += "=? ["; result += this->getPathFormula().toString(); result += "]"; return result; diff --git a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h index 3d3971b35..308a6e80e 100644 --- a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h +++ b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h @@ -148,7 +148,7 @@ public: // 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 b(mayBeStatesSetBitCount); - this->getModel().getTransitionMatrix()->getConstrainedRowCountVector(maybeStates, statesWithProbability1, &b); + this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, statesWithProbability1, &b); // Solve the corresponding system of linear equations. this->solveLinearEquationSystem(*gmmxxMatrix, x, b); @@ -346,9 +346,10 @@ public: */ static void putOptions(boost::program_options::options_description* desc) { desc->add_options()("lemethod", boost::program_options::value()->default_value("bicgstab")->notifier(&validateLeMethod), "Sets the method used for linear equation solving. Must be in {bicgstab, qmr}."); - desc->add_options()("lemaxiter", boost::program_options::value()->default_value(10000), "Sets the maximal number of iterations for iterative linear equation solving."); - desc->add_options()("precision", boost::program_options::value()->default_value(10e-6), "Sets the precision for iterative linear equation solving."); + desc->add_options()("maxiter", boost::program_options::value()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving."); + desc->add_options()("precision", boost::program_options::value()->default_value(1e-6), "Sets the precision for iterative equation solving."); desc->add_options()("precond", boost::program_options::value()->default_value("ilu")->notifier(&validatePreconditioner), "Sets the preconditioning technique for linear equation solving. Must be in {ilu, diagonal, ildlt, none}."); + desc->add_options()("relative", boost::program_options::value()->default_value(true), "Sets whether the relative or absolute error is considered for detecting convergence."); } /*! @@ -388,7 +389,7 @@ private: // Prepare an iteration object that determines the accuracy, maximum number of iterations // and the like. - gmm::iteration iter(s->get("precision"), 0, s->get("lemaxiter")); + gmm::iteration iter(s->get("precision"), 0, s->get("maxiter")); // Now do the actual solving. LOG4CPLUS_INFO(logger, "Starting iterative solver."); diff --git a/src/modelchecker/GmmxxMdpPrctlModelChecker.h b/src/modelchecker/GmmxxMdpPrctlModelChecker.h index d0d940eba..41efbd145 100644 --- a/src/modelchecker/GmmxxMdpPrctlModelChecker.h +++ b/src/modelchecker/GmmxxMdpPrctlModelChecker.h @@ -5,8 +5,8 @@ * Author: Christian Dehnert */ -#ifndef STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_ -#define STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_ +#ifndef STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_ +#define STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_ #include @@ -36,12 +36,12 @@ namespace modelChecker { * A model checking engine that makes use of the gmm++ backend. */ template -class GmmxxDtmcPrctlModelChecker : public MdpPrctlModelChecker { +class GmmxxMdpPrctlModelChecker : public MdpPrctlModelChecker { public: - explicit GmmxxDtmcPrctlModelChecker(storm::models::Mdp& mdp) : MdpPrctlModelChecker(mdp) { } + explicit GmmxxMdpPrctlModelChecker(storm::models::Mdp& mdp) : MdpPrctlModelChecker(mdp) { } - virtual ~GmmxxDtmcPrctlModelChecker() { } + virtual ~GmmxxMdpPrctlModelChecker() { } virtual std::vector* checkBoundedUntil(const storm::formula::BoundedUntil& formula) const { // First, we need to compute the states that satisfy the sub-formulas of the until-formula. @@ -67,13 +67,13 @@ public: // Create vector for result of multiplication, which is reduced to the result vector after // each multiplication. - std::vector* multiplyResult = new std::vector(this->getModel().getTransitionMatrix().getRowCount()); + std::vector* multiplyResult = new std::vector(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()) { + if (this->minimumOperatorStack.top()) { storm::utility::reduceVectorMin(*multiplyResult, result, *nondeterministicChoiceIndices); } else { storm::utility::reduceVectorMax(*multiplyResult, result, *nondeterministicChoiceIndices); @@ -103,7 +103,7 @@ public: delete nextStates; // Create resulting vector. - std::vector* temporaryResult = new std::vector(this->getModel().getTransitionMatrix().getRowCount()); + std::vector* temporaryResult = new std::vector(this->getModel().getTransitionMatrix()->getRowCount()); // Perform the actual computation, namely matrix-vector multiplication. gmm::mult(*gmmxxMatrix, *result, *temporaryResult); @@ -112,7 +112,7 @@ public: // vector properly. std::shared_ptr> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices(); - if (minimumOperatorStack.top()) { + if (this->minimumOperatorStack.top()) { storm::utility::reduceVectorMin(*temporaryResult, result, *nondeterministicChoiceIndices); } else { storm::utility::reduceVectorMax(*temporaryResult, result, *nondeterministicChoiceIndices); @@ -133,7 +133,7 @@ public: // 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()) { + if (this->minimumOperatorStack.top()) { storm::utility::GraphAnalyzer::performProb01Min(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1); } else { storm::utility::GraphAnalyzer::performProb01Max(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1); @@ -152,32 +152,36 @@ public: std::vector* result = new std::vector(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* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates); + uint_fast64_t maybeStatesSetBitCount = maybeStates.getNumberOfSetBits(); + if (maybeStatesSetBitCount > 0) { + // First, we can eliminate the rows and columns from the original transition probability matrix for states + // whose probabilities are already known. + storm::storage::SparseMatrix* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices()); - // Transform the submatrix to the gmm++ format to use its solvers. + // Get the "new" nondeterministic choice indices for the submatrix. + std::shared_ptr> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates); + + // Transform the submatrix to the gmm++ format to use its capabilities. gmm::csr_matrix* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*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 x(mayBeStatesSetBitCount, Type(0.5)); + // Create vector for results for maybe states. + std::vector* x = new std::vector(maybeStatesSetBitCount); // 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 b(mayBeStatesSetBitCount); - this->getModel().getTransitionMatrix()->getConstrainedRowCountVector(maybeStates, statesWithProbability1, &b); + std::vector b(submatrix->getRowCount()); + this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, *this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, &b); + delete submatrix; - // Solve the corresponding system of linear equations. - this->solveEquationSystem(*gmmxxMatrix, x, b); + // Solve the corresponding system of equations. + this->solveEquationSystem(*gmmxxMatrix, x, b, *subNondeterministicChoiceIndices); // Set values of resulting vector according to result. - storm::utility::setVectorValues(result, maybeStates, x); + storm::utility::setVectorValues(result, maybeStates, *x); - // Delete temporary matrix. + // Delete temporary matrix and vector. delete gmmxxMatrix; + delete x; } // Set values of resulting vector that are known exactly. @@ -272,7 +276,8 @@ public: // 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); + // TODO: just commented out to make it compile + //storm::utility::GraphAnalyzer::performProb1(this->getModel(), trueStates, *targetStates, &infinityStates); infinityStates.complement(); // Create resulting vector. @@ -325,7 +330,8 @@ public: } // Solve the corresponding system of linear equations. - this->solveLinearEquationSystem(*gmmxxMatrix, x, *b); + // TODO: just commented out to make it compile + // this->solveEquationSystem(*gmmxxMatrix, x, *b); // Set values of resulting vector according to result. storm::utility::setVectorValues(result, maybeStates, x); @@ -344,31 +350,6 @@ public: 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 getOptionTrigger() { - return std::pair("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()->default_value(10000), "Sets the maximal number of iterations for iterative equation solving."); - desc->add_options()("precision", boost::program_options::value()->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. @@ -380,25 +361,70 @@ private: * @return The solution of the system of linear equations in form of the elements of the vector * x. */ - void solveLinearEquationSystem(gmm::csr_matrix const& A, std::vector& x, std::vector const& b) const { - // Get the settings object to customize linear solving. + void solveEquationSystem(gmm::csr_matrix const& A, 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 precison = s->get("precision"); - unsigned maxIterations = s->get("lemaxiter"); + // 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* temporaryResult = new std::vector(b.size()); + std::vector* newX = new std::vector(x->size()); + std::vector* swap = nullptr; + bool converged = false; + uint_fast64_t iterations = 0; + + // 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. + gmm::mult(A, *x, *temporaryResult); + gmm::add(b, *temporaryResult); + + // Reduce the vector x' by applying min/max for all non-deterministic choices. + if (this->minimumOperatorStack.top()) { + storm::utility::reduceVectorMin(*temporaryResult, newX, nondeterministicChoiceIndices); + } else { + storm::utility::reduceVectorMax(*temporaryResult, newX, nondeterministicChoiceIndices); + } + // Determine whether the method converged. + converged = storm::utility::equalModuloPrecision(*x, *newX, precision, relative); + // Update environment variables. + swap = x; + x = newX; + newX = swap; + ++iterations; + } + delete temporaryResult; // Check if the solver converged and issue a warning otherwise. - if (iter.converged()) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iter.get_iteration() << " iterations."); + 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(); + std::shared_ptr> subNondeterministicChoiceIndices(new std::vector(constraint.getNumberOfSetBits() + 1)); + uint_fast64_t currentRowCount = 0; + uint_fast64_t currentIndexCount = 1; + (*subNondeterministicChoiceIndices)[0] = 0; + for (auto index : constraint) { + (*subNondeterministicChoiceIndices)[currentIndexCount] = currentRowCount + (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index]; + currentRowCount += (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index]; + ++currentIndexCount; + } + (*subNondeterministicChoiceIndices)[constraint.getNumberOfSetBits()] = currentRowCount; + return subNondeterministicChoiceIndices; } }; @@ -406,4 +432,4 @@ private: } //namespace storm -#endif /* STORM_MODELCHECKER_GMMXXDTMCPRCTLMODELCHECKER_H_ */ +#endif /* STORM_MODELCHECKER_GMMXXMDPPRCTLMODELCHECKER_H_ */ diff --git a/src/modelchecker/MdpPrctlModelChecker.h b/src/modelchecker/MdpPrctlModelChecker.h index e3139062f..7dfb88825 100644 --- a/src/modelchecker/MdpPrctlModelChecker.h +++ b/src/modelchecker/MdpPrctlModelChecker.h @@ -18,6 +18,7 @@ #include "src/utility/Vector.h" #include "src/modelchecker/AbstractModelChecker.h" #include +#include #include "log4cplus/logger.h" #include "log4cplus/loggingmacros.h" @@ -361,11 +362,11 @@ public: */ virtual std::vector* checkReachabilityReward(const storm::formula::ReachabilityReward& formula) const = 0; -protected: - std::stack minimumOperatorStack; - private: storm::models::Mdp& model; + +protected: + mutable std::stack minimumOperatorStack; }; } //namespace modelChecker diff --git a/src/models/GraphTransitions.h b/src/models/GraphTransitions.h index 6f4da36f1..3af0524b4 100644 --- a/src/models/GraphTransitions.h +++ b/src/models/GraphTransitions.h @@ -216,7 +216,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 (uint_fast64_t j = choiceIndices->at(i); j < choiceIndices->at(i + 1); ++j) { + for (uint_fast64_t j = (*choiceIndices)[i]; j < (*choiceIndices)[i + 1]; ++j) { for (auto rowIt = transitionMatrix->beginConstColumnIterator(j); rowIt != transitionMatrix->endConstColumnIterator(j); ++rowIt) { this->successorList[nextIndicesList[*rowIt]++] = i; } diff --git a/src/parser/NondeterministicSparseTransitionParser.cpp b/src/parser/NondeterministicSparseTransitionParser.cpp index c5118eee4..1cd08604b 100644 --- a/src/parser/NondeterministicSparseTransitionParser.cpp +++ b/src/parser/NondeterministicSparseTransitionParser.cpp @@ -131,9 +131,17 @@ uint_fast64_t NondeterministicSparseTransitionParser::firstPass(char* buf, uint_ */ nonzero++; + // The PRISM output format lists the name of the transition in the fourth column, + // but omits the fourth column if it is an internal action. In either case, however, the third column + // is followed by a space. We need to skip over that space first (instead of trimming whitespaces), + // before we can skip to the line end, because trimming the white spaces will proceed to the next line + // in case there is no action label in the fourth column. + ++buf; + /* * Proceed to beginning of next line. */ + buf += strcspn(buf, " \t\n\r"); buf = trimWhitespaces(buf); } @@ -277,10 +285,12 @@ NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(s /* * Proceed to beginning of next line in file and next row in matrix. */ + ++buf; + buf += strcspn(buf, " \t\n\r"); buf = trimWhitespaces(buf); } - this->rowMapping->at(maxnode+1) = curRow; + this->rowMapping->at(maxnode+1) = curRow + 1; if (!fixDeadlocks && hadDeadlocks) 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/SparseMatrix.h b/src/storage/SparseMatrix.h index a35df8b47..32ddfb7a4 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -668,7 +668,6 @@ public: uint_fast64_t rowEnd = rowIndications[row + 1]; if (rowStart >= rowEnd) { - this->print(); LOG4CPLUS_ERROR(logger, "Cannot make row absorbing, because there is no entry in this row."); throw storm::exceptions::InvalidStateException("Cannot make row absorbing, because there is no entry in this row."); } @@ -722,13 +721,30 @@ public: * @param resultVector A pointer to the resulting vector that has at least * as many elements as there are bits set to true in the constraint. */ - void getConstrainedRowCountVector(const storm::storage::BitVector& rowConstraint, const storm::storage::BitVector& columnConstraint, std::vector* resultVector) const { + void getConstrainedRowSumVector(const storm::storage::BitVector& rowConstraint, const storm::storage::BitVector& columnConstraint, std::vector* resultVector) const { uint_fast64_t currentRowCount = 0; for (auto row : rowConstraint) { (*resultVector)[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); } } + /*! + * Computes a vector in which each element is the sum of those elements in the + * corresponding row whose column bits are set to one in the given constraint. + * @param rowConstraint A bit vector that indicates for which rows to perform summation. + * @param columnConstraint A bit vector that indicates which columns to add. + * @param resultVector A pointer to the resulting vector that has at least + * as many elements as there are bits set to true in the constraint. + */ + void getConstrainedRowSumVector(const storm::storage::BitVector& rowConstraint, std::vector const& nondeterministicChoiceIndices, const storm::storage::BitVector& columnConstraint, std::vector* resultVector) const { + uint_fast64_t currentRowCount = 0; + for (auto index : rowConstraint) { + for (uint_fast64_t row = nondeterministicChoiceIndices[index]; row < nondeterministicChoiceIndices[index + 1]; ++row) { + (*resultVector)[currentRowCount++] = getConstrainedRowSum(row, columnConstraint); + } + } + } + /*! * Creates a sub-matrix of the current matrix by dropping all rows and * columns whose bits are not set to one in the given bit vector. @@ -761,7 +777,7 @@ public: // Create a temporary array that stores for each index whose bit is set // to true the number of bits that were set before that particular index. - uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[rowCount]; + uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[colCount]; uint_fast64_t lastIndex = 0; uint_fast64_t currentNumberOfSetBits = 0; for (auto index : constraint) { @@ -792,6 +808,70 @@ public: return result; } + SparseMatrix* getSubmatrix(storm::storage::BitVector& constraint, std::vector const& rowIndices) const { + LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size)."); + + // Check for valid constraint. + if (constraint.getNumberOfSetBits() == 0) { + LOG4CPLUS_ERROR(logger, "Trying to create a sub-matrix of size 0."); + throw storm::exceptions::InvalidArgumentException("Trying to create a sub-matrix of size 0."); + } + + // First, we need to determine the number of non-zero entries and the number of rows of the + // sub-matrix. + uint_fast64_t subNonZeroEntries = 0; + uint_fast64_t subRowCount = 0; + for (auto index : constraint) { + subRowCount += rowIndices[index + 1] - rowIndices[index]; + for (uint_fast64_t i = rowIndices[index]; i < rowIndices[index + 1]; ++i) { + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { + if (constraint.get(columnIndications[j])) { + ++subNonZeroEntries; + } + } + } + } + + LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << constraint.getNumberOfSetBits() << "."); + + // Create and initialize resulting matrix. + SparseMatrix* result = new SparseMatrix(subRowCount, constraint.getNumberOfSetBits()); + result->initialize(subNonZeroEntries); + + // Create a temporary array that stores for each index whose bit is set + // to true the number of bits that were set before that particular index. + uint_fast64_t* bitsSetBeforeIndex = new uint_fast64_t[colCount]; + uint_fast64_t lastIndex = 0; + uint_fast64_t currentNumberOfSetBits = 0; + for (auto index : constraint) { + while (lastIndex <= index) { + bitsSetBeforeIndex[lastIndex++] = currentNumberOfSetBits; + } + ++currentNumberOfSetBits; + } + + // Copy over selected entries. + uint_fast64_t rowCount = 0; + for (auto index : constraint) { + for (uint_fast64_t i = rowIndices[index]; i < rowIndices[index + 1]; ++i) { + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { + if (constraint.get(columnIndications[j])) { + result->addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]); + } + } + ++rowCount; + } + } + + // Dispose of the temporary array. + delete[] bitsSetBeforeIndex; + + // Finalize sub-matrix and return result. + result->finalize(); + LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); + return result; + } + void convertToEquationSystem() { invertDiagonal(); negateAllNonDiagonalElements(); @@ -1005,11 +1085,82 @@ public: return true; } - void print() const { - std::cout << "entries in (" << rowCount << "x" << colCount << " matrix):" << std::endl; - std::cout << rowIndications << std::endl; - std::cout << columnIndications << std::endl; - std::cout << valueStorage << std::endl; + /*! + * Retrieves a compressed string representation of the matrix. + * @return a compressed string representation of the matrix. + */ + std::string toStringCompressed() const { + std::stringstream result; + result << rowIndications << std::endl; + result << columnIndications << std::endl; + result << valueStorage << std::endl; + return result.str(); + } + + /*! + * Retrieves a (non-compressed) string representation of the matrix. + * Note: the matrix is presented densely. That is, all zeros are filled in and are part of the string + * representation. + * @param nondeterministicChoiceIndices A vector indicating which rows belong together. If given, rows belonging + * 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::stringstream result; + uint_fast64_t currentNondeterministicChoiceIndex = 0; + + // Print column numbers in header. + result << "\t\t"; + for (uint_fast64_t i = 0; i < colCount; ++i) { + result << i << "\t"; + } + result << std::endl; + + // Iterate over all rows. + for (uint_fast64_t i = 0; i < rowCount; ++i) { + uint_fast64_t nextIndex = rowIndications[i]; + + // If we need to group of rows, print a dashed line in case we have moved to the next group of rows. + if (nondeterministicChoiceIndices != nullptr) { + if (i == (*nondeterministicChoiceIndices)[currentNondeterministicChoiceIndex]) { + if (i != 0) { + result << "\t(\t"; + for (uint_fast64_t j = 0; j < colCount - 2; ++j) { + result << "----"; + if (j == 1) { + result << "\t" << currentNondeterministicChoiceIndex << "\t"; + } + } + result << "\t)" << std::endl; + } + ++currentNondeterministicChoiceIndex; + } + } + + // Print the actual row. + result << i << "\t(\t"; + uint_fast64_t currentRealIndex = 0; + while (currentRealIndex < colCount) { + if (currentRealIndex == columnIndications[nextIndex] && nextIndex < rowIndications[i + 1]) { + result << valueStorage[nextIndex] << "\t"; + ++nextIndex; + } else { + result << "0\t"; + } + ++currentRealIndex; + } + result << "\t)\t" << i << std::endl; + } + + // Print column numbers in footer. + result << "\t\t"; + for (uint_fast64_t i = 0; i < colCount; ++i) { + result << i << "\t"; + } + result << std::endl; + + // Return final result. + return result.str(); } private: diff --git a/src/storm.cpp b/src/storm.cpp index 326b86e96..62152c31c 100644 --- a/src/storm.cpp +++ b/src/storm.cpp @@ -24,6 +24,7 @@ #include "src/models/AtomicPropositionsLabeling.h" #include "src/modelchecker/EigenDtmcPrctlModelChecker.h" #include "src/modelchecker/GmmxxDtmcPrctlModelChecker.h" +#include "src/modelchecker/GmmxxMdpPrctlModelChecker.h" #include "src/parser/AutoParser.h" #include "src/parser/PrctlParser.h" #include "src/utility/Settings.h" @@ -125,6 +126,10 @@ bool parseOptions(const int argc, const char* argv[]) { return true; } +void setUp() { + std::cout.precision(10); +} + /*! * Function to perform some cleanup. */ @@ -206,18 +211,36 @@ void testCheckingSynchronousLeader(storm::models::Dtmc& dtmc, uint_fast6 } void testCheckingDice(storm::models::Mdp& mdp) { - storm::storage::BitVector* targetStates = mdp.getLabeledStates(std::string("two")); - *targetStates |= *mdp.getLabeledStates(std::string("three")); + storm::formula::Ap* threeFormula = new storm::formula::Ap("three"); + storm::formula::Eventually* eventuallyFormula = new storm::formula::Eventually(threeFormula); + storm::formula::ProbabilisticNoBoundOperator* probFormula = new storm::formula::ProbabilisticNoBoundOperator(eventuallyFormula, false); + + storm::modelChecker::GmmxxMdpPrctlModelChecker* mc = new storm::modelChecker::GmmxxMdpPrctlModelChecker(mdp); + + mc->check(*probFormula); + delete probFormula; + + delete mc; +} + +void testCheckingAsynchLeader(storm::models::Mdp& mdp) { + storm::formula::Ap* electedFormula = new storm::formula::Ap("elected"); + storm::formula::Eventually* eventuallyFormula = new storm::formula::Eventually(electedFormula); + storm::formula::ProbabilisticNoBoundOperator* probMaxFormula = new storm::formula::ProbabilisticNoBoundOperator(eventuallyFormula, false); - storm::storage::BitVector* trueStates = new storm::storage::BitVector(mdp.getNumberOfStates(), true); - storm::storage::BitVector* statesWithProbability0 = new storm::storage::BitVector(mdp.getNumberOfStates()); - storm::storage::BitVector* statesWithProbability1 = new storm::storage::BitVector(mdp.getNumberOfStates()); + storm::modelChecker::GmmxxMdpPrctlModelChecker* mc = new storm::modelChecker::GmmxxMdpPrctlModelChecker(mdp); - storm::utility::GraphAnalyzer::performProb01Max(mdp, *trueStates, *targetStates, statesWithProbability0, statesWithProbability1); + mc->check(*probMaxFormula); + delete probMaxFormula; - delete statesWithProbability0; - delete statesWithProbability1; - delete trueStates; + electedFormula = new storm::formula::Ap("elected"); + eventuallyFormula = new storm::formula::Eventually(electedFormula); + storm::formula::ProbabilisticNoBoundOperator* probMinFormula = new storm::formula::ProbabilisticNoBoundOperator(eventuallyFormula, true); + + mc->check(*probMinFormula); + delete probMinFormula; + + delete mc; } /*! @@ -240,6 +263,7 @@ void testChecking() { mdp->printModelInformationToStream(std::cout); // testCheckingDice(*mdp); + // testCheckingAsynchLeader(*mdp); } else { std::cout << "Input is neither a DTMC nor an MDP." << std::endl; } @@ -253,6 +277,7 @@ int main(const int argc, const char* argv[]) { if (!parseOptions(argc, argv)) { return 0; } + setUp(); LOG4CPLUS_INFO(logger, "StoRM was invoked."); printHeader(argc, argv); diff --git a/src/utility/GraphAnalyzer.h b/src/utility/GraphAnalyzer.h index f104465b0..c7325befd 100644 --- a/src/utility/GraphAnalyzer.h +++ b/src/utility/GraphAnalyzer.h @@ -232,13 +232,12 @@ public: for(auto it = backwardTransitions.beginStateSuccessorsIterator(currentState); it != backwardTransitions.endStateSuccessorsIterator(currentState); ++it) { if (phiStates.get(*it) && !nextStates->get(*it)) { - - // Check whether the predecessor has at only successors in the current state set for one of the + // Check whether the predecessor has only successors in the current state set for one of the // nondeterminstic choices. - for (auto rowIt = nondeterministicChoiceIndices->begin() + *it; rowIt != nondeterministicChoiceIndices->begin() + *it + 1; ++rowIt) { + for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) { bool allSuccessorsInCurrentStates = true; - for (auto colIt = transitionMatrix->beginConstColumnIterator(*rowIt); colIt != transitionMatrix->endConstColumnIterator(*rowIt); ++colIt) { - if (currentStates->get(*colIt)) { + for (auto colIt = transitionMatrix->beginConstColumnIterator(row); colIt != transitionMatrix->endConstColumnIterator(row); ++colIt) { + if (!currentStates->get(*colIt)) { allSuccessorsInCurrentStates = false; break; } @@ -316,7 +315,6 @@ public: for(auto it = backwardTransitions.beginStateSuccessorsIterator(currentState); it != backwardTransitions.endStateSuccessorsIterator(currentState); ++it) { if (phiStates.get(*it) && !statesWithProbability0->get(*it)) { - // Check whether the predecessor has at least one successor in the current state // set for every nondeterministic choice. bool addToStatesWithProbability0 = true; @@ -349,9 +347,70 @@ public: template static void performProb1A(storm::models::AbstractNondeterministicModel& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::storage::BitVector* statesWithProbability1) { - // This result is a rough guess and does not compute all states with probability 1. - // TODO: Check whether it makes sense to implement the precise but complicated algorithm here. - *statesWithProbability1 = psiStates; + // Check for valid parameters. + if (statesWithProbability1 == nullptr) { + LOG4CPLUS_ERROR(logger, "Parameter 'statesWithProbability1' must not be null."); + throw storm::exceptions::InvalidArgumentException("Parameter 'statesWithProbability1' must not be null."); + } + + // Get some temporaries for convenience. + std::shared_ptr> transitionMatrix = model.getTransitionMatrix(); + 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::storage::BitVector* currentStates = new storm::storage::BitVector(model.getNumberOfStates(), true); + + std::vector stack; + stack.reserve(model.getNumberOfStates()); + + bool done = false; + while (!done) { + stack.clear(); + storm::storage::BitVector* nextStates = new storm::storage::BitVector(psiStates); + psiStates.addSetIndicesToList(stack); + + while (!stack.empty()) { + uint_fast64_t currentState = stack.back(); + stack.pop_back(); + + for(auto it = backwardTransitions.beginStateSuccessorsIterator(currentState); it != backwardTransitions.endStateSuccessorsIterator(currentState); ++it) { + if (phiStates.get(*it) && !nextStates->get(*it)) { + // Check whether the predecessor has only successors in the current state set for all of the + // nondeterminstic choices. + bool allSuccessorsInCurrentStatesForAllChoices = true; + for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) { + for (auto colIt = transitionMatrix->beginConstColumnIterator(row); colIt != transitionMatrix->endConstColumnIterator(row); ++colIt) { + if (!currentStates->get(*colIt)) { + allSuccessorsInCurrentStatesForAllChoices = false; + goto afterCheckLoop; + } + } + } + + afterCheckLoop: + // If all successors for all nondeterministic choices are in the current state set, we + // add it to the set of states for the next iteration and perform a backward search from + // that state. + if (allSuccessorsInCurrentStatesForAllChoices) { + nextStates->set(*it, true); + stack.push_back(*it); + } + } + } + } + + // Check whether we need to perform an additional iteration. + if (*currentStates == *nextStates) { + done = true; + } else { + *currentStates = *nextStates; + } + } + + *statesWithProbability1 = *currentStates; + delete currentStates; } }; diff --git a/src/utility/Vector.h b/src/utility/Vector.h index 680990e38..7e5224af4 100644 --- a/src/utility/Vector.h +++ b/src/utility/Vector.h @@ -12,6 +12,11 @@ #include "ConstTemplates.h" #include +#include "log4cplus/logger.h" +#include "log4cplus/loggingmacros.h" + +extern log4cplus::Logger logger; + namespace storm { namespace utility { @@ -56,20 +61,19 @@ void subtractFromConstantOneVector(std::vector* vector) { template void reduceVectorMin(std::vector const& source, std::vector* target, std::vector 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) { + uint_fast64_t currentTargetRow = -1; + for (auto it = source.cbegin(); it != source.cend(); ++it, ++currentSourceRow) { + // Check whether we have considered all source rows for the current target row. + if (filter[currentTargetRow + 1] <= currentSourceRow || currentSourceRow == 0) { ++currentTargetRow; - (*target)[currentTargetRow] = (*source)[currentSourceRow]; + (*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; + (*target)[currentTargetRow] = *it; } } } @@ -77,22 +81,39 @@ void reduceVectorMin(std::vector const& source, std::vector* target, std:: template void reduceVectorMax(std::vector const& source, std::vector* target, std::vector 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) { + uint_fast64_t currentTargetRow = -1; + for (auto it = source.cbegin(); it != source.cend(); ++it, ++currentSourceRow) { + // Check whether we have considered all source rows for the current target row. + if (filter[currentTargetRow + 1] <= currentSourceRow || currentSourceRow == 0) { ++currentTargetRow; - (*target)[currentTargetRow] = (*source)[currentSourceRow]; + (*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; + (*target)[currentTargetRow] = *it; + } + } +} + +template +bool equalModuloPrecision(std::vector const& vectorLeft, std::vector const& vectorRight, T precision, bool relativeError) { + if (vectorLeft.size() != vectorRight.size()) { + LOG4CPLUS_ERROR(logger, "Length 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 i = 0; i < vectorLeft.size(); ++i) { + if (relativeError) { + if (std::abs(vectorLeft[i] - vectorRight[i])/vectorRight[i] > precision) return false; + } else { + if (std::abs(vectorLeft[i] - vectorRight[i]) > precision) return false; } } + + return true; } } //namespace utility