From 9ab33528b40f7df0955a5c0aea373ad6867f0259 Mon Sep 17 00:00:00 2001 From: dehnert Date: Sun, 26 Jun 2016 23:10:35 +0200 Subject: [PATCH] started to fill value iteration implementation in new general min-max solver Former-commit-id: e54cb8a0f900f1ee48a13b3f6e7869385b6cc29b --- .../prctl/helper/SparseMdpPrctlHelper.cpp | 1 - src/settings/SettingsManager.cpp | 6 +- src/settings/modules/CoreSettings.cpp | 14 -- src/settings/modules/CoreSettings.h | 21 +-- .../modules/MinMaxEquationSolverSettings.cpp | 72 ++++++++++ .../modules/MinMaxEquationSolverSettings.h | 91 ++++++++++++ src/solver/EigenLinearEquationSolver.cpp | 31 +++- src/solver/EigenLinearEquationSolver.h | 4 +- .../EliminationLinearEquationSolver.cpp | 35 +++-- src/solver/EliminationLinearEquationSolver.h | 2 +- src/solver/GmmxxLinearEquationSolver.cpp | 8 ++ src/solver/GmmxxLinearEquationSolver.h | 4 +- src/solver/LinearEquationSolver.h | 14 +- src/solver/MinMaxLinearEquationSolver.cpp | 8 +- src/solver/NativeLinearEquationSolver.cpp | 17 +++ src/solver/NativeLinearEquationSolver.h | 2 +- src/solver/SolverSelectionOptions.cpp | 8 +- src/solver/SolverSelectionOptions.h | 2 +- .../StandardMinMaxLinearEquationSolver.cpp | 133 ++++++++++++++++-- .../StandardMinMaxLinearEquationSolver.h | 20 ++- .../solver/MinMaxTechniqueSelectionTest.cpp | 6 +- 21 files changed, 422 insertions(+), 77 deletions(-) create mode 100644 src/settings/modules/MinMaxEquationSolverSettings.cpp create mode 100644 src/settings/modules/MinMaxEquationSolverSettings.h diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index d2613fca8..7879ff845 100644 --- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -39,7 +39,6 @@ namespace storm { maybeStates &= ~psiStates; STORM_LOG_INFO("Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states."); - if (!maybeStates.empty()) { // We can eliminate the rows and columns from the original transition probability matrix that have probability 0. storm::storage::SparseMatrix submatrix = transitionMatrix.getSubmatrix(true, maybeStates, maybeStates, false); diff --git a/src/settings/SettingsManager.cpp b/src/settings/SettingsManager.cpp index c88f7017f..3cd1bbd2d 100644 --- a/src/settings/SettingsManager.cpp +++ b/src/settings/SettingsManager.cpp @@ -23,11 +23,12 @@ #include "src/settings/modules/EigenEquationSolverSettings.h" #include "src/settings/modules/GmmxxEquationSolverSettings.h" #include "src/settings/modules/NativeEquationSolverSettings.h" +#include "src/settings/modules/EliminationSettings.h" +#include "src/settings/modules/MinMaxEquationSolverSettings.h" #include "src/settings/modules/BisimulationSettings.h" #include "src/settings/modules/GlpkSettings.h" #include "src/settings/modules/GurobiSettings.h" #include "src/settings/modules/ParametricSettings.h" -#include "src/settings/modules/EliminationSettings.h" #include "src/settings/modules/TopologicalValueIterationEquationSolverSettings.h" #include "src/settings/modules/ExplorationSettings.h" #include "src/utility/macros.h" @@ -511,12 +512,13 @@ namespace storm { storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); + storm::settings::addModule(); + storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); storm::settings::addModule(); - storm::settings::addModule(); storm::settings::addModule(); } diff --git a/src/settings/modules/CoreSettings.cpp b/src/settings/modules/CoreSettings.cpp index 105d67d7f..5617d959c 100644 --- a/src/settings/modules/CoreSettings.cpp +++ b/src/settings/modules/CoreSettings.cpp @@ -142,20 +142,6 @@ namespace storm { this->engine = newEngine; } - storm::solver::MinMaxTechnique CoreSettings::getMinMaxEquationSolvingTechnique() const { - std::string minMaxEquationSolvingTechnique = this->getOption(minMaxEquationSolvingTechniqueOptionName).getArgumentByName("name").getValueAsString(); - if (minMaxEquationSolvingTechnique == "value-iteration" || minMaxEquationSolvingTechnique == "vi") { - return storm::solver::MinMaxTechnique::ValueIteration; - } else if (minMaxEquationSolvingTechnique == "policy-iteration" || minMaxEquationSolvingTechnique == "pi") { - return storm::solver::MinMaxTechnique::PolicyIteration; - } - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'."); - } - - bool CoreSettings::isMinMaxEquationSolvingTechniqueSet() const { - return this->getOption(minMaxEquationSolvingTechniqueOptionName).getHasOptionBeenSet(); - } - void CoreSettings::finalize() { // Finalize engine. std::string engineStr = this->getOption(engineOptionName).getArgumentByName("name").getValueAsString(); diff --git a/src/settings/modules/CoreSettings.h b/src/settings/modules/CoreSettings.h index 204081dd2..9b7189e2c 100644 --- a/src/settings/modules/CoreSettings.h +++ b/src/settings/modules/CoreSettings.h @@ -10,7 +10,7 @@ namespace storm { namespace solver { enum class EquationSolverType; enum class LpSolverType; - enum class MinMaxTechnique; + enum class MinMaxMethod; enum class SmtSolverType; } @@ -22,7 +22,7 @@ namespace storm { namespace modules { /*! - * This class represents the markov chain settings. + * This class represents the core settings. */ class CoreSettings : public ModuleSettings { public: @@ -32,7 +32,7 @@ namespace storm { }; /*! - * Creates a new set of markov chain settings. + * Creates a new set of core settings. */ CoreSettings(); @@ -128,21 +128,6 @@ namespace storm { */ void setEngine(Engine); - /*! - * 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. - */ - storm::solver::MinMaxTechnique getMinMaxEquationSolvingTechnique() const; - - bool check() const override; void finalize() override; diff --git a/src/settings/modules/MinMaxEquationSolverSettings.cpp b/src/settings/modules/MinMaxEquationSolverSettings.cpp new file mode 100644 index 000000000..8a7067094 --- /dev/null +++ b/src/settings/modules/MinMaxEquationSolverSettings.cpp @@ -0,0 +1,72 @@ +#include "src/settings/modules/MinMaxEquationSolverSettings.h" + +#include "src/settings/Option.h" +#include "src/settings/ArgumentBuilder.h" +#include "src/settings/OptionBuilder.h" + +#include "src/utility/macros.h" +#include "src/exceptions/IllegalArgumentValueException.h" + +namespace storm { + namespace settings { + namespace modules { + + const std::string MinMaxEquationSolverSettings::moduleName = "minmax"; + const std::string MinMaxEquationSolverSettings::solvingMethodOptionName = "method"; + const std::string MinMaxEquationSolverSettings::maximalIterationsOptionName = "maxiter"; + const std::string MinMaxEquationSolverSettings::maximalIterationsOptionShortName = "i"; + const std::string MinMaxEquationSolverSettings::precisionOptionName = "precision"; + + MinMaxEquationSolverSettings::MinMaxEquationSolverSettings() : ModuleSettings(moduleName) { + std::vector minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration"}; + this->addOption(storm::settings::OptionBuilder(moduleName, solvingMethodOptionName, 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: value-iteration (vi) and policy-iteration (pi).").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(minMaxSolvingTechniques)).setDefaultValueString("vi").build()).build()); + + this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build()); + + this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + + this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build()); + } + + storm::solver::MinMaxMethod MinMaxEquationSolverSettings::getMinMaxEquationSolvingMethod() const { + std::string minMaxEquationSolvingTechnique = this->getOption(solvingMethodOptionName).getArgumentByName("name").getValueAsString(); + if (minMaxEquationSolvingTechnique == "value-iteration" || minMaxEquationSolvingTechnique == "vi") { + return storm::solver::MinMaxMethod::ValueIteration; + } else if (minMaxEquationSolvingTechnique == "policy-iteration" || minMaxEquationSolvingTechnique == "pi") { + return storm::solver::MinMaxMethod::PolicyIteration; + } + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'."); + } + + bool MinMaxEquationSolverSettings::isMinMaxEquationSolvingMethodSet() const { + return this->getOption(solvingMethodOptionName).getHasOptionBeenSet(); + } + + bool MinMaxEquationSolverSettings::isMaximalIterationCountSet() const { + return this->getOption(maximalIterationsOptionName).getHasOptionBeenSet(); + } + + uint_fast64_t MinMaxEquationSolverSettings::getMaximalIterationCount() const { + return this->getOption(maximalIterationsOptionName).getArgumentByName("count").getValueAsUnsignedInteger(); + } + + bool MinMaxEquationSolverSettings::isPrecisionSet() const { + return this->getOption(precisionOptionName).getHasOptionBeenSet(); + } + + double MinMaxEquationSolverSettings::getPrecision() const { + return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble(); + } + + bool MinMaxEquationSolverSettings::isConvergenceCriterionSet() const { + return this->getOption(absoluteOptionName).getHasOptionBeenSet(); + } + + MinMaxEquationSolverSettings::ConvergenceCriterion MinMaxEquationSolverSettings::getConvergenceCriterion() const { + return this->getOption(absoluteOptionName).getHasOptionBeenSet() ? MinMaxEquationSolverSettings::ConvergenceCriterion::Absolute : MinMaxEquationSolverSettings::ConvergenceCriterion::Relative; + } + + } + } +} \ No newline at end of file diff --git a/src/settings/modules/MinMaxEquationSolverSettings.h b/src/settings/modules/MinMaxEquationSolverSettings.h new file mode 100644 index 000000000..33e65a8ca --- /dev/null +++ b/src/settings/modules/MinMaxEquationSolverSettings.h @@ -0,0 +1,91 @@ +#pragma once + +#include "storm-config.h" +#include "src/settings/modules/ModuleSettings.h" + +#include "src/solver/SolverSelectionOptions.h" + +namespace storm { + namespace settings { + namespace modules { + + /*! + * This class represents the min/max solver settings. + */ + class MinMaxEquationSolverSettings : public ModuleSettings { + public: + // An enumeration of all available convergence criteria. + enum class ConvergenceCriterion { Absolute, Relative }; + + MinMaxEquationSolverSettings(); + + /*! + * Retrieves whether a min/max equation solving technique has been set. + * + * @return True iff an equation solving technique has been set. + */ + bool isMinMaxEquationSolvingMethodSet() const; + + /*! + * Retrieves the selected min/max equation solving technique. + * + * @return The selected min/max equation solving technique. + */ + storm::solver::MinMaxMethod getMinMaxEquationSolvingMethod() const; + + /*! + * Retrieves whether the maximal iteration count has been set. + * + * @return True iff the maximal iteration count has been set. + */ + bool isMaximalIterationCountSet() const; + + /*! + * Retrieves the maximal number of iterations to perform until giving up on converging. + * + * @return The maximal iteration count. + */ + uint_fast64_t getMaximalIterationCount() const; + + /*! + * Retrieves whether the precision has been set. + * + * @return True iff the precision has been set. + */ + bool isPrecisionSet() const; + + /*! + * Retrieves the precision that is used for detecting convergence. + * + * @return The precision to use for detecting convergence. + */ + double getPrecision() const; + + /*! + * Retrieves whether the convergence criterion has been set. + * + * @return True iff the convergence criterion has been set. + */ + bool isConvergenceCriterionSet() const; + + /*! + * Retrieves the selected convergence criterion. + * + * @return The selected convergence criterion. + */ + ConvergenceCriterion getConvergenceCriterion() const; + + // The name of the module. + static const std::string moduleName; + + private: + static const std::string solvingMethodOptionName; + static const std::string maximalIterationsOptionName; + static const std::string maximalIterationsOptionShortName; + static const std::string precisionOptionName; + static const std::string absoluteOptionName; + }; + + } + } +} \ No newline at end of file diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp index 91fc9cad4..cac22e645 100644 --- a/src/solver/EigenLinearEquationSolver.cpp +++ b/src/solver/EigenLinearEquationSolver.cpp @@ -196,6 +196,34 @@ namespace storm { } } + template + void EigenLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b) const { + // Typedef the map-type so we don't have to spell it out. + typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; + + auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); + auto eigenResult = Eigen::Matrix::Map(result.data(), result.size()); + + std::unique_ptr eigenB; + if (b != nullptr) { + eigenB = std::make_unique(Eigen::Matrix::Map(b->data(), b->size())); + } + + if (&x != &result) { + if (b != nullptr) { + eigenResult.noalias() = *eigenA * eigenX + *eigenB; + } else { + eigenResult.noalias() = *eigenA * eigenX; + } + } else { + if (b != nullptr) { + eigenResult = *eigenA * eigenX + *eigenB; + } else { + eigenResult = *eigenA * 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. @@ -223,7 +251,6 @@ namespace storm { } 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. @@ -286,6 +313,7 @@ namespace storm { } 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. @@ -338,6 +366,7 @@ namespace storm { } 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. diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h index 68e2969ab..7dbd5d519 100644 --- a/src/solver/EigenLinearEquationSolver.h +++ b/src/solver/EigenLinearEquationSolver.h @@ -62,9 +62,9 @@ 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& result, std::vector const* b = nullptr) const override; + EigenLinearEquationSolverSettings& getSettings(); EigenLinearEquationSolverSettings const& getSettings() const; diff --git a/src/solver/EliminationLinearEquationSolver.cpp b/src/solver/EliminationLinearEquationSolver.cpp index abc6dc173..36c1144a3 100644 --- a/src/solver/EliminationLinearEquationSolver.cpp +++ b/src/solver/EliminationLinearEquationSolver.cpp @@ -32,17 +32,17 @@ namespace storm { storm::settings::modules::EliminationSettings::EliminationOrder EliminationLinearEquationSolverSettings::getEliminationOrder() const { return order; } - + template EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A, EliminationLinearEquationSolverSettings const& settings) : localA(nullptr), A(A), settings(settings) { // Intentionally left empty. } - + template EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings) : localA(std::make_unique>(std::move(A))), A(*localA), settings(settings) { // Intentionally left empty. } - + template void EliminationLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver."); @@ -65,7 +65,7 @@ namespace storm { storm::storage::SparseMatrix const& transitionMatrix = localA ? *localA : locallyConvertedMatrix; storm::storage::SparseMatrix backwardTransitions = transitionMatrix.transpose(); - + // Initialize the solution to the right-hand side of the equation system. x = b; @@ -136,7 +136,24 @@ namespace storm { delete copyX; } } - + + template + void EliminationLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b) const { + if (&x != &result) { + A.multiplyWithVector(x, result); + if (b != nullptr) { + storm::utility::vector::addVectors(result, *b, result); + } + } else { + // If the two vectors are aliases, we need to create a temporary. + std::vector tmp(result.size()); + A.multiplyWithVector(x, tmp); + if (b != nullptr) { + storm::utility::vector::addVectors(tmp, *b, result); + } + } + } + template EliminationLinearEquationSolverSettings& EliminationLinearEquationSolver::getSettings() { return settings; @@ -161,12 +178,12 @@ namespace storm { EliminationLinearEquationSolverSettings& EliminationLinearEquationSolverFactory::getSettings() { return settings; } - + template EliminationLinearEquationSolverSettings const& EliminationLinearEquationSolverFactory::getSettings() const { return settings; } - + template std::unique_ptr> EliminationLinearEquationSolverFactory::clone() const { return std::make_unique>(*this); @@ -175,11 +192,11 @@ namespace storm { template class EliminationLinearEquationSolverSettings; template class EliminationLinearEquationSolverSettings; template class EliminationLinearEquationSolverSettings; - + template class EliminationLinearEquationSolver; template class EliminationLinearEquationSolver; template class EliminationLinearEquationSolver; - + template class EliminationLinearEquationSolverFactory; template class EliminationLinearEquationSolverFactory; template class EliminationLinearEquationSolverFactory; diff --git a/src/solver/EliminationLinearEquationSolver.h b/src/solver/EliminationLinearEquationSolver.h index 5cd5d7ea2..43d93c42e 100644 --- a/src/solver/EliminationLinearEquationSolver.h +++ b/src/solver/EliminationLinearEquationSolver.h @@ -30,8 +30,8 @@ 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(); EliminationLinearEquationSolverSettings const& getSettings() const; diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index a043c03c7..d9d331f9f 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -214,6 +214,14 @@ namespace storm { } } + 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); + } + } + template uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult) const { // Get a Jacobi decomposition of the matrix A. diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 423a10e71..ca14791f2 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -90,9 +90,9 @@ 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(); GmmxxLinearEquationSolverSettings const& getSettings() const; diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h index 342020eba..8585f2136 100644 --- a/src/solver/LinearEquationSolver.h +++ b/src/solver/LinearEquationSolver.h @@ -40,13 +40,25 @@ namespace storm { * matrix A has to be given upon construction time of the solver object. * * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal - * to the number of rows of A. + * to the number of columns of A. * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal * to the number of rows of A. * @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 = 0; + + /*! + * Performs on matrix-vector multiplication x' = A*x + b. + * + * @param x The input vector with which to multiply the matrix. Its length must be equal + * to the number of columns of A. + * @param result The target vector into which to write the multiplication result. Its length must be equal + * to the number of rows of A. + * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal + * to the number of rows of A. + */ + virtual void performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b = nullptr) const = 0; }; template diff --git a/src/solver/MinMaxLinearEquationSolver.cpp b/src/solver/MinMaxLinearEquationSolver.cpp index 887b60a7b..22cefe139 100644 --- a/src/solver/MinMaxLinearEquationSolver.cpp +++ b/src/solver/MinMaxLinearEquationSolver.cpp @@ -7,7 +7,7 @@ #include "src/solver/TopologicalMinMaxLinearEquationSolver.h" #include "src/settings/SettingsManager.h" -#include "src/settings/modules/CoreSettings.h" +#include "src/settings/modules/MinMaxEquationSolverSettings.h" #include "src/utility/macros.h" #include "src/exceptions/NotImplementedException.h" @@ -118,10 +118,10 @@ namespace storm { template std::unique_ptr> GeneralMinMaxLinearEquationSolverFactory::selectSolver(MatrixType&& matrix) const { std::unique_ptr> result; - auto technique = storm::settings::getModule().getMinMaxEquationSolvingTechnique(); - if (technique == MinMaxTechnique::ValueIteration || technique == MinMaxTechnique::PolicyIteration) { + auto method = storm::settings::getModule().getMinMaxEquationSolvingMethod(); + if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration) { result = std::make_unique>(std::forward(matrix), std::make_unique>()); - } else if (technique == MinMaxTechnique::Topological) { + } else if (method == MinMaxMethod::Topological) { result = std::make_unique>(std::forward(matrix)); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index 16de87bb1..5c73d4b5e 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -219,6 +219,23 @@ namespace storm { } } + template + void NativeLinearEquationSolver::performMatrixVectorMultiplication(std::vector& x, std::vector& result, std::vector const* b) const { + if (&x != &result) { + A.multiplyWithVector(x, result); + if (b != nullptr) { + storm::utility::vector::addVectors(result, *b, result); + } + } else { + // If the two vectors are aliases, we need to create a temporary. + std::vector tmp(result.size()); + A.multiplyWithVector(x, tmp); + if (b != nullptr) { + storm::utility::vector::addVectors(tmp, *b, result); + } + } + } + template NativeLinearEquationSolverSettings& NativeLinearEquationSolver::getSettings() { return settings; diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index 9c315d8ca..c76a45b4e 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -47,8 +47,8 @@ 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(); NativeLinearEquationSolverSettings const& getSettings() const; diff --git a/src/solver/SolverSelectionOptions.cpp b/src/solver/SolverSelectionOptions.cpp index f50803761..0578d2d0b 100644 --- a/src/solver/SolverSelectionOptions.cpp +++ b/src/solver/SolverSelectionOptions.cpp @@ -2,13 +2,13 @@ namespace storm { namespace solver { - std::string toString(MinMaxTechnique m) { + std::string toString(MinMaxMethod m) { switch(m) { - case MinMaxTechnique::PolicyIteration: + case MinMaxMethod::PolicyIteration: return "policy"; - case MinMaxTechnique::ValueIteration: + case MinMaxMethod::ValueIteration: return "value"; - case MinMaxTechnique::Topological: + case MinMaxMethod::Topological: return "topological"; } diff --git a/src/solver/SolverSelectionOptions.h b/src/solver/SolverSelectionOptions.h index 2cffc4ae0..eb1448630 100644 --- a/src/solver/SolverSelectionOptions.h +++ b/src/solver/SolverSelectionOptions.h @@ -6,7 +6,7 @@ namespace storm { namespace solver { - ExtendEnumsWithSelectionField(MinMaxTechnique, PolicyIteration, ValueIteration, Topological) + ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, Topological) ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk) ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination) diff --git a/src/solver/StandardMinMaxLinearEquationSolver.cpp b/src/solver/StandardMinMaxLinearEquationSolver.cpp index aec71080e..d6213ddd7 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp @@ -1,13 +1,14 @@ #include "src/solver/StandardMinMaxLinearEquationSolver.h" #include "src/settings/SettingsManager.h" -#include "src/settings/modules/CoreSettings.h" +#include "src/settings/modules/MinMaxEquationSolverSettings.h" #include "src/solver/GmmxxLinearEquationSolver.h" #include "src/solver/EigenLinearEquationSolver.h" #include "src/solver/NativeLinearEquationSolver.h" #include "src/solver/EliminationLinearEquationSolver.h" +#include "src/utility/vector.h" #include "src/utility/macros.h" #include "src/exceptions/InvalidSettingsException.h" @@ -15,10 +16,15 @@ namespace storm { namespace solver { StandardMinMaxLinearEquationSolverSettings::StandardMinMaxLinearEquationSolverSettings() { - auto technique = storm::settings::getModule().getMinMaxEquationSolvingTechnique(); - switch (technique) { - case MinMaxTechnique::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break; - case MinMaxTechnique::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break; + // Get the settings object to customize linear solving. + storm::settings::modules::MinMaxEquationSolverSettings const& settings = storm::settings::getModule(); + + maximalNumberOfIterations = settings.getMaximalIterationCount(); + + auto method = settings.getMinMaxEquationSolvingMethod(); + switch (method) { + case MinMaxMethod::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break; + case MinMaxMethod::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break; default: STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); } @@ -28,10 +34,34 @@ namespace storm { this->solutionMethod = solutionMethod; } + void StandardMinMaxLinearEquationSolverSettings::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) { + this->maximalNumberOfIterations = maximalNumberOfIterations; + } + + void StandardMinMaxLinearEquationSolverSettings::setRelativeTerminationCriterion(bool value) { + this->relative = value; + } + + void StandardMinMaxLinearEquationSolverSettings::setPrecision(double precision) { + this->precision = precision; + } + StandardMinMaxLinearEquationSolverSettings::SolutionMethod const& StandardMinMaxLinearEquationSolverSettings::getSolutionMethod() const { return solutionMethod; } + uint64_t StandardMinMaxLinearEquationSolverSettings::getMaximalNumberOfIterations() const { + return maximalNumberOfIterations; + } + + double StandardMinMaxLinearEquationSolverSettings::getPrecision() const { + return precision; + } + + 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) { // Intentionally left empty. @@ -43,13 +73,98 @@ namespace storm { } template - void StandardMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection d, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { - //FIXME: implement + 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; + } } template - void StandardMinMaxLinearEquationSolver::performMatrixVectorMultiplication(OptimizationDirection d, std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { - //FIXME: implement + void StandardMinMaxLinearEquationSolver::solveEquationSystemPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + // FIXME. + } + + template + void StandardMinMaxLinearEquationSolver::solveEquationSystemValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + // Create scratch memory if none was provided. + bool multiplyResultMemoryProvided = multiplyResult != nullptr; + if (multiplyResult == nullptr) { + multiplyResult = new std::vector(this->A.getRowCount()); + } + std::vector* currentX = &x; + bool xMemoryProvided = newX != nullptr; + if (newX == nullptr) { + newX = new std::vector(x.size()); + } + + // Keep track of which of the vectors for x is the auxiliary copy. + std::vector* copyX = newX; + + uint64_t iterations = 0; + bool converged = false; + + // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. + while (!converged && iterations < this->getSettings().getMaximalNumberOfIterations() && (!this->hasCustomTerminationCondition() || this->getTerminationCondition().terminateNow(*currentX))) { + // Compute x' = A*x + b. + this->A.multiplyWithVector(*currentX, *multiplyResult); + storm::utility::vector::addVectors(*multiplyResult, b, *multiplyResult); + + // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost + // element of the min/max operator stack. + storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, *newX, this->A.getRowGroupIndices()); + + // Determine whether the method converged. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()); + + // Update environment variables. + std::swap(currentX, newX); + ++iterations; + } + + // Check if the solver converged and issue a warning otherwise. + if (converged) { + STORM_LOG_INFO("Iterative solver converged after " << iterations << " iterations."); + } else { + STORM_LOG_WARN("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 + void StandardMinMaxLinearEquationSolver::performMatrixVectorMultiplication(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { + std::unique_ptr> solver = linearEquationSolverFactory->create(A); + + // If scratch memory was not provided, we create it. + bool multiplyResultMemoryProvided = multiplyResult != nullptr; + if (!multiplyResult) { + multiplyResult = new std::vector(this->A.getRowCount()); + } + + for (uint64_t i = 0; i < n; ++i) { + solver->performMatrixVectorMultiplication(x, *multiplyResult, b); + + // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost + // element of the min/max operator stack. + storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices()); + } + + if (!multiplyResultMemoryProvided) { + delete multiplyResult; + } } template diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h index c1781648e..50b4fa88b 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/solver/StandardMinMaxLinearEquationSolver.h @@ -15,11 +15,20 @@ namespace storm { }; void setSolutionMethod(SolutionMethod const& solutionMethod); - + void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations); + void setRelativeTerminationCriterion(bool value); + void setPrecision(double precision); + SolutionMethod const& getSolutionMethod() const; - + uint64_t getMaximalNumberOfIterations() const; + double getPrecision() const; + bool getRelativeTerminationCriterion() const; + private: SolutionMethod solutionMethod; + uint64_t maximalNumberOfIterations; + double precision; + bool relative; }; template @@ -28,13 +37,16 @@ namespace storm { 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 d, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override; - virtual void performMatrixVectorMultiplication(OptimizationDirection d, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + 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(); private: + void solveEquationSystemPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const; + void solveEquationSystemValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const; + /// The settings of this solver. StandardMinMaxLinearEquationSolverSettings settings; diff --git a/test/functional/solver/MinMaxTechniqueSelectionTest.cpp b/test/functional/solver/MinMaxTechniqueSelectionTest.cpp index 78cc8696e..102cc5519 100644 --- a/test/functional/solver/MinMaxTechniqueSelectionTest.cpp +++ b/test/functional/solver/MinMaxTechniqueSelectionTest.cpp @@ -2,9 +2,9 @@ #include "src/solver/MinMaxLinearEquationSolver.h" -TEST( MinMaxTechnique, Simple ) { - storm::solver::MinMaxTechniqueSelection ts = storm::solver::MinMaxTechniqueSelection::PolicyIteration; - storm::solver::MinMaxTechnique t = storm::solver::MinMaxTechnique::PolicyIteration; +TEST( MinMaxMethod, Simple ) { + storm::solver::MinMaxMethodSelection ts = storm::solver::MinMaxMethodSelection::PolicyIteration; + storm::solver::MinMaxMethod t = storm::solver::MinMaxMethod::PolicyIteration; ASSERT_EQ(convert(ts), t);