From 287393abc437b31498e2e62352b9105490209fa9 Mon Sep 17 00:00:00 2001 From: PBerger Date: Sat, 30 May 2015 16:06:00 +0200 Subject: [PATCH] Added Policy Iteration to the NativeMinMaxLinearEquationSolver. Added a test. Former-commit-id: 087934eb473d9283bdf1c6f5070cebeab69b1daf --- .../GmmxxMinMaxLinearEquationSolver.cpp | 3 +- src/solver/NativeLinearEquationSolver.cpp | 3 +- .../NativeMinMaxLinearEquationSolver.cpp | 213 ++++++++++++------ src/solver/NativeMinMaxLinearEquationSolver.h | 5 +- .../NativeMinMaxLinearEquationSolverTest.cpp | 21 ++ 5 files changed, 178 insertions(+), 67 deletions(-) diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp index e19a37548..6c5756c1a 100644 --- a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp +++ b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp @@ -95,7 +95,6 @@ namespace storm { } } 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); @@ -150,7 +149,7 @@ namespace storm { } // Determine whether the method converged. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative); + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, static_cast(this->precision), this->relative); // Update environment variables. std::swap(currentX, newX); diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index ab4a0355b..a2480ddc5 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -63,7 +63,7 @@ namespace storm { std::swap(nextX, currentX); // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative); + converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, static_cast(precision), relative); // Increase iteration count so we can abort if convergence is too slow. ++iterationCount; @@ -128,5 +128,6 @@ namespace storm { // Explicitly instantiate the linear equation solver. template class NativeLinearEquationSolver; + template class NativeLinearEquationSolver; } } \ No newline at end of file diff --git a/src/solver/NativeMinMaxLinearEquationSolver.cpp b/src/solver/NativeMinMaxLinearEquationSolver.cpp index 51da3e692..a33b7e869 100644 --- a/src/solver/NativeMinMaxLinearEquationSolver.cpp +++ b/src/solver/NativeMinMaxLinearEquationSolver.cpp @@ -4,6 +4,7 @@ #include "src/settings/SettingsManager.h" #include "src/utility/vector.h" +#include "src/solver/NativeLinearEquationSolver.h" namespace storm { namespace solver { @@ -12,82 +13,168 @@ namespace storm { NativeMinMaxLinearEquationSolver::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A) : A(A) { // Get the settings object to customize solving. storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings(); + storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings(); // Get appropriate settings. maximalNumberOfIterations = settings.getMaximalIterationCount(); precision = settings.getPrecision(); relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative; + useValueIteration = (generalSettings.getMinMaxEquationSolvingTechnique() == storm::settings::modules::GeneralSettings::MinMaxTechnique::ValueIteration); } template - NativeMinMaxLinearEquationSolver::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) { + NativeMinMaxLinearEquationSolver::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative) : A(A), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), useValueIteration(useValueIteration), relative(relative) { // Intentionally left empty. } template void NativeMinMaxLinearEquationSolver::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(A.getRowCount()); - 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. - 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. - if (minimize) { - storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, A.getRowGroupIndices()); - } else { - storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, A.getRowGroupIndices()); - } - - // Determine whether the method converged. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, static_cast(precision), 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(A.getRowCount()); + 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. + 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. + if (minimize) { + storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, A.getRowGroupIndices()); + } else { + storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, A.getRowGroupIndices()); + } + + // Determine whether the method converged. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, static_cast(precision), 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(A.getRowGroupIndices().size() - 1); + + // Create our own multiplyResult for solving the deterministic instances. + std::vector deterministicMultiplyResult(A.getRowGroupIndices().size() - 1); + std::vector subB(A.getRowGroupIndices().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 = A.selectRowsFromRowGroups(choiceVector, true); + submatrix.convertToEquationSystem(); + + NativeLinearEquationSolver nativeLinearEquationSolver(submatrix); + storm::utility::vector::selectVectorValues(subB, choiceVector, A.getRowGroupIndices(), b); + + // Copy X since we will overwrite it + std::copy(currentX->begin(), currentX->end(), newX->begin()); + + // Solve the resulting linear equation system + nativeLinearEquationSolver.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. + A.multiplyWithVector(*newX, *multiplyResult); + storm::utility::vector::addVectors(*multiplyResult, b, *multiplyResult); + + // Reduce the vector x by applying min/max over all nondeterministic choices. + if (minimize) { + storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, A.getRowGroupIndices(), &choiceVector); + } else { + storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, A.getRowGroupIndices(), &choiceVector); + } + + // Determine whether the method converged. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, static_cast(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/NativeMinMaxLinearEquationSolver.h b/src/solver/NativeMinMaxLinearEquationSolver.h index 9044f5ca8..682f79cca 100644 --- a/src/solver/NativeMinMaxLinearEquationSolver.h +++ b/src/solver/NativeMinMaxLinearEquationSolver.h @@ -29,7 +29,7 @@ namespace storm { * @param relative If set, the relative error rather than the absolute error is considered for convergence * detection. */ - NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true); + NativeMinMaxLinearEquationSolver(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* newX = nullptr) const override; @@ -47,6 +47,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/NativeMinMaxLinearEquationSolverTest.cpp b/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp index d2e014dd5..10c630395 100644 --- a/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp +++ b/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp @@ -62,4 +62,25 @@ TEST(NativeMinMaxLinearEquationSolver, MatrixVectorMultiplication) { x = {0, 1, 0}; ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, x, nullptr, 20)); ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::nativeEquationSolverSettings().getPrecision()); +} + +TEST(NativeMinMaxLinearEquationSolver, 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::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings(); + storm::solver::NativeMinMaxLinearEquationSolver solver(A, settings.getPrecision(), settings.getMaximalIterationCount(), false, settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative); + + ASSERT_NO_THROW(solver.solveEquationSystem(true, x, b)); + ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::nativeEquationSolverSettings().getPrecision()); + + ASSERT_NO_THROW(solver.solveEquationSystem(false, x, b)); + ASSERT_LT(std::abs(x[0] - 0.99), storm::settings::nativeEquationSolverSettings().getPrecision()); } \ No newline at end of file