diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h index c48e3b6b3..e35c9bca1 100644 --- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h +++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h @@ -2,6 +2,7 @@ #define STORM_MODELCHECKER_CSL_SPARSEMARKOVAUTOMATONCSLMODELCHECKER_H_ #include +#include #include "src/modelchecker/csl/AbstractModelChecker.h" #include "src/models/MarkovAutomaton.h" @@ -69,35 +70,105 @@ namespace storm { return result; } - std::vector checkTimeBoundedEventually(bool min, storm::storage::BitVector const& goalStates) const { + std::vector checkTimeBoundedEventually(bool min, storm::storage::BitVector const& goalStates, uint_fast64_t lowerBound, uint_fast64_t upperBound) const { if (!this->getModel().isClosed()) { throw storm::exceptions::InvalidArgumentException() << "Unable to compute time-bounded reachability on non-closed Markov automaton."; } // (1) Compute the number of steps we need to take. + ValueType lambda = this->getModel().getMaximalExitRate(); + ValueType delta = (2 * storm::settings::Settings::getInstance()->getOptionByLongName("precision").getArgument(0).getValueAsDouble()) / (upperBound * lambda * lambda); + LOG4CPLUS_INFO(logger, "Determined delta to be " << delta << "."); + + // Get some data fields for convenient access. + std::vector const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices(); + std::vector const& exitRates = this->getModel().getExitRates(); + typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); // (2) Compute four sparse matrices: - // * a matrix A_MSwG with all (discretized!) transitions from Markovian non-goal states to *all other* non-goal states. Note: this matrix has more columns than rows. - // * a matrix A_PSwg with all (non-discretized) transitions from probabilistic non-goal states to other probabilistic non-goal states. This matrix has more rows than columns. + // * a matrix A_MS with all (discretized) transitions from Markovian non-goal states to all Markovian non-goal states. + // * a matrix A_MStoPS with all (discretized) transitions from Markovian non-goal states to all probabilistic non-goal states. + // * a matrix A_PS with all (non-discretized) transitions from probabilistic non-goal states to other probabilistic non-goal states. This matrix has more rows than columns. // * a matrix A_PStoMS with all (non-discretized) transitions from probabilistic non-goal states to all Markovian non-goal states. This matrix may have any shape. + storm::storage::BitVector const& markovianNonGoalStates = this->getModel().getMarkovianStates() & ~goalStates; + storm::storage::BitVector const& probabilisticNonGoalStates = ~this->getModel().getMarkovianStates() & ~goalStates; + + typename storm::storage::SparseMatrix aMarkovian = this->getModel().getTransitionMatrix().getSubmatrix(markovianNonGoalStates, this->getModel().getNondeterministicChoiceIndices(), true); + typename storm::storage::SparseMatrix aMarkovianToProbabilistic = this->getModel().getTransitionMatrix().getSubmatrix(markovianNonGoalStates, probabilisticNonGoalStates, nondeterministicChoiceIndices); + typename storm::storage::SparseMatrix aProbabilistic = this->getModel().getTransitionMatrix().getSubmatrix(probabilisticNonGoalStates, this->getModel().getNondeterministicChoiceIndices()); + typename storm::storage::SparseMatrix aProbabilisticToMarkovian = this->getModel().getTransitionMatrix().getSubmatrix(probabilisticNonGoalStates, markovianNonGoalStates, nondeterministicChoiceIndices); + + // The matrices with transitions from Markovian states need to be digitized. + + // Digitize aMarkovian. Based on whether the transition is a self-loop or not, we apply the two digitization rules. + uint_fast64_t rowIndex = 0; + for (auto state : markovianNonGoalStates) { + typename storm::storage::SparseMatrix::MutableRows row = aMarkovian.getRow(rowIndex); + for (auto element : row) { + ValueType eTerm = std::exp(-exitRates[state] * delta); + if (element.column() == rowIndex) { + element.value() = (storm::utility::constGetOne() - eTerm) * element.value() + eTerm; + } else { + element.value() = (storm::utility::constGetOne() - eTerm) * element.value(); + } + } + ++rowIndex; + } - // (3) Initialize three used vectors: + // Digitize aMarkovianToProbabilistic. As there are no self-loops in this case, we only need to apply the digitization formula for regular successors. + rowIndex = 0; + for (auto state : markovianNonGoalStates) { + typename storm::storage::SparseMatrix::MutableRows row = aMarkovianToProbabilistic.getRow(rowIndex); + for (auto element : row) { + element.value() = (1 - std::exp(-exitRates[state] * delta)) * element.value(); + } + ++rowIndex; + } + + // (3) Initialize two vectors: // * v_PS holds the probability values of probabilistic non-goal states. // * v_MS holds the probability values of Markovian non-goal states. - // * v_all holds the probability values of all non-goal states. This vector is needed for the Markov step. + std::vector vProbabilistic(probabilisticNonGoalStates.getNumberOfSetBits()); + std::vector vMarkovian(markovianNonGoalStates.getNumberOfSetBits()); + + // (4) Compute the two fixed right-hand side vectors, one for Markovian states and one for the probabilistic ones. + std::vector bProbabilistic = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(probabilisticNonGoalStates, nondeterministicChoiceIndices, ~this->getModel().getMarkovianStates() & goalStates, aProbabilistic.getRowCount()); + storm::storage::BitVector probabilisticGoalStates = ~this->getModel().getMarkovianStates() & goalStates; + std::vector bMarkovian; + bMarkovian.reserve(markovianNonGoalStates.getNumberOfSetBits()); + for (auto state : markovianNonGoalStates) { + bMarkovian.push_back(storm::utility::constGetZero()); + + typename storm::storage::SparseMatrix::Rows row = transitionMatrix.getRow(rowIndex); + for (auto element : row) { + bMarkovian.back() += (1 - std::exp(-exitRates[state] * delta)) * element.value(); + } + } + + std::cout << delta << std::endl; + std::cout << markovianNonGoalStates.toString() << std::endl; + std::cout << aMarkovian.toString() << std::endl; + std::cout << aMarkovianToProbabilistic.toString() << std::endl; + std::cout << probabilisticNonGoalStates.toString() << std::endl; + std::cout << aProbabilistic.toString() << std::endl; + std::cout << aProbabilisticToMarkovian.toString() << std::endl; + std::cout << bProbabilistic << std::endl; + std::cout << bMarkovian << std::endl; // (3) Perform value iteration - // * initialize the vectors v_PS, v_MS and v_all. // * loop until the step bound has been reached // * in the loop: - // * perform value iteration using A_PSwG, v_PS and the vector x where x = x_PS + x_MS, x_PS = (A * 1_G)|PS with x_MS = A_PStoMS * v_MS + // * perform value iteration using A_PSwG, v_PS and the vector b where b = (A * 1_G)|PS + A_PStoMS * v_MS // and 1_G being the characteristic vector for all goal states. - // * copy the values from v_PS to the correct positions into v_all - // * perform one timed-step using A_MSwG, v_all and x_MS = (A * 1_G)|MS and obtain v_MS - // * copy the values from v_MS to the correct positions into v_all + // * perform one timed-step using v_MS := A_MSwG * v_MS + A_MStoPS * v_PS + (A * 1_G)|MS // // After the loop, perform one more step of the value iteration for PS states. - // Finally, create the result vector out of 1_G and v_all. + + + + + + // (4) Finally, create the result vector out of 1_G and v_all. // Return dummy vector for the time being. return std::vector(); diff --git a/src/models/MarkovAutomaton.h b/src/models/MarkovAutomaton.h index 67e2a8acf..2d53745f2 100644 --- a/src/models/MarkovAutomaton.h +++ b/src/models/MarkovAutomaton.h @@ -100,6 +100,14 @@ namespace storm { return this->exitRates[state]; } + T const& getMaximalExitRate() const { + T result = storm::utility::constGetZero(); + for (auto markovianState : this->markovianStates) { + result = std::max(result, this->exitRates[markovianState]); + } + return result; + } + storm::storage::BitVector const& getMarkovianStates() const { return this->markovianStates; } diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index b52a63269..db277459c 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -172,6 +173,89 @@ public: uint_fast64_t const* columnPtr; }; + /*! + * A class representing an iterator over a continuous number of rows of the matrix. + */ + class Iterator { + public: + /*! + * Constructs an iterator from the given parameters. + * + * @param valuePtr A pointer to the value of the first element that is to be iterated over. + * @param columnPtr A pointer to the column of the first element that is to be iterated over. + */ + Iterator(T* valuePtr, uint_fast64_t* columnPtr) : valuePtr(valuePtr), columnPtr(columnPtr) { + // Intentionally left empty. + } + + /*! + * Moves the iterator to the next non-zero element. + * + * @return A reference to itself. + */ + Iterator& operator++() { + ++valuePtr; + ++columnPtr; + return *this; + } + + /*! + * Dereferences the iterator by returning a reference to itself. This is needed for making use of the range-based + * for loop over transitions. + * + * @return A reference to itself. + */ + Iterator& operator*() { + return *this; + } + + /*! + * Compares the two iterators for inequality. + * + * @return True iff the given iterator points to a different index as the current iterator. + */ + bool operator!=(Iterator const& other) const { + return this->valuePtr != other.valuePtr; + } + + /*! + * Assigns the position of the given iterator to the current iterator. + * + * @return A reference to itself. + */ + Iterator& operator=(Iterator const& other) { + this->valuePtr = other.valuePtr; + this->columnPtr = other.columnPtr; + return *this; + } + + /*! + * Retrieves the column that is associated with the current non-zero element to which this iterator + * points. + * + * @return The column of the current non-zero element to which this iterator points. + */ + uint_fast64_t column() { + return *columnPtr; + } + + /*! + * Retrieves the value of the current non-zero element to which this iterator points. + * + * @return The value of the current non-zero element to which this iterator points. + */ + T& value() { + return *valuePtr; + } + + private: + // A pointer to the value of the current non-zero element. + T* valuePtr; + + // A pointer to the column of the current non-zero element. + uint_fast64_t* columnPtr; + }; + /*! * This class represents a number of consecutive rows of the matrix. */ @@ -216,6 +300,51 @@ public: // The number of non-zero entries in the rows. uint_fast64_t entryCount; }; + + /*! + * This class represents a number of consecutive rows of the matrix. + */ + class MutableRows { + public: + /*! + * Constructs a row from the given parameters. + * + * @param valuePtr A pointer to the value of the first non-zero element of the rows. + * @param columnPtr A pointer to the column of the first non-zero element of the rows. + * @param entryCount The number of non-zero elements of the rows. + */ + MutableRows(T* valuePtr, uint_fast64_t* columnPtr, uint_fast64_t entryCount) : valuePtr(valuePtr), columnPtr(columnPtr), entryCount(entryCount) { + // Intentionally left empty. + } + + /*! + * Retrieves an iterator that points to the beginning of the rows. + * + * @return An iterator that points to the beginning of the rows. + */ + Iterator begin() { + return Iterator(valuePtr, columnPtr); + } + + /*! + * Retrieves an iterator that points past the last element of the rows. + * + * @return An iterator that points past the last element of the rows. + */ + Iterator end() { + return Iterator(valuePtr + entryCount, columnPtr + entryCount); + } + + private: + // The pointer to the value of the first element. + T* valuePtr; + + // The pointer to the column of the first element. + uint_fast64_t* columnPtr; + + // The number of non-zero entries in the rows. + uint_fast64_t entryCount; + }; /*! * This class represents an iterator to all rows of the matrix. @@ -919,29 +1048,53 @@ public: * @returns A matrix corresponding to a submatrix of the current matrix in which only row groups * and columns given by the row group constraint are kept and all others are dropped. */ - SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices) const { + SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector const& rowGroupIndices, bool insertDiagonalEntries = false) const { + return getSubmatrix(rowGroupConstraint, rowGroupConstraint, rowGroupIndices, insertDiagonalEntries); + } + + /*! + * Creates a submatrix of the current matrix by keeping only row groups and columns in the given + * row group and column constraint, respectively. + * + * @param rowGroupConstraint A bit vector indicating which row groups to keep. + * @param columnConstraint A bit vector indicating which columns to keep. + * @param rowGroupIndices A vector indicating which rows belong to a given row group. + * @returns A matrix corresponding to a submatrix of the current matrix in which only row groups + * and columns given by the row group constraint are kept and all others are dropped. + */ + SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector const& rowGroupIndices, bool insertDiagonalEntries = false) const { LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size)."); - + // 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 : rowGroupConstraint) { subRowCount += rowGroupIndices[index + 1] - rowGroupIndices[index]; for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { + bool foundDiagonalElement = false; + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { - if (rowGroupConstraint.get(columnIndications[j])) { + if (columnConstraint.get(columnIndications[j])) { ++subNonZeroEntries; + + if (index == columnIndications[j]) { + foundDiagonalElement = true; + } } } + + if (insertDiagonalEntries && !foundDiagonalElement) { + ++subNonZeroEntries; + } } } - + LOG4CPLUS_DEBUG(logger, "Determined size of submatrix to be " << subRowCount << "x" << rowGroupConstraint.getNumberOfSetBits() << "."); - + // Create and initialize resulting matrix. SparseMatrix result(subRowCount, rowGroupConstraint.getNumberOfSetBits()); result.initialize(subNonZeroEntries); - + // Create a temporary vector that stores for each index whose bit is set // to true the number of bits that were set before that particular index. std::vector bitsSetBeforeIndex; @@ -950,27 +1103,45 @@ public: // Compute the information to fill this vector. uint_fast64_t lastIndex = 0; uint_fast64_t currentNumberOfSetBits = 0; - for (auto index : rowGroupConstraint) { + + // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows + storm::storage::BitVector columnBitCountConstraint = columnConstraint; + if (insertDiagonalEntries) { + columnBitCountConstraint |= rowGroupConstraint; + } + for (auto index : columnBitCountConstraint) { while (lastIndex <= index) { bitsSetBeforeIndex.push_back(currentNumberOfSetBits); ++lastIndex; } ++currentNumberOfSetBits; } - + // Copy over selected entries. uint_fast64_t rowCount = 0; for (auto index : rowGroupConstraint) { for (uint_fast64_t i = rowGroupIndices[index]; i < rowGroupIndices[index + 1]; ++i) { + bool insertedDiagonalElement = false; + for (uint_fast64_t j = rowIndications[i]; j < rowIndications[i + 1]; ++j) { - if (rowGroupConstraint.get(columnIndications[j])) { + if (columnConstraint.get(columnIndications[j])) { + if (index == columnIndications[j]) { + insertedDiagonalElement = true; + } else if (insertDiagonalEntries && !insertedDiagonalElement && columnIndications[j] > index) { + result.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::constGetZero()); + insertedDiagonalElement = true; + } result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]); } } + if (insertDiagonalEntries && !insertedDiagonalElement) { + result.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::constGetZero()); + } + ++rowCount; } } - + // Finalize sub-matrix and return result. result.finalize(); LOG4CPLUS_DEBUG(logger, "Done creating sub-matrix."); @@ -980,8 +1151,7 @@ public: SparseMatrix getSubmatrix(std::vector const& rowGroupToRowIndexMapping, std::vector const& rowGroupIndices, bool insertDiagonalEntries = true) const { LOG4CPLUS_DEBUG(logger, "Creating a sub-matrix (of unknown size)."); - // First, we need to count how many non-zero entries the resulting matrix will have and reserve space for diagonal - // entries. + // First, we need to count how many non-zero entries the resulting matrix will have and reserve space for diagonal entries. uint_fast64_t subNonZeroEntries = 0; for (uint_fast64_t rowGroupIndex = 0, rowGroupIndexEnd = rowGroupToRowIndexMapping.size(); rowGroupIndex < rowGroupIndexEnd; ++rowGroupIndex) { // Determine which row we need to select from the current row group. @@ -1227,6 +1397,27 @@ public: return getRows(row, row); } + /*! + * Returns an object representing the consecutive rows given by the parameters. + * + * @param startRow The starting row. + * @param endRow The ending row (which is included in the result). + * @return An object representing the consecutive rows given by the parameters. + */ + MutableRows getRows(uint_fast64_t startRow, uint_fast64_t endRow) { + return MutableRows(this->valueStorage.data() + this->rowIndications[startRow], this->columnIndications.data() + this->rowIndications[startRow], this->rowIndications[endRow + 1] - this->rowIndications[startRow]); + } + + /*! + * Returns an object representing the given row. + * + * @param row The chosen row. + * @return An object representing the given row. + */ + MutableRows getRow(uint_fast64_t row) { + return getRows(row, row); + } + /*! * Returns a const iterator to the rows of the matrix. * @@ -1435,7 +1626,7 @@ public: uint_fast64_t currentRealIndex = 0; while (currentRealIndex < colCount) { if (nextIndex < rowIndications[i + 1] && currentRealIndex == columnIndications[nextIndex]) { - result << valueStorage[nextIndex] << "\t"; + result << std::setprecision(8) << valueStorage[nextIndex] << "\t"; ++nextIndex; } else { result << "0\t"; diff --git a/src/storm.cpp b/src/storm.cpp index 8668f01f5..97e5c1b6b 100644 --- a/src/storm.cpp +++ b/src/storm.cpp @@ -478,6 +478,8 @@ int main(const int argc, const char* argv[]) { std::cout << mc.checkExpectedTime(false, markovAutomaton->getLabeledStates("goal")) << std::endl; std::cout << mc.checkLongRunAverage(true, markovAutomaton->getLabeledStates("goal")) << std::endl; std::cout << mc.checkLongRunAverage(false, markovAutomaton->getLabeledStates("goal")) << std::endl; + std::cout << mc.checkTimeBoundedEventually(true, markovAutomaton->getLabeledStates("goal"), 0, 1) << std::endl; + std::cout << mc.checkTimeBoundedEventually(false, markovAutomaton->getLabeledStates("goal"), 0, 1) << std::endl; break; } case storm::models::Unknown: