Parameter lifting for mdps, some fixes

src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.cpp

@@ -195,7 +195,20 @@ namespace storm {
                 }
                 return std::make_unique<storm::modelchecker::ExplicitQuantitativeCheckResult<ConstantType>>(std::move(result));
             }
-
+            
+                
+            template <typename SparseModelType, typename ConstantType>
+            void SparseDtmcParameterLiftingModelChecker<SparseModelType, ConstantType>::reset() {
+                maybeStates.resize(0);
+                resultsForNonMaybeStates.clear();
+                stepBound = boost::none;
+                parameterLifter = nullptr;
+                minSched = boost::none;
+                maxSched = boost::none;
+                x.clear();
+                lowerResultBound = boost::none;
+                upperResultBound = boost::none;
+            }
     
             template class SparseDtmcParameterLiftingModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>, double>;
 

src/storm/modelchecker/parametric/SparseDtmcParameterLiftingModelChecker.h

@@ -6,6 +6,7 @@
 
 #include "storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h"
 #include "storm/storage/BitVector.h"
+#include "storm/storage/TotalScheduler.h"
 #include "storm/solver/MinMaxLinearEquationSolver.h"
 #include "storm/transformer/ParameterLifter.h"
 
@@ -28,7 +29,8 @@ namespace storm {
                 
                 virtual std::unique_ptr<CheckResult> computeQuantitativeValues(ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dirForParameters) override;
                 
-
+                virtual void reset() override;
+                
             private:
                 storm::storage::BitVector maybeStates;
                 std::vector<ConstantType> resultsForNonMaybeStates;

src/storm/modelchecker/parametric/SparseMdpParameterLiftingModelChecker.cpp

@@ -0,0 +1,273 @@
+#include "SparseMdpParameterLiftingModelChecker.h"
+
+#include "storm/adapters/CarlAdapter.h"
+#include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h"
+#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/StandardRewardModel.h"
+#include "storm/utility/vector.h"
+#include "storm/utility/graph.h"
+#include "storm/solver/GameSolver.h"
+
+#include "storm/exceptions/InvalidArgumentException.h"
+#include "storm/exceptions/InvalidPropertyException.h"
+#include "storm/exceptions/NotSupportedException.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace parametric {
+            
+            template <typename SparseModelType, typename ConstantType>
+            SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::SparseMdpParameterLiftingModelChecker(SparseModelType const& parametricModel) : SparseParameterLiftingModelChecker<SparseModelType, ConstantType>(parametricModel), solverFactory(std::make_unique<storm::utility::solver::GameSolverFactory<ConstantType>>()) {
+            }
+            
+            template <typename SparseModelType, typename ConstantType>
+            SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::SparseMdpParameterLiftingModelChecker(SparseModelType const& parametricModel, std::unique_ptr<storm::utility::solver::GameSolverFactory<ConstantType>>&& solverFactory) : SparseParameterLiftingModelChecker<SparseModelType, ConstantType>(parametricModel), solverFactory(std::move(solverFactory)) {
+            }
+    
+            template <typename SparseModelType, typename ConstantType>
+            void SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyBoundedUntilFormula(CheckTask<logic::BoundedUntilFormula, ConstantType> const& checkTask) {
+                
+                // get the step bound
+                STORM_LOG_THROW(!checkTask.getFormula().hasLowerBound(), storm::exceptions::NotSupportedException, "Lower step bounds are not supported.");
+                STORM_LOG_THROW(checkTask.getFormula().hasUpperBound(), storm::exceptions::NotSupportedException, "Expected a bounded until formula with an upper bound.");
+                STORM_LOG_THROW(checkTask.getFormula().isStepBounded(), storm::exceptions::NotSupportedException, "Expected a bounded until formula with step bounds.");
+                stepBound = checkTask.getFormula().getUpperBound().evaluateAsInt();
+                STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Can not apply parameter lifting on step bounded formula: The step bound has to be positive.");
+                if (checkTask.getFormula().isUpperBoundStrict()) {
+                    STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Expected a strict upper step bound that is greater than zero.");
+                    --(*stepBound);
+                }
+                STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Can not apply parameter lifting on step bounded formula: The step bound has to be positive.");
+
+                // get the results for the subformulas
+                storm::modelchecker::SparsePropositionalModelChecker<SparseModelType> propositionalChecker(this->parametricModel);
+                STORM_LOG_THROW(propositionalChecker.canHandle(checkTask.getFormula().getLeftSubformula()) && propositionalChecker.canHandle(checkTask.getFormula().getRightSubformula()), storm::exceptions::NotSupportedException, "Parameter lifting with non-propositional subformulas is not supported");
+                storm::storage::BitVector phiStates = std::move(propositionalChecker.check(checkTask.getFormula().getLeftSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+                storm::storage::BitVector psiStates = std::move(propositionalChecker.check(checkTask.getFormula().getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+                
+                // get the maybeStates
+                maybeStates = storm::solver::minimize(checkTask.getOptimizationDirection()) ?
+                              storm::utility::graph::performProbGreater0A(this->parametricModel.getTransitionMatrix(), this->parametricModel.getTransitionMatrix().getRowGroupIndices(), this->parametricModel.getBackwardTransitions(), phiStates, psiStates, true, *stepBound) :
+                              storm::utility::graph::performProbGreater0E(this->parametricModel.getBackwardTransitions(), phiStates, psiStates, true, *stepBound);
+                maybeStates &= ~psiStates;
+                
+                // set the result for all non-maybe states
+                resultsForNonMaybeStates = std::vector<ConstantType>(this->parametricModel.getNumberOfStates(), storm::utility::zero<ConstantType>());
+                storm::utility::vector::setVectorValues(resultsForNonMaybeStates, psiStates, storm::utility::one<ConstantType>());
+                
+                // if there are maybestates, create the parameterLifter
+                if (!maybeStates.empty()) {
+                    // Create the vector of one-step probabilities to go to target states.
+                    std::vector<typename SparseModelType::ValueType> b = this->parametricModel.getTransitionMatrix().getConstrainedRowSumVector(storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getRowCount(), true), psiStates);
+                    
+                    parameterLifter = std::make_unique<storm::transformer::ParameterLifter<typename SparseModelType::ValueType, ConstantType>>(this->parametricModel.getTransitionMatrix(), b, this->parametricModel.getTransitionMatrix().getRowIndicesOfRowGroups(maybeStates), maybeStates);
+                    computePlayer1Matrix();
+                    
+                    applyPreviousResultAsHint = false;
+                }
+                
+                // We know some bounds for the results
+                lowerResultBound = storm::utility::zero<ConstantType>();
+                upperResultBound = storm::utility::one<ConstantType>();
+            }
+    
+            template <typename SparseModelType, typename ConstantType>
+            void SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyUntilFormula(CheckTask<logic::UntilFormula, ConstantType> const& checkTask) {
+                
+                // get the results for the subformulas
+                storm::modelchecker::SparsePropositionalModelChecker<SparseModelType> propositionalChecker(this->parametricModel);
+                STORM_LOG_THROW(propositionalChecker.canHandle(checkTask.getFormula().getLeftSubformula()) && propositionalChecker.canHandle(checkTask.getFormula().getRightSubformula()), storm::exceptions::NotSupportedException, "Parameter lifting with non-propositional subformulas is not supported");
+                storm::storage::BitVector phiStates = std::move(propositionalChecker.check(checkTask.getFormula().getLeftSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+                storm::storage::BitVector psiStates = std::move(propositionalChecker.check(checkTask.getFormula().getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+                
+                // get the maybeStates
+                std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::solver::minimize(checkTask.getOptimizationDirection()) ?
+                                                                                                          storm::utility::graph::performProb01Min(this->parametricModel.getTransitionMatrix(), this->parametricModel.getTransitionMatrix().getRowGroupIndices(), this->parametricModel.getBackwardTransitions(), phiStates, psiStates) :
+                                                                                                          storm::utility::graph::performProb01Max(this->parametricModel.getTransitionMatrix(), this->parametricModel.getTransitionMatrix().getRowGroupIndices(), this->parametricModel.getBackwardTransitions(), phiStates, psiStates);
+                maybeStates = ~(statesWithProbability01.first | statesWithProbability01.second);
+                
+                // set the result for all non-maybe states
+                resultsForNonMaybeStates = std::vector<ConstantType>(this->parametricModel.getNumberOfStates(), storm::utility::zero<ConstantType>());
+                storm::utility::vector::setVectorValues(resultsForNonMaybeStates, statesWithProbability01.second, storm::utility::one<ConstantType>());
+                
+                // if there are maybestates, create the parameterLifter
+                if (!maybeStates.empty()) {
+                    // Create the vector of one-step probabilities to go to target states.
+                    std::vector<typename SparseModelType::ValueType> b = this->parametricModel.getTransitionMatrix().getConstrainedRowSumVector(storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getRowCount(), true), psiStates);
+                    
+                    parameterLifter = std::make_unique<storm::transformer::ParameterLifter<typename SparseModelType::ValueType, ConstantType>>(this->parametricModel.getTransitionMatrix(), b, this->parametricModel.getTransitionMatrix().getRowIndicesOfRowGroups(maybeStates), maybeStates);
+                    computePlayer1Matrix();
+                    
+                    // Check whether there is an EC consisting of maybestates
+                    applyPreviousResultAsHint = storm::solver::minimize(checkTask.getOptimizationDirection()) || // when minimizing, there can not be an EC within the maybestates
+                                                storm::utility::graph::performProb1A(this->parametricModel.getTransitionMatrix(), this->parametricModel.getTransitionMatrix().getRowGroupIndices(), this->parametricModel.getBackwardTransitions(), maybeStates, ~maybeStates).full();
+                }
+                                
+                // We know some bounds for the results
+                lowerResultBound = storm::utility::zero<ConstantType>();
+                upperResultBound = storm::utility::one<ConstantType>();
+            }
+    
+            template <typename SparseModelType, typename ConstantType>
+            void SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyReachabilityRewardFormula(CheckTask<logic::EventuallyFormula, ConstantType> const& checkTask) {
+                
+                // get the results for the subformula
+                storm::modelchecker::SparsePropositionalModelChecker<SparseModelType> propositionalChecker(this->parametricModel);
+                STORM_LOG_THROW(propositionalChecker.canHandle(checkTask.getFormula().getSubformula()), storm::exceptions::NotSupportedException, "Parameter lifting with non-propositional subformulas is not supported");
+                storm::storage::BitVector targetStates = std::move(propositionalChecker.check(checkTask.getFormula().getSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+                
+                // get the maybeStates
+                storm::storage::BitVector infinityStates = storm::solver::minimize(checkTask.getOptimizationDirection()) ?
+                                                           storm::utility::graph::performProb1E(this->parametricModel.getTransitionMatrix(), this->parametricModel.getTransitionMatrix().getRowGroupIndices(), this->parametricModel.getBackwardTransitions(), storm::storage::BitVector(this->parametricModel.getNumberOfStates(), true), targetStates) :
+                                                           storm::utility::graph::performProb1A(this->parametricModel.getTransitionMatrix(), this->parametricModel.getTransitionMatrix().getRowGroupIndices(), this->parametricModel.getBackwardTransitions(), storm::storage::BitVector(this->parametricModel.getNumberOfStates(), true), targetStates);
+                infinityStates.complement();
+                maybeStates = ~(targetStates | infinityStates);
+                
+                // set the result for all the non-maybe states
+                resultsForNonMaybeStates = std::vector<ConstantType>(this->parametricModel.getNumberOfStates(), storm::utility::zero<ConstantType>());
+                storm::utility::vector::setVectorValues(resultsForNonMaybeStates, infinityStates, storm::utility::infinity<ConstantType>());
+                
+                // if there are maybestates, create the parameterLifter
+                if (!maybeStates.empty()) {
+                    // Create the reward vector
+                    STORM_LOG_THROW((checkTask.isRewardModelSet() && this->parametricModel.hasRewardModel(checkTask.getRewardModel())) || (!checkTask.isRewardModelSet() && this->parametricModel.hasUniqueRewardModel()), storm::exceptions::InvalidPropertyException, "The reward model specified by the CheckTask is not available in the given model.");
+
+                    typename SparseModelType::RewardModelType const& rewardModel = checkTask.isRewardModelSet() ? this->parametricModel.getRewardModel(checkTask.getRewardModel()) : this->parametricModel.getUniqueRewardModel();
+
+                    std::vector<typename SparseModelType::ValueType> b = rewardModel.getTotalRewardVector(this->parametricModel.getTransitionMatrix());
+                    
+                    // We need to handle choices that lead to an infinity state.
+                    // As a maybeState does not have reward infinity, such a choice will never be picked. Hence, we can unselect the corresponding rows
+                    storm::storage::BitVector selectedRows = this->parametricModel.getTransitionMatrix().getRowIndicesOfRowGroups(maybeStates);
+                    for (uint_fast64_t row : selectedRows) {
+                        for (auto const& element : this->parametricModel.getTransitionMatrix().getRow(row)) {
+                            if (infinityStates.get(element.getColumn())) {
+                                selectedRows.set(row, false);
+                                break;
+                            }
+                        }
+                    }
+                    
+                    parameterLifter = std::make_unique<storm::transformer::ParameterLifter<typename SparseModelType::ValueType, ConstantType>>(this->parametricModel.getTransitionMatrix(), b, selectedRows, maybeStates);
+                    computePlayer1Matrix();
+                    
+                    // Check whether there is an EC consisting of maybestates
+                    applyPreviousResultAsHint = !storm::solver::minimize(checkTask.getOptimizationDirection()) || // when maximizing, there can not be an EC within the maybestates
+                                                storm::utility::graph::performProb1A(this->parametricModel.getTransitionMatrix(), this->parametricModel.getTransitionMatrix().getRowGroupIndices(), this->parametricModel.getBackwardTransitions(), maybeStates, ~maybeStates).full();
+                }
+                
+                // We only know a lower bound for the result
+                lowerResultBound = storm::utility::zero<ConstantType>();
+
+            }
+    
+            template <typename SparseModelType, typename ConstantType>
+            void SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyCumulativeRewardFormula(CheckTask<logic::CumulativeRewardFormula, ConstantType> const& checkTask) {
+                
+                // Obtain the stepBound
+                stepBound = checkTask.getFormula().getBound().evaluateAsInt();
+                if (checkTask.getFormula().isBoundStrict()) {
+                    STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Expected a strict upper step bound that is greater than zero.");
+                    --(*stepBound);
+                }
+                STORM_LOG_THROW(*stepBound > 0, storm::exceptions::NotSupportedException, "Can not apply parameter lifting on step bounded formula: The step bound has to be positive.");
+                
+                // Create the reward vector
+                STORM_LOG_THROW((checkTask.isRewardModelSet() && this->parametricModel.hasRewardModel(checkTask.getRewardModel())) || (!checkTask.isRewardModelSet() && this->parametricModel.hasUniqueRewardModel()), storm::exceptions::InvalidPropertyException, "The reward model specified by the CheckTask is not available in the given model.");
+                typename SparseModelType::RewardModelType const& rewardModel = checkTask.isRewardModelSet() ? this->parametricModel.getRewardModel(checkTask.getRewardModel()) : this->parametricModel.getUniqueRewardModel();
+                std::vector<typename SparseModelType::ValueType> b = rewardModel.getTotalRewardVector(this->parametricModel.getTransitionMatrix());
+                    
+                parameterLifter = std::make_unique<storm::transformer::ParameterLifter<typename SparseModelType::ValueType, ConstantType>>(this->parametricModel.getTransitionMatrix(), b, storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getRowCount(), true), storm::storage::BitVector(this->parametricModel.getTransitionMatrix().getColumnCount(), true));
+                
+                applyPreviousResultAsHint = false;
+                
+                // We only know a lower bound for the result
+                lowerResultBound = storm::utility::zero<ConstantType>();
+            }
+            
+            template <typename SparseModelType, typename ConstantType>
+            std::unique_ptr<CheckResult> SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::computeQuantitativeValues(ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dirForParameters) {
+                
+                if(maybeStates.empty()) {
+                    return std::make_unique<storm::modelchecker::ExplicitQuantitativeCheckResult<ConstantType>>(resultsForNonMaybeStates);
+                }
+                
+                parameterLifter->specifyRegion(region, dirForParameters);
+                
+                // Set up the solver
+                auto solver = solverFactory->create(player1Matrix, parameterLifter->getMatrix());
+                solver->setTrackScheduler(true);
+                if(lowerResultBound) solver->setLowerBound(lowerResultBound.get());
+                if(upperResultBound) solver->setUpperBound(upperResultBound.get());
+                if(applyPreviousResultAsHint) {
+                    x.resize(maybeStates.getNumberOfSetBits(), storm::utility::zero<ConstantType>());
+                    if(storm::solver::minimize(dirForParameters) && minSched && player1Sched) solver->setSchedulerHint(std::move(player1Sched.get()), std::move(minSched.get()));
+                    if(storm::solver::maximize(dirForParameters) && maxSched && player1Sched) solver->setSchedulerHint(std::move(player1Sched.get()), std::move(maxSched.get()));
+                } else {
+                    x.assign(maybeStates.getNumberOfSetBits(), storm::utility::zero<ConstantType>());
+                }
+                
+                // Invoke the solver
+                if(stepBound) {
+                    assert(*stepBound > 0);
+                    solver->repeatedMultiply(this->currentCheckTask->getOptimizationDirection(), dirForParameters, x, &parameterLifter->getVector(), *stepBound);
+                } else {
+                    solver->solveGame(this->currentCheckTask->getOptimizationDirection(), dirForParameters, x, parameterLifter->getVector());
+                    if(storm::solver::minimize(dirForParameters)) {
+                        minSched = solver->getPlayer2Scheduler();
+                    } else {
+                        maxSched = solver->getPlayer2Scheduler();
+                    }
+                    player1Sched = solver->getPlayer1Scheduler();
+                }
+                
+                // Get the result for the complete model (including maybestates)
+                std::vector<ConstantType> result = resultsForNonMaybeStates;
+                auto maybeStateResIt = x.begin();
+                for(auto const& maybeState : maybeStates) {
+                    result[maybeState] = *maybeStateResIt;
+                    ++maybeStateResIt;
+                }
+                return std::make_unique<storm::modelchecker::ExplicitQuantitativeCheckResult<ConstantType>>(std::move(result));
+            }
+            
+            template <typename SparseModelType, typename ConstantType>
+            void SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::computePlayer1Matrix() {
+                uint_fast64_t n = 0;
+                for(auto const& maybeState : maybeStates) {
+                    n += this->parametricModel.getTransitionMatrix().getRowGroupSize(maybeState);
+                }
+                // The player 1 matrix is the identity matrix of size n with the row groups as given by the original matrix
+                storm::storage::SparseMatrixBuilder<storm::storage::sparse::state_type> matrixBuilder(n, n, n, true, true, maybeStates.getNumberOfSetBits());
+                uint_fast64_t p1MatrixRow = 0;
+                for (auto maybeState : maybeStates){
+                    matrixBuilder.newRowGroup(p1MatrixRow);
+                    for (uint_fast64_t row = this->parametricModel.getTransitionMatrix().getRowGroupIndices()[maybeState]; row < this->parametricModel.getTransitionMatrix().getRowGroupIndices()[maybeState + 1]; ++row){
+                        matrixBuilder.addNextValue(p1MatrixRow, p1MatrixRow, storm::utility::one<storm::storage::sparse::state_type>());
+                        ++p1MatrixRow;
+                    }
+                }
+                player1Matrix = matrixBuilder.build();
+            }
+    
+            template <typename SparseModelType, typename ConstantType>
+            void SparseMdpParameterLiftingModelChecker<SparseModelType, ConstantType>::reset() {
+                maybeStates.resize(0);
+                resultsForNonMaybeStates.clear();
+                stepBound = boost::none;
+                player1Matrix = storm::storage::SparseMatrix<storm::storage::sparse::state_type>();
+                parameterLifter = nullptr;
+                minSched = boost::none;
+                maxSched = boost::none;
+                x.clear();
+                lowerResultBound = boost::none;
+                upperResultBound = boost::none;
+                applyPreviousResultAsHint = false;
+            }
+    
+            template class SparseMdpParameterLiftingModelChecker<storm::models::sparse::Mdp<storm::RationalFunction>, double>;
+        }
+    }
+}

src/storm/modelchecker/parametric/SparseMdpParameterLiftingModelChecker.h

@@ -0,0 +1,57 @@
+#pragma once
+
+#include <vector>
+#include <memory>
+#include <boost/optional.hpp>
+
+#include "storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h"
+#include "storm/storage/BitVector.h"
+#include "storm/storage/SparseMatrix.h"
+#include "storm/storage/sparse/StateType.h"
+#include "storm/utility/solver.h"
+#include "storm/transformer/ParameterLifter.h"
+#include "storm/storage/TotalScheduler.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace parametric {
+            
+            template <typename SparseModelType, typename ConstantType>
+            class SparseMdpParameterLiftingModelChecker : public SparseParameterLiftingModelChecker<SparseModelType, ConstantType> {
+            public:
+                SparseMdpParameterLiftingModelChecker(SparseModelType const& parametricModel);
+                SparseMdpParameterLiftingModelChecker(SparseModelType const& parametricModel, std::unique_ptr<storm::utility::solver::GameSolverFactory<ConstantType>>&& solverFactory);
+                
+            protected:
+                
+                virtual void specifyBoundedUntilFormula(CheckTask<storm::logic::BoundedUntilFormula, ConstantType> const& checkTask) override;
+                virtual void specifyUntilFormula(CheckTask<storm::logic::UntilFormula, ConstantType> const& checkTask) override;
+                virtual void specifyReachabilityRewardFormula(CheckTask<storm::logic::EventuallyFormula, ConstantType> const& checkTask) override;
+                virtual void specifyCumulativeRewardFormula(CheckTask<storm::logic::CumulativeRewardFormula, ConstantType> const& checkTask) override;
+                
+                virtual std::unique_ptr<CheckResult> computeQuantitativeValues(ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dirForParameters) override;
+                
+                virtual void reset() override;
+                
+
+            private:
+                void computePlayer1Matrix();
+                
+                storm::storage::BitVector maybeStates;
+                std::vector<ConstantType> resultsForNonMaybeStates;
+                boost::optional<uint_fast64_t> stepBound;
+                
+                storm::storage::SparseMatrix<storm::storage::sparse::state_type> player1Matrix;
+                std::unique_ptr<storm::transformer::ParameterLifter<typename SparseModelType::ValueType, ConstantType>> parameterLifter;
+                std::unique_ptr<storm::utility::solver::GameSolverFactory<ConstantType>> solverFactory;
+                
+                // Results from the most recent solver call.
+                boost::optional<storm::storage::TotalScheduler> minSched, maxSched;
+                boost::optional<storm::storage::TotalScheduler> player1Sched;
+                std::vector<ConstantType> x;
+                boost::optional<ConstantType> lowerResultBound, upperResultBound;
+                bool applyPreviousResultAsHint;
+            };
+        }
+    }
+}

src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.cpp

@@ -23,6 +23,7 @@ namespace storm {
             template <typename SparseModelType, typename ConstantType>
             void SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyFormula(storm::modelchecker::CheckTask<storm::logic::Formula, typename SparseModelType::ValueType> const& checkTask) {
                 
+                reset();
                 currentFormula = checkTask.getFormula().asSharedPointer();
                 currentCheckTask = std::make_unique<storm::modelchecker::CheckTask<storm::logic::Formula, ConstantType>>(checkTask.substituteFormula(*currentFormula).template convertValueType<ConstantType>());
                 

src/storm/modelchecker/parametric/SparseParameterLiftingModelChecker.h

@@ -42,6 +42,7 @@ namespace storm {
                 
                 virtual std::unique_ptr<CheckResult> computeQuantitativeValues(ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dirForParameters) = 0;
                 
+                virtual void reset() = 0;
                 
                 SparseModelType const& parametricModel;
                 std::unique_ptr<CheckTask<storm::logic::Formula, ConstantType>> currentCheckTask;

src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@@ -252,7 +252,7 @@ namespace storm {
                 // Initialize result to either the state rewards of the model or the null vector.
                 std::vector<ValueType> result;
                 if (rewardModel.hasStateRewards()) {
-                    result = rewardModel.getStateRewardVector();
+                    result = rewardModel.getStateRewardVector(); // Whyyy?
                 } else {
                     result.resize(transitionMatrix.getRowGroupCount());
                 }

src/storm/modelchecker/region/SparseMdpRegionModelChecker.cpp

@@ -33,6 +33,7 @@
 
 #include "storm/transformer/SparseParametricMdpSimplifier.h"
 #include "storm/modelchecker/parametric/SparseMdpInstantiationModelChecker.h"
+#include "storm/modelchecker/parametric/SparseMdpParameterLiftingModelChecker.h"
 
 
 namespace storm {
@@ -320,7 +321,22 @@ namespace storm {
                         
             template<typename ParametricSparseModelType, typename ConstantType>
             void SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::OperatorFormula const> formula) {
-                std::cout << "approx model for mdps not implemented" << std::endl;
+                STORM_LOG_DEBUG("Initializing the Approx Model....");
+                std::chrono::high_resolution_clock::time_point timeInitApproximationModelStart = std::chrono::high_resolution_clock::now();
+                this->approximationModel = std::make_shared<storm::modelchecker::parametric::SparseMdpParameterLiftingModelChecker<ParametricSparseModelType, ConstantType>>(model);
+                // replace bounds by optimization direction to obtain a quantitative formula
+                if(formula->isProbabilityOperatorFormula()) {
+                    auto quantitativeFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula const>(formula->getSubformula().asSharedPointer(), storm::logic::OperatorInformation(storm::logic::isLowerBound(formula->getComparisonType()) ? storm::solver::OptimizationDirection::Minimize : storm::solver::OptimizationDirection::Maximize));
+                    this->approximationModel->specifyFormula(*quantitativeFormula);
+                } else {
+                    auto quantitativeFormula = std::make_shared<storm::logic::RewardOperatorFormula>(formula->getSubformula().asSharedPointer(), formula->asRewardOperatorFormula().getOptionalRewardModelName(),  storm::logic::OperatorInformation(storm::logic::isLowerBound(formula->getComparisonType()) ? storm::solver::OptimizationDirection::Minimize : storm::solver::OptimizationDirection::Maximize));
+                    this->approximationModel->specifyFormula(*quantitativeFormula);
+                }
+                std::chrono::high_resolution_clock::time_point timeInitApproximationModelEnd = std::chrono::high_resolution_clock::now();
+                STORM_LOG_DEBUG("Initialized Approximation Model");
+                
+                
+                
             }
             
 #ifdef STORM_HAVE_CARL

src/storm/transformer/ParameterLifter.cpp

@@ -5,6 +5,8 @@
 
 #include "storm/utility/vector.h"
 
+#include "storm/exceptions/UnexpectedException.h"
+#include "storm/exceptions/NotSupportedException.h"
 
 namespace storm {
     namespace transformer {
@@ -203,6 +205,23 @@ namespace storm {
             return seed;
         }
     
+        template<typename ParametricType, typename ConstantType>
+        typename ParameterLifter<ParametricType, ConstantType>::AbstractValuation ParameterLifter<ParametricType, ConstantType>::AbstractValuation::getSubValuation(std::set<VariableType> const& pars) const {
+            AbstractValuation result;
+            for (auto const& p : pars) {
+                if (std::find(lowerPars.begin(), lowerPars.end(), p) != lowerPars.end()) {
+                    result.addParameterLower(p);
+                } else if (std::find(upperPars.begin(), upperPars.end(), p) != upperPars.end()) {
+                    result.addParameterUpper(p);
+                } else if (std::find(unspecifiedPars.begin(), unspecifiedPars.end(), p) != unspecifiedPars.end()) {
+                    result.addParameterUnspecified(p);
+                } else {
+                    STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Tried to obtain a subvaluation for parameters that are not specified by this valuation");
+                }
+            }
+            return result;
+        }
+            
         template<typename ParametricType, typename ConstantType>
         std::vector<storm::utility::parametric::Valuation<ParametricType>> ParameterLifter<ParametricType, ConstantType>::AbstractValuation::getConcreteValuations(storm::modelchecker::parametric::ParameterRegion<ParametricType> const& region) const {
             auto result = region.getVerticesOfRegion(unspecifiedPars);
@@ -219,11 +238,18 @@ namespace storm {
         
         template<typename ParametricType, typename ConstantType>
         ConstantType& ParameterLifter<ParametricType, ConstantType>::FunctionValuationCollector::add(ParametricType const& function, AbstractValuation const& valuation) {
-            
+            ParametricType simplifiedFunction = function;
+            storm::utility::simplify(simplifiedFunction);
+            std::set<VariableType> variablesInFunction;
+            storm::utility::parametric::gatherOccurringVariables(simplifiedFunction, variablesInFunction);
+            AbstractValuation simplifiedValuation = valuation.getSubValuation(variablesInFunction);
             // insert the function and the valuation
-            auto functionsIt = collectedFunctions.insert(std::pair<FunctionValuation, ConstantType>(FunctionValuation(function, valuation), storm::utility::one<ConstantType>())).first;
+            auto insertionRes = collectedFunctions.insert(std::pair<FunctionValuation, ConstantType>(FunctionValuation(std::move(simplifiedFunction), std::move(simplifiedValuation)), storm::utility::one<ConstantType>()));
+            if(insertionRes.second) {
+                STORM_LOG_THROW(storm::utility::parametric::isMultiLinearPolynomial(insertionRes.first->first.first), storm::exceptions::NotSupportedException, "Parameter lifting for non-multilinear polynomial " << insertionRes.first->first.first << " is not supported");
+            }
             //Note that references to elements of an unordered map remain valid after calling unordered_map::insert.
-            return functionsIt->second;
+            return insertionRes.first->second;
         }
     
         template<typename ParametricType, typename ConstantType>
@@ -247,7 +273,6 @@ namespace storm {
             }
         }
         
-        
         template class ParameterLifter<storm::RationalFunction, double>;
     }
 }

src/storm/transformer/ParameterLifter.h

@@ -71,6 +71,7 @@ namespace storm {
                 void addParameterUnspecified(VariableType const& var);
                 
                 std::size_t getHashValue() const;
+                AbstractValuation getSubValuation(std::set<VariableType> const& pars) const;
                 
                 /*!
                  * Returns the concrete valuation(s) (w.r.t. the provided region) represented by this abstract valuation.

src/storm/utility/parametric.cpp

@@ -31,6 +31,25 @@ namespace storm {
             void gatherOccurringVariables<storm::RationalFunction>(storm::RationalFunction const& function, std::set<typename VariableType<storm::RationalFunction>::type>& variableSet){
                 function.gatherVariables(variableSet);
             }
+            
+            template<>
+            bool isLinear<storm::RationalFunction>(storm::RationalFunction const& function) {
+                return storm::utility::isConstant(function.denominator()) && function.nominator().isLinear();
+            }
+            
+            template<>
+            bool isMultiLinearPolynomial<storm::RationalFunction>(storm::RationalFunction const& function) {
+                if (!storm::utility::isConstant(function.denominator())) {
+                    return false;
+                }
+                auto varInfos = function.nominator().getVarInfo<false>();
+                for (auto const& varInfo : varInfos) {
+                    if (varInfo.second.maxDegree() > 1) {
+                        return false;
+                    }
+                }
+                return true;
+            }
 #endif
         }
     }

src/storm/utility/parametric.h

@@ -48,6 +48,18 @@ namespace storm {
             template<typename FunctionType>
             void gatherOccurringVariables(FunctionType const& function, std::set<typename VariableType<FunctionType>::type>& variableSet);
             
+            /*!
+             *  Checks whether the function is linear (in one parameter)
+             */
+            template<typename FunctionType>
+            bool isLinear(FunctionType const& function);
+            
+            /*!
+             *  Checks whether the function is a multilinear polynomial, i.e., a polynomial which only considers variables with exponent at most 1
+             */
+            template<typename FunctionType>
+            bool isMultiLinearPolynomial(FunctionType const& function);
+            
         }
         
     }