refactored the quick implementation of the recent bugfix

Former-commit-id: 5c0d1fa3b9

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>();
+                ConstantType dummyNonZeroValue = storm::utility::one<ConstantType>();
                 storm::storage::SparseMatrixBuilder<ConstantType> 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<ConstantType>(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<ParametricType, CoefficientType>(storm::utility::zero<CoefficientType>());
                         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<ConstantType>(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<ConstantType>(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>();
+                ConstantType dummyNonZeroValue = storm::utility::one<ConstantType>();
                 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<VariableType, RegionBoundary>::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<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(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<double>() : storm::utility::zero<double>()));
             }
             
             template<>
             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->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<double>() : storm::utility::zero<double>()));
             }
             
                 

src/modelchecker/region/ApproximationModel.h

@@ -117,6 +117,7 @@ namespace storm {
                     storm::storage::SparseMatrix<ConstantType> matrix; //The matrix itself.
                     std::vector<std::pair<typename storm::storage::SparseMatrix<ConstantType>::iterator, ConstantType&>> assignment; // Connection of matrix entries with placeholders
                     std::vector<std::size_t> 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<ConstantType> vector; //The vector itself.

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<ConstantType>(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<storm::models::sparse::Mdp<storm::RationalFunction>, double>::invokeSolver(){
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<double>> solver = storm::solver::configureMinMaxLinearEquationSolver(this->solverData.solveGoal, storm::utility::solver::MinMaxLinearEquationSolverFactory<double>(), 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<double>() : storm::utility::zero<double>()));
             }
             
             

src/modelchecker/region/SamplingModel.h

@@ -89,6 +89,7 @@ namespace storm {
                 struct MatrixData {
                     storm::storage::SparseMatrix<ConstantType> matrix; //The matrix itself.
                     std::vector<std::pair<typename storm::storage::SparseMatrix<ConstantType>::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<ConstantType> vector; //The vector itself.

src/solver/GameSolver.cpp

@@ -19,45 +19,10 @@ namespace storm {
         }
 
         template <typename ValueType>
-        void GameSolver<ValueType>::solveGame(OptimizationDirection player1Goal, OptimizationDirection player2Goal, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<storm::storage::sparse::state_type>* initialPlayer1Policy, std::vector<storm::storage::sparse::state_type>* 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<storm::storage::sparse::state_type> 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<ValueType> eqSysMatrix = player2Matrix.selectRowsFromRowIndexSequence(selectedRows, true);
-                std::vector<ValueType> subB(numberOfPlayer1States);
-                storm::utility::vector::selectVectorValues<ValueType>(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<ValueType>();
-                }
-                for(auto const& probZeroState : prob01States.second){
-                    x[probZeroState] = storm::utility::one<ValueType>();
-                }
-                eqSysMatrix.convertToEquationSystem();
-                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = storm::utility::solver::LinearEquationSolverFactory<ValueType>().create(eqSysMatrix);
-                solver->solveEquationSystem(x, subB);
-            }
-            
+        void GameSolver<ValueType>::solveGame(OptimizationDirection player1Goal, OptimizationDirection player2Goal, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             // Set up the environment for value iteration.
             bool converged = false;
+            uint_fast64_t numberOfPlayer1States = x.size();
             std::vector<ValueType> tmpResult(numberOfPlayer1States);
             std::vector<ValueType> nondetResult(player2Matrix.getRowCount());
             std::vector<ValueType> player2Result(player2Matrix.getRowGroupCount());
@@ -145,6 +110,16 @@ namespace storm {
             }
         }
 
+        template <typename ValueType>
+        storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& GameSolver<ValueType>::getPlayer1Matrix() const {
+            return player1Matrix;
+        }
+        
+        template <typename ValueType>
+        storm::storage::SparseMatrix<ValueType> const& GameSolver<ValueType>::getPlayer2Matrix() const {
+            return player2Matrix;
+        }
+
         template class GameSolver<double>;
     }
 }

src/solver/GameSolver.h

@@ -38,18 +38,20 @@ namespace storm {
             GameSolver(storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& player1Matrix, storm::storage::SparseMatrix<ValueType> 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<ValueType>& x, std::vector<ValueType> const& b, std::vector<storm::storage::sparse::state_type>* initialPlayer1Policy = nullptr, std::vector<storm::storage::sparse::state_type>* initialPlayer2Policy = nullptr) const;
+            virtual void solveGame(OptimizationDirection player1Goal, OptimizationDirection player2Goal, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            
+            
+            storm::storage::SparseMatrix<ValueType> const& getPlayer2Matrix() const;
+            storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& getPlayer1Matrix() const;
 
         private:
             // The matrix defining the choices of player 1.

src/solver/GmmxxMinMaxLinearEquationSolver.cpp

@@ -31,7 +31,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, std::vector<storm::storage::sparse::state_type>* initialPolicy) const {
+        void GmmxxMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* 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<ValueType> submatrix = this->A.selectRowsFromRowGroups(*initialPolicy, true);
-                                    std::vector<ValueType> subB(rowGroupIndices.size() - 1);
-                                    storm::utility::vector::selectVectorValues<ValueType>(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<ValueType>();
-                                    }
-                                    for(auto const& probZeroState : prob01States.second){
-                                        (*currentX)[probZeroState] = storm::utility::one<ValueType>();
-                                    }
-                                    
-                                    submatrix.convertToEquationSystem();
-                                    GmmxxLinearEquationSolver<ValueType> 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<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
-                                } else {
-                                    this->policy = *initialPolicy;
-                                }
+                                this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
 
 				// Create our own multiplyResult for solving the deterministic sub-instances.
 				std::vector<ValueType> deterministicMultiplyResult(rowGroupIndices.size() - 1);

src/solver/GmmxxMinMaxLinearEquationSolver.h

@@ -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, std::vector<storm::storage::sparse::state_type>* initialPolicy = 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;
             
         private:
             // The (gmm++) matrix associated with this equation solver.

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<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;
+            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;
             
             /*!
              * As solveEquationSystem with an optimization-direction, but this uses the internally set direction.
@@ -137,6 +137,10 @@ namespace storm {
             void setEarlyTerminationCriterion(std::unique_ptr<AllowEarlyTerminationCondition<ValueType>> v) {
                 earlyTermination = std::move(v);
             }
+
+            storm::storage::SparseMatrix<ValueType> const& getMatrix() const {
+                return A;
+            }
             
             
         protected:

src/solver/NativeMinMaxLinearEquationSolver.cpp

@@ -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, std::vector<storm::storage::sparse::state_type>* initialPolicy) const {
+        void NativeMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* 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<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;
 
@@ -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<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
-                                } else {
-                                    this->policy = *initialPolicy;
-                                }
+                                this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
 
 				// Create our own multiplyResult for solving the deterministic sub-instances.
 				std::vector<ValueType> deterministicMultiplyResult(this->A.getRowGroupIndices().size() - 1);

src/solver/NativeMinMaxLinearEquationSolver.h

@@ -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, std::vector<storm::storage::sparse::state_type>* initialPolicy = 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;
 
         };
     } // namespace solver

src/solver/TopologicalMinMaxLinearEquationSolver.cpp

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

src/solver/TopologicalMinMaxLinearEquationSolver.h

@@ -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, std::vector<storm::storage::sparse::state_type>* initialPolicy = 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;
 		private:
 
 			bool enableCuda;

src/utility/policyguessing.cpp

@@ -0,0 +1,190 @@
+/* 
+ * File:   Regions.cpp
+ * Author: Tim Quatmann
+ * 
+ * Created on November 16, 2015,
+ */
+#include <stdint.h>
+
+#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 <typename ValueType>
+            void solveGame( storm::solver::GameSolver<ValueType>& solver,
+                            std::vector<ValueType>& x,
+                            std::vector<ValueType> const& b,
+                            OptimizationDirection player1Goal,
+                            OptimizationDirection player2Goal,
+                            std::vector<storm::storage::sparse::state_type>& player1Policy, 
+                            std::vector<storm::storage::sparse::state_type>& 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 <typename ValueType>
+            void solveInducedEquationSystem(storm::solver::GameSolver<ValueType> const& solver,
+                                            std::vector<ValueType>& x,
+                                            std::vector<ValueType> const& b,
+                                            std::vector<storm::storage::sparse::state_type> const& player1Policy, 
+                                            std::vector<storm::storage::sparse::state_type> 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<storm::storage::sparse::state_type> 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<ValueType> inducedA = solver.getPlayer2Matrix().selectRowsFromRowIndexSequence(selectedRows, false);
+                std::vector<ValueType> inducedB(numberOfPlayer1States);
+                storm::utility::vector::selectVectorValues<ValueType>(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<ValueType> subA = inducedA.getSubmatrix(true, probGreater0States, probGreater0States, true);
+                subA.convertToEquationSystem();
+                std::vector<ValueType> subX(probGreater0States.getNumberOfSetBits());
+                storm::utility::vector::selectVectorValues(subX, probGreater0States, x);
+                std::vector<ValueType> subB(probGreater0States.getNumberOfSetBits());
+                storm::utility::vector::selectVectorValues(subB, probGreater0States, inducedB);
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory<ValueType>().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 <typename ValueType>
+            void solveMinMaxLinearEquationSystem( storm::solver::MinMaxLinearEquationSolver<ValueType>& solver,
+                                    std::vector<ValueType>& x,
+                                    std::vector<ValueType> const& b,
+                                    OptimizationDirection goal,
+                                    std::vector<storm::storage::sparse::state_type>& 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 <typename ValueType>
+            void solveInducedEquationSystem(storm::solver::MinMaxLinearEquationSolver<ValueType> const& solver,
+                                            std::vector<ValueType>& x,
+                                            std::vector<ValueType> const& b,
+                                            std::vector<storm::storage::sparse::state_type> 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<ValueType> inducedA = solver.getMatrix().selectRowsFromRowGroups(policy, false);
+                std::vector<ValueType> inducedB(numberOfStates);
+                storm::utility::vector::selectVectorValues<ValueType>(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<ValueType> subA = inducedA.getSubmatrix(true, probGreater0States, probGreater0States, true);
+                subA.convertToEquationSystem();
+                std::vector<ValueType> subX(probGreater0States.getNumberOfSetBits());
+                storm::utility::vector::selectVectorValues(subX, probGreater0States, x);
+                std::vector<ValueType> subB(probGreater0States.getNumberOfSetBits());
+                storm::utility::vector::selectVectorValues(subB, probGreater0States, inducedB);
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linEqSysSolver = storm::utility::solver::LinearEquationSolverFactory<ValueType>().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<double>( storm::solver::GameSolver<double>& solver,
+                                std::vector<double>& x,
+                                std::vector<double> const& b,
+                                OptimizationDirection player1Goal,
+                                OptimizationDirection player2Goal,
+                                std::vector<storm::storage::sparse::state_type>& player1Policy, 
+                                std::vector<storm::storage::sparse::state_type>& player2Policy,
+                                storm::storage::BitVector const& targetChoices,
+                                double const& prob0Value
+                            );
+            
+            template void solveInducedEquationSystem<double>(storm::solver::GameSolver<double> const& solver,
+                                std::vector<double>& x,
+                                std::vector<double> const& b,
+                                std::vector<storm::storage::sparse::state_type> const& player1Policy, 
+                                std::vector<storm::storage::sparse::state_type> const& player2Policy,
+                                storm::storage::BitVector const& targetChoices,
+                                double const& prob0Value
+                            );
+            
+            template void solveMinMaxLinearEquationSystem<double>( storm::solver::MinMaxLinearEquationSolver<double>& solver,
+                                std::vector<double>& x,
+                                std::vector<double> const& b,
+                                OptimizationDirection goal,
+                                std::vector<storm::storage::sparse::state_type>& policy,
+                                storm::storage::BitVector const& targetChoices,
+                                double const& prob0Value
+                            );
+            
+            template void solveInducedEquationSystem<double>(storm::solver::MinMaxLinearEquationSolver<double> const& solver,
+                                std::vector<double>& x,
+                                std::vector<double> const& b,
+                                std::vector<storm::storage::sparse::state_type> const& policy,
+                                storm::storage::BitVector const& targetChoices,
+                                double const& prob0Value
+                            );
+                
+        }
+    }
+}

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<typename ValueType>
+            void solveGame( storm::solver::GameSolver<ValueType>& solver,
+                                    std::vector<ValueType>& x,
+                                    std::vector<ValueType> const& b,
+                                    OptimizationDirection player1Goal,
+                                    OptimizationDirection player2Goal,
+                                    std::vector<storm::storage::sparse::state_type>& player1Policy,
+                                    std::vector<storm::storage::sparse::state_type>& 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<typename ValueType>
+            void solveInducedEquationSystem(storm::solver::GameSolver<ValueType> const& solver,
+                                                    std::vector<ValueType>& x,
+                                                    std::vector<ValueType> const& b,
+                                                    std::vector<storm::storage::sparse::state_type> const& player1Policy,
+                                                    std::vector<storm::storage::sparse::state_type> 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<typename ValueType>
+            void solveMinMaxLinearEquationSystem( storm::solver::MinMaxLinearEquationSolver<ValueType>& solver,
+                                    std::vector<ValueType>& x,
+                                    std::vector<ValueType> const& b,
+                                    OptimizationDirection goal,
+                                    std::vector<storm::storage::sparse::state_type>& 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<typename ValueType>
+            void solveInducedEquationSystem(storm::solver::MinMaxLinearEquationSolver<ValueType> const& solver,
+                                            std::vector<ValueType>& x,
+                                            std::vector<ValueType> const& b,
+                                            std::vector<storm::storage::sparse::state_type> const& policy,
+                                            storm::storage::BitVector const& targetChoices,
+                                            ValueType const& prob0Value
+                                        );
+            
+        }
+    }
+}
+
+
+#endif	/* STORM_UTILITY_REGIONS_H */
+