From e2cfa54d5b9f101e1b091e693d67c427e17f465b Mon Sep 17 00:00:00 2001 From: TimQu Date: Wed, 28 Jun 2017 18:37:36 +0200 Subject: [PATCH 1/2] Fixed an issue when computing expected rewards of Markov Automata --- .../helper/SparseMarkovAutomatonCslHelper.cpp | 36 +++++++++++-------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index c555fd0e8..78f38109a 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -200,7 +200,7 @@ namespace storm { template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { - std::vector stateRewardWeights(transitionMatrix.getRowGroupCount()); + std::vector stateRewardWeights(transitionMatrix.getRowGroupCount(), storm::utility::zero()); for (auto const markovianState : markovianStates) { stateRewardWeights[markovianState] = storm::utility::one() / exitRateVector[markovianState]; } @@ -382,10 +382,10 @@ namespace storm { uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount(); - // First, we need to check which states have infinite expected time (by definition). + // First, we need to check which states have infinite expected reward (by definition). storm::storage::BitVector infinityStates; if (dir == OptimizationDirection::Minimize) { - // If we need to compute the minimum expected times, we have to set the values of those states to infinity that, under all schedulers, + // If we need to compute the minimum expected rewards, we have to set the values of those states to infinity that, under all schedulers, // reach a bottom SCC without a goal state. // So we start by computing all bottom SCCs without goal states. @@ -406,8 +406,7 @@ namespace storm { infinityStates = storm::utility::graph::performProbGreater0A(transitionMatrix, transitionMatrix.getRowGroupIndices(), backwardTransitions, - storm::storage::BitVector( - numberOfStates, true), + ~goalStates, unionOfNonGoalBSccs); } else { // Otherwise, we have no infinity states. @@ -432,8 +431,7 @@ namespace storm { if (!unionOfNonGoalMaximalEndComponents.empty()) { // Now we need to check for which states there exists a scheduler that reaches one of the previously computed states. infinityStates = storm::utility::graph::performProbGreater0E(backwardTransitions, - storm::storage::BitVector( - numberOfStates, true), + ~goalStates, unionOfNonGoalMaximalEndComponents); } else { // Otherwise, we have no infinity states. @@ -447,17 +445,27 @@ namespace storm { std::vector result(numberOfStates); if (!maybeStates.empty()) { + // Then, we can eliminate the rows and columns for all states whose values are already known. - std::vector x(maybeStates.getNumberOfSetBits()); - storm::storage::SparseMatrix submatrix = transitionMatrix.getSubmatrix(true, maybeStates, - maybeStates); - - // Finally, prepare the actual right-hand side. - std::vector b(submatrix.getRowCount()); - storm::utility::vector::selectVectorValues(b, maybeStates, + storm::storage::SparseMatrix submatrix; + std::vector b; + if (infinityStates.empty()) { + submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates); + b.resize(submatrix.getRowCount()); + storm::utility::vector::selectVectorValues(b, maybeStates, transitionMatrix.getRowGroupIndices(), stateActionRewardVector); + } else { + // If there are infinity states, we also have to eliminate the choices that lead from a maybe state to an infinity state + storm::storage::BitVector selectedChoices = transitionMatrix.getRowFilter(maybeStates, ~infinityStates); + submatrix = transitionMatrix.getSubmatrix(false, selectedChoices, maybeStates); + b.resize(submatrix.getRowCount()); + storm::utility::vector::selectVectorValues(b, selectedChoices, stateActionRewardVector); + } + // Initialize the solution vector + std::vector x(maybeStates.getNumberOfSetBits()); + // Solve the corresponding system of equations. std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create( submatrix); From 752af135cd8127105f14582beb8a9feb712c97a0 Mon Sep 17 00:00:00 2001 From: TimQu Date: Thu, 29 Jun 2017 09:38:24 +0200 Subject: [PATCH 2/2] When computing expected rewards for Markov Automata, we now invoke the MDP implementation (instead of the rather inefficient MA implementation) --- .../helper/SparseMarkovAutomatonCslHelper.cpp | 125 ++---------------- .../helper/SparseMarkovAutomatonCslHelper.h | 16 --- 2 files changed, 10 insertions(+), 131 deletions(-) diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index 78f38109a..50e26ddae 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -200,12 +200,16 @@ namespace storm { template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + + // Get a reward model where the state rewards are scaled accordingly std::vector stateRewardWeights(transitionMatrix.getRowGroupCount(), storm::utility::zero()); for (auto const markovianState : markovianStates) { stateRewardWeights[markovianState] = storm::utility::one() / exitRateVector[markovianState]; } std::vector totalRewardVector = rewardModel.getTotalActionRewardVector(transitionMatrix, stateRewardWeights); - return computeExpectedRewards(dir, transitionMatrix, backwardTransitions, psiStates, totalRewardVector, minMaxLinearEquationSolverFactory); + RewardModelType scaledRewardModel(boost::none, std::move(totalRewardVector)); + + return SparseMdpPrctlHelper::computeReachabilityRewards(dir, transitionMatrix, backwardTransitions, scaledRewardModel, psiStates, false, false, minMaxLinearEquationSolverFactory).values; } template @@ -365,120 +369,15 @@ namespace storm { template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + + // Get a reward model representing expected sojourn times std::vector rewardValues(transitionMatrix.getRowCount(), storm::utility::zero()); for (auto const markovianState : markovianStates) { rewardValues[transitionMatrix.getRowGroupIndices()[markovianState]] = storm::utility::one() / exitRateVector[markovianState]; } - return computeExpectedRewards(dir, transitionMatrix, backwardTransitions, psiStates, rewardValues, minMaxLinearEquationSolverFactory); - } - - template - std::vector SparseMarkovAutomatonCslHelper::computeExpectedRewards(OptimizationDirection dir, - storm::storage::SparseMatrix const &transitionMatrix, - storm::storage::SparseMatrix const &backwardTransitions, - storm::storage::BitVector const &goalStates, - std::vector const &stateActionRewardVector, - storm::solver::MinMaxLinearEquationSolverFactory const &minMaxLinearEquationSolverFactory) { - - uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount(); - - // First, we need to check which states have infinite expected reward (by definition). - storm::storage::BitVector infinityStates; - if (dir == OptimizationDirection::Minimize) { - // If we need to compute the minimum expected rewards, we have to set the values of those states to infinity that, under all schedulers, - // reach a bottom SCC without a goal state. - - // So we start by computing all bottom SCCs without goal states. - storm::storage::StronglyConnectedComponentDecomposition sccDecomposition(transitionMatrix, - ~goalStates, true, - true); - - // Now form the union of all these SCCs. - storm::storage::BitVector unionOfNonGoalBSccs(numberOfStates); - for (auto const &scc : sccDecomposition) { - for (auto state : scc) { - unionOfNonGoalBSccs.set(state); - } - } - - // Finally, if this union is non-empty, compute the states such that all schedulers reach some state of the union. - if (!unionOfNonGoalBSccs.empty()) { - infinityStates = storm::utility::graph::performProbGreater0A(transitionMatrix, - transitionMatrix.getRowGroupIndices(), - backwardTransitions, - ~goalStates, - unionOfNonGoalBSccs); - } else { - // Otherwise, we have no infinity states. - infinityStates = storm::storage::BitVector(numberOfStates); - } - } else { - // If we maximize the property, the expected time of a state is infinite, if an end-component without any goal state is reachable. - - // So we start by computing all MECs that have no goal state. - storm::storage::MaximalEndComponentDecomposition mecDecomposition(transitionMatrix, - backwardTransitions, - ~goalStates); - - // Now we form the union of all states in these end components. - storm::storage::BitVector unionOfNonGoalMaximalEndComponents(numberOfStates); - for (auto const &mec : mecDecomposition) { - for (auto const &stateActionPair : mec) { - unionOfNonGoalMaximalEndComponents.set(stateActionPair.first); - } - } - - if (!unionOfNonGoalMaximalEndComponents.empty()) { - // Now we need to check for which states there exists a scheduler that reaches one of the previously computed states. - infinityStates = storm::utility::graph::performProbGreater0E(backwardTransitions, - ~goalStates, - unionOfNonGoalMaximalEndComponents); - } else { - // Otherwise, we have no infinity states. - infinityStates = storm::storage::BitVector(numberOfStates); - } - } - // Now we identify the states for which values need to be computed. - storm::storage::BitVector maybeStates = ~(goalStates | infinityStates); - - // Create resulting vector. - std::vector result(numberOfStates); - - if (!maybeStates.empty()) { - - // Then, we can eliminate the rows and columns for all states whose values are already known. - storm::storage::SparseMatrix submatrix; - std::vector b; - if (infinityStates.empty()) { - submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates); - b.resize(submatrix.getRowCount()); - storm::utility::vector::selectVectorValues(b, maybeStates, - transitionMatrix.getRowGroupIndices(), - stateActionRewardVector); - } else { - // If there are infinity states, we also have to eliminate the choices that lead from a maybe state to an infinity state - storm::storage::BitVector selectedChoices = transitionMatrix.getRowFilter(maybeStates, ~infinityStates); - submatrix = transitionMatrix.getSubmatrix(false, selectedChoices, maybeStates); - b.resize(submatrix.getRowCount()); - storm::utility::vector::selectVectorValues(b, selectedChoices, stateActionRewardVector); - } - - // Initialize the solution vector - std::vector x(maybeStates.getNumberOfSetBits()); - - // Solve the corresponding system of equations. - std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create( - submatrix); - solver->solveEquations(dir, x, b); - - // Set values of resulting vector according to previous result and return the result. - storm::utility::vector::setVectorValues(result, maybeStates, x); - } - - storm::utility::vector::setVectorValues(result, goalStates, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::infinity()); - - return result; + storm::models::sparse::StandardRewardModel rewardModel(boost::none, std::move(rewardValues)); + + return SparseMdpPrctlHelper::computeReachabilityRewards(dir, transitionMatrix, backwardTransitions, rewardModel, psiStates, false, false, minMaxLinearEquationSolverFactory).values; } template @@ -558,8 +457,6 @@ namespace storm { template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec); - template std::vector SparseMarkovAutomatonCslHelper::computeExpectedRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& goalStates, std::vector const& stateRewards, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template std::vector SparseMarkovAutomatonCslHelper::computeBoundedUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair const& boundsPair, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template std::vector SparseMarkovAutomatonCslHelper::computeUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); @@ -574,8 +471,6 @@ namespace storm { template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec); - template std::vector SparseMarkovAutomatonCslHelper::computeExpectedRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& goalStates, std::vector const& stateRewards, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - } } } diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index d61baf0c2..20d6c9598 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -55,22 +55,6 @@ namespace storm { */ template static ValueType computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::MaximalEndComponent const& mec); - - /*! - * Computes the expected reward that is gained from each state before entering any of the goal states. - * - * @param dir Indicates whether minimal or maximal rewards are to be computed. - * @param transitionMatrix The transition matrix of the underlying Markov automaton. - * @param backwardTransitions The reversed transition relation of the underlying Markov automaton. - * @param goalStates The goal states that define until which point rewards are gained. - * @param stateRewards A vector that defines the reward gained in each state. For probabilistic states, - * this is an instantaneous reward that is fully gained and for Markovian states the actually gained - * reward is dependent on the expected time to stay in the state, i.e. it is gouverned by the exit rate - * of the state. - * @return A vector that contains the expected reward for each state of the model. - */ - template - static std::vector computeExpectedRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& goalStates, std::vector const& stateRewards, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); }; }