step towards steady-state for CTMCs

Former-commit-id: 4ab4d6b8b6

src/modelchecker/csl/SparseCtmcCslModelChecker.cpp

@@ -455,6 +455,15 @@ namespace storm {
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(probabilityMatrix, modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative)));
         }
         
+        template<class ValueType>
+        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();
+            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)));
+        }
+        
         // Explicitly instantiate the model checker.
         template class SparseCtmcCslModelChecker<double>;
         

src/modelchecker/csl/SparseCtmcCslModelChecker.h

@@ -29,6 +29,7 @@ namespace storm {
             virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
             virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
 
         protected:
             storm::models::sparse::Ctmc<ValueType> const& getModel() const override;

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp

@@ -304,23 +304,22 @@ namespace storm {
         }
         
 		template<typename ValueType>
-		std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const {
+        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeLongRunAverageHelper(storm::models::sparse::DeterministicModel<ValueType> const& model, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
 			// If there are no goal states, we avoid the computation and directly return zero.
-			auto numOfStates = this->getModel().getNumberOfStates();
+			uint_fast64_t numOfStates = transitionMatrix.getRowCount();
 			if (psiStates.empty()) {
 				return std::vector<ValueType>(numOfStates, storm::utility::zero<ValueType>());
 			}
 
-			// Likewise, if all bits are set, we can avoid the computation and set.
+			// Likewise, if all bits are set, we can avoid the computation.
-			if ((~psiStates).empty()) {
+			if (psiStates.full()) {
 				return std::vector<ValueType>(numOfStates, storm::utility::one<ValueType>());
 			}
 
 			// Start by decomposing the DTMC into its BSCCs.
-			storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(this->getModel(), false, true);
+			storm::storage::StronglyConnectedComponentDecomposition<double> bsccDecomposition(model, false, true);
 
 			// Get some data members for convenience.
-			typename storm::storage::SparseMatrix<ValueType> const& transitionMatrix = this->getModel().getTransitionMatrix();
 			ValueType one = storm::utility::one<ValueType>();
 			ValueType zero = storm::utility::zero<ValueType>();
 
@@ -341,10 +340,11 @@ namespace storm {
 
 			storm::storage::BitVector statesNotInBsccs = ~statesInBsccs;
 
-			// calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of the transposed transition matrix for the bsccs
+			// 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
+			// Subtract identity matrix.
 			for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) {
 				for (auto& entry : bsccEquationSystem.getRow(row)) {
 					if (entry.getColumn() == row) {
@@ -352,18 +352,21 @@ namespace storm {
 					}
 				}
 			}
-			//now transpose, this internally removes all explicit zeros from the matrix that where introduced when subtracting the identity matrix
+			// Now transpose the matrix. This internally removes all explicit zeros from the matrix that were.
+            // introduced when subtracting the identity matrix.
 			bsccEquationSystem = bsccEquationSystem.transpose();
-
+            
 			std::vector<ValueType> bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero);
 			std::vector<ValueType> bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one);
 			{
-				auto solver = this->linearEquationSolverFactory->create(bsccEquationSystem);
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(bsccEquationSystem);
 				solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide);
 			}
 
-			//calculate LRA Value for each BSCC from steady state distribution in BSCCs
+			// Calculate LRA Value for each BSCC from steady state distribution in BSCCs.
-			// we have to do some scaling, as the probabilities for each BSCC have to sum up to one, which they don't necessarily do in the solution of the equation system
+			// We have to scale the results, as the probabilities for each BSCC have to sum up to one, which they don't
+            // necessarily do in the solution of the equation system.
 			std::vector<ValueType> bsccLra(bsccDecomposition.size(), zero);
 			std::vector<ValueType> bsccTotalValue(bsccDecomposition.size(), zero);
 			size_t i = 0;
@@ -376,50 +379,53 @@ namespace storm {
 			for (i = 0; i < bsccLra.size(); ++i) {
 				bsccLra[i] /= bsccTotalValue[i];
 			}
-
+            std::vector<ValueType> rewardSolution;
-			//calculate LRA for states not in bsccs as expected reachability rewards
+            if (!statesNotInBsccs.empty()) {
-			//target states are states in bsccs, transition reward is the lra of the bscc for each transition into a bscc and 0 otherwise
+                // Calculate LRA for states not in bsccs as expected reachability rewards.
-			//this corresponds to sum of LRAs in BSCC weighted by the reachability probability of the BSCC
+                // Target states are states in bsccs, transition reward is the lra of the bscc for each transition into a
-
+                // bscc and 0 otherwise. This corresponds to the sum of LRAs in BSCC weighted by the reachability probability
-			std::vector<ValueType> rewardRightSide;
+                // of the BSCC.
-			rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
+                std::vector<ValueType> rewardRightSide;
-
+                rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits());
-			for (auto state : statesNotInBsccs) {
+                for (auto state : statesNotInBsccs) {
-				ValueType reward = zero;
+                    ValueType reward = zero;
-				for (auto entry : transitionMatrix.getRow(state)) {
+                    for (auto entry : transitionMatrix.getRow(state)) {
-					if (statesInBsccs.get(entry.getColumn())) {
+                        if (statesInBsccs.get(entry.getColumn())) {
-						reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]];
+                            reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]];
-					}
+                        }
-				}
+                    }
-				rewardRightSide.push_back(reward);
+                    rewardRightSide.push_back(reward);
-			}
+                }
-
+                storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
-			storm::storage::SparseMatrix<ValueType> rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true);
+                rewardEquationSystemMatrix.convertToEquationSystem();
-			rewardEquationSystemMatrix.convertToEquationSystem();
+                rewardSolution = std::vector<ValueType>(rewardEquationSystemMatrix.getColumnCount(), one);
-
+                {
-			std::vector<ValueType> rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one);
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix);
-
+                    solver->solveEquationSystem(rewardSolution, rewardRightSide);
-			{
+                }
-				auto solver = this->linearEquationSolverFactory->create(rewardEquationSystemMatrix);
+            }
-				solver->solveEquationSystem(rewardSolution, rewardRightSide);
+            // Fill the result vector.
-			}
+            std::vector<ValueType> result(numOfStates);
-
+            auto rewardSolutionIter = rewardSolution.begin();
-			// now fill the result vector
+            for (size_t state = 0; state < numOfStates; ++state) {
-			std::vector<ValueType> result(numOfStates);
+                if (statesInBsccs.get(state)) {
-
+                    //assign the value of the bscc the state is in
-			auto rewardSolutionIter = rewardSolution.begin();
+                    result[state] = bsccLra[stateToBsccIndexMap[state]];
-			for (size_t state = 0; state < numOfStates; ++state) {
+                } else {
-				if (statesInBsccs.get(state)) {
+                    STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution.");
-					//assign the value of the bscc the state is in
+                    //take the value from the reward computation
-					result[state] = bsccLra[stateToBsccIndexMap[state]];
+                    //since the n-th state not in any bscc is the n-th entry in rewardSolution we can just take the next value from the iterator
-				} else {
+                    result[state] = *rewardSolutionIter;
-					assert(rewardSolutionIter != rewardSolution.end());
+                    ++rewardSolutionIter;
-					//take the value from the reward computation
+                }
-					//since the n-th state not in any bscc is the n-th entry in rewardSolution we can just take the next value from the iterator
+            }
-					result[state] = *rewardSolutionIter;
-					++rewardSolutionIter;
-				}
-			}
 
 			return result;
 
