diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp index 919c4c9f4..e28b82901 100644 --- a/src/modelchecker/region/ApproximationModel.cpp +++ b/src/modelchecker/region/ApproximationModel.cpp @@ -18,6 +18,7 @@ #include "src/utility/region.h" #include "src/utility/solver.h" #include "src/utility/vector.h" +#include "src/utility/policyguessing.h" #include "src/exceptions/UnexpectedException.h" #include "src/exceptions/InvalidArgumentException.h" #include "exceptions/NotImplementedException.h" @@ -83,7 +84,7 @@ namespace storm { * each rowgroup containing 2^#par rows, where #par is the number of parameters that occur in the original row. * We also store the substitution that needs to be applied for each row. */ - ConstantType dummyValue = storm::utility::one(); + ConstantType dummyNonZeroValue = storm::utility::one(); storm::storage::SparseMatrixBuilder matrixBuilder(0, //Unknown number of rows this->maybeStates.getNumberOfSetBits(), //columns 0, //Unknown number of entries @@ -125,7 +126,7 @@ namespace storm { //Note that this is still executed once, even if no parameters occur. for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ if(this->maybeStates.get(oldEntry.getColumn())){ - matrixBuilder.addNextValue(curRow, newIndices[oldEntry.getColumn()], dummyValue); + matrixBuilder.addNextValue(curRow, newIndices[oldEntry.getColumn()], dummyNonZeroValue); } } ++curRow; @@ -137,6 +138,7 @@ namespace storm { //Now run again through both matrices to get the remaining ingredients of the matrixData and vectorData this->matrixData.assignment.reserve(this->matrixData.matrix.getEntryCount()); + this->matrixData.targetChoices = storm::storage::BitVector(this->matrixData.matrix.getRowCount(), false); this->vectorData.vector = std::vector(this->matrixData.matrix.getRowCount()); //Important to initialize here since iterators have to remain valid auto vectorIt = this->vectorData.vector.begin(); this->vectorData.assignment.reserve(vectorData.vector.size()); @@ -145,6 +147,7 @@ namespace storm { for (std::size_t oldRow = parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup]; oldRow < parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup+1]; ++oldRow){ ParametricType targetProbability = storm::utility::region::getNewFunction(storm::utility::zero()); if(!this->computeRewards){ + //Compute the target probability to insert it in every new row for(auto const& oldEntry : parametricModel.getTransitionMatrix().getRow(oldRow)){ if(this->targetStates.get(oldEntry.getColumn())){ targetProbability += oldEntry.getValue(); @@ -161,20 +164,24 @@ namespace storm { if(storm::utility::isConstant(oldEntry.getValue())){ eqSysMatrixEntry->setValue(storm::utility::region::convertNumber(storm::utility::region::getConstantPart(oldEntry.getValue()))); } else { - auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(oldEntry.getValue(), this->matrixData.rowSubstitutions[curRow]), dummyValue)).first; + auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(oldEntry.getValue(), this->matrixData.rowSubstitutions[curRow]), dummyNonZeroValue)).first; this->matrixData.assignment.emplace_back(eqSysMatrixEntry, functionsIt->second); //Note that references to elements of an unordered map remain valid after calling unordered_map::insert. } ++eqSysMatrixEntry; } + if(this->targetStates.get(oldEntry.getColumn())){ + //Store that this row has a transition to target + this->matrixData.targetChoices.set(curRow); + } } if(!this->computeRewards){ if(storm::utility::isConstant(storm::utility::simplify(targetProbability))){ *vectorIt = storm::utility::region::convertNumber(storm::utility::region::getConstantPart(targetProbability)); } else { - auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(targetProbability, this->matrixData.rowSubstitutions[curRow]), dummyValue)).first; + auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(targetProbability, this->matrixData.rowSubstitutions[curRow]), dummyNonZeroValue)).first; this->vectorData.assignment.emplace_back(vectorIt, functionsIt->second); - *vectorIt = dummyValue; + *vectorIt = dummyNonZeroValue; } } ++vectorIt; @@ -198,7 +205,7 @@ namespace storm { // run through the state reward vector of the parametric model. // Constant entries can be set directly. // For Parametric entries we set a dummy value and insert the corresponding function and the assignment - ConstantType dummyValue = storm::utility::one(); + ConstantType dummyNonZeroValue = storm::utility::one(); auto vectorIt = this->vectorData.vector.begin(); for(auto oldState : this->maybeStates){ if(storm::utility::isConstant(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState])){ @@ -220,10 +227,10 @@ namespace storm { substitution.insert(typename std::map::value_type(rewardVar, RegionBoundary::UNSPECIFIED)); } // insert the FunctionSubstitution - auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState], this->matrixData.rowSubstitutions[matrixRow]), dummyValue)).first; + auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState], this->matrixData.rowSubstitutions[matrixRow]), dummyNonZeroValue)).first; //insert assignment and dummy data this->vectorData.assignment.emplace_back(vectorIt, functionsIt->second); - *vectorIt = dummyValue; + *vectorIt = dummyNonZeroValue; ++vectorIt; } } @@ -392,19 +399,22 @@ namespace storm { 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(goal.direction(), this->solverData.result, this->vectorData.vector, nullptr, nullptr, &policy); - policy = solver->getPolicy(); + storm::utility::policyguessing::solveMinMaxLinearEquationSystem(*solver, + this->solverData.result, this->vectorData.vector, + goal.direction(), + policy, + this->matrixData.targetChoices, (this->computeRewards ? storm::utility::infinity() : storm::utility::zero())); } template<> void ApproximationModel, double>::invokeSolver(bool computeLowerBounds, Policy& policy){ storm::solver::SolveGoal player2Goal(computeLowerBounds); std::unique_ptr> solver = storm::utility::solver::GameSolverFactory().create(this->solverData.player1Matrix, this->matrixData.matrix); - solver->setPolicyTracking(); - solver->solveGame(this->solverData.player1Goal.direction(), player2Goal.direction(), this->solverData.result, this->vectorData.vector, &this->solverData.lastPlayer1Policy, &policy); - this->solverData.lastPlayer1Policy = solver->getPlayer1Policy(); - policy = solver->getPlayer2Policy(); + storm::utility::policyguessing::solveGame(*solver, + this->solverData.result, this->vectorData.vector, + this->solverData.player1Goal.direction(), player2Goal.direction(), + this->solverData.lastPlayer1Policy, policy, + this->matrixData.targetChoices, (this->computeRewards ? storm::utility::infinity() : storm::utility::zero())); } diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h index ac6a77765..b7af5aed6 100644 --- a/src/modelchecker/region/ApproximationModel.h +++ b/src/modelchecker/region/ApproximationModel.h @@ -117,6 +117,7 @@ namespace storm { storm::storage::SparseMatrix matrix; //The matrix itself. std::vector::iterator, ConstantType&>> assignment; // Connection of matrix entries with placeholders std::vector rowSubstitutions; //used to obtain which row corresponds to which substitution (used to retrieve information from a scheduler) + storm::storage::BitVector targetChoices; //indicate which rows of the matrix have a positive value to a target state } matrixData; struct VectorData { std::vector vector; //The vector itself. diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp index f1fa8ef52..5ecc26b91 100644 --- a/src/modelchecker/region/SamplingModel.cpp +++ b/src/modelchecker/region/SamplingModel.cpp @@ -17,6 +17,7 @@ #include "src/utility/region.h" #include "src/utility/solver.h" #include "src/utility/vector.h" +#include "src/utility/policyguessing.h" #include "src/exceptions/UnexpectedException.h" #include "src/exceptions/InvalidArgumentException.h" #include "storage/dd/CuddBdd.h" @@ -117,6 +118,7 @@ namespace storm { //Now run again through both matrices to get the remaining ingredients of the matrixData and vectorData. //Note that we need the matrix (I-P) in case of a dtmc. this->matrixData.assignment.reserve(this->matrixData.matrix.getEntryCount()); + this->matrixData.targetChoices = storm::storage::BitVector(this->matrixData.matrix.getRowCount(), false); this->vectorData.vector = std::vector(this->matrixData.matrix.getRowCount()); //Important to initialize here since iterators have to remain valid auto vectorIt = this->vectorData.vector.begin(); this->vectorData.assignment.reserve(vectorData.vector.size()); @@ -159,8 +161,11 @@ namespace storm { } ++eqSysMatrixEntry; } - else if(!this->computeRewards && this->targetStates.get(oldEntry.getColumn())){ - targetProbability += oldEntry.getValue(); + else if(this->targetStates.get(oldEntry.getColumn())){ + if(!this->computeRewards){ + targetProbability += oldEntry.getValue(); + } + this->matrixData.targetChoices.set(curRow); } } if(!this->computeRewards){ @@ -254,9 +259,11 @@ namespace storm { template<> void SamplingModel, double>::invokeSolver(){ std::unique_ptr> solver = storm::solver::configureMinMaxLinearEquationSolver(this->solverData.solveGoal, storm::utility::solver::MinMaxLinearEquationSolverFactory(), this->matrixData.matrix); - solver->setPolicyTracking(); - solver->solveEquationSystem(this->solverData.solveGoal.direction(), this->solverData.result, this->vectorData.vector, nullptr, nullptr, &this->solverData.lastPolicy); - this->solverData.lastPolicy = solver->getPolicy(); + storm::utility::policyguessing::solveMinMaxLinearEquationSystem(*solver, + this->solverData.result, this->vectorData.vector, + this->solverData.solveGoal.direction(), + this->solverData.lastPolicy, + this->matrixData.targetChoices, (this->computeRewards ? storm::utility::infinity() : storm::utility::zero())); } diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h index e13757ce7..9d0ddb2f3 100644 --- a/src/modelchecker/region/SamplingModel.h +++ b/src/modelchecker/region/SamplingModel.h @@ -89,6 +89,7 @@ namespace storm { struct MatrixData { storm::storage::SparseMatrix matrix; //The matrix itself. std::vector::iterator, ConstantType*>> assignment; // Connection of matrix entries with placeholders + storm::storage::BitVector targetChoices; //indicate which rows of the matrix have a positive value to a target state } matrixData; struct VectorData { std::vector vector; //The vector itself. diff --git a/src/solver/GameSolver.cpp b/src/solver/GameSolver.cpp index 1e12f1dd5..516cfa48c 100644 --- a/src/solver/GameSolver.cpp +++ b/src/solver/GameSolver.cpp @@ -19,45 +19,10 @@ namespace storm { } template - void GameSolver::solveGame(OptimizationDirection player1Goal, OptimizationDirection player2Goal, std::vector& x, std::vector const& b, std::vector* initialPlayer1Policy, std::vector* initialPlayer2Policy) const { - uint_fast64_t numberOfPlayer1States = x.size(); - - if(initialPlayer1Policy != nullptr || initialPlayer2Policy != nullptr){ - //Either we work with both policies or none of them - assert(initialPlayer1Policy != nullptr && initialPlayer2Policy != nullptr); - //The policies select certain rows in the matrix of player2. - //However, note that rows can be selected more then once and in an arbitrary order. - std::vector selectedRows(numberOfPlayer1States); - for (uint_fast64_t pl1State = 0; pl1State < numberOfPlayer1States; ++pl1State){ - auto const& pl1Choice = player1Matrix.getRow(player1Matrix.getRowGroupIndices()[pl1State] + (*initialPlayer1Policy)[pl1State]); - assert(pl1Choice.getNumberOfEntries()==1); - uint_fast64_t pl2State = pl1Choice.begin()->getColumn(); - selectedRows[pl1State] = player2Matrix.getRowGroupIndices()[pl2State] + (*initialPlayer2Policy)[pl2State]; - } - storm::storage::SparseMatrix eqSysMatrix = player2Matrix.selectRowsFromRowIndexSequence(selectedRows, true); - std::vector subB(numberOfPlayer1States); - storm::utility::vector::selectVectorValues(subB, selectedRows, b); - //depending on the choices, qualitative properties might have changed. - storm::storage::BitVector targetStates(subB.size(), false); - for(std::size_t index = 0; index < targetStates.size(); ++index){ - if(!storm::utility::isZero(subB[index])){ - targetStates.set(index); - } - } - auto prob01States = storm::utility::graph::performProb01(eqSysMatrix.transpose(), storm::storage::BitVector(targetStates.size(), true), targetStates); - for(auto const& probZeroState : prob01States.first){ - x[probZeroState] = storm::utility::zero(); - } - for(auto const& probZeroState : prob01States.second){ - x[probZeroState] = storm::utility::one(); - } - eqSysMatrix.convertToEquationSystem(); - std::unique_ptr> solver = storm::utility::solver::LinearEquationSolverFactory().create(eqSysMatrix); - solver->solveEquationSystem(x, subB); - } - + void GameSolver::solveGame(OptimizationDirection player1Goal, OptimizationDirection player2Goal, std::vector& x, std::vector const& b) const { // Set up the environment for value iteration. bool converged = false; + uint_fast64_t numberOfPlayer1States = x.size(); std::vector tmpResult(numberOfPlayer1States); std::vector nondetResult(player2Matrix.getRowCount()); std::vector player2Result(player2Matrix.getRowGroupCount()); @@ -145,6 +110,16 @@ namespace storm { } } + template + storm::storage::SparseMatrix const& GameSolver::getPlayer1Matrix() const { + return player1Matrix; + } + + template + storm::storage::SparseMatrix const& GameSolver::getPlayer2Matrix() const { + return player2Matrix; + } + template class GameSolver; } } \ No newline at end of file diff --git a/src/solver/GameSolver.h b/src/solver/GameSolver.h index 06a906836..699533ee3 100644 --- a/src/solver/GameSolver.h +++ b/src/solver/GameSolver.h @@ -38,18 +38,20 @@ namespace storm { GameSolver(storm::storage::SparseMatrix const& player1Matrix, storm::storage::SparseMatrix const& player2Matrix, double precision, uint_fast64_t maximalNumberOfIterations, bool relative); /*! - * Solves the equation system defined by the game matrix. Note that the game matrix has to be given upon + * Solves the equation system defined by the game matrices. Note that the game matrices have to be given upon * construction time of the solver object. * * @param player1Goal Sets whether player 1 wants to minimize or maximize. * @param player2Goal Sets whether player 2 wants to minimize or maximize. - * @param x The initial guess of the solution. + * @param x The initial guess of the solution. For correctness, the guess has to be less (or equal) to the final solution (unless both players minimize) * @param b The vector to add after matrix-vector multiplication. - * @param initialPlayer1Policy A policy that selects rows in every rowgroup of player1. This will be used as an initial guess - * @param initialPlayer2Policy A policy that selects rows in every rowgroup of player2. This will be used as an initial guess * @return The solution vector in the for of the vector x. */ - virtual void solveGame(OptimizationDirection player1Goal, OptimizationDirection player2Goal, std::vector& x, std::vector const& b, std::vector* initialPlayer1Policy = nullptr, std::vector* initialPlayer2Policy = nullptr) const; + virtual void solveGame(OptimizationDirection player1Goal, OptimizationDirection player2Goal, std::vector& x, std::vector const& b) const; + + + storm::storage::SparseMatrix const& getPlayer2Matrix() const; + storm::storage::SparseMatrix const& getPlayer1Matrix() const; private: // The matrix defining the choices of player 1. diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp index a572f3e16..61e2f9f1b 100644 --- a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp +++ b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp @@ -31,7 +31,7 @@ namespace storm { template - void GmmxxMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX, std::vector* initialPolicy) const { + void GmmxxMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) 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; @@ -47,32 +47,6 @@ namespace storm { 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); - std::vector subB(rowGroupIndices.size() - 1); - storm::utility::vector::selectVectorValues(subB, *initialPolicy, rowGroupIndices, b); - //depending on the choice, qualitative properties might have changed. - storm::storage::BitVector targetStates(subB.size(), false); - for(std::size_t index = 0; index < targetStates.size(); ++index){ - if(!storm::utility::isZero(subB[index])){ - targetStates.set(index); - } - } - auto prob01States = storm::utility::graph::performProb01(submatrix.transpose(), storm::storage::BitVector(targetStates.size(), true), targetStates); - for(auto const& probZeroState : prob01States.first){ - (*currentX)[probZeroState] = storm::utility::zero(); - } - for(auto const& probZeroState : prob01States.second){ - (*currentX)[probZeroState] = storm::utility::one(); - } - - submatrix.convertToEquationSystem(); - GmmxxLinearEquationSolver gmmLinearEquationSolver(submatrix); - // Solve the resulting linear equation system - gmmLinearEquationSolver.solveEquationSystem(*currentX, subB); - } - uint_fast64_t iterations = 0; bool converged = false; @@ -127,11 +101,7 @@ namespace storm { } 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. - if(initialPolicy == nullptr){ - this->policy = std::vector(this->A.getRowGroupIndices().size() - 1); - } else { - this->policy = *initialPolicy; - } + this->policy = std::vector(this->A.getRowGroupIndices().size() - 1); // 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 8bf92086d..ed91278bc 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, std::vector* initialPolicy = 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; private: // The (gmm++) matrix associated with this equation solver. diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h index 934aae3f2..7d3bbe287 100644 --- a/src/solver/MinMaxLinearEquationSolver.h +++ b/src/solver/MinMaxLinearEquationSolver.h @@ -97,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, std::vector* initialPolicy = nullptr) const = 0; + virtual void solveEquationSystem(OptimizationDirection d, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const = 0; /*! * As solveEquationSystem with an optimization-direction, but this uses the internally set direction. @@ -137,6 +137,10 @@ namespace storm { void setEarlyTerminationCriterion(std::unique_ptr> v) { earlyTermination = std::move(v); } + + storm::storage::SparseMatrix const& getMatrix() const { + return A; + } protected: diff --git a/src/solver/NativeMinMaxLinearEquationSolver.cpp b/src/solver/NativeMinMaxLinearEquationSolver.cpp index 44a8e7ba2..e2eaa267b 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, std::vector* initialPolicy) const { + void NativeMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) 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; @@ -46,18 +46,6 @@ namespace storm { 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; @@ -112,11 +100,7 @@ namespace storm { } 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. - if(initialPolicy == nullptr){ - this->policy = std::vector(this->A.getRowGroupIndices().size() - 1); - } else { - this->policy = *initialPolicy; - } + this->policy = std::vector(this->A.getRowGroupIndices().size() - 1); // 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 ceaf72fcd..5fcd76aa9 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, std::vector* initialPolicy = 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; }; } // namespace solver diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp index 9fa506bc0..f16c92a07 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, std::vector* initialPolicy) const { + void TopologicalMinMaxLinearEquationSolver::solveEquationSystem(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { #ifdef GPU_USE_FLOAT #define __FORCE_FLOAT_CALCULATION true diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.h b/src/solver/TopologicalMinMaxLinearEquationSolver.h index 314c8ac58..25c325f75 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, std::vector* initialPolicy = 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; private: bool enableCuda; diff --git a/src/utility/policyguessing.cpp b/src/utility/policyguessing.cpp new file mode 100644 index 000000000..69ad63b19 --- /dev/null +++ b/src/utility/policyguessing.cpp @@ -0,0 +1,190 @@ +/* + * File: Regions.cpp + * Author: Tim Quatmann + * + * Created on November 16, 2015, + */ +#include + +#include "src/utility/policyguessing.h" + +#include "src/storage/SparseMatrix.h" +#include "src/utility/macros.h" +#include "src/utility/solver.h" +#include "src/solver/LinearEquationSolver.h" +#include "src/solver/GameSolver.h" +#include "graph.h" + +namespace storm { + namespace utility{ + namespace policyguessing { + + template + void solveGame( storm::solver::GameSolver& solver, + std::vector& x, + std::vector const& b, + OptimizationDirection player1Goal, + OptimizationDirection player2Goal, + std::vector& player1Policy, + std::vector& player2Policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value + ){ + + solveInducedEquationSystem(solver, x, b, player1Policy, player2Policy, targetChoices, prob0Value); + solver.setPolicyTracking(); + solver.solveGame(player1Goal, player2Goal, x, b); + player1Policy = solver.getPlayer1Policy(); + player2Policy = solver.getPlayer2Policy(); + } + + + template + void solveInducedEquationSystem(storm::solver::GameSolver const& solver, + std::vector& x, + std::vector const& b, + std::vector const& player1Policy, + std::vector const& player2Policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value){ + uint_fast64_t numberOfPlayer1States = x.size(); + + //Get the rows of the player2matrix that are selected by the policies + //Note that rows can be selected more then once and in an arbitrary order. + std::vector selectedRows(numberOfPlayer1States); + for (uint_fast64_t pl1State = 0; pl1State < numberOfPlayer1States; ++pl1State){ + auto const& pl1Row = solver.getPlayer1Matrix().getRow(solver.getPlayer1Matrix().getRowGroupIndices()[pl1State] + player1Policy[pl1State]); + STORM_LOG_ASSERT(pl1Row.getNumberOfEntries()==1, ""); + uint_fast64_t pl2State = pl1Row.begin()->getColumn(); + selectedRows[pl1State] = solver.getPlayer2Matrix().getRowGroupIndices()[pl2State] + player2Policy[pl2State]; + } + + //Get the matrix A, vector b, and the targetStates induced by this selection + storm::storage::SparseMatrix inducedA = solver.getPlayer2Matrix().selectRowsFromRowIndexSequence(selectedRows, false); + std::vector inducedB(numberOfPlayer1States); + storm::utility::vector::selectVectorValues(inducedB, selectedRows, b); + storm::storage::BitVector inducedTarget(numberOfPlayer1States, false); + for (uint_fast64_t pl1State = 0; pl1State < numberOfPlayer1States; ++pl1State){ + if(targetChoices.get(selectedRows[pl1State])){ + inducedTarget.set(pl1State); + } + } + + //Find the states from which no target state is reachable. + //Note that depending on the policies, qualitative properties might have changed which makes this step necessary. + storm::storage::BitVector probGreater0States = storm::utility::graph::performProbGreater0(inducedA.transpose(), storm::storage::BitVector(numberOfPlayer1States, true), inducedTarget); + + //Get the final A,x, and b and invoke linear equation solver + storm::storage::SparseMatrix subA = inducedA.getSubmatrix(true, probGreater0States, probGreater0States, true); + subA.convertToEquationSystem(); + std::vector subX(probGreater0States.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(subX, probGreater0States, x); + std::vector subB(probGreater0States.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(subB, probGreater0States, inducedB); + std::unique_ptr> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory().create(subA); + linEqSysSolver->solveEquationSystem(subX, subB); + + //fill in the result + storm::utility::vector::setVectorValues(x, probGreater0States, subX); + storm::utility::vector::setVectorValues(x, (~probGreater0States), prob0Value); + } + + + template + void solveMinMaxLinearEquationSystem( storm::solver::MinMaxLinearEquationSolver& solver, + std::vector& x, + std::vector const& b, + OptimizationDirection goal, + std::vector& policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value + ){ + + solveInducedEquationSystem(solver, x, b, policy, targetChoices, prob0Value); + solver.setPolicyTracking(); + solver.solveEquationSystem(goal, x, b); + policy = solver.getPolicy(); + } + + + template + void solveInducedEquationSystem(storm::solver::MinMaxLinearEquationSolver const& solver, + std::vector& x, + std::vector const& b, + std::vector const& policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value + ){ + uint_fast64_t numberOfStates = x.size(); + + //Get the matrix A, vector b, and the targetStates induced by the policy + storm::storage::SparseMatrix inducedA = solver.getMatrix().selectRowsFromRowGroups(policy, false); + std::vector inducedB(numberOfStates); + storm::utility::vector::selectVectorValues(inducedB, policy, solver.getMatrix().getRowGroupIndices(), b); + storm::storage::BitVector inducedTarget(numberOfStates, false); + for (uint_fast64_t state = 0; state < numberOfStates; ++state){ + if(targetChoices.get(solver.getMatrix().getRowGroupIndices()[state] + policy[state])){ + inducedTarget.set(state); + } + } + + //Find the states from which no target state is reachable. + //Note that depending on the policies, qualitative properties might have changed which makes this step necessary. + storm::storage::BitVector probGreater0States = storm::utility::graph::performProbGreater0(inducedA.transpose(), storm::storage::BitVector(numberOfStates, true), inducedTarget); + + //Get the final A,x, and b and invoke linear equation solver + storm::storage::SparseMatrix subA = inducedA.getSubmatrix(true, probGreater0States, probGreater0States, true); + subA.convertToEquationSystem(); + std::vector subX(probGreater0States.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(subX, probGreater0States, x); + std::vector subB(probGreater0States.getNumberOfSetBits()); + storm::utility::vector::selectVectorValues(subB, probGreater0States, inducedB); + std::unique_ptr> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory().create(subA); + linEqSysSolver->solveEquationSystem(subX, subB); + + //fill in the result + storm::utility::vector::setVectorValues(x, probGreater0States, subX); + storm::utility::vector::setVectorValues(x, (~probGreater0States), prob0Value); + } + + + template void solveGame( storm::solver::GameSolver& solver, + std::vector& x, + std::vector const& b, + OptimizationDirection player1Goal, + OptimizationDirection player2Goal, + std::vector& player1Policy, + std::vector& player2Policy, + storm::storage::BitVector const& targetChoices, + double const& prob0Value + ); + + template void solveInducedEquationSystem(storm::solver::GameSolver const& solver, + std::vector& x, + std::vector const& b, + std::vector const& player1Policy, + std::vector const& player2Policy, + storm::storage::BitVector const& targetChoices, + double const& prob0Value + ); + + template void solveMinMaxLinearEquationSystem( storm::solver::MinMaxLinearEquationSolver& solver, + std::vector& x, + std::vector const& b, + OptimizationDirection goal, + std::vector& policy, + storm::storage::BitVector const& targetChoices, + double const& prob0Value + ); + + template void solveInducedEquationSystem(storm::solver::MinMaxLinearEquationSolver const& solver, + std::vector& x, + std::vector const& b, + std::vector const& policy, + storm::storage::BitVector const& targetChoices, + double const& prob0Value + ); + + } + } +} diff --git a/src/utility/policyguessing.h b/src/utility/policyguessing.h new file mode 100644 index 000000000..a66cb19ca --- /dev/null +++ b/src/utility/policyguessing.h @@ -0,0 +1,154 @@ +/* + * File: Regions.h + * Author: Tim Quatmann + * + * Created on November 16, 2015, + * + * This file provides functions to apply a solver on a nondeterministic model or a two player game. + * However, policies are used to compute an initial guess which will (hopefully) speed up the value iteration techniques. + */ + + +#ifndef STORM_UTILITY_POLICYGUESSING_H +#define STORM_UTILITY_POLICYGUESSING_H + +#include "src/solver/GameSolver.h" +#include "src/solver/MinMaxLinearEquationSolver.h" +#include "src/solver/OptimizationDirection.h" +#include "src/utility/vector.h" +#include "src/storage/BitVector.h" +#include "src/storage/sparse/StateType.h" + +namespace storm { + namespace utility{ + namespace policyguessing { + + /*! + * invokes the given game solver. + * + * The given policies for player 1 and player 2 will serve as initial guess. + * A linear equation system defined by the induced Matrix A and vector b is solved before + * solving the actual game. + * Note that, depending on the policies, the qualitative properties of the graph defined by A + * might be different to the original graph of the game. + * To ensure a unique solution, we need to filter out the "prob0"-states. + * To identify these states and set the result for them correctly, it is necessary to know whether rewards or probabilities are to be computed + * + * @param solver the solver to be invoked + * @param player1Goal Sets whether player 1 wants to minimize or maximize. + * @param player2Goal Sets whether player 2 wants to minimize or maximize. + * @param x The initial guess of the solution. + * @param b The vector to add after matrix-vector multiplication. + * @param player1Policy A policy that selects rows in every rowgroup of player1. This will be used as an initial guess + * @param player2Policy A policy that selects rows in every rowgroup of player2. This will be used as an initial guess + * @param targetChoices marks the choices in the player2 matrix that have a positive probability to lead to a target state + * @param prob0Value the value that, after policy instantiation, is assigned to the states that have probability zero to reach a target + * @return The solution vector in the form of the vector x as well as the two policies. + */ + template + void solveGame( storm::solver::GameSolver& solver, + std::vector& x, + std::vector const& b, + OptimizationDirection player1Goal, + OptimizationDirection player2Goal, + std::vector& player1Policy, + std::vector& player2Policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value + ); + + + /*! + * Solves the equation system defined by the matrix A and vector b' that result from applying + * the given policies to the matrices of the two players and the given b. + * + * Note that, depending on the policies, the qualitative properties of the graph defined by A + * might be different to the original graph of the game. + * To ensure a unique solution, we need to filter out the "prob0"-states. + * (Notice that new "prob1"-states do not harm as actual target states are always excluded, independent of the choice) + * + * @param solver the solver that contains the two player matrices + * @param x The initial guess of the solution. + * @param b The vector in which to select the entries of the right hand side + * @param player1Policy A policy that selects rows in every rowgroup of player1. + * @param player2Policy A policy that selects rows in every rowgroup of player2. + * @param targetChoices marks the choices in the player2 matrix that have a positive probability to lead to a target state + * @param prob0Value the value that is assigned to the states that have probability zero to reach a target + * @return The solution vector in the form of the vector x. + */ + template + void solveInducedEquationSystem(storm::solver::GameSolver const& solver, + std::vector& x, + std::vector const& b, + std::vector const& player1Policy, + std::vector const& player2Policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value + ); + + + + /*! + * invokes the given MinMaxLinearEquationSolver. + * + * The given policy will serve as an initial guess. + * A linear equation system defined by the induced Matrix A and vector b is solved before + * solving the actual MinMax equation system. + * Note that, depending on the policy, the qualitative properties of the graph defined by A + * might be different to the original graph. + * To ensure a unique solution, we need to filter out the "prob0"-states. + * To identify these states and set the result for them correctly, it is necessary to know whether rewards or probabilities are to be computed + * + + * @param solver the solver that contains the two player matrices + * @param x The initial guess of the solution. + * @param b The vector to add after matrix-vector multiplication. + * @param goal Sets whether we want to minimize or maximize. + * @param policy A policy that selects rows in every rowgroup. + * @param targetChoices marks the rows in the matrix that have a positive probability to lead to a target state + * @param prob0Value the value that is assigned to the states that have probability zero to reach a target + * @return The solution vector in the form of the vector x. + */ + template + void solveMinMaxLinearEquationSystem( storm::solver::MinMaxLinearEquationSolver& solver, + std::vector& x, + std::vector const& b, + OptimizationDirection goal, + std::vector& policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value + ); + + /*! + * Solves the equation system defined by the matrix A and vector b' that result from applying + * the given policy to the matrices of the two players and the given b. + * + * Note that, depending on the policy, the qualitative properties of the graph defined by A + * might be different to the original graph + * To ensure a unique solution, we need to filter out the "prob0"-states. + * (Notice that new "prob1"-states do not harm as actual target states are always excluded, independent of the choice) + * + * @param solverthe solver that contains the two player matrices + * @param x The initial guess of the solution. + * @param b The vector in which to select the entries of the right hand side + * @param policy A policy that selects rows in every rowgroup. + * @param targetChoices marks the rows in the matrix that have a positive probability to lead to a target state + * @param prob0Value the value that is assigned to the states that have probability zero to reach a target + * @return The solution vector in the form of the vector x. + */ + template + void solveInducedEquationSystem(storm::solver::MinMaxLinearEquationSolver const& solver, + std::vector& x, + std::vector const& b, + std::vector const& policy, + storm::storage::BitVector const& targetChoices, + ValueType const& prob0Value + ); + + } + } +} + + +#endif /* STORM_UTILITY_REGIONS_H */ +