diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp index 005a33997..82e9c8e6c 100644 --- a/src/modelchecker/region/ApproximationModel.cpp +++ b/src/modelchecker/region/ApproximationModel.cpp @@ -269,26 +269,42 @@ namespace storm { template ConstantType ApproximationModel::computeInitialStateValue(ParameterRegion const& region, bool computeLowerBounds) { instantiate(region, computeLowerBounds); - std::vector policy; - invokeSolver(computeLowerBounds, policy); + std::shared_ptr policy; + if(computeLowerBounds && !this->minimizingPolicies.empty()){ + policy = std::make_shared(**(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(**(this->maximizingPolicies.begin())); //we need a fresh policy, initialized with the values of the old one + } else { + policy = std::make_shared(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 fixedVariables; + std::map> 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::value_type(sub.first, RegionBoundary::UNSPECIFIED)); + VarCount.insert(typename std::map>::value_type(sub.first, std::pair(0,0))); } } } //Now erase the parameters that are mapped to different boundaries. for(std::size_t rowGroup=0; rowGroupmatrixData.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 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() : storm::utility::zero(); //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( 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, double>::invokeSolver(bool computeLowerBounds, std::vector& policy){ + void ApproximationModel, double>::invokeSolver(bool computeLowerBounds, Policy& policy){ storm::solver::SolveGoal goal(computeLowerBounds); std::unique_ptr> solver = storm::solver::configureMinMaxLinearEquationSolver(goal, storm::utility::solver::MinMaxLinearEquationSolverFactory(), 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, double>::invokeSolver(bool computeLowerBounds, std::vector& policy){ + void ApproximationModel, double>::invokeSolver(bool computeLowerBounds, Policy& policy){ storm::solver::SolveGoal player2Goal(computeLowerBounds); std::unique_ptr> solver = storm::utility::solver::GameSolverFactory().create(this->player1Matrix, this->matrixData.matrix); solver->solveGame(this->player1Goal.direction(), player2Goal.direction(), this->eqSysResult, this->vectorData.vector); diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h index 787d1f33e..c4687d274 100644 --- a/src/modelchecker/region/ApproximationModel.h +++ b/src/modelchecker/region/ApproximationModel.h @@ -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 FunctionSubstitution; + typedef std::vector 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::value_type FunctionEntry; void initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector const& newIndices); void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector const& newIndices); void initializePlayer1Matrix(ParametricSparseModelType const& parametricModel); void instantiate(ParameterRegion const& region, bool computeLowerBounds); - void invokeSolver(bool computeLowerBounds, std::vector& 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 eqSysResult; //The index which represents the result for the initial state in the eqSysResult vector std::size_t eqSysInitIndex; + + std::unordered_set> minimizingPolicies; + std::unordered_set> 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) diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp index 6c82ab1cf..e65ceed42 100644 --- a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp +++ b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp @@ -30,7 +30,7 @@ namespace storm { template - void GmmxxMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void GmmxxMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX, std::vector* 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(x.size()); xMemoryProvided = false; } + + if(initialPolicy != nullptr){ + //Get initial values for x like it is done for policy iteration. + storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(*initialPolicy, true); + submatrix.convertToEquationSystem(); + GmmxxLinearEquationSolver gmmLinearEquationSolver(submatrix); + std::vector subB(rowGroupIndices.size() - 1); + storm::utility::vector::selectVectorValues(subB, *initialPolicy, rowGroupIndices, b); + // Solve the resulting linear equation system + std::vector 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(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(this->A.getRowGroupIndices().size() - 1); + if(initialPolicy == nullptr){ + this->policy = std::vector(this->A.getRowGroupIndices().size() - 1); + } else { + this->policy = *initialPolicy; + } // Create our own multiplyResult for solving the deterministic sub-instances. std::vector deterministicMultiplyResult(rowGroupIndices.size() - 1); diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.h b/src/solver/GmmxxMinMaxLinearEquationSolver.h index ed91278bc..8bf92086d 100644 --- a/src/solver/GmmxxMinMaxLinearEquationSolver.h +++ b/src/solver/GmmxxMinMaxLinearEquationSolver.h @@ -35,7 +35,7 @@ namespace storm { 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 d, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override; + virtual void solveEquationSystem(OptimizationDirection d, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr, std::vector* initialPolicy = nullptr) const override; private: // The (gmm++) matrix associated with this equation solver. diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h index 8514ee69a..934aae3f2 100644 --- a/src/solver/MinMaxLinearEquationSolver.h +++ b/src/solver/MinMaxLinearEquationSolver.h @@ -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& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const = 0; + virtual void solveEquationSystem(OptimizationDirection d, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr, std::vector* 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 computePolicy(std::vector& x, std::vector const& b) const{ - std::vector xPrime(this->A.getRowCount()); - this->A.multiplyWithVector(x, xPrime); - storm::utility::vector::addVectors(xPrime, b, xPrime); - std::vector policy(x.size()); - std::vector reduced(x.size()); - storm::utility::vector::reduceVectorMinOrMax(convert(this->direction), xPrime, reduced, this->A.getRowGroupIndices(), &(policy)); - return policy; - } - storm::storage::SparseMatrix const& A; std::unique_ptr> earlyTermination; diff --git a/src/solver/NativeMinMaxLinearEquationSolver.cpp b/src/solver/NativeMinMaxLinearEquationSolver.cpp index 4835de497..44a8e7ba2 100644 --- a/src/solver/NativeMinMaxLinearEquationSolver.cpp +++ b/src/solver/NativeMinMaxLinearEquationSolver.cpp @@ -31,7 +31,7 @@ namespace storm { } template - void NativeMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void NativeMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX, std::vector* 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(x.size()); xMemoryProvided = false; } + + if(initialPolicy != nullptr){ + //Get initial values for x like it is done for policy iteration. + storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(*initialPolicy, true); + submatrix.convertToEquationSystem(); + NativeLinearEquationSolver nativeLinearEquationSolver(submatrix); + std::vector subB(this->A.getRowGroupIndices().size() - 1); + storm::utility::vector::selectVectorValues(subB, *initialPolicy, this->A.getRowGroupIndices(), b); + // Solve the resulting linear equation system + std::vector 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(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(this->A.getRowGroupIndices().size() - 1); + if(initialPolicy == nullptr){ + this->policy = std::vector(this->A.getRowGroupIndices().size() - 1); + } else { + this->policy = *initialPolicy; + } // Create our own multiplyResult for solving the deterministic sub-instances. std::vector deterministicMultiplyResult(this->A.getRowGroupIndices().size() - 1); diff --git a/src/solver/NativeMinMaxLinearEquationSolver.h b/src/solver/NativeMinMaxLinearEquationSolver.h index 5fcd76aa9..ceaf72fcd 100644 --- a/src/solver/NativeMinMaxLinearEquationSolver.h +++ b/src/solver/NativeMinMaxLinearEquationSolver.h @@ -34,7 +34,7 @@ namespace storm { virtual void performMatrixVectorMultiplication(OptimizationDirection dir, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* newX = 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 solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr, std::vector* initialPolicy = nullptr) const override; }; } // namespace solver diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp index f16c92a07..9fa506bc0 100644 --- a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp +++ b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp @@ -46,7 +46,7 @@ namespace storm { } template - void TopologicalMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void TopologicalMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX, std::vector* initialPolicy) const { #ifdef GPU_USE_FLOAT #define __FORCE_FLOAT_CALCULATION true diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.h b/src/solver/TopologicalMinMaxLinearEquationSolver.h index 25c325f75..314c8ac58 100644 --- a/src/solver/TopologicalMinMaxLinearEquationSolver.h +++ b/src/solver/TopologicalMinMaxLinearEquationSolver.h @@ -43,7 +43,7 @@ namespace storm { */ TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true); - virtual void solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override; + virtual void solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr, std::vector* initialPolicy = nullptr) const override; private: bool enableCuda;