@@ -455,7 +461,7 @@ namespace storm {
 			std::unique_ptr<CheckResult> subResultPointer = this->check(stateFormula);
 			ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
 
-			return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(subResult.getTruthValuesVector(), qualitative)));
+			return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeLongRunAverageHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory)));
 		}
 
 

src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h

@@ -41,7 +41,7 @@ namespace storm {
             std::vector<ValueType> computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(uint_fast64_t stepBound) const;
             static std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
-			std::vector<ValueType> computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const;
+			static std::vector<ValueType> computeLongRunAverageHelper(storm::models::sparse::DeterministicModel<ValueType> const& model, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
 
 			static ValueType computeLraForBSCC(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
 

src/solver/GmmxxLinearEquationSolver.cpp

@@ -55,7 +55,7 @@ namespace storm {
         
         template<typename ValueType>
         void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner " << preconditionerToString() << "'.");
+            LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner '" << preconditionerToString() << "' (max. " << maximalNumberOfIterations << " iterations).");
             if (method == SolutionMethod::Jacobi && preconditioner != Preconditioner::None) {
                 LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored.");
             }

src/storage/expressions/ExprtkExpressionEvaluator.cpp

@@ -49,7 +49,7 @@ namespace storm {
             CompiledExpressionType& compiledExpression = result.first->second;
             compiledExpression.register_symbol_table(symbolTable);
             bool parsingOk = parser.compile(ToExprtkStringVisitor().toString(expression), compiledExpression);
-            STORM_LOG_ASSERT(parsingOk, "Expression was not properly parsed by ExprTk.");
+            STORM_LOG_ASSERT(parsingOk, "Expression was not properly parsed by ExprTk: " << *expression);
             return compiledExpression;
         }