From f63e5fc8731d82bc912bfcaf6f5c4d055194bfb9 Mon Sep 17 00:00:00 2001 From: PBerger Date: Fri, 29 May 2015 22:24:06 +0200 Subject: [PATCH] Implemented Policy Iteration inside the GmmxxMinMaxLinearEquationSolver. Added an option for selecting Value- or Policy Iteration in the GeneralSettings. Former-commit-id: 6d12f10f60418ee64347ed9e0875552938443c21 --- src/settings/modules/GeneralSettings.cpp | 19 ++ src/settings/modules/GeneralSettings.h | 22 +- .../GmmxxMinMaxLinearEquationSolver.cpp | 215 +++++++++++++----- src/solver/GmmxxMinMaxLinearEquationSolver.h | 8 +- .../GmmxxMinMaxLinearEquationSolverTest.cpp | 20 ++ 5 files changed, 218 insertions(+), 66 deletions(-) diff --git a/src/settings/modules/GeneralSettings.cpp b/src/settings/modules/GeneralSettings.cpp index bd2e03a47..807823f67 100644 --- a/src/settings/modules/GeneralSettings.cpp +++ b/src/settings/modules/GeneralSettings.cpp @@ -46,6 +46,7 @@ namespace storm { const std::string GeneralSettings::cudaOptionName = "cuda"; const std::string GeneralSettings::prismCompatibilityOptionName = "prismcompat"; const std::string GeneralSettings::prismCompatibilityOptionShortName = "pc"; + const std::string GeneralSettings::minMaxEquationSolvingTechniqueOptionName = "minMaxEquationSolvingTechnique"; #ifdef STORM_HAVE_CARL const std::string GeneralSettings::parametricOptionName = "parametric"; @@ -100,6 +101,10 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, statisticsOptionName, false, "Sets whether to display statistics if available.").setShortName(statisticsOptionShortName).build()); this->addOption(storm::settings::OptionBuilder(moduleName, cudaOptionName, false, "Sets whether to use CUDA to speed up computation time.").build()); + std::vector minMaxSolvingTechniques = {"policyIteration", "valueIteration"}; + this->addOption(storm::settings::OptionBuilder(moduleName, minMaxEquationSolvingTechniqueOptionName, false, "Sets which min/max linear equation solving technique is preferred.") + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a min/max linear equation solving technique. Available are: valueIteration and policyIteration.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(minMaxSolvingTechniques)).setDefaultValueString("valueIteration").build()).build()); + #ifdef STORM_HAVE_CARL this->addOption(storm::settings::OptionBuilder(moduleName, parametricOptionName, false, "Sets whether to use the parametric engine.").build()); #endif @@ -285,6 +290,20 @@ namespace storm { bool GeneralSettings::isPrismCompatibilityEnabled() const { return this->getOption(prismCompatibilityOptionName).getHasOptionBeenSet(); } + + GeneralSettings::MinMaxTechnique GeneralSettings::getMinMaxEquationSolvingTechnique() const { + std::string minMaxEquationSolvingTechnique = this->getOption(minMaxEquationSolvingTechniqueOptionName).getArgumentByName("name").getValueAsString(); + if (minMaxEquationSolvingTechnique == "valueIteration") { + return GeneralSettings::MinMaxTechnique::ValueIteration; + } else if (minMaxEquationSolvingTechnique == "policyIteration") { + return GeneralSettings::MinMaxTechnique::PolicyIteration; + } + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'."); + } + + bool GeneralSettings::isMinMaxEquationSolvingTechniqueSet() const { + return this->getOption(minMaxEquationSolvingTechniqueOptionName).getHasOptionBeenSet(); + } #ifdef STORM_HAVE_CARL bool GeneralSettings::isParametricSet() const { diff --git a/src/settings/modules/GeneralSettings.h b/src/settings/modules/GeneralSettings.h index e4528f64f..f2454a28f 100644 --- a/src/settings/modules/GeneralSettings.h +++ b/src/settings/modules/GeneralSettings.h @@ -21,6 +21,9 @@ namespace storm { // An enumeration of all available equation solvers. enum class EquationSolver { Gmmxx, Native }; + + // An enumeration of all available Min/Max equation solver techniques. + enum class MinMaxTechnique { ValueIteration, PolicyIteration }; /*! * Creates a new set of general settings that is managed by the given manager. @@ -316,7 +319,7 @@ namespace storm { /*! * Retrieves the selected engine. * - * @return The selecte engine. + * @return The selected engine. */ Engine getEngine() const; @@ -335,6 +338,20 @@ namespace storm { */ bool isParametricSet() const; #endif + /*! + * Retrieves whether a min/max equation solving technique has been set. + * + * @return True iff an equation solving technique has been set. + */ + bool isMinMaxEquationSolvingTechniqueSet() const; + + /*! + * Retrieves the selected min/max equation solving technique. + * + * @return The selected min/max equation solving technique. + */ + MinMaxTechnique getMinMaxEquationSolvingTechnique() const; + bool check() const override; @@ -380,7 +397,8 @@ namespace storm { static const std::string cudaOptionName; static const std::string prismCompatibilityOptionName; static const std::string prismCompatibilityOptionShortName; - + static const std::string minMaxEquationSolvingTechniqueOptionName; + #ifdef STORM_HAVE_CARL static const std::string parametricOptionName; #endif diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp index 3b349d2ff..e19a37548 100644 --- a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp +++ b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp @@ -4,91 +4,180 @@ #include "src/settings/SettingsManager.h" #include "src/adapters/GmmxxAdapter.h" +#include "src/solver/GmmxxLinearEquationSolver.h" #include "src/utility/vector.h" namespace storm { namespace solver { template - GmmxxMinMaxLinearEquationSolver::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A)), rowGroupIndices(A.getRowGroupIndices()) { + GmmxxMinMaxLinearEquationSolver::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A)), stormMatrix(A), rowGroupIndices(A.getRowGroupIndices()) { // Get the settings object to customize solving. storm::settings::modules::GmmxxEquationSolverSettings const& settings = storm::settings::gmmxxEquationSolverSettings(); + storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings(); // Get appropriate settings. maximalNumberOfIterations = settings.getMaximalIterationCount(); precision = settings.getPrecision(); relative = settings.getConvergenceCriterion() == storm::settings::modules::GmmxxEquationSolverSettings::ConvergenceCriterion::Relative; + useValueIteration = (generalSettings.getMinMaxEquationSolvingTechnique() == storm::settings::modules::GeneralSettings::MinMaxTechnique::ValueIteration); } template - GmmxxMinMaxLinearEquationSolver::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A)), rowGroupIndices(A.getRowGroupIndices()), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) { + GmmxxMinMaxLinearEquationSolver::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A)), stormMatrix(A), rowGroupIndices(A.getRowGroupIndices()), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), useValueIteration(useValueIteration), relative(relative) { // Intentionally left empty. } template void GmmxxMinMaxLinearEquationSolver::solveEquationSystem(bool minimize, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { - // Set up the environment for the power method. If scratch memory was not provided, we need to create it. - bool multiplyResultMemoryProvided = true; - if (multiplyResult == nullptr) { - multiplyResult = new std::vector(b.size()); - multiplyResultMemoryProvided = false; - } - - std::vector* currentX = &x; - bool xMemoryProvided = true; - if (newX == nullptr) { - newX = new std::vector(x.size()); - xMemoryProvided = false; - } - uint_fast64_t iterations = 0; - bool converged = false; - - // Keep track of which of the vectors for x is the auxiliary copy. - std::vector* copyX = newX; - - // Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number - // of iterations. - while (!converged && iterations < maximalNumberOfIterations) { - // Compute x' = A*x + b. - gmm::mult(*gmmxxMatrix, *currentX, *multiplyResult); - gmm::add(b, *multiplyResult); - - // Reduce the vector x by applying min/max over all nondeterministic choices. - if (minimize) { - storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, rowGroupIndices); - } else { - storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, rowGroupIndices); - } - - // Determine whether the method converged. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative); - - // Update environment variables. - std::swap(currentX, newX); - ++iterations; - } - - // Check if the solver converged and issue a warning otherwise. - if (converged) { - LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); - } else { - LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations."); - } - - // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result - // is currently stored in currentX, but x is the output vector. - if (currentX == copyX) { - std::swap(x, *currentX); - } - - if (!xMemoryProvided) { - delete copyX; - } - - if (!multiplyResultMemoryProvided) { - delete multiplyResult; - } + if (useValueIteration) { + // Set up the environment for the power method. If scratch memory was not provided, we need to create it. + bool multiplyResultMemoryProvided = true; + if (multiplyResult == nullptr) { + multiplyResult = new std::vector(b.size()); + multiplyResultMemoryProvided = false; + } + + std::vector* currentX = &x; + bool xMemoryProvided = true; + if (newX == nullptr) { + newX = new std::vector(x.size()); + xMemoryProvided = false; + } + uint_fast64_t iterations = 0; + bool converged = false; + + // Keep track of which of the vectors for x is the auxiliary copy. + std::vector* copyX = newX; + + // Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number + // of iterations. + while (!converged && iterations < maximalNumberOfIterations) { + // Compute x' = A*x + b. + gmm::mult(*gmmxxMatrix, *currentX, *multiplyResult); + gmm::add(b, *multiplyResult); + + // Reduce the vector x by applying min/max over all nondeterministic choices. + if (minimize) { + storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, rowGroupIndices); + } else { + storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, rowGroupIndices); + } + + // Determine whether the method converged. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative); + + // Update environment variables. + std::swap(currentX, newX); + ++iterations; + } + + // Check if the solver converged and issue a warning otherwise. + if (converged) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations."); + } + + // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result + // is currently stored in currentX, but x is the output vector. + if (currentX == copyX) { + std::swap(x, *currentX); + } + + if (!xMemoryProvided) { + delete copyX; + } + + if (!multiplyResultMemoryProvided) { + delete multiplyResult; + } + } else { + // We will use Policy Iteration to solve the given system. + + // We first define an initial choice resolution which will be refined after each iteration. + std::vector::index_type> choiceVector(rowGroupIndices.size() - 1); + + // Create our own multiplyResult for solving the deterministic instances. + std::vector deterministicMultiplyResult(rowGroupIndices.size() - 1); + std::vector subB(rowGroupIndices.size() - 1); + + bool multiplyResultMemoryProvided = true; + if (multiplyResult == nullptr) { + multiplyResult = new std::vector(b.size()); + multiplyResultMemoryProvided = false; + } + + std::vector* currentX = &x; + bool xMemoryProvided = true; + if (newX == nullptr) { + newX = new std::vector(x.size()); + xMemoryProvided = false; + } + + uint_fast64_t iterations = 0; + bool converged = false; + + // Keep track of which of the vectors for x is the auxiliary copy. + std::vector* copyX = newX; + + // Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number + // of iterations. + while (!converged && iterations < maximalNumberOfIterations) { + // Take the sub-matrix according to the current choices + storm::storage::SparseMatrix submatrix = stormMatrix.selectRowsFromRowGroups(choiceVector, true); + submatrix.convertToEquationSystem(); + + GmmxxLinearEquationSolver gmmLinearEquationSolver(submatrix, storm::solver::GmmxxLinearEquationSolver::SolutionMethod::Gmres, precision, maximalNumberOfIterations, storm::solver::GmmxxLinearEquationSolver::Preconditioner::None, relative, 50); + storm::utility::vector::selectVectorValues(subB, choiceVector, rowGroupIndices, b); + + // Copy X since we will overwrite it + std::copy(currentX->begin(), currentX->end(), newX->begin()); + + // Solve the resulting linear equation system + gmmLinearEquationSolver.solveEquationSystem(*newX, subB, &deterministicMultiplyResult); + + // Compute x' = A*x + b. This step is necessary to allow the choosing of the optimal policy for the next iteration. + gmm::mult(*gmmxxMatrix, *newX, *multiplyResult); + gmm::add(b, *multiplyResult); + + // Reduce the vector x by applying min/max over all nondeterministic choices. + if (minimize) { + storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, rowGroupIndices, &choiceVector); + } else { + storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, rowGroupIndices, &choiceVector); + } + + // Determine whether the method converged. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative); + + // Update environment variables. + std::swap(currentX, newX); + ++iterations; + } + + // Check if the solver converged and issue a warning otherwise. + if (converged) { + LOG4CPLUS_INFO(logger, "Iterative solver converged after " << iterations << " iterations."); + } else { + LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations."); + } + + // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result + // is currently stored in currentX, but x is the output vector. + if (currentX == copyX) { + std::swap(x, *currentX); + } + + if (!xMemoryProvided) { + delete copyX; + } + + if (!multiplyResultMemoryProvided) { + delete multiplyResult; + } + } } template diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.h b/src/solver/GmmxxMinMaxLinearEquationSolver.h index 951ae2578..a6bdb8281 100644 --- a/src/solver/GmmxxMinMaxLinearEquationSolver.h +++ b/src/solver/GmmxxMinMaxLinearEquationSolver.h @@ -31,7 +31,7 @@ namespace storm { * @param relative If set, the relative error rather than the absolute error is considered for convergence * detection. */ - GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true); + GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative = true); virtual void performMatrixVectorMultiplication(bool minimize, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; @@ -40,6 +40,9 @@ namespace storm { private: // The (gmm++) matrix associated with this equation solver. std::unique_ptr> gmmxxMatrix; + + // A reference to the original sparse matrix. + storm::storage::SparseMatrix const& stormMatrix; // A reference to the row group indices of the original matrix. std::vector const& rowGroupIndices; @@ -52,6 +55,9 @@ namespace storm { // The maximal number of iterations to do before iteration is aborted. uint_fast64_t maximalNumberOfIterations; + + // Whether value iteration or policy iteration is to be used. + bool useValueIteration; }; } // namespace solver } // namespace storm diff --git a/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp index 56357f546..bd402d915 100644 --- a/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp +++ b/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp @@ -62,4 +62,24 @@ TEST(GmmxxMinMaxLinearEquationSolver, MatrixVectorMultiplication) { x = {0, 1, 0}; ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, x, nullptr, 20)); ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::gmmxxEquationSolverSettings().getPrecision()); +} + +TEST(GmmxxMinMaxLinearEquationSolver, SolveWithPolicyIteration) { + storm::storage::SparseMatrixBuilder builder(0, 0, 0, false, true); + ASSERT_NO_THROW(builder.newRowGroup(0)); + ASSERT_NO_THROW(builder.addNextValue(0, 0, 0.9)); + + storm::storage::SparseMatrix A; + ASSERT_NO_THROW(A = builder.build(2)); + + std::vector x(1); + std::vector b = { 0.099, 0.5 }; + + storm::settings::modules::GmmxxEquationSolverSettings const& settings = storm::settings::gmmxxEquationSolverSettings(); + storm::solver::GmmxxMinMaxLinearEquationSolver solver(A, settings.getPrecision(), settings.getMaximalIterationCount(), false, settings.getConvergenceCriterion() == storm::settings::modules::GmmxxEquationSolverSettings::ConvergenceCriterion::Relative); + ASSERT_NO_THROW(solver.solveEquationSystem(true, x, b)); + ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::gmmxxEquationSolverSettings().getPrecision()); + + ASSERT_NO_THROW(solver.solveEquationSystem(false, x, b)); + ASSERT_LT(std::abs(x[0] - 0.99), storm::settings::gmmxxEquationSolverSettings().getPrecision()); } \ No newline at end of file