From 65bf06dd50b664fb77c6acf83536b2f35688d810 Mon Sep 17 00:00:00 2001 From: dehnert Date: Wed, 18 Mar 2015 19:48:04 +0100 Subject: [PATCH] Further steps towards CTMC model checking. Former-commit-id: f057eeb17efea0079ac322198b9c03293fb16049 --- src/modelchecker/AbstractModelChecker.cpp | 2 + .../csl/SparseCtmcCslModelChecker.cpp | 130 ++++++++++++++++-- .../csl/SparseCtmcCslModelChecker.h | 22 ++- .../prctl/SparseDtmcPrctlModelChecker.h | 2 + .../SparsePropositionalModelChecker.cpp | 2 + src/models/sparse/Ctmc.cpp | 19 ++- src/models/sparse/Ctmc.h | 19 ++- src/utility/cli.h | 6 + src/utility/constants.cpp | 4 + src/utility/constants.h | 2 + 10 files changed, 190 insertions(+), 18 deletions(-) diff --git a/src/modelchecker/AbstractModelChecker.cpp b/src/modelchecker/AbstractModelChecker.cpp index c7ac1b18b..7ba1de099 100644 --- a/src/modelchecker/AbstractModelChecker.cpp +++ b/src/modelchecker/AbstractModelChecker.cpp @@ -32,6 +32,8 @@ namespace storm { return this->computeGloballyProbabilities(pathFormula.asGloballyFormula(), qualitative, optimalityType); } else if (pathFormula.isUntilFormula()) { return this->computeUntilProbabilities(pathFormula.asUntilFormula(), qualitative, optimalityType); + } else if (pathFormula.isNextFormula()) { + return this->computeNextProbabilities(pathFormula.asNextFormula(), qualitative, optimalityType); } STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The given formula '" << pathFormula << "' is invalid."); } diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp index 480afd57d..4edac0a7a 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -6,6 +6,7 @@ #include "src/utility/vector.h" #include "src/utility/graph.h" #include "src/utility/solver.h" +#include "src/utility/numerical.h" #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" @@ -42,7 +43,7 @@ namespace storm { ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();; ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); std::pair const& intervalBounds = pathFormula.getIntervalBounds(); - std::unique_ptr result = std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, intervalBounds.first, intervalBounds.second))); + std::unique_ptr result = std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, intervalBounds.first, intervalBounds.second))); return result; } @@ -51,7 +52,7 @@ namespace storm { std::unique_ptr SparseCtmcCslModelChecker::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(pathFormula.getSubformula()); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - std::vector result = SparseDtmcPrctlModelChecker::computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolver); + std::vector result = SparseDtmcPrctlModelChecker::computeNextProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector(), *this->linearEquationSolver); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); } @@ -61,7 +62,7 @@ namespace storm { std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver))); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeUntilProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver))); } template @@ -74,7 +75,7 @@ namespace storm { // If the time bounds are [0, inf], we rather call untimed reachability. storm::utility::ConstantsComparator comparator; if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) { - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, qualitative, *this->linearEquationSolver))); + return this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, qualitative, *this->linearEquationSolver); } // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0 @@ -85,7 +86,8 @@ namespace storm { // If we identify the states that have probability 0 of reaching the target states, we can exclude them from the // further computations. - storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(this->getModel(), this->getModel().getBackwardTransitions(), phiStates, psiStates); + storm::storage::SparseMatrix backwardTransitions = this->getModel().getBackwardTransitions(); + storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(backwardTransitions, phiStates, psiStates); STORM_LOG_INFO("Found " << statesWithProbabilityGreater0.getNumberOfSetBits() << " states with probability greater 0."); storm::storage::BitVector statesWithProbabilityGreater0NonPsi = statesWithProbabilityGreater0 & ~psiStates; STORM_LOG_INFO("Found " << statesWithProbabilityGreater0NonPsi.getNumberOfSetBits() << " 'maybe' states."); @@ -95,26 +97,72 @@ namespace storm { // In this case, the interval is of the form [0, 0]. storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); } else { - // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. - ValueType uniformizationRate = 0; - for (auto const& state : statesWithProbabilityGreater0NonPsi) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); - } - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - if (comparator.isZero(lowerBound)) { // In this case, the interval is of the form [0, t]. // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier. - // Compute the uniformized matrix. - storm::storage::SparseMatrix uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, psiStates, uniformizationRate, exitRates); + // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. + ValueType uniformizationRate = 0; + for (auto const& state : statesWithProbabilityGreater0NonPsi) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + // Compute the uniformized matrix. + storm::storage::SparseMatrix uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), psiStates, uniformizationRate, exitRates); + + // Finally compute the transient probabilities. + std::vector psiStateValues(statesWithProbabilityGreater0.getNumberOfSetBits(), storm::utility::zero()); + storm::utility::vector::setVectorValues(psiStateValues, psiStates % statesWithProbabilityGreater0, storm::utility::one()); + std::vector subresult = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * upperBound, psiStateValues, *this->linearEquationSolver); + storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, subresult); } else if (comparator.isInfinity(upperBound)) { // In this case, the interval is of the form [t, inf] with t != 0. + // Start by computing the (unbounded) reachability probabilities of reaching psi states while + // staying in phi states. + result = this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), backwardTransitions, phiStates, psiStates, qualitative, *this->linearEquationSolver); + + ValueType uniformizationRate = 0; + for (auto const& state : statesWithProbabilityGreater0) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + + // Compute the uniformized matrix. + storm::storage::SparseMatrix uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, storm::storage::BitVector(statesWithProbabilityGreater0.getNumberOfSetBits()), uniformizationRate, exitRates); + + // Finally compute the transient probabilities. + result = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * lowerBound, result, *this->linearEquationSolver); } else { // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf. + // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. + ValueType uniformizationRate = 0; + for (auto const& state : statesWithProbabilityGreater0NonPsi) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Compute the (first) uniformized matrix. + storm::storage::SparseMatrix uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, psiStates, uniformizationRate, exitRates); + + // Start by computing the transient probabilities of reaching a psi state in time t' - t. + std::vector psiStateValues(statesWithProbabilityGreater0.getNumberOfSetBits(), storm::utility::zero()); + storm::utility::vector::setVectorValues(psiStateValues, psiStates % statesWithProbabilityGreater0, storm::utility::one()); + std::vector subresult = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * (upperBound - lowerBound), psiStateValues, *this->linearEquationSolver); + storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, subresult); + + // Then compute the transient probabilities of being in such a state after t time units. For this, + // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix. + storm::storage::BitVector absorbingStates = ~phiStates; + uniformizationRate = 0; + for (auto const& state : ~absorbingStates) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + + // Finally, we compute the second set of transient probabilities. + uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), absorbingStates, uniformizationRate, exitRates); + result = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * lowerBound, result, *this->linearEquationSolver); } } } @@ -151,10 +199,64 @@ namespace storm { } } + template + std::vector SparseCtmcCslModelChecker::computeTransientProbabilities(storm::storage::SparseMatrix const& uniformizedMatrix, ValueType const& lambda, std::vector values, storm::solver::LinearEquationSolver const& linearEquationSolver) const { + // Use Fox-Glynn to get the truncation points and the weights. + std::tuple> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e-300, 1e+300, storm::settings::generalSettings().getPrecision() / 8.0); + + // Scale the weights so they add up to one. + for (auto& element : std::get<3>(foxGlynnResult)) { + element /= std::get<2>(foxGlynnResult); + } + + // Initialize result. + std::vector result; + if (std::get<0>(foxGlynnResult) == 0) { + result = values; + storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]); + } else { + result = std::vector(values.size()); + } + std::vector multiplicationResult(result.size()); + + // Perform the matrix-vector multiplications (without adding) + if (std::get<0>(foxGlynnResult) > 1) { + linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, &result, std::get<0>(foxGlynnResult) - 1, &multiplicationResult); + } + + // For the indices that fall in between the truncation points, we need to perform the matrix-vector + // multiplication, scale and add the result. + ValueType weight = 0; + std::function addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; }; + for (uint_fast64_t index = std::get<0>(foxGlynnResult); index <= std::get<1>(foxGlynnResult); ++index) { + linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, nullptr, 1, &multiplicationResult); + + ValueType weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)]; + storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale); + } + + return result; + } + template std::vector SparseCtmcCslModelChecker::computeUntilProbabilitiesHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolver const& linearEquationSolver) { return SparseDtmcPrctlModelChecker::computeUntilProbabilitiesHelper(transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolver); } + template + storm::storage::SparseMatrix SparseCtmcCslModelChecker::computeProbabilityMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates) { + // Turn the rates into probabilities by scaling each row with the exit rate of the state. + storm::storage::SparseMatrix result(rateMatrix); + for (uint_fast64_t row = 0; row < result.getRowCount(); ++row) { + for (auto& entry : result.getRow(row)) { + entry.setValue(entry.getValue() / exitRates[row]); + } + } + return result; + } + + // Explicitly instantiate the model checker. + template class SparseCtmcCslModelChecker; + } // namespace modelchecker } // namespace storm \ No newline at end of file diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h index a88c36173..3e0966e39 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h @@ -36,8 +36,28 @@ namespace storm { * @param absorbingStates The states that need to be made absorbing. * @param uniformizationRate The rate to be used for uniformization. * @param exitRates The exit rates of all states. + * @return The uniformized matrix. */ - storm::storage::SparseMatrix computeUniformizedMatrix(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& absorbingStates, ValueType uniformizationRate, std::vector const& exitRates); + static storm::storage::SparseMatrix computeUniformizedMatrix(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& absorbingStates, ValueType uniformizationRate, std::vector const& exitRates); + + /*! + * Computes the transient probabilities for lambda time steps. + * + * @param uniformizedMatrix The uniformized transition matrix. + * @param lambda The number of time steps. + * @param values A vector mapping each state to an initial probability. + * @param linearEquationSolver The linear equation solver to use. + * @return The vector of transient probabilities. + */ + std::vector computeTransientProbabilities(storm::storage::SparseMatrix const& uniformizedMatrix, ValueType const& lambda, std::vector values, storm::solver::LinearEquationSolver const& linearEquationSolver) const; + + /*! + * Converts the given rate-matrix into a time-abstract probability matrix. + * + * @param rateMatrix The rate matrix. + * @param exitRates The exit rate vector. + */ + static storm::storage::SparseMatrix computeProbabilityMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates); // An object that is used for solving linear equations and performing matrix-vector multiplication. std::unique_ptr> linearEquationSolver; diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index 7133c12c8..5913b373c 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -14,6 +14,8 @@ namespace storm { template class SparseDtmcPrctlModelChecker : public SparsePropositionalModelChecker { public: + friend class SparseCtmcCslModelChecker; + explicit SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc const& model); explicit SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc const& model, std::unique_ptr>&& linearEquationSolver); diff --git a/src/modelchecker/propositional/SparsePropositionalModelChecker.cpp b/src/modelchecker/propositional/SparsePropositionalModelChecker.cpp index 0e8169d5d..c4c8dc236 100644 --- a/src/modelchecker/propositional/SparsePropositionalModelChecker.cpp +++ b/src/modelchecker/propositional/SparsePropositionalModelChecker.cpp @@ -3,6 +3,7 @@ #include "src/adapters/CarlAdapter.h" #include "src/models/sparse/Dtmc.h" +#include "src/models/sparse/Ctmc.h" #include "src/models/sparse/Mdp.h" #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" @@ -50,6 +51,7 @@ namespace storm { // Explicitly instantiate the template class. template storm::models::sparse::Dtmc const& SparsePropositionalModelChecker::getModelAs() const; + template storm::models::sparse::Ctmc const& SparsePropositionalModelChecker::getModelAs() const; template storm::models::sparse::Mdp const& SparsePropositionalModelChecker::getModelAs() const; template class SparsePropositionalModelChecker; diff --git a/src/models/sparse/Ctmc.cpp b/src/models/sparse/Ctmc.cpp index e6d55cdc8..e495221b1 100644 --- a/src/models/sparse/Ctmc.cpp +++ b/src/models/sparse/Ctmc.cpp @@ -12,7 +12,7 @@ namespace storm { boost::optional> const& optionalTransitionRewardMatrix, boost::optional>> const& optionalChoiceLabeling) : DeterministicModel(storm::models::ModelType::Ctmc, rateMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) { - // Intentionally left empty. + exitRates = createExitRateVector(this->getTransitionMatrix()); } template @@ -21,7 +21,22 @@ namespace storm { boost::optional>&& optionalTransitionRewardMatrix, boost::optional>>&& optionalChoiceLabeling) : DeterministicModel(storm::models::ModelType::Ctmc, std::move(rateMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix), std::move(optionalChoiceLabeling)) { - // Intentionally left empty. + // It is important to refer to the transition matrix here, because the given rate matrix has been move elsewhere. + exitRates = createExitRateVector(this->getTransitionMatrix()); + } + + template + std::vector const& Ctmc::getExitRateVector() const { + return exitRates; + } + + template + std::vector Ctmc::createExitRateVector(storm::storage::SparseMatrix const& rateMatrix) { + std::vector exitRates(rateMatrix.getRowCount()); + for (uint_fast64_t row = 0; row < rateMatrix.getRowCount(); ++row) { + exitRates[row] = rateMatrix.getRowSum(row); + } + return exitRates; } template class Ctmc; diff --git a/src/models/sparse/Ctmc.h b/src/models/sparse/Ctmc.h index abfe0aba5..53be5048e 100644 --- a/src/models/sparse/Ctmc.h +++ b/src/models/sparse/Ctmc.h @@ -52,7 +52,24 @@ namespace storm { Ctmc(Ctmc&& ctmc) = default; Ctmc& operator=(Ctmc&& ctmc) = default; #endif - + /*! + * Retrieves the vector of exit rates of the model. + * + * @return The exit rate vector. + */ + std::vector const& getExitRateVector() const; + + private: + /*! + * Computes the exit rate vector based on the given rate matrix. + * + * @param rateMatrix The rate matrix. + * @return The exit rate vector. + */ + static std::vector createExitRateVector(storm::storage::SparseMatrix const& rateMatrix); + + // A vector containing the exit rates of all states. + std::vector exitRates; }; } // namespace sparse diff --git a/src/utility/cli.h b/src/utility/cli.h index f86acabf8..07e520f68 100644 --- a/src/utility/cli.h +++ b/src/utility/cli.h @@ -67,6 +67,7 @@ log4cplus::Logger printer; #include "src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h" #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" #include "src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h" +#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h" // Headers for counterexample generation. #include "src/counterexamples/MILPMinimalLabelSetGenerator.h" @@ -441,6 +442,11 @@ namespace storm { storm::modelchecker::SparseMdpPrctlModelChecker modelchecker(*mdp); result = modelchecker.check(*formula.get()); #endif + } else if (model->getType() == storm::models::ModelType::Ctmc) { + std::shared_ptr> ctmc = sparseModel->template as>(); + + storm::modelchecker::SparseCtmcCslModelChecker modelchecker(*ctmc); + result = modelchecker.check(*formula.get()); } if (result) { diff --git a/src/utility/constants.cpp b/src/utility/constants.cpp index a69967449..d4625b578 100644 --- a/src/utility/constants.cpp +++ b/src/utility/constants.cpp @@ -110,6 +110,10 @@ namespace storm { return std::abs(value) <= precision; } + bool ConstantsComparator::isInfinity(double const& value) const { + return value == infinity(); + } + bool ConstantsComparator::isEqual(double const& value1, double const& value2) const { return std::abs(value1 - value2) <= precision; } diff --git a/src/utility/constants.h b/src/utility/constants.h index 14da9ffa5..e8573130c 100644 --- a/src/utility/constants.h +++ b/src/utility/constants.h @@ -93,6 +93,8 @@ namespace storm { bool isZero(double const& value) const; + bool isInfinity(double const& value) const; + bool isEqual(double const& value1, double const& value2) const; bool isConstant(double const& value) const;