From 9d138d86f7b0db84ba23688202f537a214f7d8af Mon Sep 17 00:00:00 2001 From: dehnert Date: Fri, 14 Aug 2015 21:52:01 +0200 Subject: [PATCH] further work on creating helper classes for model checking tasks Former-commit-id: 12cd17fb267a4e30118f15120f7123737749a084 --- src/builder/ExplicitPrismModelBuilder.cpp | 6 +- .../csl/SparseCtmcCslModelChecker.cpp | 639 +---------------- .../csl/SparseCtmcCslModelChecker.h | 61 +- .../SparseMarkovAutomatonCslModelChecker.cpp | 2 +- .../csl/helper/HybridCtmcCslHelper.cpp | 0 .../csl/helper/HybridCtmcCslHelper.h | 0 .../csl/helper/SparseCtmcCslHelper.cpp | 648 ++++++++++++++++++ .../csl/helper/SparseCtmcCslHelper.h | 63 +- .../prctl/HybridDtmcPrctlModelChecker.cpp | 254 +------ .../prctl/HybridDtmcPrctlModelChecker.h | 14 +- .../prctl/HybridMdpPrctlModelChecker.cpp | 21 +- .../prctl/HybridMdpPrctlModelChecker.h | 2 + .../prctl/SparseDtmcPrctlModelChecker.cpp | 3 +- .../prctl/SparseDtmcPrctlModelChecker.h | 7 +- .../prctl/SparseMdpPrctlModelChecker.cpp | 33 +- .../prctl/helper/HybridDtmcPrctlHelper.cpp | 241 +++++++ .../prctl/helper/HybridDtmcPrctlHelper.h | 38 + .../prctl/helper/HybridMdpPrctlHelper.cpp | 24 +- .../prctl/helper/HybridMdpPrctlHelper.h | 16 +- .../prctl/helper/SparseDtmcPrctlHelper.cpp | 20 +- .../prctl/helper/SparseDtmcPrctlHelper.h | 5 + .../prctl/helper/SparseMdpPrctlHelper.cpp | 4 +- .../prctl/helper/SparseMdpPrctlHelper.h | 4 +- 23 files changed, 1118 insertions(+), 987 deletions(-) create mode 100644 src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp create mode 100644 src/modelchecker/csl/helper/HybridCtmcCslHelper.h diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp index 469520364..601bbf991 100644 --- a/src/builder/ExplicitPrismModelBuilder.cpp +++ b/src/builder/ExplicitPrismModelBuilder.cpp @@ -173,9 +173,9 @@ namespace storm { } else { rewardModelName = ""; } - rewardModels.emplace(rewardModelName, rewardModel.hasStateRewards() ? std::move(modelComponents.stateRewards) : boost::optional>(), - boost::optional>(), - rewardModel.hasTransitionRewards() ? std::move(modelComponents.transitionRewardMatrix) : boost::optional>()); + rewardModels.emplace(rewardModelName, storm::models::sparse::StandardRewardModel(rewardModel.hasStateRewards() ? std::move(modelComponents.stateRewards) : boost::optional>(), + boost::optional>(), + rewardModel.hasTransitionRewards() ? std::move(modelComponents.transitionRewardMatrix) : boost::optional>())); } switch (program.getModelType()) { diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp index 1dba53698..0597644e0 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -1,19 +1,16 @@ #include "src/modelchecker/csl/SparseCtmcCslModelChecker.h" -#include +#include "src/modelchecker/csl/helper/SparseCtmcCslHelper.h" +#include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h" #include "src/utility/macros.h" #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" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" -#include "src/storage/StronglyConnectedComponentDecomposition.h" - #include "src/exceptions/InvalidStateException.h" #include "src/exceptions/InvalidPropertyException.h" #include "src/exceptions/NotImplementedException.h" @@ -51,15 +48,16 @@ namespace storm { upperBound = pathFormula.getDiscreteTimeBound(); } - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, lowerBound, upperBound))); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeBoundedUntilProbabilities(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, lowerBound, upperBound, *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template 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->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeNextProbabilities(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), subResult.getTruthValuesVector(), *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template @@ -68,404 +66,29 @@ 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->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory))); - } - - template - std::vector SparseCtmcCslModelChecker::computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& exitRates, bool qualitative, double lowerBound, double upperBound) const { - // If the time bounds are [0, inf], we rather call untimed reachability. - storm::utility::ConstantsComparator comparator; - if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) { - return this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, qualitative, *this->linearEquationSolverFactory); - } - - // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0 - // or t' != inf. - - // Create the result vector. - std::vector result; - - // If we identify the states that have probability 0 of reaching the target states, we can exclude them from the - // further computations. - 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."); - - if (!statesWithProbabilityGreater0NonPsi.empty()) { - if (comparator.isZero(upperBound)) { - // In this case, the interval is of the form [0, 0]. - result = std::vector(this->getModel().getNumberOfStates(), storm::utility::zero()); - storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); - } else { - 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. - - // 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]); - } - uniformizationRate *= 1.02; - 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(), statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates); - - // Compute the vector that is to be added as a compensation for removing the absorbing states. - std::vector b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(statesWithProbabilityGreater0NonPsi, psiStates); - for (auto& element : b) { - element /= uniformizationRate; - } - - // Finally compute the transient probabilities. - std::vector values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero()); - std::vector subresult = this->computeTransientProbabilities(uniformizedMatrix, &b, upperBound, uniformizationRate, values, *this->linearEquationSolverFactory); - result = std::vector(this->getModel().getNumberOfStates(), storm::utility::zero()); - - storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0NonPsi, subresult); - storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); - } 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->linearEquationSolverFactory); - - // Determine the set of states that must be considered further. - storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; - std::vector subResult(relevantStates.getNumberOfSetBits()); - storm::utility::vector::selectVectorValues(subResult, relevantStates, result); - - ValueType uniformizationRate = 0; - for (auto const& state : relevantStates) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); - } - uniformizationRate *= 1.02; - 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(), relevantStates, uniformizationRate, exitRates); - - // Compute the transient probabilities. - subResult = this->computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, subResult, *this->linearEquationSolverFactory); - - // Fill in the correct values. - storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, relevantStates, subResult); - } else { - // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf. - - if (lowerBound != upperBound) { - // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'. - - // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. - ValueType uniformizationRate = storm::utility::zero(); - for (auto const& state : statesWithProbabilityGreater0NonPsi) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); - } - uniformizationRate *= 1.02; - 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(), statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates); - - // Compute the vector that is to be added as a compensation for removing the absorbing states. - std::vector b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(statesWithProbabilityGreater0NonPsi, psiStates); - for (auto& element : b) { - element /= uniformizationRate; - } - - // Start by computing the transient probabilities of reaching a psi state in time t' - t. - std::vector values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero()); - std::vector subresult = computeTransientProbabilities(uniformizedMatrix, &b, upperBound - lowerBound, uniformizationRate, values, *this->linearEquationSolverFactory); - - storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; - std::vector newSubresult = std::vector(relevantStates.getNumberOfSetBits()); - storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult); - storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one()); - - // 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. - uniformizationRate = storm::utility::zero(); - for (auto const& state : relevantStates) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); - } - uniformizationRate *= 1.02; - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - - // Finally, we compute the second set of transient probabilities. - uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), relevantStates, uniformizationRate, exitRates); - newSubresult = computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory); - - // Fill in the correct values. - result = std::vector(this->getModel().getNumberOfStates(), storm::utility::zero()); - storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, relevantStates, newSubresult); - } else { - // In this case, the interval is of the form [t, t] with t != 0, t != inf. - - std::vector newSubresult = std::vector(statesWithProbabilityGreater0.getNumberOfSetBits()); - storm::utility::vector::setVectorValues(newSubresult, psiStates % statesWithProbabilityGreater0, storm::utility::one()); - - // 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. - ValueType uniformizationRate = storm::utility::zero(); - for (auto const& state : statesWithProbabilityGreater0) { - uniformizationRate = std::max(uniformizationRate, exitRates[state]); - } - uniformizationRate *= 1.02; - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); - - // Finally, we compute the second set of transient probabilities. - storm::storage::SparseMatrix uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, uniformizationRate, exitRates); - newSubresult = computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory); - - // Fill in the correct values. - result = std::vector(this->getModel().getNumberOfStates(), storm::utility::zero()); - storm::utility::vector::setVectorValues(result, ~statesWithProbabilityGreater0, storm::utility::zero()); - storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, newSubresult); - } - } - } - } - - return result; - } - - template - storm::storage::SparseMatrix SparseCtmcCslModelChecker::computeUniformizedMatrix(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates, ValueType uniformizationRate, std::vector const& exitRates) { - STORM_LOG_DEBUG("Computing uniformized matrix using uniformization rate " << uniformizationRate << "."); - STORM_LOG_DEBUG("Keeping " << maybeStates.getNumberOfSetBits() << " rows."); - - // Create the submatrix that only contains the states with a positive probability (including the - // psi states) and reserve space for elements on the diagonal. - storm::storage::SparseMatrix uniformizedMatrix = transitionMatrix.getSubmatrix(false, maybeStates, maybeStates, true); - - // Now we need to perform the actual uniformization. That is, all entries need to be divided by - // the uniformization rate, and the diagonal needs to be set to the negative exit rate of the - // state plus the self-loop rate and then increased by one. - uint_fast64_t currentRow = 0; - for (auto const& state : maybeStates) { - for (auto& element : uniformizedMatrix.getRow(currentRow)) { - if (element.getColumn() == currentRow) { - element.setValue((element.getValue() - exitRates[state]) / uniformizationRate + storm::utility::one()); - } else { - element.setValue(element.getValue() / uniformizationRate); - } - } - ++currentRow; - } - - return uniformizedMatrix; - } - - template - template - std::vector SparseCtmcCslModelChecker::computeTransientProbabilities(storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - - ValueType lambda = timeBound * uniformizationRate; - - // If no time can pass, the current values are the result. - if (lambda == storm::utility::zero()) { - return values; - } - - // 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); - STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult)); - - // Scale the weights so they add up to one. - for (auto& element : std::get<3>(foxGlynnResult)) { - element /= std::get<2>(foxGlynnResult); - } - - // If the cumulative reward is to be computed, we need to adjust the weights. - if (computeCumulativeReward) { - ValueType sum = storm::utility::zero(); - - for (auto& element : std::get<3>(foxGlynnResult)) { - sum += element; - element = (1 - sum) / uniformizationRate; - } - } - - STORM_LOG_DEBUG("Starting iterations with " << uniformizedMatrix.getRowCount() << " x " << uniformizedMatrix.getColumnCount() << " matrix."); - - // Initialize result. - std::vector result; - uint_fast64_t startingIteration = std::get<0>(foxGlynnResult); - if (startingIteration == 0) { - result = values; - storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]); - std::function addAndScale = [&foxGlynnResult] (ValueType const& a, ValueType const& b) { return a + std::get<3>(foxGlynnResult)[0] * b; }; - if (addVector != nullptr) { - storm::utility::vector::applyPointwise(result, *addVector, result, addAndScale); - } - ++startingIteration; - } else { - if (computeCumulativeReward) { - result = std::vector(values.size()); - std::function scaleWithUniformizationRate = [&uniformizationRate] (ValueType const& a) -> ValueType { return a / uniformizationRate; }; - storm::utility::vector::applyPointwise(values, result, scaleWithUniformizationRate); - } else { - result = std::vector(values.size()); - } - } - std::vector multiplicationResult(result.size()); - - std::unique_ptr> solver = linearEquationSolverFactory.create(uniformizedMatrix); - - if (!computeCumulativeReward && std::get<0>(foxGlynnResult) > 1) { - // Perform the matrix-vector multiplications (without adding). - solver->performMatrixVectorMultiplication(values, addVector, std::get<0>(foxGlynnResult) - 1, &multiplicationResult); - } else if (computeCumulativeReward) { - std::function addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; }; - - // For the iterations below the left truncation point, we need to add and scale the result with the uniformization rate. - for (uint_fast64_t index = 1; index < startingIteration; ++index) { - solver->performMatrixVectorMultiplication(values, nullptr, 1, &multiplicationResult); - storm::utility::vector::applyPointwise(result, values, result, addAndScale); - } - } - - // 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 = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) { - solver->performMatrixVectorMultiplication(values, addVector, 1, &multiplicationResult); - - weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)]; - storm::utility::vector::applyPointwise(result, values, result, 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::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - return SparseDtmcPrctlModelChecker::computeUntilProbabilitiesHelper(transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolverFactory); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeUntilProbabilities(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getExitRateVector(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template std::unique_ptr SparseCtmcCslModelChecker::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeInstantaneousRewardsHelper(rewardPathFormula.getContinuousTimeBound()))); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeInstantaneousRewards(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getContinuousTimeBound(), *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } - - template - std::vector SparseCtmcCslModelChecker::computeInstantaneousRewardsHelper(double timeBound) const { - // Only compute the result if the model has a state-based reward this->getModel(). - STORM_LOG_THROW(this->getModel().hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); - - // Initialize result to state rewards of the this->getModel(). - std::vector result(this->getModel().getStateRewardVector()); - - // If the time-bound is not zero, we need to perform a transient analysis. - if (timeBound > 0) { - 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 uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector()); - result = this->computeTransientProbabilities(uniformizedMatrix, nullptr, timeBound, uniformizationRate, result, *this->linearEquationSolverFactory); - } - - return result; - } - template std::unique_ptr SparseCtmcCslModelChecker::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeCumulativeRewardsHelper(rewardPathFormula.getContinuousTimeBound()))); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeCumulativeRewards(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getContinuousTimeBound(), *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } - - template - std::vector SparseCtmcCslModelChecker::computeCumulativeRewardsHelper(double timeBound) const { - // Only compute the result if the model has a state-based reward this->getModel(). - STORM_LOG_THROW(this->getModel().hasStateRewards() || this->getModel().hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); - - // If the time bound is zero, the result is the constant zero vector. - if (timeBound == 0) { - return std::vector(this->getModel().getNumberOfStates(), storm::utility::zero()); - } - - // Otherwise, we need to perform some computations. - - // Start with the uniformization. - 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 uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector()); - - // Compute the total state reward vector. - std::vector totalRewardVector; - if (this->getModel().hasTransitionRewards()) { - totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix()); - if (this->getModel().hasStateRewards()) { - storm::utility::vector::addVectors(totalRewardVector, this->getModel().getStateRewardVector(), totalRewardVector); - } - } else { - totalRewardVector = std::vector(this->getModel().getStateRewardVector()); - } - - // Finally, compute the transient probabilities. - return this->computeTransientProbabilities(uniformizedMatrix, nullptr, timeBound, uniformizationRate, totalRewardVector, *this->linearEquationSolverFactory); - } - - 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; - } - - template - storm::storage::SparseMatrix SparseCtmcCslModelChecker::computeGeneratorMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates) { - storm::storage::SparseMatrix generatorMatrix(rateMatrix, true); - - // Place the negative exit rate on the diagonal. - for (uint_fast64_t row = 0; row < generatorMatrix.getRowCount(); ++row) { - for (auto& entry : generatorMatrix.getRow(row)) { - if (entry.getColumn() == row) { - entry.setValue(-exitRates[row]); - } - } - } - - return generatorMatrix; - } - + template std::unique_ptr SparseCtmcCslModelChecker::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(rewardPathFormula.getSubformula()); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - storm::storage::SparseMatrix probabilityMatrix = computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); - - boost::optional> modifiedStateRewardVector; - if (this->getModel().hasStateRewards()) { - modifiedStateRewardVector = std::vector(this->getModel().getStateRewardVector()); - - typename std::vector::const_iterator it2 = this->getModel().getExitRateVector().begin(); - for (typename std::vector::iterator it1 = modifiedStateRewardVector.get().begin(), ite1 = modifiedStateRewardVector.get().end(); it1 != ite1; ++it1, ++it2) { - *it1 /= *it2; - } - } - return std::unique_ptr(new ExplicitQuantitativeCheckResult(SparseDtmcPrctlModelChecker::computeReachabilityRewardsHelper(probabilityMatrix, modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative))); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getExitRateVector(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template @@ -473,234 +96,10 @@ namespace storm { std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - storm::storage::SparseMatrix probabilityMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(computeLongRunAverageHelper(probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory))); - } - - template - std::vector SparseCtmcCslModelChecker::computeLongRunAverageHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - // If there are no goal states, we avoid the computation and directly return zero. - 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. - if (psiStates.full()) { - return std::vector(numOfStates, storm::utility::one()); - } - - // Start by decomposing the DTMC into its BSCCs. - storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(transitionMatrix, storm::storage::BitVector(transitionMatrix.getRowCount(), true), false, true); - - STORM_LOG_DEBUG("Found " << bsccDecomposition.size() << " BSCCs."); - - // Get some data members for convenience. - ValueType one = storm::utility::one(); - ValueType zero = storm::utility::zero(); - - // Prepare the vector holding the LRA values for each of the BSCCs. - std::vector bsccLra(bsccDecomposition.size(), zero); - - // First we check which states are in BSCCs. - storm::storage::BitVector statesInBsccs(numOfStates); - storm::storage::BitVector firstStatesInBsccs(numOfStates); - - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - - // Gather information for later use. - bool first = true; - for (auto const& state : bscc) { - statesInBsccs.set(state); - if (first) { - firstStatesInBsccs.set(state); - } - first = false; - } - } - storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; - - STORM_LOG_DEBUG("Found " << statesInBsccs.getNumberOfSetBits() << " states in BSCCs."); - - // Prepare a vector holding the index within all states that are in BSCCs for every state. - std::vector indexInStatesInBsccs; - - // Prepare a vector that maps the index within the set of all states in BSCCs to the index of the containing BSCC. - std::vector stateToBsccIndexMap; - - if (!statesInBsccs.empty()) { - firstStatesInBsccs = firstStatesInBsccs % statesInBsccs; - - // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs. - storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); - - // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a - // multiplication from the right by transposing the system. - bsccEquationSystem = bsccEquationSystem.transpose(false, true); - - // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states. - uint_fast64_t lastIndex = 0; - uint_fast64_t currentNumberOfSetBits = 0; - indexInStatesInBsccs.reserve(transitionMatrix.getRowCount()); - for (auto index : statesInBsccs) { - while (lastIndex <= index) { - indexInStatesInBsccs.push_back(currentNumberOfSetBits); - ++lastIndex; - } - ++currentNumberOfSetBits; - } - - stateToBsccIndexMap.resize(statesInBsccs.getNumberOfSetBits()); - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - for (auto const& state : bscc) { - stateToBsccIndexMap[indexInStatesInBsccs[state]] = currentBsccIndex; - } - } - - // Now build the final equation system matrix, the initial guess and the right-hand side in one go. - std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); - storm::storage::SparseMatrixBuilder builder; - for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) { - - // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the - // values for states of this BSCC must sum to one. However, in order to have a non-zero value on the - // diagonal, we add the constraint of the BSCC that produces a 1 on the diagonal. - if (firstStatesInBsccs.get(row)) { - uint_fast64_t requiredBscc = stateToBsccIndexMap[row]; - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[requiredBscc]; - - for (auto const& state : bscc) { - builder.addNextValue(row, indexInStatesInBsccs[state], one); - } - - bsccEquationSystemRightSide[row] = one; - - } else { - // Otherwise, we copy the row, and subtract 1 from the diagonal. - for (auto& entry : bsccEquationSystem.getRow(row)) { - if (entry.getColumn() == row) { - builder.addNextValue(row, entry.getColumn(), entry.getValue() - one); - } else { - builder.addNextValue(row, entry.getColumn(), entry.getValue()); - } - } - } - - } - - // Create the initial guess for the LRAs. We take a uniform distribution over all states in a BSCC. - std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero); - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; - - for (auto const& state : bscc) { - bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size(); - } - } - - bsccEquationSystem = builder.build(); - - { - std::unique_ptr> solver = linearEquationSolverFactory.create(bsccEquationSystem); - solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); - } - - // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour. - if (exitRateVector != nullptr) { - std::vector bsccTotalValue(bsccDecomposition.size(), zero); - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { - bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]] += bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter]); - } - - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { - bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] = (bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]]; - } - } - // Calculate LRA Value for each BSCC from steady state distribution in BSCCs. - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; - - for (auto const& state : bscc) { - if (psiStates.get(state)) { - bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[state]]] += bsccEquationSystemSolution[indexInStatesInBsccs[state]]; - } - } - } - - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); - } - } else { - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; - - // At this point, all BSCCs are known to contain exactly one state, which is why we can set all values - // directly (based on whether or not the contained state is a psi state). - if (psiStates.get(*bscc.begin())) { - bsccLra[bsccIndex] = 1; - } - } - - for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); - } - } - - 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[indexInStatesInBsccs[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 (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; - - for (auto const& state : bscc) { - result[state] = bsccLra[bsccIndex]; - } - } - for (auto state : statesNotInBsccs) { - 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; - } - + storm::storage::SparseMatrix probabilityMatrix = storm::modelchecker::helper::SparseCtmcCslHelper::computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverage(probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); + } // Explicitly instantiate the model checker. template class SparseCtmcCslModelChecker>; diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h index d5fadf772..7cc6289bc 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h @@ -2,20 +2,14 @@ #define STORM_MODELCHECKER_SPARSECTMCCSLMODELCHECKER_H_ #include "src/modelchecker/propositional/SparsePropositionalModelChecker.h" + #include "src/models/sparse/Ctmc.h" -#include "src/storage/dd/DdType.h" + #include "src/utility/solver.h" -#include "src/solver/LinearEquationSolver.h" namespace storm { namespace modelchecker { - template - class HybridCtmcCslModelChecker; - - template - class SparseDtmcPrctlModelChecker; - template class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker { public: @@ -36,57 +30,6 @@ namespace storm { virtual std::unique_ptr computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative = false, boost::optional const& optimalityType = boost::optional()) override; private: - std::vector computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& exitRates, bool qualitative, double lowerBound, double upperBound) const; - static std::vector computeUntilProbabilitiesHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - std::vector computeInstantaneousRewardsHelper(RewardModelType const& rewardModel, double timeBound) const; - std::vector computeCumulativeRewardsHelper(RewardModelType const& rewardModel, double timeBound) const; - static std::vector computeLongRunAverageHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - /*! - * Computes the matrix representing the transitions of the uniformized CTMC. - * - * @param transitionMatrix The matrix to uniformize. - * @param maybeStates The states that need to be considered. - * @param uniformizationRate The rate to be used for uniformization. - * @param exitRates The exit rates of all states. - * @return The uniformized matrix. - */ - static storm::storage::SparseMatrix computeUniformizedMatrix(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates, ValueType uniformizationRate, std::vector const& exitRates); - - /*! - * Computes the transient probabilities for lambda time steps. - * - * @param uniformizedMatrix The uniformized transition matrix. - * @param addVector A vector that is added in each step as a possible compensation for removing absorbing states - * with a non-zero initial value. If this is not supposed to be used, it can be set to nullptr. - * @param timeBound The time bound to use. - * @param uniformizationRate The used uniformization rate. - * @param values A vector mapping each state to an initial probability. - * @param linearEquationSolverFactory The factory to use when instantiating new linear equation solvers. - * @param useMixedPoissonProbabilities If set to true, instead of taking the poisson probabilities, mixed - * poisson probabilities are used. - * @return The vector of transient probabilities. - */ - template - static std::vector computeTransientProbabilities(storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - - /*! - * Converts the given rate-matrix into a time-abstract probability matrix. - * - * @param rateMatrix The rate matrix. - * @param exitRates The exit rate vector. - * @return The †ransition matrix of the embedded DTMC. - */ - static storm::storage::SparseMatrix computeProbabilityMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates); - - /*! - * Converts the given rate-matrix into the generator matrix. - * - * @param rateMatrix The rate matrix. - * @param exitRates The exit rate vector. - * @return The generator matrix. - */ - static storm::storage::SparseMatrix computeGeneratorMatrix(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> linearEquationSolverFactory; }; diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp index e3149901f..99c92ee85 100644 --- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp @@ -570,6 +570,6 @@ namespace storm { return result; } - template class SparseMarkovAutomatonCslModelChecker; + template class SparseMarkovAutomatonCslModelChecker>; } } \ No newline at end of file diff --git a/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp b/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp new file mode 100644 index 000000000..e69de29bb diff --git a/src/modelchecker/csl/helper/HybridCtmcCslHelper.h b/src/modelchecker/csl/helper/HybridCtmcCslHelper.h new file mode 100644 index 000000000..e69de29bb diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index e69de29bb..013df87b4 100644 --- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -0,0 +1,648 @@ +#include "src/modelchecker/csl/helper/SparseCtmcCslHelper.h" + +#include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h" + +#include "src/storage/StronglyConnectedComponentDecomposition.h" + +#include "src/utility/macros.h" +#include "src/utility/vector.h" +#include "src/utility/graph.h" + +#include "src/utility/numerical.h" + +#include "src/exceptions/InvalidStateException.h" +#include "src/exceptions/InvalidPropertyException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + template + std::vector SparseCtmcCslHelper::computeBoundedUntilProbabilities(storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& exitRates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + + uint_fast64_t numberOfStates = rateMatrix.getRowCount(); + + // If the time bounds are [0, inf], we rather call untimed reachability. + storm::utility::ConstantsComparator comparator; + if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) { + return computeUntilProbabilities(rateMatrix, backwardTransitions, exitRates, phiStates, psiStates, qualitative, linearEquationSolverFactory); + } + + // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0 + // or t' != inf. + + // Create the result vector. + std::vector result; + + // 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(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."); + + if (!statesWithProbabilityGreater0NonPsi.empty()) { + if (comparator.isZero(upperBound)) { + // In this case, the interval is of the form [0, 0]. + result = std::vector(numberOfStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); + } else { + 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. + + // 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]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Compute the uniformized matrix. + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates); + + // Compute the vector that is to be added as a compensation for removing the absorbing states. + std::vector b = rateMatrix.getConstrainedRowSumVector(statesWithProbabilityGreater0NonPsi, psiStates); + for (auto& element : b) { + element /= uniformizationRate; + } + + // Finally compute the transient probabilities. + std::vector values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero()); + std::vector subresult = computeTransientProbabilities(uniformizedMatrix, &b, upperBound, uniformizationRate, values, linearEquationSolverFactory); + result = std::vector(numberOfStates, storm::utility::zero()); + + storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0NonPsi, subresult); + storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); + } 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 = computeUntilProbabilities(rateMatrix, backwardTransitions, exitRates, phiStates, psiStates, qualitative, linearEquationSolverFactory); + + // Determine the set of states that must be considered further. + storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; + std::vector subResult(relevantStates.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(subResult, relevantStates, result); + + ValueType uniformizationRate = 0; + for (auto const& state : relevantStates) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Compute the uniformized matrix. + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, relevantStates, uniformizationRate, exitRates); + + // Compute the transient probabilities. + subResult = computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, subResult, linearEquationSolverFactory); + + // Fill in the correct values. + storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, relevantStates, subResult); + } else { + // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf. + + if (lowerBound != upperBound) { + // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'. + + // Find the maximal rate of all 'maybe' states to take it as the uniformization rate. + ValueType uniformizationRate = storm::utility::zero(); + for (auto const& state : statesWithProbabilityGreater0NonPsi) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Compute the (first) uniformized matrix. + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates); + + // Compute the vector that is to be added as a compensation for removing the absorbing states. + std::vector b = rateMatrix.getConstrainedRowSumVector(statesWithProbabilityGreater0NonPsi, psiStates); + for (auto& element : b) { + element /= uniformizationRate; + } + + // Start by computing the transient probabilities of reaching a psi state in time t' - t. + std::vector values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero()); + std::vector subresult = computeTransientProbabilities(uniformizedMatrix, &b, upperBound - lowerBound, uniformizationRate, values, linearEquationSolverFactory); + + storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates; + std::vector newSubresult = std::vector(relevantStates.getNumberOfSetBits()); + storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult); + storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one()); + + // 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. + uniformizationRate = storm::utility::zero(); + for (auto const& state : relevantStates) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Finally, we compute the second set of transient probabilities. + uniformizedMatrix = computeUniformizedMatrix(rateMatrix, relevantStates, uniformizationRate, exitRates); + newSubresult = computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, linearEquationSolverFactory); + + // Fill in the correct values. + result = std::vector(numberOfStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, relevantStates, newSubresult); + } else { + // In this case, the interval is of the form [t, t] with t != 0, t != inf. + + std::vector newSubresult = std::vector(statesWithProbabilityGreater0.getNumberOfSetBits()); + storm::utility::vector::setVectorValues(newSubresult, psiStates % statesWithProbabilityGreater0, storm::utility::one()); + + // 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. + ValueType uniformizationRate = storm::utility::zero(); + for (auto const& state : statesWithProbabilityGreater0) { + uniformizationRate = std::max(uniformizationRate, exitRates[state]); + } + uniformizationRate *= 1.02; + STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); + + // Finally, we compute the second set of transient probabilities. + storm::storage::SparseMatrix uniformizedMatrix = computeUniformizedMatrix(rateMatrix, statesWithProbabilityGreater0, uniformizationRate, exitRates); + newSubresult = computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, linearEquationSolverFactory); + + // Fill in the correct values. + result = std::vector(numberOfStates, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, ~statesWithProbabilityGreater0, storm::utility::zero()); + storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, newSubresult); + } + } + } + } + + return result; + } + + template + std::vector SparseCtmcCslHelper::computeUntilProbabilities(storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + return SparseDtmcPrctlHelper::computeUntilProbabilities(computeProbabilityMatrix(rateMatrix, exitRateVector), backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolverFactory); + } + + template + std::vector SparseCtmcCslHelper::computeNextProbabilities(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + return SparseDtmcPrctlHelper::computeNextProbabilities(computeProbabilityMatrix(rateMatrix, exitRateVector), nextStates, linearEquationSolverFactory); + } + + template + storm::storage::SparseMatrix SparseCtmcCslHelper::computeUniformizedMatrix(storm::storage::SparseMatrix const& rateMatrix, storm::storage::BitVector const& maybeStates, ValueType uniformizationRate, std::vector const& exitRates) { + STORM_LOG_DEBUG("Computing uniformized matrix using uniformization rate " << uniformizationRate << "."); + STORM_LOG_DEBUG("Keeping " << maybeStates.getNumberOfSetBits() << " rows."); + + // Create the submatrix that only contains the states with a positive probability (including the + // psi states) and reserve space for elements on the diagonal. + storm::storage::SparseMatrix uniformizedMatrix = rateMatrix.getSubmatrix(false, maybeStates, maybeStates, true); + + // Now we need to perform the actual uniformization. That is, all entries need to be divided by + // the uniformization rate, and the diagonal needs to be set to the negative exit rate of the + // state plus the self-loop rate and then increased by one. + uint_fast64_t currentRow = 0; + for (auto const& state : maybeStates) { + for (auto& element : uniformizedMatrix.getRow(currentRow)) { + if (element.getColumn() == currentRow) { + element.setValue((element.getValue() - exitRates[state]) / uniformizationRate + storm::utility::one()); + } else { + element.setValue(element.getValue() / uniformizationRate); + } + } + ++currentRow; + } + + return uniformizedMatrix; + } + + template + template + std::vector SparseCtmcCslHelper::computeTransientProbabilities(storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + + ValueType lambda = timeBound * uniformizationRate; + + // If no time can pass, the current values are the result. + if (lambda == storm::utility::zero()) { + return values; + } + + // 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); + STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult)); + + // Scale the weights so they add up to one. + for (auto& element : std::get<3>(foxGlynnResult)) { + element /= std::get<2>(foxGlynnResult); + } + + // If the cumulative reward is to be computed, we need to adjust the weights. + if (computeCumulativeReward) { + ValueType sum = storm::utility::zero(); + + for (auto& element : std::get<3>(foxGlynnResult)) { + sum += element; + element = (1 - sum) / uniformizationRate; + } + } + + STORM_LOG_DEBUG("Starting iterations with " << uniformizedMatrix.getRowCount() << " x " << uniformizedMatrix.getColumnCount() << " matrix."); + + // Initialize result. + std::vector result; + uint_fast64_t startingIteration = std::get<0>(foxGlynnResult); + if (startingIteration == 0) { + result = values; + storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]); + std::function addAndScale = [&foxGlynnResult] (ValueType const& a, ValueType const& b) { return a + std::get<3>(foxGlynnResult)[0] * b; }; + if (addVector != nullptr) { + storm::utility::vector::applyPointwise(result, *addVector, result, addAndScale); + } + ++startingIteration; + } else { + if (computeCumulativeReward) { + result = std::vector(values.size()); + std::function scaleWithUniformizationRate = [&uniformizationRate] (ValueType const& a) -> ValueType { return a / uniformizationRate; }; + storm::utility::vector::applyPointwise(values, result, scaleWithUniformizationRate); + } else { + result = std::vector(values.size()); + } + } + std::vector multiplicationResult(result.size()); + + std::unique_ptr> solver = linearEquationSolverFactory.create(uniformizedMatrix); + + if (!computeCumulativeReward && std::get<0>(foxGlynnResult) > 1) { + // Perform the matrix-vector multiplications (without adding). + solver->performMatrixVectorMultiplication(values, addVector, std::get<0>(foxGlynnResult) - 1, &multiplicationResult); + } else if (computeCumulativeReward) { + std::function addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; }; + + // For the iterations below the left truncation point, we need to add and scale the result with the uniformization rate. + for (uint_fast64_t index = 1; index < startingIteration; ++index) { + solver->performMatrixVectorMultiplication(values, nullptr, 1, &multiplicationResult); + storm::utility::vector::applyPointwise(result, values, result, addAndScale); + } + } + + // 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 = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) { + solver->performMatrixVectorMultiplication(values, addVector, 1, &multiplicationResult); + + weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)]; + storm::utility::vector::applyPointwise(result, values, result, addAndScale); + } + + return result; + } + + template + std::vector SparseCtmcCslHelper::computeInstantaneousRewards(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // Only compute the result if the model has a state-based reward this->getModel(). + STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + uint_fast64_t numberOfStates = rateMatrix.getRowCount(); + + // Initialize result to state rewards of the this->getModel(). + std::vector result(rewardModel.getStateRewardVector()); + + // If the time-bound is not zero, we need to perform a transient analysis. + if (timeBound > 0) { + ValueType uniformizationRate = 0; + for (auto const& rate : exitRateVector) { + 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 uniformizedMatrix = computeUniformizedMatrix(rateMatrix, storm::storage::BitVector(numberOfStates, true), uniformizationRate, exitRateVector); + result = computeTransientProbabilities(uniformizedMatrix, nullptr, timeBound, uniformizationRate, result, linearEquationSolverFactory); + } + + return result; + } + + template + std::vector SparseCtmcCslHelper::computeCumulativeRewards(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // Only compute the result if the model has a state-based reward this->getModel(). + STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + uint_fast64_t numberOfStates = rateMatrix.getRowCount(); + + // If the time bound is zero, the result is the constant zero vector. + if (timeBound == 0) { + return std::vector(numberOfStates, storm::utility::zero()); + } + + // Otherwise, we need to perform some computations. + + // Start with the uniformization. + ValueType uniformizationRate = 0; + for (auto const& rate : exitRateVector) { + 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 uniformizedMatrix = computeUniformizedMatrix(rateMatrix, storm::storage::BitVector(numberOfStates, true), uniformizationRate, exitRateVector); + + // Compute the total state reward vector. + std::vector totalRewardVector = rewardModel.getTotalRewardVector(uniformizedMatrix); + + // Finally, compute the transient probabilities. + return computeTransientProbabilities(uniformizedMatrix, nullptr, timeBound, uniformizationRate, totalRewardVector, linearEquationSolverFactory); + } + + template + std::vector SparseCtmcCslHelper::computeReachabilityRewards(storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + storm::storage::SparseMatrix probabilityMatrix = computeProbabilityMatrix(rateMatrix, exitRateVector); + + std::vector totalRewardVector; + if (rewardModel.hasStateRewards()) { + totalRewardVector = rewardModel.getStateRewardVector(); + typename std::vector::const_iterator it2 = exitRateVector.begin(); + for (typename std::vector::iterator it1 = totalRewardVector.begin(), ite1 = totalRewardVector.end(); it1 != ite1; ++it1, ++it2) { + *it1 /= *it2; + } + if (rewardModel.hasStateActionRewards()) { + storm::utility::vector::addVectors(totalRewardVector, rewardModel.getStateActionRewardVector(), totalRewardVector); + } + if (rewardModel.hasTransitionRewards()) { + storm::utility::vector::addVectors(totalRewardVector, probabilityMatrix.getPointwiseProductRowSumVector(rewardModel.getTransitionRewardMatrix()), totalRewardVector); + } + } else if (rewardModel.hasTransitionRewards()) { + totalRewardVector = probabilityMatrix.getPointwiseProductRowSumVector(rewardModel.getTransitionRewardMatrix()); + if (rewardModel.hasStateActionRewards()) { + storm::utility::vector::addVectors(totalRewardVector, rewardModel.getStateActionRewardVector(), totalRewardVector); + } + } else { + totalRewardVector = rewardModel.getStateActionRewardVector(); + } + + return storm::modelchecker::helper::SparseDtmcPrctlHelper::computeReachabilityRewards(probabilityMatrix, backwardTransitions, totalRewardVector, targetStates, qualitative, linearEquationSolverFactory); + } + + template + storm::storage::SparseMatrix SparseCtmcCslHelper::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; + } + + template + storm::storage::SparseMatrix SparseCtmcCslHelper::computeGeneratorMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates) { + storm::storage::SparseMatrix generatorMatrix(rateMatrix, true); + + // Place the negative exit rate on the diagonal. + for (uint_fast64_t row = 0; row < generatorMatrix.getRowCount(); ++row) { + for (auto& entry : generatorMatrix.getRow(row)) { + if (entry.getColumn() == row) { + entry.setValue(-exitRates[row]); + } + } + } + + return generatorMatrix; + } + + template + std::vector SparseCtmcCslHelper::computeLongRunAverage(storm::storage::SparseMatrix const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // If there are no goal states, we avoid the computation and directly return zero. + uint_fast64_t numberOfStates = probabilityMatrix.getRowCount(); + if (psiStates.empty()) { + return std::vector(numberOfStates, storm::utility::zero()); + } + + // Likewise, if all bits are set, we can avoid the computation. + if (psiStates.full()) { + return std::vector(numberOfStates, storm::utility::one()); + } + + // Start by decomposing the DTMC into its BSCCs. + storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(probabilityMatrix, storm::storage::BitVector(probabilityMatrix.getRowCount(), true), false, true); + + STORM_LOG_DEBUG("Found " << bsccDecomposition.size() << " BSCCs."); + + // Get some data members for convenience. + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); + + // Prepare the vector holding the LRA values for each of the BSCCs. + std::vector bsccLra(bsccDecomposition.size(), zero); + + // First we check which states are in BSCCs. + storm::storage::BitVector statesInBsccs(numberOfStates); + storm::storage::BitVector firstStatesInBsccs(numberOfStates); + + for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; + + // Gather information for later use. + bool first = true; + for (auto const& state : bscc) { + statesInBsccs.set(state); + if (first) { + firstStatesInBsccs.set(state); + } + first = false; + } + } + storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; + + STORM_LOG_DEBUG("Found " << statesInBsccs.getNumberOfSetBits() << " states in BSCCs."); + + // Prepare a vector holding the index within all states that are in BSCCs for every state. + std::vector indexInStatesInBsccs; + + // Prepare a vector that maps the index within the set of all states in BSCCs to the index of the containing BSCC. + std::vector stateToBsccIndexMap; + + if (!statesInBsccs.empty()) { + firstStatesInBsccs = firstStatesInBsccs % statesInBsccs; + + // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs. + storm::storage::SparseMatrix bsccEquationSystem = probabilityMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); + + // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a + // multiplication from the right by transposing the system. + bsccEquationSystem = bsccEquationSystem.transpose(false, true); + + // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states. + uint_fast64_t lastIndex = 0; + uint_fast64_t currentNumberOfSetBits = 0; + indexInStatesInBsccs.reserve(probabilityMatrix.getRowCount()); + for (auto index : statesInBsccs) { + while (lastIndex <= index) { + indexInStatesInBsccs.push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; + } + + stateToBsccIndexMap.resize(statesInBsccs.getNumberOfSetBits()); + for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; + for (auto const& state : bscc) { + stateToBsccIndexMap[indexInStatesInBsccs[state]] = currentBsccIndex; + } + } + + // Now build the final equation system matrix, the initial guess and the right-hand side in one go. + std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); + storm::storage::SparseMatrixBuilder builder; + for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) { + + // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the + // values for states of this BSCC must sum to one. However, in order to have a non-zero value on the + // diagonal, we add the constraint of the BSCC that produces a 1 on the diagonal. + if (firstStatesInBsccs.get(row)) { + uint_fast64_t requiredBscc = stateToBsccIndexMap[row]; + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[requiredBscc]; + + for (auto const& state : bscc) { + builder.addNextValue(row, indexInStatesInBsccs[state], one); + } + + bsccEquationSystemRightSide[row] = one; + + } else { + // Otherwise, we copy the row, and subtract 1 from the diagonal. + for (auto& entry : bsccEquationSystem.getRow(row)) { + if (entry.getColumn() == row) { + builder.addNextValue(row, entry.getColumn(), entry.getValue() - one); + } else { + builder.addNextValue(row, entry.getColumn(), entry.getValue()); + } + } + } + + } + + // Create the initial guess for the LRAs. We take a uniform distribution over all states in a BSCC. + std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero); + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + for (auto const& state : bscc) { + bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size(); + } + } + + bsccEquationSystem = builder.build(); + + { + std::unique_ptr> solver = linearEquationSolverFactory.create(bsccEquationSystem); + solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); + } + + // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour. + if (exitRateVector != nullptr) { + std::vector bsccTotalValue(bsccDecomposition.size(), zero); + for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { + bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]] += bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter]); + } + + for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++stateIter) { + bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] = (bsccEquationSystemSolution[indexInStatesInBsccs[*stateIter]] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[indexInStatesInBsccs[*stateIter]]]; + } + } + // Calculate LRA Value for each BSCC from steady state distribution in BSCCs. + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + for (auto const& state : bscc) { + if (psiStates.get(state)) { + bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[state]]] += bsccEquationSystemSolution[indexInStatesInBsccs[state]]; + } + } + } + + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); + } + } else { + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + // At this point, all BSCCs are known to contain exactly one state, which is why we can set all values + // directly (based on whether or not the contained state is a psi state). + if (psiStates.get(*bscc.begin())) { + bsccLra[bsccIndex] = 1; + } + } + + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + STORM_LOG_DEBUG("Found LRA " << bsccLra[bsccIndex] << " for BSCC " << bsccIndex << "."); + } + } + + 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 : probabilityMatrix.getRow(state)) { + if (statesInBsccs.get(entry.getColumn())) { + reward += entry.getValue() * bsccLra[stateToBsccIndexMap[indexInStatesInBsccs[entry.getColumn()]]]; + } + } + rewardRightSide.push_back(reward); + } + + storm::storage::SparseMatrix rewardEquationSystemMatrix = probabilityMatrix.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(numberOfStates); + auto rewardSolutionIter = rewardSolution.begin(); + + for (uint_fast64_t bsccIndex = 0; bsccIndex < bsccDecomposition.size(); ++bsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[bsccIndex]; + + for (auto const& state : bscc) { + result[state] = bsccLra[bsccIndex]; + } + } + for (auto state : statesNotInBsccs) { + 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; + } + + template class SparseCtmcCslHelper; + } + } +} \ No newline at end of file diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.h b/src/modelchecker/csl/helper/SparseCtmcCslHelper.h index d2edcbd18..55a143c8a 100644 --- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.h +++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.h @@ -1,11 +1,8 @@ #ifndef STORM_MODELCHECKER_SPARSE_CTMC_CSL_MODELCHECKER_HELPER_H_ #define STORM_MODELCHECKER_SPARSE_CTMC_CSL_MODELCHECKER_HELPER_H_ -#include - #include "src/models/sparse/StandardRewardModel.h" -#include "src/storage/SparseMatrix.h" #include "src/storage/BitVector.h" #include "src/utility/solver.h" @@ -16,8 +13,66 @@ namespace storm { template > class SparseCtmcCslHelper { public: + static std::vector computeBoundedUntilProbabilities(storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector const& exitRates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - } + static std::vector computeNextProbabilities(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::vector computeUntilProbabilities(storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::vector computeInstantaneousRewards(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::vector computeCumulativeRewards(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::vector computeReachabilityRewards(storm::storage::SparseMatrix const& rateMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::vector computeLongRunAverage(storm::storage::SparseMatrix const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + /*! + * Computes the matrix representing the transitions of the uniformized CTMC. + * + * @param transitionMatrix The matrix to uniformize. + * @param maybeStates The states that need to be considered. + * @param uniformizationRate The rate to be used for uniformization. + * @param exitRates The exit rates of all states. + * @return The uniformized matrix. + */ + static storm::storage::SparseMatrix computeUniformizedMatrix(storm::storage::SparseMatrix const& rateMatrix, storm::storage::BitVector const& maybeStates, ValueType uniformizationRate, std::vector const& exitRates); + + /*! + * Computes the transient probabilities for lambda time steps. + * + * @param uniformizedMatrix The uniformized transition matrix. + * @param addVector A vector that is added in each step as a possible compensation for removing absorbing states + * with a non-zero initial value. If this is not supposed to be used, it can be set to nullptr. + * @param timeBound The time bound to use. + * @param uniformizationRate The used uniformization rate. + * @param values A vector mapping each state to an initial probability. + * @param linearEquationSolverFactory The factory to use when instantiating new linear equation solvers. + * @param useMixedPoissonProbabilities If set to true, instead of taking the poisson probabilities, mixed + * poisson probabilities are used. + * @return The vector of transient probabilities. + */ + template + static std::vector computeTransientProbabilities(storm::storage::SparseMatrix const& uniformizedMatrix, std::vector const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector values, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + /*! + * Converts the given rate-matrix into a time-abstract probability matrix. + * + * @param rateMatrix The rate matrix. + * @param exitRates The exit rate vector. + * @return The †ransition matrix of the embedded DTMC. + */ + static storm::storage::SparseMatrix computeProbabilityMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates); + + /*! + * Converts the given rate-matrix into the generator matrix. + * + * @param rateMatrix The rate matrix. + * @param exitRates The exit rate vector. + * @return The generator matrix. + */ + static storm::storage::SparseMatrix computeGeneratorMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates); + }; } } } diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp index 3b61ded53..29f39f50b 100644 --- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp @@ -1,5 +1,7 @@ #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h" -#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" + +#include "src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h" +#include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h" #include "src/storage/dd/CuddOdd.h" @@ -29,88 +31,23 @@ namespace storm { bool HybridDtmcPrctlModelChecker::canHandle(storm::logic::Formula const& formula) const { return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula(); } - - template - std::unique_ptr HybridDtmcPrctlModelChecker::computeUntilProbabilitiesHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - // We need to identify the states which have to be taken out of the matrix, i.e. all states that have - // probability 0 and 1 of satisfying the until-formula. - std::pair, storm::dd::Bdd> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates); - storm::dd::Bdd maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates(); - - // Perform some logging. - STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states."); - STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states."); - STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states."); - - // Check whether we need to compute exact probabilities for some states. - if (qualitative) { - // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1. - return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5))); - } else { - // If there are maybe states, we need to solve an equation system. - if (!maybeStates.isZero()) { - // Create the ODD for the translation between symbolic and explicit storage. - storm::dd::Odd odd(maybeStates); - - // Create the matrix and the vector for the equation system. - storm::dd::Add maybeStatesAdd = maybeStates.toAdd(); - - // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting - // non-maybe states in the matrix. - storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; - - // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all - // maybe states. - storm::dd::Add prob1StatesAsColumn = statesWithProbability01.second.toAdd(); - prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs()); - storm::dd::Add subvector = submatrix * prob1StatesAsColumn; - subvector = subvector.sumAbstract(model.getColumnVariables()); - - // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed - // for solving the equation system (i.e. compute (I-A)). - submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); - submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix; - - // Create the solution vector. - std::vector x(maybeStates.getNonZeroCount(), ValueType(0.5)); - - // Translate the symbolic matrix/vector to their explicit representations and solve the equation system. - storm::storage::SparseMatrix explicitSubmatrix = submatrix.toMatrix(odd, odd); - std::vector b = subvector.template toVector(odd); - - std::unique_ptr> solver = linearEquationSolverFactory.create(explicitSubmatrix); - solver->solveEquationSystem(x, b); - - // Return a hybrid check result that stores the numerical values explicitly. - return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, statesWithProbability01.second.toAdd(), maybeStates, odd, x)); - } else { - return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), statesWithProbability01.second.toAdd())); - } - } - } - + template std::unique_ptr HybridDtmcPrctlModelChecker::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr leftResultPointer = this->check(pathFormula.getLeftSubformula()); std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); SymbolicQualitativeCheckResult const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult(); SymbolicQualitativeCheckResult const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult(); - return this->computeUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); + return storm::modelchecker::helper::HybridDtmcPrctlHelper::computeUntilProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); } template std::unique_ptr HybridDtmcPrctlModelChecker::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(pathFormula.getSubformula()); SymbolicQualitativeCheckResult const& subResult = subResultPointer->asSymbolicQualitativeCheckResult(); - return std::unique_ptr(new SymbolicQuantitativeCheckResult(this->getModel().getReachableStates(), this->computeNextProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector()))); - } - - template - storm::dd::Add HybridDtmcPrctlModelChecker::computeNextProbabilitiesHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates) { - storm::dd::Add result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd(); - return result.sumAbstract(model.getColumnVariables()); + return storm::modelchecker::helper::HybridDtmcPrctlHelper::computeNextProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector()); } - + template std::unique_ptr HybridDtmcPrctlModelChecker::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); @@ -118,185 +55,28 @@ namespace storm { std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); SymbolicQualitativeCheckResult const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult(); SymbolicQualitativeCheckResult const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult(); - return this->computeBoundedUntilProbabilitiesHelper(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); + return storm::modelchecker::helper::HybridDtmcPrctlHelper::computeBoundedUntilProbabilities(this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); } - - template - std::unique_ptr HybridDtmcPrctlModelChecker::computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - // We need to identify the states which have to be taken out of the matrix, i.e. all states that have - // probability 0 or 1 of satisfying the until-formula. - storm::dd::Bdd statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound); - storm::dd::Bdd maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates(); - - // If there are maybe states, we need to perform matrix-vector multiplications. - if (!maybeStates.isZero()) { - // Create the ODD for the translation between symbolic and explicit storage. - storm::dd::Odd odd(maybeStates); - - // Create the matrix and the vector for the equation system. - storm::dd::Add maybeStatesAdd = maybeStates.toAdd(); - - // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting - // non-maybe states in the matrix. - storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; - - // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all - // maybe states. - storm::dd::Add prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs()); - storm::dd::Add subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables()); - - // Finally cut away all columns targeting non-maybe states. - submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); - // Create the solution vector. - std::vector x(maybeStates.getNonZeroCount(), storm::utility::zero()); - - // Translate the symbolic matrix/vector to their explicit representations. - storm::storage::SparseMatrix explicitSubmatrix = submatrix.toMatrix(odd, odd); - std::vector b = subvector.template toVector(odd); - - std::unique_ptr> solver = linearEquationSolverFactory.create(explicitSubmatrix); - solver->performMatrixVectorMultiplication(x, &b, stepBound); - - // Return a hybrid check result that stores the numerical values explicitly. - return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.toAdd(), maybeStates, odd, x)); - } else { - return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), psiStates.toAdd())); - } - } - template - std::unique_ptr HybridDtmcPrctlModelChecker::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional const& optimalityType) { + std::unique_ptr HybridDtmcPrctlModelChecker::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); - return this->computeCumulativeRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); - } - - template - std::unique_ptr HybridDtmcPrctlModelChecker::computeCumulativeRewardsHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - // Only compute the result if the model has at least one reward this->getModel(). - STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); - - // Compute the reward vector to add in each step based on the available reward models. - storm::dd::Add totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero(); - if (model.hasTransitionRewards()) { - totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables()); - } - - // Create the ODD for the translation between symbolic and explicit storage. - storm::dd::Odd odd(model.getReachableStates()); - - // Create the solution vector. - std::vector x(model.getNumberOfStates(), storm::utility::zero()); - - // Translate the symbolic matrix/vector to their explicit representations. - storm::storage::SparseMatrix explicitMatrix = transitionMatrix.toMatrix(odd, odd); - std::vector b = totalRewardVector.template toVector(odd); - - // Perform the matrix-vector multiplication. - std::unique_ptr> solver = linearEquationSolverFactory.create(explicitMatrix); - solver->performMatrixVectorMultiplication(x, &b, stepBound); - - // Return a hybrid check result that stores the numerical values explicitly. - return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x)); + return storm::modelchecker::helper::HybridDtmcPrctlHelper::computeCumulativeRewards(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); } - + template - std::unique_ptr HybridDtmcPrctlModelChecker::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional const& optimalityType) { + std::unique_ptr HybridDtmcPrctlModelChecker::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); - return this->computeInstantaneousRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); - } - - template - std::unique_ptr HybridDtmcPrctlModelChecker::computeInstantaneousRewardsHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - // Only compute the result if the model has at least one reward this->getModel(). - STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); - - // Create the ODD for the translation between symbolic and explicit storage. - storm::dd::Odd odd(model.getReachableStates()); - - // Create the solution vector (and initialize it to the state rewards of the model). - std::vector x = model.getStateRewardVector().template toVector(odd); - - // Translate the symbolic matrix to its explicit representations. - storm::storage::SparseMatrix explicitMatrix = transitionMatrix.toMatrix(odd, odd); - - // Perform the matrix-vector multiplication. - std::unique_ptr> solver = linearEquationSolverFactory.create(explicitMatrix); - solver->performMatrixVectorMultiplication(x, nullptr, stepBound); - - // Return a hybrid check result that stores the numerical values explicitly. - return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x)); + return storm::modelchecker::helper::HybridDtmcPrctlHelper::computeInstantaneousRewards(this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); } template - std::unique_ptr HybridDtmcPrctlModelChecker::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional const& optimalityType) { + std::unique_ptr HybridDtmcPrctlModelChecker::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(rewardPathFormula.getSubformula()); SymbolicQualitativeCheckResult const& subResult = subResultPointer->asSymbolicQualitativeCheckResult(); - return this->computeReachabilityRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative); + return storm::modelchecker::helper::HybridDtmcPrctlHelper::computeReachabilityRewards(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); } - - template - std::unique_ptr HybridDtmcPrctlModelChecker::computeReachabilityRewardsHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::dd::Bdd const& targetStates, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, bool qualitative) { - - // Only compute the result if there is at least one reward model. - STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); - - // Determine which states have a reward of infinity by definition. - storm::dd::Bdd infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates); - infinityStates = !infinityStates && model.getReachableStates(); - storm::dd::Bdd maybeStates = (!targetStates && !infinityStates) && model.getReachableStates(); - STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states."); - STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states."); - STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states."); - - // Check whether we need to compute exact rewards for some states. - if (qualitative) { - // Set the values for all maybe-states to 1 to indicate that their reward values - // are neither 0 nor infinity. - return std::unique_ptr(new SymbolicQuantitativeCheckResult(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one()))); - } else { - // If there are maybe states, we need to solve an equation system. - if (!maybeStates.isZero()) { - // Create the ODD for the translation between symbolic and explicit storage. - storm::dd::Odd odd(maybeStates); - - // Create the matrix and the vector for the equation system. - storm::dd::Add maybeStatesAdd = maybeStates.toAdd(); - - // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting - // non-maybe states in the matrix. - storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; - - // Then compute the state reward vector to use in the computation. - storm::dd::Add subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero(); - if (transitionRewardMatrix) { - subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables()); - } - - // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed - // for solving the equation system (i.e. compute (I-A)). - submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); - submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix; - - // Create the solution vector. - std::vector x(maybeStates.getNonZeroCount(), ValueType(0.5)); - - // Translate the symbolic matrix/vector to their explicit representations. - storm::storage::SparseMatrix explicitSubmatrix = submatrix.toMatrix(odd, odd); - std::vector b = subvector.template toVector(odd); - - // Now solve the resulting equation system. - std::unique_ptr> solver = linearEquationSolverFactory.create(explicitSubmatrix); - solver->solveEquationSystem(x, b); - - // Return a hybrid check result that stores the numerical values explicitly. - return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity()), maybeStates, odd, x)); - } else { - return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity()))); - } - } - } - + template storm::models::symbolic::Dtmc const& HybridDtmcPrctlModelChecker::getModel() const { return this->template getModelAs>(); @@ -312,7 +92,7 @@ namespace storm { storm::storage::SparseMatrix explicitProbabilityMatrix = this->getModel().getTransitionMatrix().toMatrix(odd, odd); - std::vector result = SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), qualitative, *this->linearEquationSolverFactory); + std::vector result = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeLongRunAverage(explicitProbabilityMatrix, subResult.getTruthValuesVector().toVector(odd), qualitative, *this->linearEquationSolverFactory); return std::unique_ptr(new HybridQuantitativeCheckResult(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), std::move(odd), std::move(result))); } diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h index 016e2898a..36e52eebb 100644 --- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h @@ -2,19 +2,17 @@ #define STORM_MODELCHECKER_HYBRIDDTMCPRCTLMODELCHECKER_H_ #include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h" + #include "src/models/symbolic/Dtmc.h" + #include "src/utility/solver.h" namespace storm { namespace modelchecker { - template - class HybridCtmcCslModelChecker; template class HybridDtmcPrctlModelChecker : public SymbolicPropositionalModelChecker { public: - friend class HybridCtmcCslModelChecker; - explicit HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc const& model); explicit HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc const& model, std::unique_ptr>&& linearEquationSolverFactory); @@ -32,14 +30,6 @@ namespace storm { storm::models::symbolic::Dtmc const& getModel() const override; private: - // The methods that perform the actual checking. - static std::unique_ptr computeBoundedUntilProbabilitiesHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - static storm::dd::Add computeNextProbabilitiesHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates); - static std::unique_ptr computeUntilProbabilitiesHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - static std::unique_ptr computeCumulativeRewardsHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - static std::unique_ptr computeInstantaneousRewardsHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - static std::unique_ptr computeReachabilityRewardsHelper(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::dd::Bdd const& targetStates, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, bool qualitative); - // An object that is used for retrieving linear equation solvers. std::unique_ptr> linearEquationSolverFactory; }; diff --git a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp index a8af6665c..e76551809 100644 --- a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp @@ -1,13 +1,14 @@ #include "src/modelchecker/prctl/HybridMdpPrctlModelChecker.h" -#include "src/storage/dd/CuddOdd.h" +#include "src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h" -#include "src/utility/macros.h" -#include "src/utility/graph.h" +#include "src/storage/dd/CuddOdd.h" #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h" #include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h" -#include "src/modelchecker/results/HybridQuantitativeCheckResult.h" + +#include "src/utility/macros.h" +#include "src/utility/graph.h" #include "src/exceptions/InvalidStateException.h" #include "src/exceptions/InvalidPropertyException.h" @@ -36,7 +37,7 @@ namespace storm { std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); SymbolicQualitativeCheckResult const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult(); SymbolicQualitativeCheckResult const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult(); - return this->computeUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); + return storm::modelchecker::helper::HybridMdpPrctlHelper::computeUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); } template @@ -44,7 +45,7 @@ namespace storm { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); std::unique_ptr subResultPointer = this->check(pathFormula.getSubformula()); SymbolicQualitativeCheckResult const& subResult = subResultPointer->asSymbolicQualitativeCheckResult(); - return std::unique_ptr(new SymbolicQuantitativeCheckResult(this->getModel().getReachableStates(), this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector()))); + return storm::modelchecker::helper::HybridMdpPrctlHelper::computeNextProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector()); } template @@ -55,21 +56,21 @@ namespace storm { std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); SymbolicQualitativeCheckResult const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult(); SymbolicQualitativeCheckResult const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult(); - return this->computeBoundedUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); + return storm::modelchecker::helper::HybridMdpPrctlHelper::computeBoundedUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); } template std::unique_ptr HybridMdpPrctlModelChecker::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); - return this->computeCumulativeRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); + return storm::modelchecker::helper::HybridMdpPrctlHelper::computeCumulativeRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); } template std::unique_ptr HybridMdpPrctlModelChecker::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); - return this->computeInstantaneousRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); + return storm::modelchecker::helper::HybridMdpPrctlHelper::computeInstantaneousRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory); } template @@ -77,7 +78,7 @@ namespace storm { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); std::unique_ptr subResultPointer = this->check(rewardPathFormula.getSubformula()); SymbolicQualitativeCheckResult const& subResult = subResultPointer->asSymbolicQualitativeCheckResult(); - return this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative); + return storm::modelchecker::helper::HybridMdpPrctlHelper::computeReachabilityRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory); } template diff --git a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h index ee22ff766..bf5a7a48a 100644 --- a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h +++ b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h @@ -2,7 +2,9 @@ #define STORM_MODELCHECKER_HYBRIDMDPPRCTLMODELCHECKER_H_ #include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h" + #include "src/models/symbolic/Mdp.h" + #include "src/utility/solver.h" namespace storm { diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index ec9743ffe..62eb91acc 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -9,6 +9,7 @@ #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h" +#include "src/modelchecker/csl/helper/SparseCtmcCslHelper.h" namespace storm { namespace modelchecker { @@ -83,7 +84,7 @@ namespace storm { std::unique_ptr SparseDtmcPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverage(this->getModel().getTransitionMatrix(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory); + std::vector numericResult = storm::modelchecker::helper::SparseCtmcCslHelper::computeLongRunAverage(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), nullptr, qualitative, *linearEquationSolverFactory); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index cc574d289..1d9b7efba 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -9,17 +9,12 @@ namespace storm { namespace modelchecker { - template - class HybridDtmcPrctlModelChecker; - template class SparseDtmcPrctlModelChecker : public SparsePropositionalModelChecker { public: typedef typename SparseDtmcModelType::ValueType ValueType; typedef typename SparseDtmcModelType::RewardModelType RewardModelType; - - friend class HybridDtmcPrctlModelChecker; - + explicit SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model); explicit SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model, std::unique_ptr>&& linearEquationSolverFactory); diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index adceed6f8..c51721160 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -1,8 +1,5 @@ #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" -#include -#include - #include "src/utility/constants.h" #include "src/utility/macros.h" #include "src/utility/vector.h" @@ -11,6 +8,8 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h" + #include "src/solver/LpSolver.h" #include "src/exceptions/InvalidStateException.h" @@ -22,12 +21,12 @@ namespace storm { namespace modelchecker { template - SparseMdpPrctlModelChecker::SparseMdpPrctlModelChecker(SparseMdpModelType const& model) : SparsePropositionalModelChecker(model), MinMaxLinearEquationSolverFactory(new storm::utility::solver::MinMaxLinearEquationSolverFactory()) { + SparseMdpPrctlModelChecker::SparseMdpPrctlModelChecker(SparseMdpModelType const& model) : SparsePropositionalModelChecker(model), minMaxLinearEquationSolverFactory(new storm::utility::solver::MinMaxLinearEquationSolverFactory()) { // Intentionally left empty. } template - SparseMdpPrctlModelChecker::SparseMdpPrctlModelChecker(SparseMdpModelType const& model, std::unique_ptr>&& MinMaxLinearEquationSolverFactory) : SparsePropositionalModelChecker(model), MinMaxLinearEquationSolverFactory(std::move(MinMaxLinearEquationSolverFactory)) { + SparseMdpPrctlModelChecker::SparseMdpPrctlModelChecker(SparseMdpModelType const& model, std::unique_ptr>&& minMaxLinearEquationSolverFactory) : SparsePropositionalModelChecker(model), minMaxLinearEquationSolverFactory(std::move(minMaxLinearEquationSolverFactory)) { // Intentionally left empty. } @@ -43,8 +42,8 @@ namespace storm { std::unique_ptr rightResultPointer = this->check(pathFormula.getRightSubformula()); ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); - std::unique_ptr result = std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeBoundedUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound()))); - return result; + std::vector numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeBoundedUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template @@ -52,7 +51,8 @@ namespace storm { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); std::unique_ptr subResultPointer = this->check(pathFormula.getSubformula()); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector()))); + std::vector numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeNextProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template @@ -62,21 +62,24 @@ 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(SparseMdpPrctlModelChecker::computeUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), *MinMaxLinearEquationSolverFactory, qualitative))); + std::vector numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeUntilProbabilities(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template std::unique_ptr SparseMdpPrctlModelChecker::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeCumulativeRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getDiscreteTimeBound()))); + std::vector numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeCumulativeRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound(), *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template std::unique_ptr SparseMdpPrctlModelChecker::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, boost::optional const& rewardModelName, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound."); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeInstantaneousRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getDiscreteTimeBound()))); + std::vector numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeInstantaneousRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), rewardPathFormula.getDiscreteTimeBound(), *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template @@ -84,17 +87,17 @@ namespace storm { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); std::unique_ptr subResultPointer = this->check(rewardPathFormula.getSubformula()); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeReachabilityRewardsHelper(rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->MinMaxLinearEquationSolverFactory, qualitative))); + std::vector numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeReachabilityRewards(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), rewardModelName ? this->getModel().getRewardModel(rewardModelName.get()) : this->getModel().getRewardModel(""), subResult.getTruthValuesVector(), qualitative, *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template std::unique_ptr SparseMdpPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); - std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector(), qualitative))); + std::vector numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeLongRunAverage(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), qualitative, *minMaxLinearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } template class SparseMdpPrctlModelChecker>; diff --git a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp index e69de29bb..189a20006 100644 --- a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp @@ -0,0 +1,241 @@ +#include "src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h" + +#include "src/utility/graph.h" + +#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h" +#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h" +#include "src/modelchecker/results/HybridQuantitativeCheckResult.h" + +#include "src/exceptions/InvalidPropertyException.h" + +namespace storm { + namespace modelchecker { + namespace helper { + + template + std::unique_ptr HybridDtmcPrctlHelper::computeUntilProbabilities(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // We need to identify the states which have to be taken out of the matrix, i.e. all states that have + // probability 0 and 1 of satisfying the until-formula. + std::pair, storm::dd::Bdd> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates); + storm::dd::Bdd maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates(); + + // Perform some logging. + STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states."); + STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states."); + STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states."); + + // Check whether we need to compute exact probabilities for some states. + if (qualitative) { + // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1. + return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5))); + } else { + // If there are maybe states, we need to solve an equation system. + if (!maybeStates.isZero()) { + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd odd(maybeStates); + + // Create the matrix and the vector for the equation system. + storm::dd::Add maybeStatesAdd = maybeStates.toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all + // maybe states. + storm::dd::Add prob1StatesAsColumn = statesWithProbability01.second.toAdd(); + prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs()); + storm::dd::Add subvector = submatrix * prob1StatesAsColumn; + subvector = subvector.sumAbstract(model.getColumnVariables()); + + // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed + // for solving the equation system (i.e. compute (I-A)). + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix; + + // Create the solution vector. + std::vector x(maybeStates.getNonZeroCount(), ValueType(0.5)); + + // Translate the symbolic matrix/vector to their explicit representations and solve the equation system. + storm::storage::SparseMatrix explicitSubmatrix = submatrix.toMatrix(odd, odd); + std::vector b = subvector.template toVector(odd); + + std::unique_ptr> solver = linearEquationSolverFactory.create(explicitSubmatrix); + solver->solveEquationSystem(x, b); + + // Return a hybrid check result that stores the numerical values explicitly. + return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, statesWithProbability01.second.toAdd(), maybeStates, odd, x)); + } else { + return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), statesWithProbability01.second.toAdd())); + } + } + } + + template + std::unique_ptr HybridDtmcPrctlHelper::computeNextProbabilities(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates) { + storm::dd::Add result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd(); + return std::unique_ptr(new SymbolicQuantitativeCheckResult(model.getReachableStates(), result.sumAbstract(model.getColumnVariables()))); + } + + template + std::unique_ptr HybridDtmcPrctlHelper::computeBoundedUntilProbabilities(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // We need to identify the states which have to be taken out of the matrix, i.e. all states that have + // probability 0 or 1 of satisfying the until-formula. + storm::dd::Bdd statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound); + storm::dd::Bdd maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates(); + + // If there are maybe states, we need to perform matrix-vector multiplications. + if (!maybeStates.isZero()) { + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd odd(maybeStates); + + // Create the matrix and the vector for the equation system. + storm::dd::Add maybeStatesAdd = maybeStates.toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all + // maybe states. + storm::dd::Add prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs()); + storm::dd::Add subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables()); + + // Finally cut away all columns targeting non-maybe states. + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + + // Create the solution vector. + std::vector x(maybeStates.getNonZeroCount(), storm::utility::zero()); + + // Translate the symbolic matrix/vector to their explicit representations. + storm::storage::SparseMatrix explicitSubmatrix = submatrix.toMatrix(odd, odd); + std::vector b = subvector.template toVector(odd); + + std::unique_ptr> solver = linearEquationSolverFactory.create(explicitSubmatrix); + solver->performMatrixVectorMultiplication(x, &b, stepBound); + + // Return a hybrid check result that stores the numerical values explicitly. + return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.toAdd(), maybeStates, odd, x)); + } else { + return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), psiStates.toAdd())); + } + } + + + template + std::unique_ptr HybridDtmcPrctlHelper::computeInstantaneousRewards(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // Only compute the result if the model has at least one reward this->getModel(). + STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd odd(model.getReachableStates()); + + // Create the solution vector (and initialize it to the state rewards of the model). + std::vector x = model.getStateRewardVector().template toVector(odd); + + // Translate the symbolic matrix to its explicit representations. + storm::storage::SparseMatrix explicitMatrix = transitionMatrix.toMatrix(odd, odd); + + // Perform the matrix-vector multiplication. + std::unique_ptr> solver = linearEquationSolverFactory.create(explicitMatrix); + solver->performMatrixVectorMultiplication(x, nullptr, stepBound); + + // Return a hybrid check result that stores the numerical values explicitly. + return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x)); + } + + template + std::unique_ptr HybridDtmcPrctlHelper::computeCumulativeRewards(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // Only compute the result if the model has at least one reward this->getModel(). + STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + // Compute the reward vector to add in each step based on the available reward models. + storm::dd::Add totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero(); + if (model.hasTransitionRewards()) { + totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables()); + } + + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd odd(model.getReachableStates()); + + // Create the solution vector. + std::vector x(model.getNumberOfStates(), storm::utility::zero()); + + // Translate the symbolic matrix/vector to their explicit representations. + storm::storage::SparseMatrix explicitMatrix = transitionMatrix.toMatrix(odd, odd); + std::vector b = totalRewardVector.template toVector(odd); + + // Perform the matrix-vector multiplication. + std::unique_ptr> solver = linearEquationSolverFactory.create(explicitMatrix); + solver->performMatrixVectorMultiplication(x, &b, stepBound); + + // Return a hybrid check result that stores the numerical values explicitly. + return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x)); + } + + template + std::unique_ptr HybridDtmcPrctlHelper::computeReachabilityRewards(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::dd::Bdd const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + + // Only compute the result if there is at least one reward model. + STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); + + // Determine which states have a reward of infinity by definition. + storm::dd::Bdd infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates); + infinityStates = !infinityStates && model.getReachableStates(); + storm::dd::Bdd maybeStates = (!targetStates && !infinityStates) && model.getReachableStates(); + STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states."); + STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states."); + STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states."); + + // Check whether we need to compute exact rewards for some states. + if (qualitative) { + // Set the values for all maybe-states to 1 to indicate that their reward values + // are neither 0 nor infinity. + return std::unique_ptr(new SymbolicQuantitativeCheckResult(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one()))); + } else { + // If there are maybe states, we need to solve an equation system. + if (!maybeStates.isZero()) { + // Create the ODD for the translation between symbolic and explicit storage. + storm::dd::Odd odd(maybeStates); + + // Create the matrix and the vector for the equation system. + storm::dd::Add maybeStatesAdd = maybeStates.toAdd(); + + // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting + // non-maybe states in the matrix. + storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd; + + // Then compute the state reward vector to use in the computation. + storm::dd::Add subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero(); + if (transitionRewardMatrix) { + subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables()); + } + + // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed + // for solving the equation system (i.e. compute (I-A)). + submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); + submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix; + + // Create the solution vector. + std::vector x(maybeStates.getNonZeroCount(), ValueType(0.5)); + + // Translate the symbolic matrix/vector to their explicit representations. + storm::storage::SparseMatrix explicitSubmatrix = submatrix.toMatrix(odd, odd); + std::vector b = subvector.template toVector(odd); + + // Now solve the resulting equation system. + std::unique_ptr> solver = linearEquationSolverFactory.create(explicitSubmatrix); + solver->solveEquationSystem(x, b); + + // Return a hybrid check result that stores the numerical values explicitly. + return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity()), maybeStates, odd, x)); + } else { + return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity()))); + } + } + } + + template class HybridDtmcPrctlHelper; + } + } +} \ No newline at end of file diff --git a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h index e69de29bb..2c572bc96 100644 --- a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h +++ b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h @@ -0,0 +1,38 @@ +#ifndef STORM_MODELCHECKER_HYBRID_DTMC_PRCTL_MODELCHECKER_HELPER_H_ +#define STORM_MODELCHECKER_HYBRID_DTMC_PRCTL_MODELCHECKER_HELPER_H_ + +#include "src/models/symbolic/NondeterministicModel.h" + +#include "src/storage/dd/Add.h" +#include "src/storage/dd/Bdd.h" + +#include "src/utility/solver.h" + +namespace storm { + namespace modelchecker { + // Forward-declare result class. + class CheckResult; + + namespace helper { + + template + class HybridDtmcPrctlHelper { + public: + static std::unique_ptr computeBoundedUntilProbabilities(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::unique_ptr computeNextProbabilities(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates); + + static std::unique_ptr computeUntilProbabilities(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::unique_ptr computeCumulativeRewards(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::unique_ptr computeInstantaneousRewards(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + static std::unique_ptr computeReachabilityRewards(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::dd::Bdd const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + }; + + } + } +} + +#endif /* STORM_MODELCHECKER_HYBRID_DTMC_PRCTL_MODELCHECKER_HELPER_H_ */ diff --git a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp index 23cebc26a..c9e0ae290 100644 --- a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp @@ -1,11 +1,19 @@ #include "src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h" +#include "src/utility/graph.h" + +#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h" +#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h" +#include "src/modelchecker/results/HybridQuantitativeCheckResult.h" + +#include "src/exceptions/InvalidPropertyException.h" + namespace storm { namespace modelchecker { namespace helper { template - std::unique_ptr HybridMdpPrctlModelChecker::computeUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { + std::unique_ptr HybridMdpPrctlHelper::computeUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { // We need to identify the states which have to be taken out of the matrix, i.e. all states that have // probability 0 and 1 of satisfying the until-formula. std::pair, storm::dd::Bdd> statesWithProbability01; @@ -69,13 +77,13 @@ namespace storm { } template - storm::dd::Add HybridMdpPrctlModelChecker::computeNextProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates) { + std::unique_ptr HybridMdpPrctlHelper::computeNextProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates) { storm::dd::Add result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd(); - return result.sumAbstract(model.getColumnVariables()); + return std::unique_ptr(new storm::modelchecker::SymbolicQuantitativeCheckResult(model.getReachableStates(), result.sumAbstract(model.getColumnVariables()))); } template - std::unique_ptr HybridMdpPrctlModelChecker::computeBoundedUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { + std::unique_ptr HybridMdpPrctlHelper::computeBoundedUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { // We need to identify the states which have to be taken out of the matrix, i.e. all states that have // probability 0 or 1 of satisfying the until-formula. storm::dd::Bdd statesWithProbabilityGreater0; @@ -126,7 +134,7 @@ namespace storm { } template - std::unique_ptr HybridMdpPrctlModelChecker::computeInstantaneousRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { + std::unique_ptr HybridMdpPrctlHelper::computeInstantaneousRewards(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { // Only compute the result if the model has at least one reward this->getModel(). STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); @@ -148,7 +156,7 @@ namespace storm { } template - std::unique_ptr HybridMdpPrctlModelChecker::computeCumulativeRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { + std::unique_ptr HybridMdpPrctlHelper::computeCumulativeRewards(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { // Only compute the result if the model has at least one reward this->getModel(). STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); @@ -177,7 +185,7 @@ namespace storm { } template - std::unique_ptr HybridMdpPrctlModelChecker::computeReachabilityRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::dd::Bdd const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory, bool qualitative) { + std::unique_ptr HybridMdpPrctlHelper::computeReachabilityRewards(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::dd::Bdd const& targetStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory) { // Only compute the result if there is at least one reward model. STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); @@ -243,6 +251,8 @@ namespace storm { } } } + + template class HybridMdpPrctlHelper; } } } \ No newline at end of file diff --git a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h index e55b5762b..0984366aa 100644 --- a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h +++ b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.h @@ -1,7 +1,7 @@ -#ifndef STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_ -#define STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_ +#ifndef STORM_MODELCHECKER_HYBRID_MDP_PRCTL_MODELCHECKER_HELPER_H_ +#define STORM_MODELCHECKER_HYBRID_MDP_PRCTL_MODELCHECKER_HELPER_H_ -#include "src/models/symbolic/Model.h" +#include "src/models/symbolic/NondeterministicModel.h" #include "src/storage/dd/Add.h" #include "src/storage/dd/Bdd.h" @@ -10,6 +10,9 @@ namespace storm { namespace modelchecker { + // Forward-declare result class. + class CheckResult; + namespace helper { template @@ -17,7 +20,7 @@ namespace storm { public: static std::unique_ptr computeBoundedUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory); - static storm::dd::Add computeNextProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates); + static std::unique_ptr computeNextProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& nextStates); static std::unique_ptr computeUntilProbabilities(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, storm::dd::Bdd const& phiStates, storm::dd::Bdd const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory); @@ -26,9 +29,10 @@ namespace storm { static std::unique_ptr computeInstantaneousRewards(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory); static std::unique_ptr computeReachabilityRewards(bool minimize, storm::models::symbolic::NondeterministicModel const& model, storm::dd::Add const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::dd::Bdd const& targetStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory const& linearEquationSolverFactory); - } + }; + } } } -#endif /* STORM_MODELCHECKER_SPARSE_MDP_PRCTL_MODELCHECKER_HELPER_H_ */ +#endif /* STORM_MODELCHECKER_HYBRID_MDP_PRCTL_MODELCHECKER_HELPER_H_ */ diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 0cb95666f..07af96b33 100644 --- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -142,6 +142,22 @@ namespace storm { template std::vector SparseDtmcPrctlHelper::computeReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + return computeReachabilityRewards(transitionMatrix, backwardTransitions, [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates) { return rewardModel.getTotalRewardVector(numberOfRows, transitionMatrix, maybeStates); }, targetStates, qualitative, linearEquationSolverFactory); + } + + template + std::vector SparseDtmcPrctlHelper::computeReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& totalStateRewardVector, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + return computeReachabilityRewards(transitionMatrix, backwardTransitions, + [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& maybeStates) { + std::vector result(numberOfRows); + storm::utility::vector::selectVectorValues(result, maybeStates, totalStateRewardVector); + return result; + }, + targetStates, qualitative, linearEquationSolverFactory); + } + + template + std::vector SparseDtmcPrctlHelper::computeReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const&(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { // Determine which states have a reward of infinity by definition. storm::storage::BitVector trueStates(transitionMatrix.getRowCount(), true); @@ -174,7 +190,7 @@ namespace storm { std::vector x(submatrix.getColumnCount(), storm::utility::one()); // Prepare the right-hand side of the equation system. - std::vector b = rewardModel.getTotalRewardVector(submatrix.getRowCount(), transitionMatrix, maybeStates); + std::vector b = totalStateRewardVectorGetter(submatrix.getRowCount(), transitionMatrix, maybeStates); // Now solve the resulting equation system. std::unique_ptr> solver = linearEquationSolverFactory.create(submatrix); @@ -189,7 +205,7 @@ namespace storm { return result; } - + template class SparseDtmcPrctlHelper; } } diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h index 57d2a9de7..b01d0cec1 100644 --- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h +++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h @@ -29,7 +29,12 @@ namespace storm { static std::vector computeReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + static std::vector computeReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& totalStateRewardVector, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + static std::vector computeLongRunAverage(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); + + private: + static std::vector computeReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::function const&(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); }; } diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index ff278e20e..8203396aa 100644 --- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -115,7 +115,7 @@ namespace storm { } template - std::vector SparseMdpPrctlHelper::computeInstantaneousRewards(storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepCount, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + std::vector SparseMdpPrctlHelper::computeInstantaneousRewards(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { // Only compute the result if the model has a state-based reward this->getModel(). STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); @@ -129,7 +129,7 @@ namespace storm { } template - std::vector SparseMdpPrctlHelper::computeCumulativeRewards(storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + std::vector SparseMdpPrctlHelper::computeCumulativeRewards(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { // Only compute the result if the model has at least one reward this->getModel(). STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula."); diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h index 2c3aa99e4..1a41aa155 100644 --- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h +++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.h @@ -26,9 +26,9 @@ namespace storm { static std::vector computeUntilProbabilities(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::utility::solver::MinMaxLinearEquationSolverFactory const& MinMaxLinearEquationSolverFactory, bool qualitative); - static std::vector computeInstantaneousRewards(storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepCount, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + static std::vector computeInstantaneousRewards(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - static std::vector computeCumulativeRewards(storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, bool minimize, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + static std::vector computeCumulativeRewards(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); static std::vector computeReachabilityRewards(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory);