diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp
index 5ef19802d..4f1ddeb3c 100644
--- a/src/modelchecker/region/ApproximationModel.cpp
+++ b/src/modelchecker/region/ApproximationModel.cpp
@@ -6,10 +6,11 @@
  */
 
 #include "src/modelchecker/region/ApproximationModel.h"
-#include "src/modelchecker/region/ParameterRegion.h"
-#include "modelchecker/prctl/SparseMdpPrctlModelChecker.h"
+#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
-#include "exceptions/UnexpectedException.h"
+#include "src/utility/vector.h"
+#include "src/utility/regions.h"
+#include "src/exceptions/UnexpectedException.h"
 
 namespace storm {
     namespace modelchecker {
@@ -17,31 +18,98 @@ namespace storm {
         
         template<typename ParametricType, typename ConstantType>
         SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::ApproximationModel(storm::models::sparse::Dtmc<ParametricType> const& parametricModel, bool computeRewards) : computeRewards(computeRewards){
-            // Run through the rows of the original model and obtain
-            // (1) the different substitutions (this->substitutions) and the substitution used for every row
-            std::vector<std::size_t> rowSubstitutions;
-            this->substitutions.emplace_back(std::map<VariableType, TypeOfBound>()); //we want that the empty substitution is always the first one
-            // (2) the set of distinct pairs <Function, Substitution> 
-            std::set<std::pair<ParametricType, std::size_t>> distinctFuncSubs;
-            // (3) the MDP Probability matrix with some dummy entries
-            storm::storage::SparseMatrixBuilder<ConstantType> probabilityMatrixBuilder(0, parametricModel.getNumberOfStates(), 0, true, true, parametricModel.getNumberOfStates());
-            typename storm::storage::SparseMatrix<ConstantType>::index_type probMatrixRow=0;
-            std::size_t numOfNonConstProbEntries=0;
+            //Start with the probabilities
+            storm::storage::SparseMatrix<ConstantType> probabilityMatrix;
+            std::vector<std::size_t> rowSubstitutions;// the substitution used in every row (required if rewards are computed)
+            std::vector<std::size_t> matrixEntryToEvalTableMapping;// This vector will get an entry for every probability-matrix entry
+            //for the corresponding matrix entry, it stores the corresponding entry in the probability evaluation table (more precisely: the position in that table).
+            //We can later transform this mapping into the desired mapping with iterators
+            const std::size_t constantEntryIndex=std::numeric_limits<std::size_t>::max(); //this value is stored in the matrixEntrytoEvalTableMapping for every constant matrix entry. (also used for rewards later)
+            initializeProbabilities(parametricModel, probabilityMatrix, rowSubstitutions, matrixEntryToEvalTableMapping, constantEntryIndex);
+            
+            //Now consider rewards
+            boost::optional<std::vector<ConstantType>> stateRewards;
+            boost::optional<storm::storage::SparseMatrix<ConstantType>> transitionRewards;
+            std::vector<std::size_t> stateRewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping
+            std::vector<std::size_t> transitionRewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping
+            if(this->computeRewards){
+                initializeRewards(parametricModel, probabilityMatrix, rowSubstitutions, stateRewards, transitionRewards, stateRewardEntryToEvalTableMapping, transitionRewardEntryToEvalTableMapping, constantEntryIndex);
+            }
+            //Obtain other model ingredients and the approximation model itself
+            storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
+            boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;  
+            this->model=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(stateRewards), std::move(transitionRewards), std::move(noChoiceLabeling));
+            
+            //translate the matrixEntryToEvalTableMapping into the actual probability mapping
+            typename storm::storage::SparseMatrix<ConstantType>::index_type matrixRow=0;
+            auto tableIndex=matrixEntryToEvalTableMapping.begin();
             for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
-                probabilityMatrixBuilder.newRowGroup(probMatrixRow);
-                // For (1): Find the different substitutions, i.e., mappings from Variables that occur in this row to {lower, upper}
-                std::set<VariableType> occurringProbVariables;
+                for (; matrixRow<this->model->getTransitionMatrix().getRowGroupIndices()[row+1]; ++matrixRow){
+                    auto approxModelEntry = this->model->getTransitionMatrix().getRow(matrixRow).begin();
+                    for(auto const& parEntry : parametricModel.getTransitionMatrix().getRow(row)){
+                        if(*tableIndex == constantEntryIndex){
+                            approxModelEntry->setValue(storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parEntry.getValue())));
+                        } else {
+                            this->probabilityMapping.emplace_back(std::make_pair(&(std::get<2>(this->probabilityEvaluationTable[*tableIndex])), approxModelEntry));
+                        }
+                        ++approxModelEntry;
+                        ++tableIndex;
+                    }
+                }
+            }
+            if(this->computeRewards){
+                //the same for state and transition rewards
+                auto approxModelStateRewardEntry = this->model->getStateRewardVector().begin();
+                for (std::size_t const& tableIndex : stateRewardEntryToEvalTableMapping){
+                    if(tableIndex != constantEntryIndex){
+                        this->stateRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(this->rewardEvaluationTable[tableIndex])), &(std::get<3>(this->rewardEvaluationTable[tableIndex])), approxModelStateRewardEntry));
+                    }
+                    ++approxModelStateRewardEntry;
+                }
+                auto approxModelTransitionRewardEntry = this->model->getTransitionRewardMatrix().begin();
+                for (std::size_t const& tableIndex : transitionRewardEntryToEvalTableMapping){
+                    this->transitionRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(this->rewardEvaluationTable[tableIndex])), &(std::get<3>(this->rewardEvaluationTable[tableIndex])), approxModelTransitionRewardEntry));
+                    ++approxModelTransitionRewardEntry;
+                }
+                //Get the sets of reward parameters that map to "CHOSEOPTIMAL".
+                this->choseOptimalRewardPars = std::vector<std::set<VariableType>>(this->rewardSubstitutions.size());
+                for(std::size_t index = 0; index<this->rewardSubstitutions.size(); ++index){
+                    for(auto const& sub : this->rewardSubstitutions[index]){
+                        if(sub.second==TypeOfBound::CHOSEOPTIMAL){
+                            this->choseOptimalRewardPars[index].insert(sub.first);
+                        }
+                    }
+                }
+            }
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::initializeProbabilities(storm::models::sparse::Dtmc<ParametricType>const& parametricModel,
+                                                                                                                     storm::storage::SparseMatrix<ConstantType>& probabilityMatrix,
+                                                                                                                     std::vector<std::size_t>& rowSubstitutions,
+                                                                                                                     std::vector<std::size_t>& matrixEntryToEvalTableMapping,
+                                                                                                                     std::size_t const& constantEntryIndex) {
+            // Run through the rows of the original model and obtain the different substitutions, the probability evaluation table,
+            // an MDP probability matrix with some dummy entries, and the mapping between the two
+            storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(0, parametricModel.getNumberOfStates(), 0, true, true, parametricModel.getNumberOfStates());
+            this->probabilitySubstitutions.emplace_back(std::map<VariableType, TypeOfBound>()); //we want that the empty substitution is always the first one
+            std::size_t numOfNonConstEntries=0;
+            typename storm::storage::SparseMatrix<ConstantType>::index_type matrixRow=0;
+            for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
+                matrixBuilder.newRowGroup(matrixRow);
+                // Find the different substitutions, i.e., mappings from Variables that occur in this row to {lower, upper}
+                std::set<VariableType> occurringVariables;
                 for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
-                    storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringProbVariables);
+                    storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringVariables);
                 }
-                std::size_t numOfSubstitutions=1ull<<occurringProbVariables.size(); //=2^(#variables). Note that there is still 1 substitution when #variables==0 (i.e.,the empty substitution)
+                std::size_t numOfSubstitutions=1ull<<occurringVariables.size(); //=2^(#variables). Note that there is still 1 substitution when #variables==0 (i.e.,the empty substitution)
                 for(uint_fast64_t substitutionId=0; substitutionId<numOfSubstitutions; ++substitutionId){
                     //compute actual substitution from substitutionId by interpreting the Id as a bit sequence.
-                    //the occurringProbVariables.size() least significant bits of substitutionId will always represent the substitution that we have to consider
+                    //the occurringVariables.size() least significant bits of substitutionId will always represent the substitution that we have to consider
                     //(00...0 = lower bounds for all parameters, 11...1 = upper bounds for all parameters)
                     std::map<VariableType, TypeOfBound> currSubstitution;
                     std::size_t parameterIndex=0;
-                    for(auto const& parameter : occurringProbVariables){
+                    for(auto const& parameter : occurringVariables){
                         if((substitutionId>>parameterIndex)%2==0){
                             currSubstitution.insert(std::make_pair(parameter, TypeOfBound::LOWER));
                         }
@@ -50,184 +118,119 @@ namespace storm {
                         }
                         ++parameterIndex;
                     }
-                    //Find the current substitution in this->substitutions (add it if we see this substitution the first time)
-                    std::size_t substitutionIndex = std::find(this->substitutions.begin(),this->substitutions.end(), currSubstitution) - this->substitutions.begin();
-                    if(substitutionIndex==this->substitutions.size()){
-                        this->substitutions.emplace_back(std::move(currSubstitution));
-                    }
+                    std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->probabilitySubstitutions, std::move(currSubstitution));
                     rowSubstitutions.push_back(substitutionIndex);
-                    //Run again through the row and...
-                    //For (2): add pair <function, substitution> for the occuring functions and the current substitution
-                    //For (3): add a row with dummy entries for the current substitution
-                    //Note that this is still executed even if no variables occur. However, we will not put any <function, substitution> pair into distninctFuncSubs if substitution is empty.
+                    //For every substitution, run again through the row and add an entry in matrixEntryToEvalTableMapping as well as dummy entries in the matrix
+                    //Note that this is still executed once, even if no parameters occur.
                     ConstantType dummyEntry=storm::utility::one<ConstantType>();
                     for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
-                        if(!this->parametricTypeComparator.isConstant(entry.getValue())){
-                            distinctFuncSubs.insert(std::make_pair(entry.getValue(), substitutionIndex));
-                            ++numOfNonConstProbEntries;
+                        if(this->parametricTypeComparator.isConstant(entry.getValue())){
+                            matrixEntryToEvalTableMapping.emplace_back(constantEntryIndex);
+                        } else {
+                            ++numOfNonConstEntries;
+                            std::size_t tableIndex=storm::utility::vector::findOrInsert(this->probabilityEvaluationTable, std::tuple<ParametricType, std::size_t, ConstantType>(entry.getValue(), substitutionIndex, storm::utility::zero<ConstantType>()));
+                            matrixEntryToEvalTableMapping.emplace_back(tableIndex);
                         }
-                        probabilityMatrixBuilder.addNextValue(probMatrixRow, entry.getColumn(), dummyEntry);
+                        matrixBuilder.addNextValue(matrixRow, entry.getColumn(), dummyEntry);
                         dummyEntry=storm::utility::zero<ConstantType>(); 
                     }
-                    ++probMatrixRow;
+                    ++matrixRow;
                 }
             }
-            storm::storage::SparseMatrix<ConstantType> probabilityMatrix(probabilityMatrixBuilder.build());
-            this->probabilityMapping.reserve(numOfNonConstProbEntries);
+            this->probabilityMapping.reserve(numOfNonConstEntries);
+            probabilityMatrix=matrixBuilder.build();
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::initializeRewards(storm::models::sparse::Dtmc<ParametricType> const& parametricModel,
+                                                                                                               storm::storage::SparseMatrix<ConstantType> const& probabilityMatrix,
+                                                                                                               std::vector<std::size_t> const& rowSubstitutions,
+                                                                                                               boost::optional<std::vector<ConstantType> >& stateRewards,
+                                                                                                               boost::optional<storm::storage::SparseMatrix<ConstantType> >& transitionRewards,
+                                                                                                               std::vector<std::size_t>& stateRewardEntryToEvalTableMapping,
+                                                                                                               std::vector<std::size_t>& transitionRewardEntryToEvalTableMapping,
+                                                                                                               std::size_t const& constantEntryIndex) {
+            // run through the state reward vector of the parametric model.
+            // Constant entries can be set directly.
+            // For Parametric entries we have two cases:
+            // (1) make state rewards if the reward is independent of parameters that occur in probability functions.
+            // (2) make transition rewards otherwise.
             
-            // Now obtain transition and staterewards (if required)
-            boost::optional<std::vector<ConstantType>> stateRewards;
-            boost::optional<storm::storage::SparseMatrix<ConstantType>> transitionRewards;
-            std::vector<std::size_t> stateRewardSubstituions;
-            std::vector<std::size_t> transitionRewardRowSubstitutions;
-            std::set<std::pair<ParametricType, std::size_t>> distinctRewardFuncSubs;
-            if(this->computeRewards){
-                this->rewardSubstitutions.emplace_back(std::map<VariableType, TypeOfBound>()); //we want that the empty substitution is always the first one
-                //stateRewards
-                std::vector<ConstantType> stateRewardsAsVector(parametricModel.getNumberOfStates());
-                stateRewardSubstituions = std::vector<std::size_t>(parametricModel.getNumberOfStates(),0); //init with empty substitution (=0)
-                std::size_t numOfNonConstStateRewEntries=0;
-                //TransitionRewards
-                storm::storage::SparseMatrixBuilder<ConstantType> rewardMatrixBuilder(probabilityMatrix.getRowCount(), probabilityMatrix.getColumnCount(), 0, true, true, probabilityMatrix.getRowGroupCount());
-                transitionRewardRowSubstitutions = std::vector<std::size_t>(rowSubstitutions.size(),0); //init with empty substitution (=0)
-                std::size_t numOfNonConstTransitonRewEntries=0;
-                for(std::size_t state=0; state<parametricModel.getNumberOfStates(); ++state){
-                    rewardMatrixBuilder.newRowGroup(probabilityMatrix.getRowGroupIndices()[state]);
-                    std::set<VariableType> occurringRewVariables;
-                    std::set<VariableType> occurringProbVariables;
-                    bool makeStateReward=true; //we make state reward whenever no probability parameter occurs in the state reward function
-                    if(this->parametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){
-                        //constant reward for this state
-                        stateRewardsAsVector[state]=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state]));
-                    } else {
-                        //reward depends on parameters. Lets find out if also probability parameters occur here.
-                        //If this is the case, the reward depends on the nondeterministic choices and should be given as transition rewards.
-                        stateRewardsAsVector[state]=storm::utility::zero<ConstantType>();
-                        storm::utility::regions::gatherOccurringVariables(parametricModel.getStateRewardVector()[state], occurringRewVariables);
-                        for(auto const& entry : parametricModel.getTransitionMatrix().getRow(state)){
-                            storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringProbVariables);
+            //stateRewards
+            std::vector<ConstantType> stateRewardsAsVector(parametricModel.getNumberOfStates());
+            std::size_t numOfNonConstStateRewEntries=0;
+            //TransitionRewards
+            storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(probabilityMatrix.getRowCount(), probabilityMatrix.getColumnCount(), 0, true, true, probabilityMatrix.getRowGroupCount());
+            std::size_t numOfNonConstTransitonRewEntries=0;
+            this->rewardSubstitutions.emplace_back(std::map<VariableType, TypeOfBound>()); //we want that the empty substitution is always the first one
+            
+            for(std::size_t state=0; state<parametricModel.getNumberOfStates(); ++state){
+                matrixBuilder.newRowGroup(probabilityMatrix.getRowGroupIndices()[state]);
+                std::set<VariableType> occurringRewVariables;
+                std::set<VariableType> occurringProbVariables;
+                bool makeStateReward=true;
+                if(this->parametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){
+                    stateRewardsAsVector[state]=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state]));
+                    stateRewardEntryToEvalTableMapping.emplace_back(constantEntryIndex);
+                } else {
+                    //reward depends on parameters. Lets find out if probability parameters occur here.
+                    //If this is the case, the reward depends on the nondeterministic choices and should be given as transition rewards.
+                    storm::utility::regions::gatherOccurringVariables(parametricModel.getStateRewardVector()[state], occurringRewVariables);
+                    for(auto const& entry : parametricModel.getTransitionMatrix().getRow(state)){
+                        storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringProbVariables);
+                    }
+                    for( auto const& rewVar : occurringRewVariables){
+                        if(occurringProbVariables.find(rewVar)!=occurringProbVariables.end()){
+                            makeStateReward=false;
+                            break;
                         }
-                        for( auto const& rewVar : occurringRewVariables){
-                            if(occurringProbVariables.find(rewVar)!=occurringProbVariables.end()){
-                                makeStateReward=false;
-                                break;
-                            }
+                    }
+                    if(makeStateReward){
+                        //Get the corresponding substitution and substitutionIndex
+                        std::map<VariableType, TypeOfBound> rewardSubstitution;
+                        for(auto const& rewardVar : occurringRewVariables){
+                            rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL));
                         }
-                        if(makeStateReward){
+                        std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->rewardSubstitutions, std::move(rewardSubstitution));
+                        //insert table entry and dummy data in the stateRewardVector
+                        std::size_t tableIndex=storm::utility::vector::findOrInsert(this->rewardEvaluationTable, std::tuple<ParametricType, std::size_t, ConstantType, ConstantType>(parametricModel.getStateRewardVector()[state], substitutionIndex, storm::utility::zero<ConstantType>(), storm::utility::zero<ConstantType>()));
+                        stateRewardEntryToEvalTableMapping.emplace_back(tableIndex);
+                        stateRewardsAsVector[state]=storm::utility::zero<ConstantType>();
+                        ++numOfNonConstStateRewEntries;
+                    } else {
+                        for(auto matrixRow=probabilityMatrix.getRowGroupIndices()[state]; matrixRow<probabilityMatrix.getRowGroupIndices()[state+1]; ++matrixRow){
+                            //Get the corresponding substitution and substitutionIndex
                             std::map<VariableType, TypeOfBound> rewardSubstitution;
                             for(auto const& rewardVar : occurringRewVariables){
-                                rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL));
-                            }
-                            //Find the substitution in this->rewardSubstitutions (add it if we see this substitution the first time)
-                            std::size_t substitutionIndex = std::find(this->rewardSubstitutions.begin(),this->rewardSubstitutions.end(), rewardSubstitution) - this->rewardSubstitutions.begin();
-                            if(substitutionIndex==this->rewardSubstitutions.size()){
-                                this->rewardSubstitutions.emplace_back(std::move(rewardSubstitution));
-                            }
-                            distinctRewardFuncSubs.insert(std::make_pair(parametricModel.getStateRewardVector()[state], substitutionIndex));
-                            ++numOfNonConstStateRewEntries;
-                            stateRewardSubstituions[state]=substitutionIndex;
-                        } else {
-                            for(auto matrixRow=probabilityMatrix.getRowGroupIndices()[state]; matrixRow<probabilityMatrix.getRowGroupIndices()[state+1]; ++matrixRow){
-                                //Reserve space for transition rewards
-                                for(auto const& matrixEntry : probabilityMatrix.getRow(matrixRow)){
-                                    rewardMatrixBuilder.addNextValue(matrixRow,matrixEntry.getColumn(),storm::utility::zero<ConstantType>());
-                                }
-                                //Get the corresponding substitution
-                                std::map<VariableType, TypeOfBound> rewardSubstitution;
-                                for(auto const& rewardVar : occurringRewVariables){
-                                    typename std::map<VariableType, TypeOfBound>::iterator const substitutionIt=this->substitutions[rowSubstitutions[matrixRow]].find(rewardVar);
-                                    if(substitutionIt==this->substitutions[rowSubstitutions[matrixRow]].end()){
-                                        rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL));
-                                    } else {
-                                        rewardSubstitution.insert(*substitutionIt);
-                                    }
-                                }
-                                //Find the substitution in this->rewardSubstitutions (add it if we see this substitution the first time)
-                                std::size_t substitutionIndex = std::find(this->rewardSubstitutions.begin(),this->rewardSubstitutions.end(), rewardSubstitution) - this->rewardSubstitutions.begin();
-                                if(substitutionIndex==this->rewardSubstitutions.size()){
-                                    this->rewardSubstitutions.emplace_back(std::move(rewardSubstitution));
+                                typename std::map<VariableType, TypeOfBound>::iterator const substitutionIt=this->probabilitySubstitutions[rowSubstitutions[matrixRow]].find(rewardVar);
+                                if(substitutionIt==this->probabilitySubstitutions[rowSubstitutions[matrixRow]].end()){
+                                    rewardSubstitution.insert(std::make_pair(rewardVar, TypeOfBound::CHOSEOPTIMAL));
+                                } else {
+                                    rewardSubstitution.insert(*substitutionIt);
                                 }
-                                distinctRewardFuncSubs.insert(std::make_pair(parametricModel.getStateRewardVector()[state], substitutionIndex));
-                                ++numOfNonConstTransitonRewEntries;
-                                transitionRewardRowSubstitutions[matrixRow]=substitutionIndex;
                             }
-                        }
-                    }
-                }
-                stateRewards=std::move(stateRewardsAsVector);
-                this->stateRewardMapping.reserve(numOfNonConstStateRewEntries);
-                transitionRewards=rewardMatrixBuilder.build();
-                this->transitionRewardMapping.reserve(numOfNonConstTransitonRewEntries);
-            }
-            //Obtain other model ingredients and the approximation model itself
-            storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
-            boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;  
-            this->model=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(stateRewards), std::move(transitionRewards), std::move(noChoiceLabeling));
-            
-            //Get the (probability) evaluation table. Note that it remains sorted due to the fact that distinctFuncSubs is sorted
-            this->evaluationTable.reserve(distinctFuncSubs.size());
-            for(auto const& funcSub : distinctFuncSubs){
-                this->evaluationTable.emplace_back(funcSub.first, funcSub.second, storm::utility::zero<ConstantType>());
-            }
-            //Fill in the entries for the probability mapping
-            probMatrixRow=0;
-            for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
-                for (; probMatrixRow<this->model->getTransitionMatrix().getRowGroupIndices()[row+1]; ++probMatrixRow){
-                    auto appEntry = this->model->getTransitionMatrix().getRow(probMatrixRow).begin();
-                    for(auto const& parEntry : parametricModel.getTransitionMatrix().getRow(row)){
-                        if(this->parametricTypeComparator.isConstant(parEntry.getValue())){
-                            appEntry->setValue(storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parEntry.getValue())));
-                        } else {
-                            std::tuple<ParametricType, std::size_t, ConstantType> searchedTuple(parEntry.getValue(), rowSubstitutions[probMatrixRow], storm::utility::zero<ConstantType>());
-                            auto const tableIt= std::find(this->evaluationTable.begin(), this->evaluationTable.end(), searchedTuple);
-                            STORM_LOG_THROW((*tableIt==searchedTuple),storm::exceptions::UnexpectedException, "Could not find the current tuple in the evaluationTable. Either the table is missing that tuple or it is not sorted properly");
-                            this->probabilityMapping.emplace_back(std::make_pair(&(std::get<2>(*tableIt)), appEntry));
-                        }
-                        ++appEntry;
-                    }
-                }
-            }
-            if(this->computeRewards){
-                //Get the (reward) evaluation table. Note that it remains sorted due to the fact that distinctRewFuncSubs is sorted
-                this->rewardEvaluationTable.reserve(distinctRewardFuncSubs.size());
-                for(auto const& funcSub : distinctRewardFuncSubs){
-                    this->rewardEvaluationTable.emplace_back(funcSub.first, funcSub.second, storm::utility::zero<ConstantType>(), storm::utility::zero<ConstantType>());
-                }
-                
-                for (std::size_t state=0; state< parametricModel.getNumberOfStates(); ++state){
-                    //check if the state reward is not constant
-                    if(stateRewardSubstituions[state]!=0){
-                        std::tuple<ParametricType, std::size_t, ConstantType, ConstantType> searchedTuple(parametricModel.getStateRewardVector()[state], stateRewardSubstituions[state], storm::utility::zero<ConstantType>(), storm::utility::zero<ConstantType>());
-                        auto const tableIt= std::find(this->rewardEvaluationTable.begin(), this->rewardEvaluationTable.end(), searchedTuple);
-                        STORM_LOG_THROW((*tableIt==searchedTuple),storm::exceptions::UnexpectedException, "Could not find the current tuple in the rewardEvaluationTable. Either the table is missing that tuple or it is not sorted properly (stateRewards)");
-                        this->stateRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(*tableIt)), &(std::get<3>(*tableIt)), &(this->model->getStateRewardVector()[state])));
-                    }
-                    //check if transitionrewards are not constant
-                    for(auto matrixRow=this->model->getTransitionRewardMatrix().getRowGroupIndices()[state]; matrixRow<this->model->getTransitionRewardMatrix().getRowGroupIndices()[state+1]; ++matrixRow){
-                        if(transitionRewardRowSubstitutions[matrixRow]!=0){
-                            //Note that all transitions of the same row will have the same reward value
-                            std::tuple<ParametricType, std::size_t, ConstantType, ConstantType> searchedTuple(parametricModel.getStateRewardVector()[state], transitionRewardRowSubstitutions[matrixRow], storm::utility::zero<ConstantType>(), storm::utility::zero<ConstantType>());
-                            auto const tableIt= std::find(this->rewardEvaluationTable.begin(), this->rewardEvaluationTable.end(), searchedTuple);
-                            STORM_LOG_THROW((*tableIt==searchedTuple),storm::exceptions::UnexpectedException, "Could not find the current tuple in the rewardEvaluationTable. Either the table is missing that tuple or it is not sorted properly (transitionRewards)");
-                            for(auto transitionMatrixEntry=this->model->getTransitionRewardMatrix().getRow(matrixRow).begin(); transitionMatrixEntry<this->model->getTransitionRewardMatrix().getRow(matrixRow).end(); ++transitionMatrixEntry){
-                                this->transitionRewardMapping.emplace_back(std::make_tuple(&(std::get<2>(*tableIt)), &(std::get<3>(*tableIt)), transitionMatrixEntry));
+                            std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->rewardSubstitutions, std::move(rewardSubstitution));
+                            //insert table entries and dummy data
+                            std::size_t tableIndex=storm::utility::vector::findOrInsert(this->rewardEvaluationTable, std::tuple<ParametricType, std::size_t, ConstantType, ConstantType>(parametricModel.getStateRewardVector()[state], substitutionIndex, storm::utility::zero<ConstantType>(), storm::utility::zero<ConstantType>()));
+                            for(auto const& matrixEntry : probabilityMatrix.getRow(matrixRow)){
+                                transitionRewardEntryToEvalTableMapping.emplace_back(tableIndex);
+                                matrixBuilder.addNextValue(matrixRow, matrixEntry.getColumn(), storm::utility::zero<ConstantType>());
+                                ++numOfNonConstTransitonRewEntries;
                             }
                         }
-                    }
-                }
-                //Get the sets of reward parameters that map to "CHOSEOPTIMAL".
-                this->choseOptimalRewardPars = std::vector<std::set<VariableType>>(this->rewardSubstitutions.size());
-                for(std::size_t index = 0; index<this->rewardSubstitutions.size(); ++index){
-                    for(auto const& sub : this->rewardSubstitutions[index]){
-                        if(sub.second==TypeOfBound::CHOSEOPTIMAL){
-                            this->choseOptimalRewardPars[index].insert(sub.first);
-                        }
+                        stateRewardsAsVector[state]=storm::utility::zero<ConstantType>();
+                        stateRewardEntryToEvalTableMapping.emplace_back(constantEntryIndex);
                     }
                 }
             }
