Started implementation of time-bounded reachability of Markov automata.

Former-commit-id: 512bb117a6

src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h

@@ -2,6 +2,7 @@
 #define STORM_MODELCHECKER_CSL_SPARSEMARKOVAUTOMATONCSLMODELCHECKER_H_
 
 #include <stack>
+#include <iomanip>
 
 #include "src/modelchecker/csl/AbstractModelChecker.h"
 #include "src/models/MarkovAutomaton.h"
@@ -69,35 +70,105 @@ namespace storm {
                     return result;
                 }
                 
-                std::vector<ValueType> checkTimeBoundedEventually(bool min, storm::storage::BitVector const& goalStates) const {
+                std::vector<ValueType> 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<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+                    std::vector<ValueType> const& exitRates = this->getModel().getExitRates();
+                    typename storm::storage::SparseMatrix<ValueType> 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<ValueType> aMarkovian = this->getModel().getTransitionMatrix().getSubmatrix(markovianNonGoalStates, this->getModel().getNondeterministicChoiceIndices(), true);
+                    typename storm::storage::SparseMatrix<ValueType> aMarkovianToProbabilistic = this->getModel().getTransitionMatrix().getSubmatrix(markovianNonGoalStates, probabilisticNonGoalStates, nondeterministicChoiceIndices);
+                    typename storm::storage::SparseMatrix<ValueType> aProbabilistic = this->getModel().getTransitionMatrix().getSubmatrix(probabilisticNonGoalStates, this->getModel().getNondeterministicChoiceIndices());
+                    typename storm::storage::SparseMatrix<ValueType> 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<ValueType>::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<ValueType>() - eTerm) * element.value() + eTerm;
+                            } else {
+                                element.value() = (storm::utility::constGetOne<ValueType>() - eTerm) * element.value();
+                            }
+                        }
+                        ++rowIndex;
+                    }
+                    
+                    // 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<ValueType>::MutableRows row = aMarkovianToProbabilistic.getRow(rowIndex);
+                        for (auto element : row) {
+                            element.value() = (1 - std::exp(-exitRates[state] * delta)) * element.value();
+                        }
+                        ++rowIndex;
+                    }
                     
-                    // (3) Initialize three used vectors:
+                    // (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<ValueType> vProbabilistic(probabilisticNonGoalStates.getNumberOfSetBits());
+                    std::vector<ValueType> vMarkovian(markovianNonGoalStates.getNumberOfSetBits());
+                    
+                    // (4) Compute the two fixed right-hand side vectors, one for Markovian states and one for the probabilistic ones.
+                    std::vector<ValueType> bProbabilistic = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(probabilisticNonGoalStates, nondeterministicChoiceIndices, ~this->getModel().getMarkovianStates() & goalStates, aProbabilistic.getRowCount());
+                    storm::storage::BitVector probabilisticGoalStates = ~this->getModel().getMarkovianStates() & goalStates;
+                    std::vector<ValueType> bMarkovian;
+                    bMarkovian.reserve(markovianNonGoalStates.getNumberOfSetBits());
+                    for (auto state : markovianNonGoalStates) {
+                        bMarkovian.push_back(storm::utility::constGetZero<ValueType>());
+                        
+                        typename storm::storage::SparseMatrix<ValueType>::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<ValueType>();

src/models/MarkovAutomaton.h

@@ -100,6 +100,14 @@ namespace storm {
                 return this->exitRates[state];
             }
             
+            T const& getMaximalExitRate() const {
+                T result = storm::utility::constGetZero<T>();
+                for (auto markovianState : this->markovianStates) {
+                    result = std::max(result, this->exitRates[markovianState]);
+                }
+                return result;
+            }
+            
             storm::storage::BitVector const& getMarkovianStates() const {
                 return this->markovianStates;
             }

src/storage/SparseMatrix.h

@@ -15,6 +15,7 @@
 #include <new>
 #include <algorithm>
 #include <iostream>
+#include <iomanip>
 #include <iterator>
 #include <set>
 #include <cstdint>
@@ -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.
      */
@@ -217,6 +301,51 @@ public:
         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,7 +1048,21 @@ 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<uint_fast64_t> const& rowGroupIndices) const {
+	SparseMatrix getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, std::vector<uint_fast64_t> 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<uint_fast64_t> 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.
@@ -928,12 +1071,22 @@ public:
 		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() << ".");
@@ -950,7 +1103,13 @@ 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;
@@ -962,11 +1121,23 @@ public:
 		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<T>());
+                            insertedDiagonalElement = true;
+                        }
 						result.addNextValue(rowCount, bitsSetBeforeIndex[columnIndications[j]], valueStorage[j]);
 					}
 				}
+                if (insertDiagonalEntries && !insertedDiagonalElement) {
+                    result.addNextValue(rowCount, bitsSetBeforeIndex[index], storm::utility::constGetZero<T>());
+                }
+                
 				++rowCount;
 			}
 		}
@@ -980,8 +1151,7 @@ public:
     SparseMatrix getSubmatrix(std::vector<uint_fast64_t> const& rowGroupToRowIndexMapping, std::vector<uint_fast64_t> 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";

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: