steady state working for CTMCs

Former-commit-id: 9b2cf09400

src/modelchecker/csl/SparseCtmcCslModelChecker.cpp

@@ -459,10 +459,18 @@ namespace storm {
         std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            std::cout << this->getModel().getTransitionMatrix() << std::endl;
-            storm::storage::SparseMatrix<ValueType> probabilityMatrix = computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector());
-                        
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(this->getModel(), probabilityMatrix, subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory)));
+
+            // Since we use the uniformized matrix to compute the steady state probabilities, we need to build it first.
+            ValueType uniformizationRate = 0;
+            for (auto const& rate : this->getModel().getExitRateVector()) {
+                uniformizationRate = std::max(uniformizationRate, rate);
+            }
+            uniformizationRate *= 1.02;
+            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+            
+            storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector());
+            
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(this->getModel(), uniformizedMatrix, subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory)));
         }
         
         // Explicitly instantiate the model checker.

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp

@@ -343,7 +343,7 @@ namespace storm {
 			// Calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of
             // the transposed transition matrix for the bsccs
 			storm::storage::SparseMatrix<ValueType> bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true);
-
+            
 			// Subtract identity matrix.
 			for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
 				for (auto& entry : bsccEquationSystem.getRow(row)) {
@@ -357,16 +357,14 @@ namespace storm {
             // introduced when subtracting the identity matrix.
 			bsccEquationSystem = bsccEquationSystem.transpose();
 
-            std::cout << bsccEquationSystem << std::endl;
-
             // Add a row to the matrix that expresses that the sum over all entries needs to be one.
             storm::storage::SparseMatrixBuilder<ValueType> builder(std::move(bsccEquationSystem));
             typename storm::storage::SparseMatrixBuilder<ValueType>::index_type row = builder.getLastRow();
             for (uint_fast64_t i = 0; i <= row; ++i) {
                 builder.addNextValue(row + 1, i, 1);
             }
+            builder.addNextValue(row + 1, row + 1, 0);
             bsccEquationSystem = builder.build();
-            std::cout << bsccEquationSystem << std::endl;
             
 			std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
             bsccEquationSystemRightSide.back() = one;
@@ -375,10 +373,6 @@ namespace storm {
                 std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(bsccEquationSystem);
 				solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide);
 			}
-            
-            for (auto const& elem : bsccEquationSystemSolution) {
-                std::cout << "sol " << elem << std::endl;
-            }
 
 			// Calculate LRA Value for each BSCC from steady state distribution in BSCCs.
 			// We have to scale the results, as the probabilities for each BSCC have to sum up to one, which they don't

src/storage/SparseMatrix.cpp

@@ -861,6 +861,22 @@ namespace storm {
         }
 #endif
 
+        template<typename ValueType>
+        void SparseMatrix<ValueType>::multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const {
+            const_iterator it = this->begin();
+            const_iterator ite;
+            std::vector<index_type>::const_iterator rowIterator = rowIndications.begin();
+            std::vector<index_type>::const_iterator rowIteratorEnd = rowIndications.end();
+            
+            uint_fast64_t currentRow = 0;
+            for (; rowIterator != rowIteratorEnd - 1; ++rowIterator) {
+                for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) {
+                    result[it->getColumn()] += it->getValue() * vector[currentRow];
+                }
+                ++currentRow;
+            }
+        }
+        
         template<typename ValueType>
         std::size_t SparseMatrix<ValueType>::getSizeInBytes() const {
             uint_fast64_t size = sizeof(*this);

src/storage/SparseMatrix.h

@@ -671,6 +671,17 @@ namespace storm {
              * @return The product of the matrix and the given vector as the content of the given result vector.
              */
             void multiplyWithVector(std::vector<value_type> const& vector, std::vector<value_type>& result) const;
+
+            /*!
+             * Multiplies the vector to the matrix from the left and writes the result to the given result vector.
+             *
+             * @param vector The vector with which the matrix is to be multiplied. This vector is interpreted as being
+             * a row vector.
+             * @param result The vector that is supposed to hold the result of the multiplication after the operation.
+             * @return The product of the matrix and the given vector as the content of the given result vector. The 
+             * result is to be interpreted as a row vector.
+             */
+            void multiplyVectorWithMatrix(std::vector<value_type> const& vector, std::vector<value_type>& result) const;
             
             /*!
              * Multiplies the matrix with the given vector in a sequential way and writes the result to the given result