sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Added Policy Iteration to the NativeMinMaxLinearEquationSolver. 

Added a test.


Former-commit-id: 087934eb47
tempestpy_adaptions
PBerger 10 years ago
parent
56e1b79085
commit
287393abc4
5 changed files with 178 additions and 67 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 3
      src/solver/GmmxxMinMaxLinearEquationSolver.cpp
  2. 3
      src/solver/NativeLinearEquationSolver.cpp
  3. 213
      src/solver/NativeMinMaxLinearEquationSolver.cpp
  4. 5
      src/solver/NativeMinMaxLinearEquationSolver.h
  5. 21
      test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp

3
src/solver/GmmxxMinMaxLinearEquationSolver.cpp
View File

@ -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<storm::storage::SparseMatrix<ValueType>::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<ValueType>(*currentX, *newX, static_cast<ValueType>(this->precision), this->relative);
// Update environment variables.
std::swap(currentX, newX);

3
src/solver/NativeLinearEquationSolver.cpp
View File

@ -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<ValueType>(*currentX, *nextX, static_cast<ValueType>(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<double>;
template class NativeLinearEquationSolver<float>;
}
}

213
src/solver/NativeMinMaxLinearEquationSolver.cpp
View File

@ -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<ValueType>::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> 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<typename ValueType>
NativeMinMaxLinearEquationSolver<ValueType>::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
NativeMinMaxLinearEquationSolver<ValueType>::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> 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<typename ValueType>
void NativeMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* 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<ValueType>(A.getRowCount());
multiplyResultMemoryProvided = false;
}
std::vector<ValueType>* currentX = &x;
bool xMemoryProvided = true;
if (newX == nullptr) {
newX = new std::vector<ValueType>(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<ValueType>* 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<ValueType>(*currentX, *newX, static_cast<ValueType>(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<ValueType>(A.getRowCount());
multiplyResultMemoryProvided = false;
}
std::vector<ValueType>* currentX = &x;
bool xMemoryProvided = true;
if (newX == nullptr) {
newX = new std::vector<ValueType>(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<ValueType>* 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<ValueType>(*currentX, *newX, static_cast<ValueType>(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<storm::storage::SparseMatrix<ValueType>::index_type> choiceVector(A.getRowGroupIndices().size() - 1);
// Create our own multiplyResult for solving the deterministic instances.
std::vector<ValueType> deterministicMultiplyResult(A.getRowGroupIndices().size() - 1);
std::vector<ValueType> subB(A.getRowGroupIndices().size() - 1);
bool multiplyResultMemoryProvided = true;
if (multiplyResult == nullptr) {
multiplyResult = new std::vector<ValueType>(b.size());
multiplyResultMemoryProvided = false;
}
std::vector<ValueType>* currentX = &x;
bool xMemoryProvided = true;
if (newX == nullptr) {
newX = new std::vector<ValueType>(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<ValueType>* 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<ValueType> submatrix = A.selectRowsFromRowGroups(choiceVector, true);
submatrix.convertToEquationSystem();
NativeLinearEquationSolver<ValueType> nativeLinearEquationSolver(submatrix);
storm::utility::vector::selectVectorValues<ValueType>(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<ValueType>(*currentX, *newX, static_cast<ValueType>(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<typename ValueType>

5
src/solver/NativeMinMaxLinearEquationSolver.h
View File

@ -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<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative = true);
virtual void performMatrixVectorMultiplication(bool minimize, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* 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

21
test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp
View File

@ -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<double> builder(0, 0, 0, false, true);
ASSERT_NO_THROW(builder.newRowGroup(0));
ASSERT_NO_THROW(builder.addNextValue(0, 0, 0.9));
storm::storage::SparseMatrix<double> A;
ASSERT_NO_THROW(A = builder.build(2));
std::vector<double> x(1);
std::vector<double> b = { 0.099, 0.5 };
storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
storm::solver::NativeMinMaxLinearEquationSolver<double> 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());
}
Write Preview
Loading…
Cancel
Save