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

MinMaxLinearEqSolvers can now use some initial policy as a first guess. 

First steps to use this for region approximation


Former-commit-id: 9a8151607f
tempestpy_adaptions
TimQu 9 years ago
parent
a723cfca14
commit
ca917a651c
9 changed files with 110 additions and 99 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. 58
      src/modelchecker/region/ApproximationModel.cpp
  2. 70
      src/modelchecker/region/ApproximationModel.h
  3. 29
      src/solver/GmmxxMinMaxLinearEquationSolver.cpp
  4. 2
      src/solver/GmmxxMinMaxLinearEquationSolver.h
  5. 14
      src/solver/MinMaxLinearEquationSolver.h
  6. 30
      src/solver/NativeMinMaxLinearEquationSolver.cpp
  7. 2
      src/solver/NativeMinMaxLinearEquationSolver.h
  8. 2
      src/solver/TopologicalMinMaxLinearEquationSolver.cpp
  9. 2
      src/solver/TopologicalMinMaxLinearEquationSolver.h

58
src/modelchecker/region/ApproximationModel.cpp
View File

@ -269,26 +269,42 @@ namespace storm {
template<typename ParametricSparseModelType, typename ConstantType>
ConstantType ApproximationModel<ParametricSparseModelType, ConstantType>::computeInitialStateValue(ParameterRegion<ParametricType> const& region, bool computeLowerBounds) {
instantiate(region, computeLowerBounds);
std::vector<std::size_t> policy;
invokeSolver(computeLowerBounds, policy);
std::shared_ptr<Policy> policy;
if(computeLowerBounds && !this->minimizingPolicies.empty()){
policy = std::make_shared<Policy>(**(this->minimizingPolicies.begin())); //we need a fresh policy, initialized with the values of the old one
} else if(!computeLowerBounds && !this->maximizingPolicies.empty()){
policy = std::make_shared<Policy>(**(this->maximizingPolicies.begin())); //we need a fresh policy, initialized with the values of the old one
} else {
policy = std::make_shared<Policy>(this->matrixData.matrix.getRowGroupCount());
}
invokeSolver(computeLowerBounds, *policy);
if(computeLowerBounds){
this->minimizingPolicies.insert(policy);
std::cout << minimizingPolicies.size() << std::endl;
} else {
this->maximizingPolicies.insert(policy);
std::cout << maximizingPolicies.size() << std::endl;
}
//TODO: policy for games.
//TODO: do this only when necessary.
//TODO: (maybe) when a few parameters are mapped to another value, build a "nicer" scheduler and check whether it induces values that are more optimal.
if(policy.empty()) return this->eqSysResult[this->eqSysInitIndex];
if(policy->empty()) return this->eqSysResult[this->eqSysInitIndex];
//Get the set of parameters which are (according to the policy) always mapped to the same region boundary.
//First, collect all (relevant) parameters, i.e., the ones that are set to the lower or upper boundary.
std::map<VariableType, RegionBoundary> fixedVariables;
std::map<VariableType, std::pair<std::size_t, std::size_t>> VarCount;
std::size_t substitutionCount =0;
for(auto const& substitution : this->funcSubData.substitutions){
for( auto const& sub : substitution){
if(sub.second!= RegionBoundary::UNSPECIFIED){
fixedVariables.insert(typename std::map<VariableType, RegionBoundary>::value_type(sub.first, RegionBoundary::UNSPECIFIED));
VarCount.insert(typename std::map<VariableType, std::pair<std::size_t, std::size_t>>::value_type(sub.first, std::pair<std::size_t, std::size_t>(0,0)));
}
}
}
//Now erase the parameters that are mapped to different boundaries.
for(std::size_t rowGroup=0; rowGroup<this->matrixData.matrix.getRowGroupCount(); ++rowGroup){
std::size_t row = this->matrixData.matrix.getRowGroupIndices()[rowGroup] + policy[rowGroup];
std::size_t row = this->matrixData.matrix.getRowGroupIndices()[rowGroup] + (*policy)[rowGroup];
for(std::pair<VariableType, RegionBoundary> const& sub : this->funcSubData.substitutions[this->matrixData.rowSubstitutions[row]]){
auto fixedVarIt = fixedVariables.find(sub.first);
if(fixedVarIt != fixedVariables.end() && fixedVarIt->second != sub.second){
@ -299,15 +315,29 @@ namespace storm {
fixedVariables.erase(fixedVarIt);
}
}
auto varcountIt = VarCount.find(sub.first);
if(sub.second==RegionBoundary::LOWER){
++(varcountIt->second.first);
} else if (sub.second==RegionBoundary::UPPER){
++(varcountIt->second.second);
}
++substitutionCount;
}
if (fixedVariables.empty()){
break;
// break;
}
}
// std::cout << "Used Approximation" << std::endl;
for (auto const& varcount : VarCount){
if(varcount.second.first > 0 && varcount.second.second > 0){
// std::cout << " Variable " << varcount.first << " has been set to lower " << varcount.second.first << " times and to upper " << varcount.second.second << " times. (total: " << substitutionCount << ")" << std::endl;
}
}
for (auto const& fixVar : fixedVariables){
//std::cout << "APPROXMODEL: variable " << fixVar.first << " is always mapped to " << fixVar.second << std::endl;
//std::cout << " APPROXMODEL: variable " << fixVar.first << " is always mapped to " << fixVar.second << std::endl;
}
// std::cout << " Result is " << this->eqSysResult[this->eqSysInitIndex] << std::endl;
return this->eqSysResult[this->eqSysInitIndex];
}
@ -342,17 +372,17 @@ namespace storm {
auto& result = functionResult.second;
result = computeLowerBounds ? storm::utility::infinity<ConstantType>() : storm::utility::zero<ConstantType>();
//Iterate over the different combinations of lower and upper bounds and update the min and max values
auto const& vertices=region.getVerticesOfRegion(unspecifiedParameters[funcSub.getSubstitution()]);
auto const& vertices=region.getVerticesOfRegion(unspecifiedParameters[funcSub.second]);
for(auto const& vertex : vertices){
//extend the substitution
for(auto const& vertexSub : vertex){
instantiatedSubs[funcSub.getSubstitution()][vertexSub.first]=vertexSub.second;
instantiatedSubs[funcSub.second][vertexSub.first]=vertexSub.second;
}
//evaluate the function
ConstantType currValue = storm::utility::region::convertNumber<ConstantType>(
storm::utility::region::evaluateFunction(
funcSub.getFunction(),
instantiatedSubs[funcSub.getSubstitution()]
funcSub.first,
instantiatedSubs[funcSub.second]
)
);
result = computeLowerBounds ? std::min(result, currValue) : std::max(result, currValue);
@ -370,16 +400,16 @@ namespace storm {
template<>
void ApproximationModel<storm::models::sparse::Dtmc<storm::RationalFunction>, double>::invokeSolver(bool computeLowerBounds, std::vector<std::size_t>& policy){
void ApproximationModel<storm::models::sparse::Dtmc<storm::RationalFunction>, double>::invokeSolver(bool computeLowerBounds, Policy& policy){
storm::solver::SolveGoal goal(computeLowerBounds);
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<double>> solver = storm::solver::configureMinMaxLinearEquationSolver(goal, storm::utility::solver::MinMaxLinearEquationSolverFactory<double>(), this->matrixData.matrix);
solver->setPolicyTracking();
solver->solveEquationSystem(this->eqSysResult, this->vectorData.vector);
solver->solveEquationSystem(goal.direction(), this->eqSysResult, this->vectorData.vector, nullptr, nullptr, &policy);
policy = solver->getPolicy();
}
template<>
void ApproximationModel<storm::models::sparse::Mdp<storm::RationalFunction>, double>::invokeSolver(bool computeLowerBounds, std::vector<std::size_t>& policy){
void ApproximationModel<storm::models::sparse::Mdp<storm::RationalFunction>, double>::invokeSolver(bool computeLowerBounds, Policy& policy){
storm::solver::SolveGoal player2Goal(computeLowerBounds);
std::unique_ptr<storm::solver::GameSolver<double>> solver = storm::utility::solver::GameSolverFactory<double>().create(this->player1Matrix, this->matrixData.matrix);
solver->solveGame(this->player1Goal.direction(), player2Goal.direction(), this->eqSysResult, this->vectorData.vector);

70
src/modelchecker/region/ApproximationModel.h
View File

@ -58,74 +58,26 @@ namespace storm {
private:
//A class that represents a function and how it should be substituted (i.e. which variables should be replaced with lower and which with upper bounds of the region)
//The substitution is given as an index in funcSubData.substitutions (allowing to instantiate the substitutions more easily).
class FunctionSubstitution {
public:
FunctionSubstitution(ParametricType const& fun, std::size_t const& sub) : hash(computeHash(fun,sub)), function(fun), substitution(sub) {
//intentionally left empty
}
FunctionSubstitution(ParametricType&& fun, std::size_t&& sub) : hash(computeHash(fun,sub)), function(std::move(fun)), substitution(std::move(sub)) {
//intentionally left empty
}
FunctionSubstitution(FunctionSubstitution const& other) : hash(other.hash), function(other.function), substitution(other.substitution){
//intentionally left empty
}
FunctionSubstitution(FunctionSubstitution&& other) : hash(std::move(other.hash)), function(std::move(other.function)), substitution(std::move(other.substitution)){
//intentionally left empty
}
FunctionSubstitution() = default;
~FunctionSubstitution() = default;
bool operator==(FunctionSubstitution const& other) const {
return this->hash==other.hash && this->substitution==other.substitution && this->function==other.function;
}
ParametricType const& getFunction() const{
return this->function;
}
std::size_t const& getSubstitution() const{
return this->substitution;
}
std::size_t const& getHash() const{
return this->hash;
}
private:
static std::size_t computeHash(ParametricType const& fun, std::size_t const& sub) {
std::size_t seed = 0;
boost::hash_combine(seed, fun);
boost::hash_combine(seed, sub);
return seed;
}
std::size_t hash;
ParametricType function;
std::size_t substitution;
};
typedef std::pair<ParametricType, std::size_t> FunctionSubstitution;
typedef std::vector<storm::storage::sparse::state_type> Policy;
class FuncSubHash{
public:
std::size_t operator()(FunctionSubstitution const& fs) const {
return fs.getHash();
// return fs.getHash();
std::size_t seed = 0;
boost::hash_combine(seed, fs.first);
boost::hash_combine(seed, fs.second);
return seed;
}
};
typedef typename std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash>::value_type FunctionEntry;
void initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices);
void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices);
void initializePlayer1Matrix(ParametricSparseModelType const& parametricModel);
void instantiate(ParameterRegion<ParametricType> const& region, bool computeLowerBounds);
void invokeSolver(bool computeLowerBounds, std::vector<std::size_t>& policy);
void invokeSolver(bool computeLowerBounds, Policy& policy);
//Some designated states in the original model
storm::storage::BitVector targetStates, maybeStates;
@ -134,6 +86,10 @@ namespace storm {
std::vector<ConstantType> eqSysResult;
//The index which represents the result for the initial state in the eqSysResult vector
std::size_t eqSysInitIndex;
std::unordered_set<std::shared_ptr<Policy>> minimizingPolicies;
std::unordered_set<std::shared_ptr<Policy>> maximizingPolicies;
//A flag that denotes whether we compute probabilities or rewards
bool computeRewards;
//Player 1 represents the nondeterminism of the given mdp (so, this is irrelevant if we approximate values of a DTMC)

29
src/solver/GmmxxMinMaxLinearEquationSolver.cpp
View File

@ -30,7 +30,7 @@ namespace storm {
template<typename ValueType>
void GmmxxMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
void GmmxxMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX, std::vector<storm::storage::sparse::state_type>* initialPolicy) const {
if (this->useValueIteration) {
// Set up the environment for the power method. If scratch memory was not provided, we need to create it.
bool multiplyResultMemoryProvided = true;
@ -45,6 +45,19 @@ namespace storm {
newX = new std::vector<ValueType>(x.size());
xMemoryProvided = false;
}
if(initialPolicy != nullptr){
//Get initial values for x like it is done for policy iteration.
storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(*initialPolicy, true);
submatrix.convertToEquationSystem();
GmmxxLinearEquationSolver<ValueType> gmmLinearEquationSolver(submatrix);
std::vector<ValueType> subB(rowGroupIndices.size() - 1);
storm::utility::vector::selectVectorValues<ValueType>(subB, *initialPolicy, rowGroupIndices, b);
// Solve the resulting linear equation system
std::vector<ValueType> deterministicMultiplyResult(rowGroupIndices.size() - 1);
gmmLinearEquationSolver.solveEquationSystem(*currentX, subB, &deterministicMultiplyResult);
}
uint_fast64_t iterations = 0;
bool converged = false;
@ -81,6 +94,11 @@ namespace storm {
if (currentX == copyX) {
std::swap(x, *currentX);
}
if(this->trackPolicy){
this->policy = std::vector<storm::storage::sparse::state_type>(x.size());
storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, rowGroupIndices, &(this->policy));
}
if (!xMemoryProvided) {
delete copyX;
@ -90,14 +108,15 @@ namespace storm {
delete multiplyResult;
}
if(this->trackPolicy){
this->policy = this->computePolicy(x,b);
}
} else {
// We will use Policy Iteration to solve the given system.
// We first guess an initial choice resolution which will be refined after each iteration.
this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
if(initialPolicy == nullptr){
this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
} else {
this->policy = *initialPolicy;
}
// Create our own multiplyResult for solving the deterministic sub-instances.
std::vector<ValueType> deterministicMultiplyResult(rowGroupIndices.size() - 1);

2
src/solver/GmmxxMinMaxLinearEquationSolver.h
View File

@ -35,7 +35,7 @@ namespace storm {
virtual void performMatrixVectorMultiplication(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
virtual void solveEquationSystem(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
virtual void solveEquationSystem(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr, std::vector<storm::storage::sparse::state_type>* initialPolicy = nullptr) const override;
private:
// The (gmm++) matrix associated with this equation solver.

14
src/solver/MinMaxLinearEquationSolver.h
View File

@ -8,7 +8,6 @@
#include "src/storage/sparse/StateType.h"
#include "AllowEarlyTerminationCondition.h"
#include "OptimizationDirection.h"
#include "src/utility/vector.h"
namespace storm {
namespace storage {
@ -98,7 +97,7 @@ namespace storm {
* vector must be equal to the length of the vector x (and thus to the number of columns of A).
* @return The solution vector x of the system of linear equations as the content of the parameter x.
*/
virtual void solveEquationSystem(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const = 0;
virtual void solveEquationSystem(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr, std::vector<storm::storage::sparse::state_type>* initialPolicy = nullptr) const = 0;
/*!
* As solveEquationSystem with an optimization-direction, but this uses the internally set direction.
@ -141,17 +140,6 @@ namespace storm {
protected:
std::vector<storm::storage::sparse::state_type> computePolicy(std::vector<ValueType>& x, std::vector<ValueType> const& b) const{
std::vector<ValueType> xPrime(this->A.getRowCount());
this->A.multiplyWithVector(x, xPrime);
storm::utility::vector::addVectors(xPrime, b, xPrime);
std::vector<storm::storage::sparse::state_type> policy(x.size());
std::vector<ValueType> reduced(x.size());
storm::utility::vector::reduceVectorMinOrMax(convert(this->direction), xPrime, reduced, this->A.getRowGroupIndices(), &(policy));
return policy;
}
storm::storage::SparseMatrix<ValueType> const& A;
std::unique_ptr<AllowEarlyTerminationCondition<ValueType>> earlyTermination;

30
src/solver/NativeMinMaxLinearEquationSolver.cpp
View File

@ -31,7 +31,7 @@ namespace storm {
}
template<typename ValueType>
void NativeMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
void NativeMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX, std::vector<storm::storage::sparse::state_type>* initialPolicy) const {
if (this->useValueIteration) {
// Set up the environment for the power method. If scratch memory was not provided, we need to create it.
bool multiplyResultMemoryProvided = true;
@ -45,6 +45,19 @@ namespace storm {
newX = new std::vector<ValueType>(x.size());
xMemoryProvided = false;
}
if(initialPolicy != nullptr){
//Get initial values for x like it is done for policy iteration.
storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(*initialPolicy, true);
submatrix.convertToEquationSystem();
NativeLinearEquationSolver<ValueType> nativeLinearEquationSolver(submatrix);
std::vector<ValueType> subB(this->A.getRowGroupIndices().size() - 1);
storm::utility::vector::selectVectorValues<ValueType>(subB, *initialPolicy, this->A.getRowGroupIndices(), b);
// Solve the resulting linear equation system
std::vector<ValueType> deterministicMultiplyResult(this->A.getRowGroupIndices().size() - 1);
nativeLinearEquationSolver.solveEquationSystem(*currentX, subB, &deterministicMultiplyResult);
}
uint_fast64_t iterations = 0;
bool converged = false;
@ -83,6 +96,11 @@ namespace storm {
std::swap(x, *currentX);
}
if(this->trackPolicy){
this->policy = std::vector<storm::storage::sparse::state_type>(x.size());
storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices(), &(this->policy));
}
if (!xMemoryProvided) {
delete copyX;
}
@ -91,14 +109,14 @@ namespace storm {
delete multiplyResult;
}
if(this->trackPolicy){
this->policy = this->computePolicy(x,b);
}
} else {
// We will use Policy Iteration to solve the given system.
// We first guess an initial choice resolution which will be refined after each iteration.
this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
if(initialPolicy == nullptr){
this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
} else {
this->policy = *initialPolicy;
}
// Create our own multiplyResult for solving the deterministic sub-instances.
std::vector<ValueType> deterministicMultiplyResult(this->A.getRowGroupIndices().size() - 1);

2
src/solver/NativeMinMaxLinearEquationSolver.h
View File

@ -34,7 +34,7 @@ namespace storm {
virtual void performMatrixVectorMultiplication(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* newX = nullptr) const override;
virtual void solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
virtual void solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr, std::vector<storm::storage::sparse::state_type>* initialPolicy = nullptr) const override;
};
} // namespace solver

2
src/solver/TopologicalMinMaxLinearEquationSolver.cpp
View File

@ -46,7 +46,7 @@ namespace storm {
}
template<typename ValueType>
void TopologicalMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
void TopologicalMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX, std::vector<storm::storage::sparse::state_type>* initialPolicy) const {
#ifdef GPU_USE_FLOAT
#define __FORCE_FLOAT_CALCULATION true

2
src/solver/TopologicalMinMaxLinearEquationSolver.h
View File

@ -43,7 +43,7 @@ namespace storm {
*/
TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
virtual void solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
virtual void solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr, std::vector<storm::storage::sparse::state_type>* initialPolicy = nullptr) const override;
private:
bool enableCuda;

Write Preview
Loading…
Cancel
Save