From b1f2c26df046fd187cdd8c77519818033df4c0cd Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 27 Jun 2016 20:58:26 +0200 Subject: [PATCH] made all instantiations to call MDP model checking with rational numbers Former-commit-id: d3f8df780412d11a13f7aafc42762acf36a54015 --- .../prctl/SparseMdpPrctlModelChecker.cpp | 1 + .../prctl/helper/SparseMdpPrctlHelper.cpp | 24 ++- src/solver/EigenLinearEquationSolver.cpp | 160 +++++------------- src/solver/EigenLinearEquationSolver.h | 2 +- .../EliminationLinearEquationSolver.cpp | 37 ---- src/solver/EliminationLinearEquationSolver.h | 1 - src/solver/GmmxxLinearEquationSolver.cpp | 44 +---- src/solver/GmmxxLinearEquationSolver.h | 1 - src/solver/LinearEquationSolver.h | 2 +- src/solver/MinMaxLinearEquationSolver.cpp | 14 +- src/solver/NativeLinearEquationSolver.cpp | 37 ---- src/solver/NativeLinearEquationSolver.h | 1 - src/solver/SolveGoal.cpp | 7 +- .../StandardMinMaxLinearEquationSolver.cpp | 81 ++++++--- .../StandardMinMaxLinearEquationSolver.h | 23 +-- src/solver/TerminationCondition.cpp | 8 +- src/solver/TerminationCondition.h | 2 +- .../MaximalEndComponentDecomposition.cpp | 3 + src/utility/constants.cpp | 14 ++ src/utility/constants.h | 3 + src/utility/graph.cpp | 4 + src/utility/storm.h | 28 +-- src/utility/vector.h | 7 +- 23 files changed, 198 insertions(+), 306 deletions(-) diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index f58539f05..629cf29eb 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -147,5 +147,6 @@ namespace storm { } template class SparseMdpPrctlModelChecker>; + template class SparseMdpPrctlModelChecker>; } } \ No newline at end of file diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 7879ff845..0b43a524b 100644 --- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -20,6 +20,7 @@ #include "src/exceptions/InvalidStateException.h" #include "src/exceptions/InvalidPropertyException.h" #include "src/exceptions/InvalidSettingsException.h" +#include "src/exceptions/IllegalFunctionCallException.h" namespace storm { namespace modelchecker { @@ -70,7 +71,6 @@ namespace storm { return result; } - template MDPSparseModelCheckingHelperReturnType SparseMdpPrctlHelper::computeUntilProbabilities(storm::solver::SolveGoal const& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, bool produceScheduler, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { @@ -107,7 +107,7 @@ namespace storm { // 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. - storm::utility::vector::setVectorValues(result, maybeStates, ValueType(0.5)); + storm::utility::vector::setVectorValues(result, maybeStates, storm::utility::convertNumber(0.5)); } else { if (!maybeStates.empty()) { // In this case we have have to compute the probabilities. @@ -270,6 +270,11 @@ namespace storm { }, \ targetStates, qualitative, minMaxLinearEquationSolverFactory); } + + template<> + std::vector SparseMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& intervalRewardModel, bool lowerBoundOfIntervals, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + STORM_LOG_THROW(false, storm::exceptions::IllegalFunctionCallException, "Computing reachability rewards is unsupported for this data type."); + } #endif template @@ -357,7 +362,7 @@ namespace storm { } // Start by decomposing the MDP into its MECs. - storm::storage::MaximalEndComponentDecomposition mecDecomposition(transitionMatrix, backwardTransitions); + storm::storage::MaximalEndComponentDecomposition mecDecomposition(transitionMatrix, backwardTransitions); // Get some data members for convenience. std::vector const& nondeterministicChoiceIndices = transitionMatrix.getRowGroupIndices(); @@ -517,12 +522,12 @@ namespace storm { ValueType r = 0; for (auto element : transitionMatrix.getRow(choice)) { - constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); + constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(storm::utility::convertNumber(element.getValue())); if (psiStates.get(element.getColumn())) { r += element.getValue(); } } - constraint = solver->getConstant(r) + constraint; + constraint = solver->getConstant(storm::utility::convertNumber(r)) + constraint; if (dir == OptimizationDirection::Minimize) { constraint = stateToVariableMap.at(state) <= constraint; @@ -534,7 +539,7 @@ namespace storm { } solver->optimize(); - return solver->getContinuousValue(lambda); + return storm::utility::convertNumber(solver->getContinuousValue(lambda)); } template @@ -644,7 +649,12 @@ namespace storm { template std::vector SparseMdpPrctlHelper::computeInstantaneousRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, uint_fast64_t stepCount, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template std::vector SparseMdpPrctlHelper::computeCumulativeRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template std::vector SparseMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - + + template class SparseMdpPrctlHelper; + template std::vector SparseMdpPrctlHelper::computeInstantaneousRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, uint_fast64_t stepCount, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template std::vector SparseMdpPrctlHelper::computeCumulativeRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::models::sparse::StandardRewardModel const& rewardModel, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template std::vector SparseMdpPrctlHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + } } } diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp index e7d4bd2d6..a910801f3 100644 --- a/src/solver/EigenLinearEquationSolver.cpp +++ b/src/solver/EigenLinearEquationSolver.cpp @@ -224,45 +224,45 @@ namespace storm { } } - template - void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { - // Typedef the map-type so we don't have to spell it out. - typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; - - bool multiplyResultProvided = multiplyResult != nullptr; - if (!multiplyResult) { - multiplyResult = new std::vector(eigenA->cols()); - } - auto eigenMultiplyResult = Eigen::Matrix::Map(multiplyResult->data(), multiplyResult->size()); - - // Map the input vectors x and b to Eigen's format. - std::unique_ptr eigenB; - if (b != nullptr) { - eigenB = std::make_unique(Eigen::Matrix::Map(b->data(), b->size())); - } - auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); - - // Perform n matrix-vector multiplications. - auto currentX = &eigenX; - auto nextX = &eigenMultiplyResult; - for (uint64_t iteration = 0; iteration < n; ++iteration) { - if (eigenB) { - nextX->noalias() = *eigenA * *currentX + *eigenB; - } else { - nextX->noalias() = *eigenA * *currentX; - } - std::swap(nextX, currentX); - } - - // If the last result we obtained is not the one in the input vector x, we swap the result there. - if (currentX != &eigenX) { - std::swap(*nextX, *currentX); - } - - if (!multiplyResultProvided) { - delete multiplyResult; - } - } +// template +// void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { +// // Typedef the map-type so we don't have to spell it out. +// typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; +// +// bool multiplyResultProvided = multiplyResult != nullptr; +// if (!multiplyResult) { +// multiplyResult = new std::vector(eigenA->cols()); +// } +// auto eigenMultiplyResult = Eigen::Matrix::Map(multiplyResult->data(), multiplyResult->size()); +// +// // Map the input vectors x and b to Eigen's format. +// std::unique_ptr eigenB; +// if (b != nullptr) { +// eigenB = std::make_unique(Eigen::Matrix::Map(b->data(), b->size())); +// } +// auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); +// +// // Perform n matrix-vector multiplications. +// auto currentX = &eigenX; +// auto nextX = &eigenMultiplyResult; +// for (uint64_t iteration = 0; iteration < n; ++iteration) { +// if (eigenB) { +// nextX->noalias() = *eigenA * *currentX + *eigenB; +// } else { +// nextX->noalias() = *eigenA * *currentX; +// } +// std::swap(nextX, currentX); +// } +// +// // If the last result we obtained is not the one in the input vector x, we swap the result there. +// if (currentX != &eigenX) { +// std::swap(*nextX, *currentX); +// } +// +// if (!multiplyResultProvided) { +// delete multiplyResult; +// } +// } template EigenLinearEquationSolverSettings& EigenLinearEquationSolver::getSettings() { @@ -287,46 +287,6 @@ namespace storm { solver._solve_impl(eigenB, eigenX); } - template<> - void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { - // Typedef the map-type so we don't have to spell it out. - typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; - - bool multiplyResultProvided = multiplyResult != nullptr; - if (!multiplyResult) { - multiplyResult = new std::vector(eigenA->cols()); - } - auto eigenMultiplyResult = Eigen::Matrix::Map(multiplyResult->data(), multiplyResult->size()); - - // Map the input vectors x and b to Eigen's format. - std::unique_ptr eigenB; - if (b != nullptr) { - eigenB = std::make_unique(Eigen::Matrix::Map(b->data(), b->size())); - } - auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); - - // Perform n matrix-vector multiplications. - auto currentX = &eigenX; - auto nextX = &eigenMultiplyResult; - for (uint64_t iteration = 0; iteration < n; ++iteration) { - if (eigenB) { - nextX->noalias() = *eigenA * *currentX + *eigenB; - } else { - nextX->noalias() = *eigenA * *currentX; - } - std::swap(nextX, currentX); - } - - // If the last result we obtained is not the one in the input vector x, we swap the result there. - if (currentX != &eigenX) { - std::swap(*nextX, *currentX); - } - - if (!multiplyResultProvided) { - delete multiplyResult; - } - } - // Specialization form storm::RationalFunction template<> @@ -339,47 +299,7 @@ namespace storm { solver.compute(*eigenA); solver._solve_impl(eigenB, eigenX); } - - template<> - void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { - // Typedef the map-type so we don't have to spell it out. - typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; - - bool multiplyResultProvided = multiplyResult != nullptr; - if (!multiplyResult) { - multiplyResult = new std::vector(eigenA->cols()); - } - auto eigenMultiplyResult = Eigen::Matrix::Map(multiplyResult->data(), multiplyResult->size()); - - // Map the input vectors x and b to Eigen's format. - std::unique_ptr eigenB; - if (b != nullptr) { - eigenB = std::make_unique(Eigen::Matrix::Map(b->data(), b->size())); - } - auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); - - // Perform n matrix-vector multiplications. - auto currentX = &eigenX; - auto nextX = &eigenMultiplyResult; - for (uint64_t iteration = 0; iteration < n; ++iteration) { - if (eigenB) { - nextX->noalias() = *eigenA * *currentX + *eigenB; - } else { - nextX->noalias() = *eigenA * *currentX; - } - std::swap(nextX, currentX); - } - - // If the last result we obtained is not the one in the input vector x, we swap the result there. - if (currentX != &eigenX) { - std::swap(*nextX, *currentX); - } - - if (!multiplyResultProvided) { - delete multiplyResult; - } - } - + template std::unique_ptr> EigenLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::make_unique>(matrix, settings); diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h index 7dbd5d519..2f3f4ca23 100644 --- a/src/solver/EigenLinearEquationSolver.h +++ b/src/solver/EigenLinearEquationSolver.h @@ -62,7 +62,7 @@ namespace storm { EigenLinearEquationSolver(storm::storage::SparseMatrix&& A, EigenLinearEquationSolverSettings const& settings = EigenLinearEquationSolverSettings()); virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; +// virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; virtual void performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; EigenLinearEquationSolverSettings& getSettings(); diff --git a/src/solver/EliminationLinearEquationSolver.cpp b/src/solver/EliminationLinearEquationSolver.cpp index 36c1144a3..a56e1e584 100644 --- a/src/solver/EliminationLinearEquationSolver.cpp +++ b/src/solver/EliminationLinearEquationSolver.cpp @@ -100,43 +100,6 @@ namespace storm { }; } - template - void EliminationLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { - // Set up some temporary variables so that we can just swap pointers instead of copying the result after - // each iteration. - std::vector* currentX = &x; - - bool multiplyResultProvided = true; - std::vector* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; - - // Now perform matrix-vector multiplication as long as we meet the bound. - for (uint_fast64_t i = 0; i < n; ++i) { - A.multiplyWithVector(*currentX, *nextX); - std::swap(nextX, currentX); - - // If requested, add an offset to the current result vector. - if (b != nullptr) { - storm::utility::vector::addVectors(*currentX, *b, *currentX); - } - } - - // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, - // because the output is supposed to be stored in the input vector x. - if (currentX == copyX) { - std::swap(x, *currentX); - } - - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; - } - } - template void EliminationLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b) const { if (&x != &result) { diff --git a/src/solver/EliminationLinearEquationSolver.h b/src/solver/EliminationLinearEquationSolver.h index 43d93c42e..580a77438 100644 --- a/src/solver/EliminationLinearEquationSolver.h +++ b/src/solver/EliminationLinearEquationSolver.h @@ -30,7 +30,6 @@ namespace storm { EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings = EliminationLinearEquationSolverSettings()); virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; virtual void performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; EliminationLinearEquationSolverSettings& getSettings(); diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index d9d331f9f..53e643741 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -177,48 +177,12 @@ namespace storm { } } - template - void GmmxxLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { - // Set up some temporary variables so that we can just swap pointers instead of copying the result after - // each iteration. - std::vector* currentX = &x; - - bool multiplyResultProvided = true; - std::vector* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; - - // Now perform matrix-vector multiplication as long as we meet the bound. - for (uint_fast64_t i = 0; i < n; ++i) { - gmm::mult(*gmmxxMatrix, *currentX, *nextX); - std::swap(nextX, currentX); - - // If requested, add an offset to the current result vector. - if (b != nullptr) { - gmm::add(*b, *currentX); - } - } - - // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, - // because the output is supposed to be stored in the input vector x. - if (currentX == copyX) { - std::swap(x, *currentX); - } - - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; - } - } - template void GmmxxLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b) const { - gmm::mult(*gmmxxMatrix, x, result); - if (b != nullptr) { - gmm::add(*b, result); + if (b) { + gmm::mult_add(*gmmxxMatrix, x, *b, result); + } else { + gmm::mult(*gmmxxMatrix, x, result); } } diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index ca14791f2..da57e8325 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -90,7 +90,6 @@ namespace storm { GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings = GmmxxLinearEquationSolverSettings()); virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; virtual void performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; GmmxxLinearEquationSolverSettings& getSettings(); diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h index a6755872e..e12fa6eb5 100644 --- a/src/solver/LinearEquationSolver.h +++ b/src/solver/LinearEquationSolver.h @@ -58,7 +58,7 @@ namespace storm { * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this * vector must be equal to the number of rows of A. */ - virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const; + void performMatrixVectorMultiplication(std::vector& x, std::vector const* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const; }; template diff --git a/src/solver/MinMaxLinearEquationSolver.cpp b/src/solver/MinMaxLinearEquationSolver.cpp index 22cefe139..f068b75aa 100644 --- a/src/solver/MinMaxLinearEquationSolver.cpp +++ b/src/solver/MinMaxLinearEquationSolver.cpp @@ -128,12 +128,24 @@ namespace storm { } return result; } - + + template<> + template + std::unique_ptr> GeneralMinMaxLinearEquationSolverFactory::selectSolver(MatrixType&& matrix) const { + std::unique_ptr> result; + auto method = storm::settings::getModule().getMinMaxEquationSolvingMethod(); + STORM_LOG_THROW(method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration, storm::exceptions::InvalidSettingsException, "For this data type only value iteration and policy iteration are available."); + return std::make_unique>(std::forward(matrix), std::make_unique>()); + } + template class MinMaxLinearEquationSolver; template class MinMaxLinearEquationSolver; + template class MinMaxLinearEquationSolver; template class MinMaxLinearEquationSolverFactory; template class GeneralMinMaxLinearEquationSolverFactory; + template class MinMaxLinearEquationSolverFactory; + template class GeneralMinMaxLinearEquationSolverFactory; } } diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index 5c73d4b5e..ef4e1c2ff 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -182,43 +182,6 @@ namespace storm { } } - template - void NativeLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { - // Set up some temporary variables so that we can just swap pointers instead of copying the result after - // each iteration. - std::vector* currentX = &x; - - bool multiplyResultProvided = true; - std::vector* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; - - // Now perform matrix-vector multiplication as long as we meet the bound. - for (uint_fast64_t i = 0; i < n; ++i) { - A.multiplyWithVector(*currentX, *nextX); - std::swap(nextX, currentX); - - // If requested, add an offset to the current result vector. - if (b != nullptr) { - storm::utility::vector::addVectors(*currentX, *b, *currentX); - } - } - - // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, - // because the output is supposed to be stored in the input vector x. - if (currentX == copyX) { - std::swap(x, *currentX); - } - - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; - } - } - template void NativeLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b) const { if (&x != &result) { diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index c76a45b4e..6d6125741 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -47,7 +47,6 @@ namespace storm { NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings = NativeLinearEquationSolverSettings()); virtual void solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void performMatrixVectorMultiplication(std::vector& x, std::vector const* b, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; virtual void performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; NativeLinearEquationSolverSettings& getSettings(); diff --git a/src/solver/SolveGoal.cpp b/src/solver/SolveGoal.cpp index ebe8d8467..2f0cd365c 100644 --- a/src/solver/SolveGoal.cpp +++ b/src/solver/SolveGoal.cpp @@ -17,7 +17,7 @@ namespace storm { std::unique_ptr> configureMinMaxLinearEquationSolver(BoundedGoal const& goal, storm::solver::MinMaxLinearEquationSolverFactory const& factory, storm::storage::SparseMatrix const& matrix) { std::unique_ptr> p = factory.create(matrix); p->setOptimizationDirection(goal.direction()); - p->setTerminationCondition(std::make_unique>(goal.relevantValues(), goal.boundIsStrict(), goal.thresholdValue(), goal.minimize())); + p->setTerminationCondition(std::make_unique>(goal.relevantValues(), goal.boundIsStrict(), goal.thresholdValue(), goal.minimize())); return p; } @@ -50,6 +50,9 @@ namespace storm { template std::unique_ptr> configureMinMaxLinearEquationSolver(SolveGoal const& goal, storm::solver::MinMaxLinearEquationSolverFactory const& factory, storm::storage::SparseMatrix const& matrix); template std::unique_ptr> configureLinearEquationSolver(BoundedGoal const& goal, storm::solver::LinearEquationSolverFactory const& factory, storm::storage::SparseMatrix const& matrix); template std::unique_ptr> configureLinearEquationSolver(SolveGoal const& goal, storm::solver::LinearEquationSolverFactory const& factory, storm::storage::SparseMatrix const& matrix); - + + template std::unique_ptr> configureMinMaxLinearEquationSolver(BoundedGoal const& goal, storm::solver::MinMaxLinearEquationSolverFactory const& factory, storm::storage::SparseMatrix const& matrix); + template std::unique_ptr> configureMinMaxLinearEquationSolver(SolveGoal const& goal, storm::solver::MinMaxLinearEquationSolverFactory const& factory, storm::storage::SparseMatrix const& matrix); + } } diff --git a/src/solver/StandardMinMaxLinearEquationSolver.cpp b/src/solver/StandardMinMaxLinearEquationSolver.cpp index 0e9791969..f485e0de2 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp @@ -15,12 +15,13 @@ namespace storm { namespace solver { - StandardMinMaxLinearEquationSolverSettings::StandardMinMaxLinearEquationSolverSettings() { + template + StandardMinMaxLinearEquationSolverSettings::StandardMinMaxLinearEquationSolverSettings() { // Get the settings object to customize linear solving. storm::settings::modules::MinMaxEquationSolverSettings const& settings = storm::settings::getModule(); maximalNumberOfIterations = settings.getMaximalIterationCount(); - precision = settings.getPrecision(); + precision = storm::utility::convertNumber(settings.getPrecision()); relative = settings.getConvergenceCriterion() == storm::settings::modules::MinMaxEquationSolverSettings::ConvergenceCriterion::Relative; auto method = settings.getMinMaxEquationSolvingMethod(); @@ -32,53 +33,61 @@ namespace storm { } } - void StandardMinMaxLinearEquationSolverSettings::setSolutionMethod(SolutionMethod const& solutionMethod) { + template + void StandardMinMaxLinearEquationSolverSettings::setSolutionMethod(SolutionMethod const& solutionMethod) { this->solutionMethod = solutionMethod; } - void StandardMinMaxLinearEquationSolverSettings::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) { + template + void StandardMinMaxLinearEquationSolverSettings::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) { this->maximalNumberOfIterations = maximalNumberOfIterations; } - void StandardMinMaxLinearEquationSolverSettings::setRelativeTerminationCriterion(bool value) { + template + void StandardMinMaxLinearEquationSolverSettings::setRelativeTerminationCriterion(bool value) { this->relative = value; } - void StandardMinMaxLinearEquationSolverSettings::setPrecision(double precision) { + template + void StandardMinMaxLinearEquationSolverSettings::setPrecision(ValueType precision) { this->precision = precision; } - StandardMinMaxLinearEquationSolverSettings::SolutionMethod const& StandardMinMaxLinearEquationSolverSettings::getSolutionMethod() const { + template + typename StandardMinMaxLinearEquationSolverSettings::SolutionMethod const& StandardMinMaxLinearEquationSolverSettings::getSolutionMethod() const { return solutionMethod; } - uint64_t StandardMinMaxLinearEquationSolverSettings::getMaximalNumberOfIterations() const { + template + uint64_t StandardMinMaxLinearEquationSolverSettings::getMaximalNumberOfIterations() const { return maximalNumberOfIterations; } - double StandardMinMaxLinearEquationSolverSettings::getPrecision() const { + template + ValueType StandardMinMaxLinearEquationSolverSettings::getPrecision() const { return precision; } - bool StandardMinMaxLinearEquationSolverSettings::getRelativeTerminationCriterion() const { + template + bool StandardMinMaxLinearEquationSolverSettings::getRelativeTerminationCriterion() const { return relative; } template - StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings) : settings(settings), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(nullptr), A(A) { + StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings) : settings(settings), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(nullptr), A(A) { // Intentionally left empty. } template - StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings) : settings(settings), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(std::make_unique>(std::move(A))), A(*localA) { + StandardMinMaxLinearEquationSolver::StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings) : settings(settings), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), localA(std::make_unique>(std::move(A))), A(*localA) { // Intentionally left empty. } template void StandardMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { switch (this->getSettings().getSolutionMethod()) { - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: solveEquationSystemValueIteration(dir, x, b, multiplyResult, newX); break; - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: solveEquationSystemPolicyIteration(dir, x, b, multiplyResult, newX); break; + case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: solveEquationSystemValueIteration(dir, x, b, multiplyResult, newX); break; + case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: solveEquationSystemPolicyIteration(dir, x, b, multiplyResult, newX); break; } } @@ -90,7 +99,7 @@ namespace storm { if (multiplyResult == nullptr) { multiplyResult = new std::vector(this->A.getRowCount()); } - + // Create the initial scheduler. std::vector scheduler(this->A.getRowGroupCount()); @@ -107,12 +116,14 @@ namespace storm { storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); submatrix.convertToEquationSystem(); storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); - + // Solve the equation system for the 'DTMC'. // FIXME: we need to remove the 0- and 1- states to make the solution unique. + // HOWEVER: if we start with a valid scheduler, then we will never get an illegal one, because staying + // within illegal MECs will never strictly improve the value. Is this true? auto solver = linearEquationSolverFactory->create(submatrix); solver->solveEquationSystem(x, subB, &deterministicMultiplyResult); - + // Go through the multiplication result and see whether we can improve any of the choices. bool schedulerImproved = false; for (uint_fast64_t group = 0; group < this->A.getRowGroupCount(); ++group) { @@ -121,7 +132,7 @@ namespace storm { if (choice - this->A.getRowGroupIndices()[group] == scheduler[group]) { continue; } - + // Create the value of the choice. ValueType choiceValue = storm::utility::zero(); for (auto const& entry : this->A.getRow(choice)) { @@ -204,7 +215,7 @@ namespace storm { storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, *newX, this->A.getRowGroupIndices()); // Determine whether the method converged. - if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion())) { + if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) { status = Status::Converged; } @@ -278,12 +289,12 @@ namespace storm { } template - StandardMinMaxLinearEquationSolverSettings const& StandardMinMaxLinearEquationSolver::getSettings() const { + StandardMinMaxLinearEquationSolverSettings const& StandardMinMaxLinearEquationSolver::getSettings() const { return settings; } template - StandardMinMaxLinearEquationSolverSettings& StandardMinMaxLinearEquationSolver::getSettings() { + StandardMinMaxLinearEquationSolverSettings& StandardMinMaxLinearEquationSolver::getSettings() { return settings; } @@ -300,10 +311,20 @@ namespace storm { template StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler) { switch (solverType) { - case EquationSolverType::Gmmxx: linearEquationSolverFactory = std::make_unique>(); break; - case EquationSolverType::Eigen: linearEquationSolverFactory = std::make_unique>(); break; - case EquationSolverType::Native: linearEquationSolverFactory = std::make_unique>(); break; - case EquationSolverType::Elimination: linearEquationSolverFactory = std::make_unique>(); break; + case EquationSolverType::Gmmxx: linearEquationSolverFactory = std::make_unique>(); break; + case EquationSolverType::Eigen: linearEquationSolverFactory = std::make_unique>(); break; + case EquationSolverType::Native: linearEquationSolverFactory = std::make_unique>(); break; + case EquationSolverType::Elimination: linearEquationSolverFactory = std::make_unique>(); break; + } + } + + template<> + StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(EquationSolverType const& solverType, bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler) { + switch (solverType) { + case EquationSolverType::Eigen: linearEquationSolverFactory = std::make_unique>(); break; + case EquationSolverType::Elimination: linearEquationSolverFactory = std::make_unique>(); break; + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Cannot create the requested solver for this data type."); } } @@ -331,12 +352,12 @@ namespace storm { } template - StandardMinMaxLinearEquationSolverSettings& StandardMinMaxLinearEquationSolverFactory::getSettings() { + StandardMinMaxLinearEquationSolverSettings& StandardMinMaxLinearEquationSolverFactory::getSettings() { return settings; } template - StandardMinMaxLinearEquationSolverSettings const& StandardMinMaxLinearEquationSolverFactory::getSettings() const { + StandardMinMaxLinearEquationSolverSettings const& StandardMinMaxLinearEquationSolverFactory::getSettings() const { return settings; } @@ -361,12 +382,16 @@ namespace storm { } template class StandardMinMaxLinearEquationSolver; - template class StandardMinMaxLinearEquationSolverFactory; template class GmmxxMinMaxLinearEquationSolverFactory; template class EigenMinMaxLinearEquationSolverFactory; template class NativeMinMaxLinearEquationSolverFactory; template class EliminationMinMaxLinearEquationSolverFactory; + template class StandardMinMaxLinearEquationSolver; + template class StandardMinMaxLinearEquationSolverFactory; + template class EigenMinMaxLinearEquationSolverFactory; + template class EliminationMinMaxLinearEquationSolverFactory; + } } \ No newline at end of file diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h index e3931a13d..0bf4349be 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/solver/StandardMinMaxLinearEquationSolver.h @@ -6,6 +6,7 @@ namespace storm { namespace solver { + template class StandardMinMaxLinearEquationSolverSettings { public: StandardMinMaxLinearEquationSolverSettings(); @@ -17,31 +18,31 @@ namespace storm { void setSolutionMethod(SolutionMethod const& solutionMethod); void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations); void setRelativeTerminationCriterion(bool value); - void setPrecision(double precision); + void setPrecision(ValueType precision); SolutionMethod const& getSolutionMethod() const; uint64_t getMaximalNumberOfIterations() const; - double getPrecision() const; + ValueType getPrecision() const; bool getRelativeTerminationCriterion() const; private: SolutionMethod solutionMethod; uint64_t maximalNumberOfIterations; - double precision; + ValueType precision; bool relative; }; template class StandardMinMaxLinearEquationSolver : public MinMaxLinearEquationSolver { public: - StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); - StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); + StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); + StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); virtual void solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override; virtual void performMatrixVectorMultiplication(OptimizationDirection dir, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; - StandardMinMaxLinearEquationSolverSettings const& getSettings() const; - StandardMinMaxLinearEquationSolverSettings& getSettings(); + StandardMinMaxLinearEquationSolverSettings const& getSettings() const; + StandardMinMaxLinearEquationSolverSettings& getSettings(); private: void solveEquationSystemPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const; @@ -57,7 +58,7 @@ namespace storm { void reportStatus(Status status, uint64_t iterations) const; /// The settings of this solver. - StandardMinMaxLinearEquationSolverSettings settings; + StandardMinMaxLinearEquationSolverSettings settings; /// The factory used to obtain linear equation solvers. std::unique_ptr> linearEquationSolverFactory; @@ -81,11 +82,11 @@ namespace storm { virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const override; virtual std::unique_ptr> create(storm::storage::SparseMatrix&& matrix) const override; - StandardMinMaxLinearEquationSolverSettings& getSettings(); - StandardMinMaxLinearEquationSolverSettings const& getSettings() const; + StandardMinMaxLinearEquationSolverSettings& getSettings(); + StandardMinMaxLinearEquationSolverSettings const& getSettings() const; private: - StandardMinMaxLinearEquationSolverSettings settings; + StandardMinMaxLinearEquationSolverSettings settings; std::unique_ptr> linearEquationSolverFactory; }; diff --git a/src/solver/TerminationCondition.cpp b/src/solver/TerminationCondition.cpp index f640ac7e0..0503723bc 100644 --- a/src/solver/TerminationCondition.cpp +++ b/src/solver/TerminationCondition.cpp @@ -1,6 +1,8 @@ #include "src/solver/TerminationCondition.h" #include "src/utility/vector.h" +#include "src/adapters/CarlAdapter.h" + #include "src/utility/macros.h" namespace storm { @@ -24,7 +26,7 @@ namespace storm { } template - TerminateIfFilteredExtremumExceedsThreshold::TerminateIfFilteredExtremumExceedsThreshold(storm::storage::BitVector const& filter, ValueType const& threshold, bool strict, bool useMinimum) : TerminateIfFilteredSumExceedsThreshold(filter, threshold, strict), useMinimum(useMinimum) { + TerminateIfFilteredExtremumExceedsThreshold::TerminateIfFilteredExtremumExceedsThreshold(storm::storage::BitVector const& filter, bool strict, ValueType const& threshold, bool useMinimum) : TerminateIfFilteredSumExceedsThreshold(filter, threshold, strict), useMinimum(useMinimum) { // Intentionally left empty. } @@ -37,5 +39,9 @@ namespace storm { template class TerminateIfFilteredSumExceedsThreshold; template class TerminateIfFilteredExtremumExceedsThreshold; + + template class TerminateIfFilteredSumExceedsThreshold; + template class TerminateIfFilteredExtremumExceedsThreshold; + } } diff --git a/src/solver/TerminationCondition.h b/src/solver/TerminationCondition.h index 42437bc68..c685d3ca6 100644 --- a/src/solver/TerminationCondition.h +++ b/src/solver/TerminationCondition.h @@ -34,7 +34,7 @@ namespace storm { template class TerminateIfFilteredExtremumExceedsThreshold : public TerminateIfFilteredSumExceedsThreshold{ public: - TerminateIfFilteredExtremumExceedsThreshold(storm::storage::BitVector const& filter, ValueType const& threshold, bool strict, bool useMinimum); + TerminateIfFilteredExtremumExceedsThreshold(storm::storage::BitVector const& filter, bool strict, ValueType const& threshold, bool useMinimum); bool terminateNow(std::vector const& currentValue) const; protected: diff --git a/src/storage/MaximalEndComponentDecomposition.cpp b/src/storage/MaximalEndComponentDecomposition.cpp index d977b54d0..d410b2546 100644 --- a/src/storage/MaximalEndComponentDecomposition.cpp +++ b/src/storage/MaximalEndComponentDecomposition.cpp @@ -193,5 +193,8 @@ namespace storm { template class MaximalEndComponentDecomposition; template MaximalEndComponentDecomposition::MaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel const& model); + template class MaximalEndComponentDecomposition; + template MaximalEndComponentDecomposition::MaximalEndComponentDecomposition(storm::models::sparse::NondeterministicModel const& model); + } } diff --git a/src/utility/constants.cpp b/src/utility/constants.cpp index 3075dc5fd..4821d1b94 100644 --- a/src/utility/constants.cpp +++ b/src/utility/constants.cpp @@ -109,6 +109,11 @@ namespace storm { return number; } + template + ValueType abs(ValueType const& number) { + return std::fabs(number); + } + template<> RationalFunction& simplify(RationalFunction& value); @@ -153,6 +158,11 @@ namespace storm { return RationalFunction(carl::rationalize(number)); } + template<> + storm::RationalNumber abs(storm::RationalNumber const& number) { + return carl::abs(number); + } + template storm::storage::MatrixEntry simplify(storm::storage::MatrixEntry matrixEntry) { simplify(matrixEntry.getValue()); @@ -188,6 +198,8 @@ namespace storm { template storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry); template storm::storage::MatrixEntry&& simplify(storm::storage::MatrixEntry&& matrixEntry); + template double abs(double const& number); + template bool isOne(float const& value); template bool isZero(float const& value); template bool isConstant(float const& value); @@ -252,6 +264,8 @@ namespace storm { template double convertNumber(storm::RationalNumber const& number); template storm::RationalNumber convertNumber(double const& number); + template storm::RationalNumber abs(storm::RationalNumber const& number); + // template storm::RationalNumber pow(storm::RationalNumber const& value, uint_fast64_t exponent); template storm::RationalNumber simplify(storm::RationalNumber value); diff --git a/src/utility/constants.h b/src/utility/constants.h index 80222b74f..4ec4e7539 100644 --- a/src/utility/constants.h +++ b/src/utility/constants.h @@ -53,6 +53,9 @@ namespace storm { template TargetType convertNumber(SourceType const& number); + + template + ValueType abs(ValueType const& number); } } diff --git a/src/utility/graph.cpp b/src/utility/graph.cpp index 58ff2efd8..8a7e1e185 100644 --- a/src/utility/graph.cpp +++ b/src/utility/graph.cpp @@ -1200,7 +1200,11 @@ namespace storm { template std::pair performProb01(storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates); + template storm::storage::PartialScheduler computeSchedulerProbGreater0E(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, boost::optional const& rowFilter); + template storm::storage::PartialScheduler computeSchedulerProb0E(storm::storage::BitVector const& prob0EStates, storm::storage::SparseMatrix const& transitionMatrix); + + template storm::storage::PartialScheduler computeSchedulerProb1E(storm::storage::BitVector const& prob1EStates, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates); template storm::storage::BitVector performProbGreater0E(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& nondeterministicChoiceIndices, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) ; diff --git a/src/utility/storm.h b/src/utility/storm.h index e85c098d0..41a17a4c7 100644 --- a/src/utility/storm.h +++ b/src/utility/storm.h @@ -294,6 +294,18 @@ namespace storm { return result; } + template + std::unique_ptr verifySparseMdp(std::shared_ptr> mdp, storm::modelchecker::CheckTask const& task) { + std::unique_ptr result; + storm::modelchecker::SparseMdpPrctlModelChecker> modelchecker(*mdp); + if (modelchecker.canHandle(task)) { + result = modelchecker.check(task); + } else { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The property " << task.getFormula() << " is not supported."); + } + return result; + } + template std::unique_ptr verifySparseModel(std::shared_ptr> model, std::shared_ptr const& formula, bool onlyInitialStatesRelevant = false) { storm::modelchecker::CheckTask task(*formula, onlyInitialStatesRelevant); @@ -302,19 +314,7 @@ namespace storm { if (model->getType() == storm::models::ModelType::Dtmc) { result = verifySparseDtmc(model->template as>(), task); } else if (model->getType() == storm::models::ModelType::Mdp) { - std::shared_ptr> mdp = model->template as>(); -#ifdef STORM_HAVE_CUDA - if (storm::settings::getModule().isCudaSet()) { - storm::modelchecker::TopologicalValueIterationMdpPrctlModelChecker modelchecker(*mdp); - result = modelchecker.check(task); - } else { - storm::modelchecker::SparseMdpPrctlModelChecker> modelchecker(*mdp); - result = modelchecker.check(task); - } -#else - storm::modelchecker::SparseMdpPrctlModelChecker> modelchecker(*mdp); - result = modelchecker.check(task); -#endif + result = verifySparseMdp(model->template as>(), task); } else if (model->getType() == storm::models::ModelType::Ctmc) { result = verifySparseCtmc(model->template as>(), task); } else if (model->getType() == storm::models::ModelType::MarkovAutomaton) { @@ -341,6 +341,8 @@ namespace storm { result = verifySparseDtmc(model->template as>(), task); } else if (model->getType() == storm::models::ModelType::Ctmc) { result = verifySparseCtmc(model->template as>(), task); + } else if (model->getType() == storm::models::ModelType::Mdp) { + result = verifySparseMdp(model->template as>(), task); } else { STORM_LOG_ASSERT(false, "Illegal model type."); } diff --git a/src/utility/vector.h b/src/utility/vector.h index 0a4da1e03..6faf88c7c 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -12,6 +12,7 @@ #include #include "src/storage/BitVector.h" +#include "src/utility/constants.h" #include "src/utility/macros.h" #include "src/solver/OptimizationDirection.h" @@ -590,13 +591,13 @@ namespace storm { bool equalModuloPrecision(T const& val1, T const& val2, T precision, bool relativeError = true) { if (relativeError) { if (val2 == 0) { - return (std::abs(val1) <= precision); + return (storm::utility::abs(val1) <= precision); } - if (std::abs((val1 - val2)/val2) > precision) { + if (storm::utility::abs((val1 - val2)/val2) > precision) { return false; } } else { - if (std::abs(val1 - val2) > precision) return false; + if (storm::utility::abs(val1 - val2) > precision) return false; } return true; }