+            stateRewards=std::move(stateRewardsAsVector);
+            this->stateRewardMapping.reserve(numOfNonConstStateRewEntries);
+            transitionRewards=matrixBuilder.build();
+            this->transitionRewardMapping.reserve(numOfNonConstTransitonRewEntries);
         }
 
 
+
         template<typename ParametricType, typename ConstantType>
         SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::~ApproximationModel() {
             //Intentionally left empty
@@ -241,9 +244,9 @@ namespace storm {
         template<typename ParametricType, typename ConstantType>
         void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::ApproximationModel::instantiate(const ParameterRegion& region) {
             //Instantiate the substitutions
-            std::vector<std::map<VariableType, CoefficientType>> instantiatedSubs(this->substitutions.size());
-            for(uint_fast64_t substitutionIndex=0; substitutionIndex<this->substitutions.size(); ++substitutionIndex){
-                for(std::pair<VariableType, TypeOfBound> const& sub : this->substitutions[substitutionIndex]){
+            std::vector<std::map<VariableType, CoefficientType>> instantiatedSubs(this->probabilitySubstitutions.size());
+            for(uint_fast64_t substitutionIndex=0; substitutionIndex<this->probabilitySubstitutions.size(); ++substitutionIndex){
+                for(std::pair<VariableType, TypeOfBound> const& sub : this->probabilitySubstitutions[substitutionIndex]){
                     switch(sub.second){
                         case TypeOfBound::LOWER:
                             instantiatedSubs[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first)));
@@ -257,7 +260,7 @@ namespace storm {
                 }
             }
             //write entries into evaluation table
-            for(auto& tableEntry : this->evaluationTable){
+            for(auto& tableEntry : this->probabilityEvaluationTable){
                 std::get<2>(tableEntry)=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(
                         storm::utility::regions::evaluateFunction(
                                 std::get<0>(tableEntry),
@@ -316,7 +319,6 @@ namespace storm {
                 }
                 //Note: the rewards are written to the model as soon as the optimality type is known (i.e. in computeValues)
             }
-            
         }
 
         template<typename ParametricType, typename ConstantType>
diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h
index 540cc75a8..e7037ddae 100644
--- a/src/modelchecker/region/ApproximationModel.h
+++ b/src/modelchecker/region/ApproximationModel.h
@@ -10,6 +10,8 @@
 
 #include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
 #include "src/models/sparse/Mdp.h"
+#include "src/modelchecker/region/ParameterRegion.h"
+#include "src/storage/SparseMatrix.h"
 
 namespace storm {
     namespace modelchecker {
@@ -54,24 +56,28 @@ namespace storm {
                 CHOSEOPTIMAL
             }; 
            
+            void initializeProbabilities(storm::models::sparse::Dtmc<ParametricType> const& parametricModel, storm::storage::SparseMatrix<ConstantType>& probabilityMatrix, std::vector<std::size_t>& rowSubstitutions, std::vector<std::size_t>& matrixEntryToEvalTableMapping, std::size_t const& constantEntryIndex);
+            void initializeRewards(storm::models::sparse::Dtmc<ParametricType> const& parametricModel, storm::storage::SparseMatrix<ConstantType> const& probabilityMatrix, std::vector<std::size_t> const& rowSubstitutions, boost::optional<std::vector<ConstantType>>& stateRewards, boost::optional<storm::storage::SparseMatrix<ConstantType>>& transitionRewards, std::vector<std::size_t>& stateRewardEntryToEvalTableMapping, std::vector<std::size_t>& transitionRewardEntryToEvalTableMapping, std::size_t const& constantEntryIndex);
+            
             //Vector has one entry for every (non-constant) matrix entry.
             //pair.first points to an entry in the evaluation table,
             //pair.second is an iterator to the corresponding matrix entry
             std::vector<std::pair<ConstantType*, typename storm::storage::SparseMatrix<ConstantType>::iterator>> probabilityMapping;
             //similar for the rewards. But now the first entry points to a minimal and the second one to a maximal value.
             //The third entry points to the state reward vector and the transitionRewardMatrix, respectively.
-            std::vector<std::tuple<ConstantType*, ConstantType*, ConstantType*>> stateRewardMapping;
+            std::vector<std::tuple<ConstantType*, ConstantType*, typename std::vector<ConstantType>::iterator>> stateRewardMapping;
             std::vector<std::tuple<ConstantType*, ConstantType*, typename storm::storage::SparseMatrix<ConstantType>::iterator>> transitionRewardMapping;
             
-            //Vector has one entry for
-            //(every distinct, non-constant function that occurs somewhere in the model) x (the required combinations of lower and upper bounds of the region)
-            //The second entry represents a substitution as an index in the substitutions vector
+            //Vector has one (unique) entry for every occurring pair of a non-constant function and 
+            // a substitution, i.e., a mapping of variables to a TypeOfBound
+            //The second entry represents the substitution as an index in the substitutions vector
             //The third entry should contain the result when evaluating the function in the first entry, regarding the substitution given by the second entry.
-            std::vector<std::tuple<ParametricType, std::size_t, ConstantType>> evaluationTable;
+            std::vector<std::tuple<ParametricType, std::size_t, ConstantType>> probabilityEvaluationTable;
+            //For rewards, we have the third entry for the minimal value and the fourth entry for the maximal value
             std::vector<std::tuple<ParametricType, std::size_t, ConstantType, ConstantType>> rewardEvaluationTable;
             
             //Vector has one entry for every required substitution (=replacement of parameters with lower/upper bounds of region)
-            std::vector<std::map<VariableType, TypeOfBound>> substitutions;
+            std::vector<std::map<VariableType, TypeOfBound>> probabilitySubstitutions;
             //Same for the different substitutions for the reward functions.
             //In addition, we store the parameters for which the correct substitution 
             //depends on the region and whether to minimize/maximize
diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp
index bdacf0693..370070392 100644
--- a/src/modelchecker/region/SamplingModel.cpp
+++ b/src/modelchecker/region/SamplingModel.cpp
@@ -6,85 +6,113 @@
  */
 
 #include "src/modelchecker/region/SamplingModel.h"
-#include "modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
+#include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
-#include "exceptions/UnexpectedException.h"
+#include "src/utility/vector.h"
+#include "src/utility/regions.h"
+#include "src/exceptions/UnexpectedException.h"
 
 namespace storm {
     namespace modelchecker {
         
         
         template<typename ParametricType, typename ConstantType>
-        SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::SamplingModel(storm::models::sparse::Dtmc<ParametricType> const& parametricModel, bool computeRewards) : probabilityMapping(), stateRewardMapping(), evaluationTable(), computeRewards(computeRewards){
-            // Run through the rows of the original model and obtain the set of distinct functions as well as a matrix with dummy entries
-            std::set<ParametricType> functionSet;
-            storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates(), parametricModel.getTransitionMatrix().getEntryCount());
-            std::size_t numOfNonConstProbEntries=0;
-            for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
-                ConstantType dummyEntry=storm::utility::one<ConstantType>();
-                for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
-                    if(!this->parametricTypeComparator.isConstant(entry.getValue())){
-                        functionSet.insert(entry.getValue());
-                        ++numOfNonConstProbEntries;
-                    }
-                    matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry);
-                    dummyEntry=storm::utility::zero<ConstantType>();
-                }
-            }
-            this->probabilityMapping.reserve(numOfNonConstProbEntries);
+        SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::SamplingModel(storm::models::sparse::Dtmc<ParametricType> const& parametricModel, bool computeRewards) : probabilityMapping(), stateRewardMapping(), probabilityEvaluationTable(), computeRewards(computeRewards){
+            //Start with the probabilities
+            storm::storage::SparseMatrix<ConstantType> probabilityMatrix;
+            // This vector will get an entry for every probability matrix entry.
+            // For the corresponding matrix entry, it stores the right entry in the probability evaluation table (more precisely: the position in that table).
+            std::vector<std::size_t> matrixEntryToEvalTableMapping;
+            const std::size_t constantEntryIndex=std::numeric_limits<std::size_t>::max(); //this value is stored in the matrixEntrytoEvalTableMapping for every constant matrix entry. (also used for rewards later)
+            initializeProbabilities(parametricModel, probabilityMatrix, matrixEntryToEvalTableMapping, constantEntryIndex);
             
-            //Obtain other model ingredients and the sampling model itself
-            storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
+            //Now consider rewards
             boost::optional<std::vector<ConstantType>> stateRewards;
+            std::vector<std::size_t> rewardEntryToEvalTableMapping; //does a similar thing as matrixEntryToEvalTableMapping
             if(this->computeRewards){
-                std::size_t numOfNonConstRewEntries=0;
-                std::vector<ConstantType> stateRewardsAsVector(parametricModel.getNumberOfStates());
-                for(std::size_t state=0; state<parametricModel.getNumberOfStates(); ++state){
-                    if(this->parametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){
-                        stateRewardsAsVector[state] = storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state]));
-                    } else {
-                        stateRewardsAsVector[state] = storm::utility::zero<ConstantType>();
-                        functionSet.insert(parametricModel.getStateRewardVector()[state]);
-                        ++numOfNonConstRewEntries;
-                    }
-                }
-                stateRewards=std::move(stateRewardsAsVector);
-                this->stateRewardMapping.reserve(numOfNonConstRewEntries);
+                initializeRewards(parametricModel, stateRewards, rewardEntryToEvalTableMapping, constantEntryIndex);
             }
+            
+            //Obtain other model ingredients and the sampling model itself
+            storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
             boost::optional<storm::storage::SparseMatrix<ConstantType>> noTransitionRewards;
             boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;  
-            this->model=std::make_shared<storm::models::sparse::Dtmc<ConstantType>>(matrixBuilder.build(), std::move(labeling), std::move(stateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
-            
-            //Get the evaluation table. Note that it remains sorted due to the fact that probFunctionSet is sorted
-            this->evaluationTable.reserve(functionSet.size());
-            for(auto const& func : functionSet){
-                this->evaluationTable.emplace_back(func, storm::utility::zero<ConstantType>());
-            }
-            
-            //Fill in the entries for the mappings
-            auto samEntry = this->model->getTransitionMatrix().begin();
-            for(auto const& parEntry : parametricModel.getTransitionMatrix()){
-                if(this->parametricTypeComparator.isConstant(parEntry.getValue())){
-                    samEntry->setValue(storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parEntry.getValue())));
+            this->model=std::make_shared<storm::models::sparse::Dtmc<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(stateRewards), std::move(noTransitionRewards), std::move(noChoiceLabeling));
+                
+            //translate the matrixEntryToEvalTableMapping into the actual probability mapping
+            auto sampleModelEntry = this->model->getTransitionMatrix().begin();
+            auto parModelEntry = parametricModel.getTransitionMatrix().begin();
+            for(std::size_t const& tableIndex : matrixEntryToEvalTableMapping){
+                STORM_LOG_THROW(sampleModelEntry->getColumn()==parModelEntry->getColumn(), storm::exceptions::UnexpectedException, "The entries of the given parametric model and the constructed sampling model do not match");
+                if(tableIndex == constantEntryIndex){
+                    sampleModelEntry->setValue(storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parModelEntry->getValue())));
+                } else {
+                    this->probabilityMapping.emplace_back(std::make_pair(&(this->probabilityEvaluationTable[tableIndex].second), sampleModelEntry));
                 }
-                else {
-                    std::pair<ParametricType, ConstantType> searchedPair(parEntry.getValue(), storm::utility::zero<ConstantType>());
-                    auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair);
-                    STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly");
-                    this->probabilityMapping.emplace_back(std::make_pair(&(tableIt->second), samEntry));
+                ++sampleModelEntry;
+                ++parModelEntry;
+            }
+            //also do this for the rewards thing
+            if(this->computeRewards){    
+                auto sampleModelStateRewardEntry = this->model->getStateRewardVector().begin();
+                for(std::size_t const& tableIndex : rewardEntryToEvalTableMapping){
+                    if(tableIndex != constantEntryIndex){
+                        this->stateRewardMapping.emplace_back(std::make_pair(&(this->rewardEvaluationTable[tableIndex].second), sampleModelStateRewardEntry));
+                    }
+                    ++sampleModelStateRewardEntry;
                 }
-                ++samEntry;
             }
-            if(this->computeRewards){
-                for(std::size_t state=0; state<parametricModel.getNumberOfStates(); ++state){
-                    if(!this->parametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){
-                        std::pair<ParametricType, ConstantType> searchedPair(parametricModel.getStateRewardVector()[state], storm::utility::zero<ConstantType>());
-                        auto const tableIt= std::lower_bound(evaluationTable.begin(), evaluationTable.end(), searchedPair);
-                        STORM_LOG_THROW((*tableIt==searchedPair), storm::exceptions::UnexpectedException, "Could not find the current pair in the evaluationTable. Either the table is missing that pair or it is not sorted properly");
-                        this->stateRewardMapping.emplace_back(std::make_pair(&(tableIt->second), &(this->model->getStateRewardVector()[state])));
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::initializeProbabilities(storm::models::sparse::Dtmc<ParametricType> const& parametricModel,
+                                                                                                                storm::storage::SparseMatrix<ConstantType>& probabilityMatrix,
+                                                                                                                std::vector<std::size_t>& matrixEntryToEvalTableMapping,
+                                                                                                                std::size_t const& constantEntryIndex) {
+            // Run through the rows of the original model and obtain the probability evaluation table, a matrix with dummy entries and the mapping between the two.
+            storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates(), parametricModel.getTransitionMatrix().getEntryCount());
+            matrixEntryToEvalTableMapping.reserve(parametricModel.getTransitionMatrix().getEntryCount());
+            std::size_t numOfNonConstEntries=0;
+            for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
+                ConstantType dummyEntry=storm::utility::one<ConstantType>();
+                for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
+                    if(this->parametricTypeComparator.isConstant(entry.getValue())){
+                        matrixEntryToEvalTableMapping.emplace_back(constantEntryIndex);
+                    } else {
+                        ++numOfNonConstEntries;
+                        std::size_t tableIndex=storm::utility::vector::findOrInsert(this->probabilityEvaluationTable, std::pair<ParametricType, ConstantType>(entry.getValue(), storm::utility::zero<ConstantType>()));
+                        matrixEntryToEvalTableMapping.emplace_back(tableIndex);
                     }
+                    matrixBuilder.addNextValue(row,entry.getColumn(), dummyEntry);
+                    dummyEntry=storm::utility::zero<ConstantType>();
+                }
+            }
+            this->probabilityMapping.reserve(numOfNonConstEntries);
+            probabilityMatrix=matrixBuilder.build();
+        }
+        
+        template<typename ParametricType, typename ConstantType>
+        void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::initializeRewards(storm::models::sparse::Dtmc<ParametricType> const& parametricModel,
+                                                                                                          boost::optional<std::vector<ConstantType>>& stateRewards,
+                                                                                                          std::vector<std::size_t>& rewardEntryToEvalTableMapping,
+                                                                                                          std::size_t const& constantEntryIndex) {
+            // run through the state reward vector of the parametric model. Constant entries can be set directly. Parametric entries are inserted into the table
+            std::vector<ConstantType> stateRewardsAsVector(parametricModel.getNumberOfStates());
+            rewardEntryToEvalTableMapping.reserve(parametricModel.getNumberOfStates());
+            std::size_t numOfNonConstEntries=0;
+            for(std::size_t state=0; state<parametricModel.getNumberOfStates(); ++state){
+                if(this->parametricTypeComparator.isConstant(parametricModel.getStateRewardVector()[state])){
+                    stateRewardsAsVector[state] = storm::utility::regions::convertNumber<CoefficientType, ConstantType>(storm::utility::regions::getConstantPart(parametricModel.getStateRewardVector()[state]));
+                    rewardEntryToEvalTableMapping.emplace_back(constantEntryIndex);
+                } else {
+                    ++numOfNonConstEntries;
+                    stateRewardsAsVector[state] = storm::utility::zero<ConstantType>();
+                    std::size_t tableIndex=storm::utility::vector::findOrInsert(this->rewardEvaluationTable, std::pair<ParametricType, ConstantType>(parametricModel.getStateRewardVector()[state], storm::utility::zero<ConstantType>()));
+                    rewardEntryToEvalTableMapping.emplace_back(tableIndex);
                 }
             }
+            this->stateRewardMapping.reserve(numOfNonConstEntries);
+            stateRewards=std::move(stateRewardsAsVector);
         }
 
         template<typename ParametricType, typename ConstantType>
@@ -99,11 +127,16 @@ namespace storm {
 
         template<typename ParametricType, typename ConstantType>
         void SparseDtmcRegionModelChecker<ParametricType, ConstantType>::SamplingModel::instantiate(std::map<VariableType, CoefficientType>const& point) {
-            //write entries into evaluation table
-            for(auto& tableEntry : this->evaluationTable){
+            //write entries into evaluation tables
+            for(auto& tableEntry : this->probabilityEvaluationTable){
                 tableEntry.second=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(
                         storm::utility::regions::evaluateFunction(tableEntry.first, point));
             }
+            for(auto& tableEntry : this->rewardEvaluationTable){
+                tableEntry.second=storm::utility::regions::convertNumber<CoefficientType, ConstantType>(
+                        storm::utility::regions::evaluateFunction(tableEntry.first, point));
+            }
+            
             //write the instantiated values to the matrix according to the mappings
             for(auto& mappingPair : this->probabilityMapping){
                 mappingPair.second->setValue(*(mappingPair.first));
diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h
index ca719acb3..9eb1c2210 100644
--- a/src/modelchecker/region/SamplingModel.h
+++ b/src/modelchecker/region/SamplingModel.h
@@ -10,6 +10,7 @@
 
 #include "src/modelchecker/region/SparseDtmcRegionModelChecker.h"
 #include "src/models/sparse/Dtmc.h"
+#include "src/storage/SparseMatrix.h"
 
 namespace storm {
     namespace modelchecker{
@@ -47,16 +48,19 @@ namespace storm {
             
         private:
             
+            void initializeProbabilities(storm::models::sparse::Dtmc<ParametricType> const& parametricModel, storm::storage::SparseMatrix<ConstantType>& probabilityMatrix, std::vector<std::size_t>& matrixEntryToEvalTableMapping, std::size_t const& constantEntryIndex);
+            void initializeRewards(storm::models::sparse::Dtmc<ParametricType> const& parametricModel, boost::optional<std::vector<ConstantType>>& stateRewards, std::vector<std::size_t>& rewardEntryToEvalTableMapping, std::size_t const& constantEntryIndex);
             
             //Vector has one entry for every (non-constant) matrix entry.
             //pair.first points to an entry in the evaluation table,
             //pair.second is an iterator to the corresponding matrix entry
             std::vector<std::pair<ConstantType*, typename storm::storage::SparseMatrix<ConstantType>::iterator>> probabilityMapping;
-            std::vector<std::pair<ConstantType*, ConstantType*>> stateRewardMapping;
+            std::vector<std::pair<ConstantType*, typename std::vector<ConstantType>::iterator>> stateRewardMapping;
             
             //Vector has one entry for every distinct, non-constant function that occurs somewhere in the model.
             //The second entry should contain the result when evaluating the function in the first entry.
-            std::vector<std::pair<ParametricType, ConstantType>> evaluationTable;
+            std::vector<std::pair<ParametricType, ConstantType>> probabilityEvaluationTable;
+            std::vector<std::pair<ParametricType, ConstantType>> rewardEvaluationTable;
             
             //The model with which we work
             std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> model;
diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
index d7af4d61e..5bc5ba4a7 100644
--- a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
+++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
@@ -91,7 +91,8 @@ namespace storm {
                     initializeApproximationModel(*this->simplifiedModel);
                 }
                 //now create the model used for Sampling
-                if(storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE || (storm::settings::regionSettings().getApproxMode()==storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST)){
+                if(storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE ||
+                        (!storm::settings::regionSettings().doSample() && storm::settings::regionSettings().getApproxMode()==storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST)){
                     initializeSamplingModel(*this->simplifiedModel);
                 }
                 //Check if the reachability function needs to be computed
@@ -295,7 +296,7 @@ namespace storm {
             labeling.addLabel("sink", std::move(sinkLabel));
             // other ingredients
             if(this->computeRewards){
-                storm::utility::vector::selectVectorValues(stateRewards.get(), subsystem, stateRewards.get());
+                storm::utility::vector::selectVectorValues(stateRewards.get(), subsystem);
                 stateRewards->push_back(storm::utility::zero<ParametricType>()); //target state
                 stateRewards->push_back(storm::utility::zero<ParametricType>()); //sink state
             }
@@ -379,6 +380,9 @@ namespace storm {
                 // first.
                 storm::utility::vector::selectVectorValues(stateRewards, maybeStates, this->model.getStateRewardVector());
             }
+            for(auto& stateReward: stateRewards){
+                storm::utility::simplify(stateReward);
+            }
         }
 
 
diff --git a/src/utility/vector.h b/src/utility/vector.h
index 634d2bc5a..d9163d5c8 100644
--- a/src/utility/vector.h
+++ b/src/utility/vector.h
@@ -21,6 +21,27 @@ namespace storm {
     namespace utility {
         namespace vector {
             
+            /*!
+             * Finds the given element in the given vector.
+             * If the vector does not contain the element, it is inserted (at the end of vector).
+             * Either way, the returned value will be the position of the element inside the vector,
+             * 
+             * @note old indices to other elements remain valid, as the vector will not be sorted.
+             * 
+             * @param vector The vector in which the element is searched and possibly insert
+             * @param element The element that will be searched for (or inserted)
+             * 
+             * @return The position at which the element is located
+             */
+            template<class T>
+            std::size_t findOrInsert(std::vector<T>& vector, T&& element){
+                std::size_t position=std::find(vector.begin(), vector.end(), element) - vector.begin();
+                if(position==vector.size()){
+                    vector.emplace_back(std::move(element));
+                }
+                return position;
+            }
+            
             /*!
              * Sets the provided values at the provided positions in the given vector.
              *
@@ -64,6 +85,25 @@ namespace storm {
                 }
             }
             
+            /*!
+             * Selects the elements from the given vector at the specified positions and writes them consecutively into the same vector.
+             * The size of the vector is reduced accordingly.
+             * @param vector The vector from which to select the elements and into which the selected elements are to be written.
+             * @param positions The positions at which to select the elements from the values vector.
+             */
+            template<class T>
+            void selectVectorValues(std::vector<T>& vector, storm::storage::BitVector const& positions) {
+                uint_fast64_t oldPosition = 0;
+                for (auto position : positions) {
+                    if(oldPosition!=position){
+                        vector[oldPosition++] = std::move(vector[position]);
+                    } else {
+                        ++oldPosition;
+                    }
+                }
+                vector.resize(oldPosition);
+            }
+            
             /*!
              * Selects groups of elements from a vector at the specified positions and writes them consecutively into another vector.
              *