From 1130efe0dcf148fd1a3c6ed639864707d7cb319a Mon Sep 17 00:00:00 2001 From: dehnert Date: Wed, 17 Jun 2015 18:11:41 +0200 Subject: [PATCH] step towards steady-state for CTMCs Former-commit-id: 4ab4d6b8b6869d066ab41646dfbc6d663ad6d9d9 --- .../csl/SparseCtmcCslModelChecker.cpp | 9 ++ .../csl/SparseCtmcCslModelChecker.h | 1 + .../prctl/SparseDtmcPrctlModelChecker.cpp | 122 +++++++++--------- .../prctl/SparseDtmcPrctlModelChecker.h | 2 +- src/solver/GmmxxLinearEquationSolver.cpp | 2 +- .../expressions/ExprtkExpressionEvaluator.cpp | 2 +- 6 files changed, 77 insertions(+), 61 deletions(-) diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp index ecd550f4c..f915a4842 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -455,6 +455,15 @@ namespace storm { return std::unique_ptr(new ExplicitQuantitativeCheckResult(SparseDtmcPrctlModelChecker::computeReachabilityRewardsHelper(probabilityMatrix, modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative))); } + template + std::unique_ptr SparseCtmcCslModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { + std::unique_ptr subResultPointer = this->check(stateFormula); + ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); + storm::storage::SparseMatrix probabilityMatrix = computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); + + return std::unique_ptr(new ExplicitQuantitativeCheckResult(SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(this->getModel(), probabilityMatrix, subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory))); + } + // Explicitly instantiate the model checker. template class SparseCtmcCslModelChecker; diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h index 5108a183f..529c23128 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h @@ -29,6 +29,7 @@ namespace storm { virtual std::unique_ptr computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; virtual std::unique_ptr computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; + virtual std::unique_ptr computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; protected: storm::models::sparse::Ctmc const& getModel() const override; diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index f16951d12..ae76c832f 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -304,23 +304,22 @@ namespace storm { } template - std::vector SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const { + std::vector SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(storm::models::sparse::DeterministicModel const& model, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory 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(numOfStates, storm::utility::zero()); } - // Likewise, if all bits are set, we can avoid the computation and set. - if ((~psiStates).empty()) { + // Likewise, if all bits are set, we can avoid the computation. + if (psiStates.full()) { return std::vector(numOfStates, storm::utility::one()); } // Start by decomposing the DTMC into its BSCCs. - storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(this->getModel(), false, true); + storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(model, false, true); // Get some data members for convenience. - typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); ValueType one = storm::utility::one(); ValueType zero = storm::utility::zero(); @@ -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 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 bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one); { - auto solver = this->linearEquationSolverFactory->create(bsccEquationSystem); + std::unique_ptr> solver = linearEquationSolverFactory.create(bsccEquationSystem); solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); } - //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 + // 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 + // necessarily do in the solution of the equation system. std::vector bsccLra(bsccDecomposition.size(), zero); std::vector 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]; } - - //calculate LRA for states not in bsccs as expected reachability rewards - //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 sum of LRAs in BSCC weighted by the reachability probability of the BSCC - - std::vector rewardRightSide; - rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); - - for (auto state : statesNotInBsccs) { - ValueType reward = zero; - for (auto entry : transitionMatrix.getRow(state)) { - if (statesInBsccs.get(entry.getColumn())) { - reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]]; - } - } - rewardRightSide.push_back(reward); - } - - storm::storage::SparseMatrix rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); - rewardEquationSystemMatrix.convertToEquationSystem(); - - std::vector rewardSolution(rewardEquationSystemMatrix.getColumnCount(), one); - - { - auto solver = this->linearEquationSolverFactory->create(rewardEquationSystemMatrix); - solver->solveEquationSystem(rewardSolution, rewardRightSide); - } - - // now fill the result vector - std::vector result(numOfStates); - - auto rewardSolutionIter = rewardSolution.begin(); - for (size_t state = 0; state < numOfStates; ++state) { - if (statesInBsccs.get(state)) { - //assign the value of the bscc the state is in - result[state] = bsccLra[stateToBsccIndexMap[state]]; - } else { - assert(rewardSolutionIter != rewardSolution.end()); - //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; - } - } + + std::vector rewardSolution; + if (!statesNotInBsccs.empty()) { + // Calculate LRA for states not in bsccs as expected reachability rewards. + // 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 + // of the BSCC. + + std::vector rewardRightSide; + rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); + + for (auto state : statesNotInBsccs) { + ValueType reward = zero; + for (auto entry : transitionMatrix.getRow(state)) { + if (statesInBsccs.get(entry.getColumn())) { + reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]]; + } + } + rewardRightSide.push_back(reward); + } + + storm::storage::SparseMatrix rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); + rewardEquationSystemMatrix.convertToEquationSystem(); + + rewardSolution = std::vector(rewardEquationSystemMatrix.getColumnCount(), one); + + { + std::unique_ptr> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix); + solver->solveEquationSystem(rewardSolution, rewardRightSide); + } + } + + // Fill the result vector. + std::vector result(numOfStates); + auto rewardSolutionIter = rewardSolution.begin(); + for (size_t state = 0; state < numOfStates; ++state) { + if (statesInBsccs.get(state)) { + //assign the value of the bscc the state is in + result[state] = bsccLra[stateToBsccIndexMap[state]]; + } else { + STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution."); + //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 subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(subResult.getTruthValuesVector(), qualitative))); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory))); } diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index b55f00faa..212b99588 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -41,7 +41,7 @@ namespace storm { std::vector computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const; std::vector computeCumulativeRewardsHelper(uint_fast64_t stepBound) const; static std::vector computeReachabilityRewardsHelper(storm::storage::SparseMatrix const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, bool qualitative); - std::vector computeLongRunAverageHelper(storm::storage::BitVector const& psiStates, bool qualitative) const; + static std::vector computeLongRunAverageHelper(storm::models::sparse::DeterministicModel const& model, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); static ValueType computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index 789ee6e55..4fb0a7ed8 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -55,7 +55,7 @@ namespace storm { template void GmmxxLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* 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."); } diff --git a/src/storage/expressions/ExprtkExpressionEvaluator.cpp b/src/storage/expressions/ExprtkExpressionEvaluator.cpp index b621b64a6..4e79a8dc0 100644 --- a/src/storage/expressions/ExprtkExpressionEvaluator.cpp +++ b/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; }